solving a second spatial derivative using ode45

3 views (last 30 days)
I am trying to solve this above system using ode45 but I am confused on how to code for solving the second spatial derivative. The way this problem is set up is a 1D cabel model with a certain set of nodes, spacing between the nodes, and specific nodes where a stimulus is applied. So the above ODE's are being solved at each node.
function dUdt = fhn(t,U,stim,numNodes,dx,stimNodes,u0)
% FitzHugh-Nagumo equations defined to be plugged into ode45
%
% INPUT: t time
% U vector that holds u and v
% numNodes number of nodes along the fiber
% dx spatial step, spacing between nodes
% stimNodes set nodes where stim is applied
%
% OUTPUT: dUdt vector containing solutions to partial derivatives of u and
% v at each node
% Other needed variables
a = 0.13;
b = 0.013;
c1 = 0.26;
c2 = 0.1;
D = 5;
% There are 2 ODE's per node
% Using no-flux boundary conditions on the two ends of the cable
% At the first and last nodes, use a special version of the diffusion term
% in which du/dx is 0
for i = 1:numNodes
% If at the first node solve the system using no flux boundary conditions
if i == 1
dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+D*((U(1)-u0)/dx^2);
dvdt = b*(U(1)-U(2));
% If at the last node then solve with no flux boundary conditions
elseif i == numNodes
dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+D*((U(1)-u0)/dx^2);
dvdt = b*(U(1)-U(2));
% If at one of the specified stimulus nodes then the stim term needs to
% be incorporated
elseif i == stimNodes
dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+ stim;
dvdt = b*(U(1)-U(2));
% Otherwise just solve the system normally without stimulus or no flux
% boundary conditions
else
dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+(U(1(i+1))/dx^2);
dvdt = b*(U(1)-U(2));
end
end
% Create a two element vector that holds the derivative equations of u and
% v
% U(1) is u and U(2) is v
dUdt = [dudt; dvdt];
Above is my code so far, but I am confused on how to solve the second spatial derivative since the equation given to solve for it on a 1D cable model uses central finite differences. The equation is ( u(n+1) + u(n-1) - 2u ) / dx^2 where u corresponds to U(1) in the code and n corresponds to the number of nodes.

Answers (2)

J. Alex Lee
J. Alex Lee on 1 Mar 2020
This comment is incorrect:
% Create a two element vector that holds the derivative equations of u and
% v
% U(1) is u and U(2) is v
dUdt = [dudt; dvdt];
You say: "So the above ODE's are being solved at each node", but that's not quite right either. You are discretizing your spatial derivatives into finite differences (you're wanting to implement a method of lines), so that you have numNodes equations for the discretized , but is still a scalar. So in total, your fhn() function should return a vector of length numNodes+1.

darova
darova on 1 Mar 2020
You difference scheme is wrong
dudt = (c1*U(1))*(U(1)-a)*(1-U(1))-(c2*U(1)*U(2))+D*((U(1)-u0)/dx^2)
Here is correct version
uu = u(i,j);
vv = v(i,j)
(u(i+1,j)-uu)/dt = c1*uu*(uu-a)*(1-uu) - c2*uu*vv + stim + D*(u(i,j+1)-2*uu+u(i,j-1))/dx^2
u(i+1,j) = LONG_EXPRESSION*dt + uu;
Then for loop will
for i = 1:m-1
for j = 2:n-1
u(i+1,j) = LONG_EXPRESSION*dt + u(i,j);
v(i+1,j) = b*(u(i,j)-v(i,j))*dt + v(i,j);
end
end
u and v - matrices
Do you have more initial/boundary conditions except of this?
% Using no-flux boundary conditions on the two ends of the cable
% At the first and last nodes, use a special version of the diffusion term
% in which du/dx is 0

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!