solving a second spatial derivative using ode45
7 views (last 30 days)
Show older comments
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.
0 Comments
Answers (2)
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.
0 Comments
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
0 Comments
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!