tricky ode system from m file

Question

Karlis 2013-1-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/58165-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.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Jan 2013-1-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/58165-tricky-ode-system-from-m-file#answer_70399

编辑：Jan 2013-1-7

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.

2 个评论
显示无隐藏无

Karlis 2013-1-7

编辑：Karlis 2013-1-7

在 MATLAB Online 中打开

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 2013-1-7

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

Thanks all. All problems gone.

请先登录，再进行评论。

Answer 2

Ryan G 2013-1-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/58165-tricky-ode-system-from-m-file#answer_70398

在 MATLAB Online 中打开

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.

2 个评论
显示无隐藏无

Karlis 2013-1-7

编辑：Karlis 2013-1-7

在 MATLAB Online 中打开

 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)