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

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!

Answer 1

[Star Strider]

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.

[Jan-Willem]

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]

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.

Answer 2

[Andrew Newell]

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];