Changing the parameters of the dynamical system in the middle of simulation, Part 2

In my previous question, I asked "suppose I want to run

f = @(t,x) [-x(1)*x(2);x(1)*x(2)-x(2);x(2)]
[t,xa]=ode45(f,[0 6], [4 0.1 0]);

this runs the system from time 0 to 6. Suppose now I want to

*run the system from time 0 to 2

*from time 2, want to run

f = @(t,x) [-0.5*x(1)*x(2);0.7*x(1)*x(2)-x(2);x(2)]

How can one write a script?" Thank you for answers.

Now I want to consider a following senario:

"Suppose I want to

*initially run the system

*but once x(2) becomes more than 1, switch to

f = @(t,x) [-0.5*x(1)*x(2);0.7*x(1)*x(2)-x(2);x(2)]

How can one write such a script?

Answer 1

[Bjorn Gustavsson]

For such conditions you might just get away with simply changing the output of the f-function. This is much easier to achieve if you write f as a proper matlab-function. Something like:

function dSIRdt = odeSIRvariableF(t,y,a,b,level)

dSIRdt = zeros(3,1);
dSIRdt(1) = -a*y(1)*y(2);
if y(2) < level
  dSIRdt(2) = a*y(1)*y(2) - b(1)*y(2);
  dSIRdt(3) = b(1)*y(2);
else
  dSIRdt(2) = a*y(1)*y(2) - b(2)*y(2);
  dSIRdt(3) = b(2)*y(2);
end

If this doesn't work nicely you should have a look at ballode - where you can see how to handl the case of a bouncing ball.

HTH

[alpedhuez]

so you mean that one deals with this issue with the definition of functions? that is "changing the intput of the f-function"?

[alpedhuez]

Just to make sure, will this definition ensure continuity of path at y(2) = level?

[Bjorn Gustavsson]

Yeah, you can define derivatives with just any functional dependency on time and values of the dependent parameter, some will be much more difficult for the odeNN to integrate (that's why there are a number of them with different strengths and applicabilities). The main problem is if you have discrete changes of the independent variables - have a look at ballode - it is a demo so you should be able to first just run this. This modification of yours should be a very minor modification, it is only a change in the derivatives not the dependent parameters, so should be OK. If you want to make sure that the ODE-solution exactly hit your transition-level you might have to use the events-trick you'll find in ballode. And yes this change will not make the solution discontinuous. You should see a discontinuity of the time-derivative when the solution moves across the transition - but that is what you aim for.

[alpedhuez]

Thank you. Then after the process I want to calculate \int_{0}^{endtime}y(1,t)dt. How can I do that? Thank you.

[Bjorn Gustavsson]

Since you get the solution as two arrays t and y you can always calculate the integral numericaly with trapz.

[alpedhuez]

But I meant

\int_{0}^{endtime}(dy(1,t)/dt)dt

This is just y(endtime).

[alpedhuez]

or

T=0:1/100:2;
trapz(T,y(1))

?

[Bjorn Gustavsson]

Yes, because you dont feed trapz the array...

You do something like this:

[t,y] = ode45(@(t,y) yourodefcn(t,y,pars),[0,6],y0);

where you'll get t and y, check the sizes of those:

whos t y

to integrate you call trapz with 2 arrays:

y_int = trapz(t,y(:,1)); % if that is the dimension you should integrate over

hth

[alpedhuez]

I see, trapz looks at the x axis and y axis and calcuate integral, Makes sense.

[alpedhuez]

Well, another question. How can I plot this time derivative (d/dt)y(1) instead of y(1)?

t=[0 1 2 3 4 5] % Vector time
y=[1 2 4 8 12 13] % Your signal
% The derivative can be approximated :
diff(y)./diff(t)

or is it recorded somewhere in the calculation?

[Bjorn Gustavsson]

That is one way to calculate the derivatives. If you look at the help for diff you will also be adviced to look at the gradient function that you can use. You should make a habit of reading the help and documentation of the functions you use and pay attention to the SEE also part, another useful tool for finding functions is lookfor.

Since you have an ode-function that calculates derivatives as a function of t and y you should be able to figure out some way to use that to your advantage.

[alpedhuez]

Thank you. I will take a look.