tricky ode system from m file

Starting from t=0 there is one equation system, then after some time, constraint is met and other system is necessary, then this 2nd system meets constrain and changes back to first and so on. I tried some ideas, but they didn't quite work.

Constraints are simple: once x(2) reaches X2 other system should turn in and once x(3) reaches X3 1st system should turn back in.

 function xp=uzdevums1(t,x)
 %parameters
   r1 = 0.1; r2 = 1; r3 = 0.2;K1=100;K2 = 100; K3 = 100;X2=25;X3=10;speedx2 = 0.02; speedx3=0.02;

 xp=zeros(3,1);

 if (x(2)>X2)&&(xp(3)>=0)
    xp(1)=r1*(1-x(1)/(x(2)+x(3)))*x(1);
    xp(2)=r2*(1-x(2)/K2)*x(2)-x(2)*x(1)*speedx2; 
    xp(3)=r3*(1-x(3)/K3)*x(3);
 elseif (x(3)>X3)&&(xp(2)>=0)
    xp(1)=r1*(1-x(1)/(x(2)+x(3)))*x(1);
    xp(2)=r2*(1-x(2)/K2)*x(2); 
    xp(3)=r3*(1-x(3)/K3)*x(3)-x(3)*x(1)*speedx3;
 end

This is one way that doesn't quite work. I hope there is some easier way to do this, because i don't see how to do it in the hard way.

Answer 1

[Jan]

And you want to integrate uzdevums1 by ODE45 or a similar function?

Do not do this. The ODE-integrators have a stepsize control, which is substantially confused by discontinuities. In consequence the number of steps and function evaluations can be increased very much, such that the accumulated rounding errors will influence the results.

Use event functions to get the switching points, stop the integrator and restart it with the former final values. See doc odeset.

[Karlis]

OK. So I should define f1 - 1st ode system; f2 - 2nd ode system; events1 - events for 1st system; events2 - events for 2nd system;

and then do something like

 tstart = 0;
 tfinal = 50;
 y0 = [10; 80; 20]; %initial conditions 

 options = odeset('Events',@events1);

 for i = 1:5 %(cycles)
  [t,y,te,ye,ie] = ode45(@f1,[tstart tfinal],y0,options);

nt = length(t);

y0 = [y(nt,1),y(nt,2),y(nt,3)];%new initial conditions 
tstart = t(nt);

options = odeset(options,'Events',@events2);
[t,y,te,ye,ie] = ode45(@f2,[tstart tfinal],y0,options);

nt = length(t);

y0 = [y(nt,1),y(nt,2),y(nt,3)];%new initial conditions 

options = odeset(options,'Events',@events1);
tstart = t(nt);
end

This seem to work, but i don't get how to correctly define event functions.

 function [value,isterminal,direction] = events1(t,y)
value = y(2);  % this part not clear for me.
isterminal = 1;   % Stop the integration
direction = -1;   % Negative direction only

How do i make it stop integrating when y(2) reaches some number, like 25?

[Karlis]

OK, figured it out. value = y(2)-25; %works

Thanks all. All problems gone.

Answer 2

[Ryan G]

What do you mean doesn't work? It sounds like you need a new state to define what condition you are currently in. This way you will have a default operating point instead of using the if-else condition every time.

Technically I would use Simulink/Stateflow to accomplish this and it would be ideal. However, within the bounds of the ode solver you could do something like this:

if(myState == true)
  %ode equations
  myState = x(3)>X3)&&(xp(2)>=0;
else
  %other ode equations
  myState = (x(2)>X2)&&(xp(3)>=0);
end

This will keep you in the current state until the condition is met to go back to the previous (or next) state.

[Karlis]

 if t==0
    myState = 1;
 end

 if (myState == 1)
    ode
    if (condition to go in 2nd state)
        myState =2;
    end
 end

 if (myState == 2)
    ode2
    if (condition to go in 1st state)
        myState =1;
    end   
end

Something like this? I thought about it, but apparently MATLAB doesn't remember value of myState in previous time step and produces some error. (I am using ode45)

[Jan]

You could store the state persistenly, see "help persistent". But as explained in my answer, I recommend to avoid this.