The program should output a time vector and separate 3D arrays for u and v. The first dimension of the u and v arrays should be time and the second and third should be the two spatial dimensions.
solving 2D FitzHugh-Nagumo system
7 views (last 30 days)
Show older comments
function [U,V,t] = fhn1d(tend,dx,stimMask,stimStr,stimDur,U0,V0,t0)
% function [U,V,t] = fhn1d(tend,dx,stimMask,stimStr,stimDur,U0,V0,t0)
% simulates propagation of a FitxHugh-Nagumo cable
%
% INPUT: tend time to integrate to
% dx spatial time step in nm
% stimMask length of this vector specifies how many nodes are in
% the cable and specifies which nodes will receive a
% stimulus current
% stimStr strength of current stimulus
% stimDur duration of current stimulus
% U0 vector the same length of stimMask that contains the
% initial condition for u on each node
% V0 vector the same length of stimMask that contains the
% initial condition for v on each node
% t0 time at which to begin integrating
%
% OUTPUT: U contains solution where it has a column for each node on
% the cable and a row for each timestep
% V contains solution where it has a column for each node on
% the cable and a row for each timestep
% t vector with the time of each timestep
% Set argument defaults
% The first three arguments are required (others are optional)
if nargin < 8 || isempty(t0)
t0 = 0;
end
if nargin < 7 || isempty(V0)
V0 = zeros(size(stimMask));
end
if nargin < 6 || isempty(U0)
U0 = zeros(size(stimMask));
end
if nargin < 5 || isempty(stimDur)
stimDur = 0;
end
if nargin < 4 || isempty(stimStr)
stimStr = 0;
end
% Set membrane parameters and cable conductivity
a = 0.13;
b = 0.013;
c1 = 0.26;
c2 = 0.1;
D = 5;
% stimMask is a column vector
stimMask = stimMask(:);
% Number of nodes
nx = length(stimMask);
% Stack up initial values for u and v into a long column vector
% All u values are first and then the v values
UV0 = [U0(:);V0(:)];
% Set time at which stimulus turns off
stimOff = stimDur + t0;
% Make a handle to an anonymous function that will evaluate the
% right-hand-sides of all the ODEs in the system
fhnHandle = @(t,UV) fhn1(t,UV,stimStr,stimOff,stimMask,dx,a,b,c1,c2,D);
% Call ode 45 to do the integration
% On output, t is a column vector of times
% UV has one column for each ODE in the system and one row per timestep
% fhn1 is written so that the first nx columns in UV contain u values going
% down the cable and the next nx columns contain the v values
[t,UV] = ode45(fhnHandle,[t0, tend],UV0);
% Unpack UV into seperate U and V arrays
U = UV(:,1:nx);
V = UV(:,nx+1:end);
plotFiber(U,V,t)
end
function dUVdt = fhn1(t,UV,stimStr,stimOff,stimMask,dx,a,b,c1,c2,D)
% function dUVdt = fhn1(t,UV,stimStr,stimOff,stimMask,dx,a,b,c1,c2,D)
% evaluates the right-hand-sides of a system of ODEs that stimulate
% propagation on an FHN cable; called by ode45()
%
% INPUT: t current time
% UV column vector with the current values of the dependent
% variables
% stimStr strength of current stimulus
% stimOff time at which the stimulus switches off
% stimMask vector that specifies which nodes receive the stimulus
% dx spatial time step in nm
% a,b,c1,c2,D model parameters
%
% OUTPUT: dUVdt column vector containing solutions
% Get number of nodes
nx = length(stimMask);
% Unpack u and v
u = UV(1:nx);
v = UV(nx+1:end);
% Compute the diffusive term for each node, d2u/dx2, and store in column
% vector L
% Preallocate L
L = zeros(nx,1);
% Special term for first node that imposes no-flux boundary
L(1) = D*(u(2)-u(1))/(dx*dx);
% Diffusive term for all interior nodes
for i = 2:nx-1
L(i) = D*(u(i-1) + u(i+1) - 2*u(i))/(dx*dx);
end
% No-flux boundary condition for last node
L(nx) = D*(u(nx-1)-u(nx))/(dx*dx);
% Make a vector that contains the current stimulus for each node
if t < stimOff
IStim = stimMask*stimStr;
else
IStim = stimMask*0;
end
% Get fhn current for each node
I = c1.*u.*(u-a).*(1-u)-c2*u.*v;
% Add all the terms together to give dUdt on each node
dUdt = L + I + IStim;
% Compute dVdt
dVdt = b*(u-v);
% Stack all of our dUdt and dVdt values into one column vector
dUVdt = [dUdt;dVdt];
end
How do I take my code for the solving the fitzhugh-nagumo system in one dimension and solve it two dimensions according to the equations below?
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/275382/image.jpeg)
The membrane parameters (a,b,c1,c2) are the same. D11 and D22 are conductivity constants in the x and y directions respectively; both can be set to five. How do you specify the number of nodes in the sheet in both directions? The spatial steps dx and dy should be input arguments to the function. How do you do the no flux boundary conditions for d2u/dx2 on the left and right edges of the sheet and the top and bottom edges for d2u/dy2? Then once the diffusive terms are figured out for this system, how do you pass them into ode45? Overall, I want to know how to convert my 1d version into what needs to be solved for the 2d version.
5 Comments
Answers (1)
darova
on 7 Mar 2020
- How do you do the no flux boundary conditions for d2u/dx2 on the left and right edges of the sheet and the top and bottom edges for d2u/dy2?
It's question to you. What are boundary conditions (temperature or heat flux) ?
- Then once the diffusive terms are figured out for this system, how do you pass them into ode45? Overall, I want to know how to convert my 1d version into what needs to be solved for the 2d version.
Once you have "2D version" you can't use ode45
I attach similar problem solved using difference scheme. See if it helps
2 Comments
darova
on 7 Mar 2020
ode45 (Ordinary Differencial Equation) is for ordinary diff equations
You have partial differential equations (function depends on 2 or more variables)
You can't solve your equations using ode45
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!