using ode45

Dear all,

I would like to solve the following system of differential equations with input with Matlab. I tried with ODE45. Unfortunately, I always get a very strange behaviour of the solution, almost insensitive to the input. The values of the parameters (taus tauf etc..) are taken from a work where integration of the same equations was succesfull, and are supposed to be meaningful. The input should be an impulse at given time. The solution should look null before the impulse, increasing afterwards, reaching zero again after about 20 sec.

It is probably some stupid mistake but I cannot get it right. Any suggestion?

Thanks in advance, Sara

%%---------------------------------

pnts = 7680;

dt = 0.0039;

x0 = zeros(1,2);

u = zeros(pnts,1);

u(768) = 1;

t = [1 pnts]*dt;

[T, SOL] = ode45(@fx, t, x0,[],u,dt);

%% ----------------------------

function dx = fx(t,x,u,dt)

ui = u(ceil(t/dt));

x = x';

taus = 0.8;

tauf = 0.4;

epsilon = 0.5;

dx(:,1) = epsilon*ui-x(:,1)./taus-(x(:,2)-1)./tauf;

dx(:,2) = x(:,2);

dx = dx';

end

%%---------------------------------------------------------------

[Andrew Newell]

It might help if you tell us what the equation is that you are trying to solve (in case you're not representing it right).

[Sara]

The system of equations describes the haemodynamic forward model for fMRI, even though I reported here only the first subblock, which is called "the balloon component". I do not know if this can help...

[Andrew Newell]

No, I mean something like

dx1/dt = ...

dx2/dt = ...

[Sara]

Oh.. sorry. Here it is:

dx1/dt = epsilon*u-x1/taus-(x2-1)/tauf;

dx2/dt = x1;

Hope I get what you mean this time..

PS: there is a typo in my original question: in my code dx(:,2) is indeed equal to x(:,1);

[Andrew Newell]

What is u? You have it equal to zero except at one point.

[Sara]

u is an impulse that is one only at one time point and zero otherwise. I did try to use eps instead of zero, but with no difference in the results. I also try with u=randn(pnts,1), but still i do not get the output I expect. I am wondering whther it is a problem of transient behaviour of the equations or ode45 is not the appropriate solver to use...

Answer 1

[Andrew Newell]

See How Do I Obtain an Impulse Response? for the right way to model impulses.

EDIT: Here is how to apply it to your problem, assuming that your impulse occurs at t=3.

First, integrate the impulse-free equation for the first three seconds. Your equations can be represented by a simple online function:

taus = 0.8;
tauf = 0.4;
fx = @(t,x) [-x(1)./taus-(x(2)-1)./tauf; x(1)];

x0 = zeros(1,2);
tspan = [0 3];
[T, SOL] = ode45(fx, tspan, x0);

Now add the impulse at t=3 and use the result as the initial condition for the next seven seconds:

epsilon = 0.5;
x0 = SOL(end,:) + [epsilon 0];

tspan = [3 10];
[T1, SOL1] = ode45(fx,tspan,x0);

Finally, combine the two segments to get your solution:

T = [T; T1]; SOL = [SOL; SOL1];

[Sara]

Thanks a lot for your detailed and elegant answer. I am going to look into that more carefully but after a quick try, it seems to me that I do still have an issue with transient behaviour before reaching the equilibrium. Do you know how I can upload an image? It would be easier to explain what my problem is...

[Andrew Newell]

See http://www.mathworks.com/matlabcentral/answers/7924-where-can-i-upload-images-and-files-for-use-on-matlab-answers.

[Sara]

I ran the code you suggested a couple of times. IN this figure (http://fileslap.com/93x/ComparisonImpulseTime) you can see the results (green is x1, blue is x2; thick is with an impulse at 3; dashed is with impulse at 10). The dashed lines look correct. So, my understanding is that the system as some kind of transient, that I have to wait before giving the impulse. I guess I can calculate this transient time from the time costants, but this is not matlab-related I would say.

Nevertheless, I do have an additional question: if I have to simulate a train of impulses as input do I have to break the calculation into segments for each input? Is by any chance the "event" feature of the ode45 related to this matter?

Thanks

Sara

[Andrew Newell]

The only difference between the dashed output and that of my code is the time intervals. Just replace 3 by 10 for the impulse time and 10 by 20 for the total time. As for a train of impulses, if you know the time of each impulse you can create a loop to piece together multiple time segments.