Basic differential solving and plotting problem for a part of an exam.

13 views (last 30 days)
I've been given following task and got stuck. Apparently I am returning a vector of length 4 instead of the 2 it should be.
The differential equation I am trying to plot is: dv/dt=v*(1-v)(v-alpha)-w+C and dw/dt=varepsilon(v-gamma*w)
function dydt = fun (v, alpha, w, C, gamma, varepsilon)
dydt = [v*(1-v)*(v-alpha)-w+C;varepsilon*(v-gamma*w)]
end
This is the Fitzhugh-Nagumo model (<http://en.wikipedia.org/wiki/FitzHugh-Nagumo_model>)
Now I want to solve this equation via ode45 (this is a prerequisite) in the following script:
tspan = [0, 6]; y0 = [1; 1];
alpha = 0.7;
gamma = 0.8;
varepsilon = 12.5;
C = 0.5;
ode = @(v,w) fun (v, alpha, w, C, gamma, varepsilon);
[v,w] = ode45(ode, tspan, y0);
plot(t,v(:,1),'r')
xlabel ('t')
ylabel('solution v & w')
title ('opdracht 1, name')
hold on
plot(t,w(:,1),'b')
shg
The goal is to plot v and w in function of t (time, defined by tspan). But I get the following error:
Error using odearguments (line 92) @(V,W)OPDRACHT1FUNCTIEV(V,ALPHA,W,C,GAMMA,VAREPSILON) returns a vector of length 4, but the length of initial conditions vector is 2. The vector returned by @(V,W)OPDRACHT1FUNCTIEV(V,ALPHA,W,C,GAMMA,VAREPSILON) and the initial conditions vector must have the same number of elements.
Can anyone tell me where these 4 solutions come from whilst there only should be 2? The handle @(v,w) contains but two values, does this not suffice? Thanks in advance!

Answers (2)

Star Strider
Star Strider on 3 May 2014
The ODE equations have to be functions of (t,y) and they have to return a column vector t and a vector or matrix of y values.
I had to change your code somewhat (incorporating some assumptions in the process), but it now integrates and plots:
% Define vw(1) = v, vw(2) = w
fun = @(t, vw, alpha, C, gamma, varepsilon) [vw(1)*(1-vw(1))*(vw(1)-alpha)-vw(2)+C; varepsilon*(vw(1)-gamma*vw(2))];
tspan = [0, 6]; y0 = [1; 1];
alpha = 0.7;
gamma = 0.8;
varepsilon = 12.5;
C = 0.5;
ode = @(t,vw) fun (t, vw, alpha, C, gamma, varepsilon);
[t,vw] = ode45(ode, tspan, y0);
plot(t,vw(:,1),'r')
xlabel ('t')
ylabel('solution v & w')
title ('opdracht 1, name')
hold on
plot(t,vw(:,2),'b')
shg
If I made incorrect guesses or assumptions and this is not producing the results you want, please clarify. We can sort this.
  2 Comments
Jan-Willem
Jan-Willem on 4 May 2014
Thanks for your help! Your script seems to work fine, but the result is not like I would expect it to be. The graph should look somewhat like this: http://en.wikipedia.org/wiki/File:Fitzhugh_Nagumo_Model_Time_Domain_Graph.png even though the function is perfectly correct, there is no form of oscillation. The values on the other hand are the same as used for the wikipedia-version (with tspan[0,200]). I personally find this very strange.
Star Strider
Star Strider on 4 May 2014
My pleasure!
However something must be wrong somewhere, possibly in your derivations. This model:
ode = @(t,vw) [vw(1) - (vw(1).^3)/3 - vw(2) + 1; (vw(1) + 0.1 - 0.2.*vw(2))/1.5];
taken directly from the Wikipedia article on the Fitzhugh-Nagamo model, oscillates normally, even with the completely arbitrary constants I chose for it.
I don’t know the details of what you’re doing or the details of your derivation, so I can’t suggest a solution.

Sign in to comment.


Andrew Newell
Andrew Newell on 3 May 2014
When you add some scalars to a vector, you get a vector, so
[v*(1-v)*(v-alpha)-w+C; varepsilon*(v-gamma*w)]
stacks 2-vectors to get a 4-vector. As a simple example, see what happens when you run this:
w=[1; 1];
[1-w; 2*w];

Products

Community Treasure Hunt

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

Start Hunting!