Switch between 3 different ODE using event functions

Hello, I have a system of 3 different equations of motion of second order which I set up as ODEs. These ODEs should toggle on and off depending on the time dependent parameters y(1) & y(3) as seen below, using 1 or more event functions. But this doesn't work with my current code. It only calculates the last of my 3 ODEs. If I delete the last one it calculates only the 2nd one and if I delete the 2nd and 3rd one it doesn't seem to work at all (which is another problem). I need it to automatically switch back and forth between the 3 ODEs. In the end I need a continous plot of my desired variables as also seen below.

Note: I only now that eqmi is the starting equation but I don't know which of the other 2 equation comes next and the time at which the event occurs. It also can switch back and forth between the 3 ODEs multiple times in one timespan.

This was my first question regarding my problem and might help understanding what I'm actually doing, just theoretically, without events. This was my 2nd question, where the answers helped me a lot and shows some of the progress I made but it still doesn't work as needed.

Can anyone help me out with this? Down below is my current Matlab-Code.

function mfc
clc
%%parameters
    m1=10; m2=10; m3=500; k3=300; k2=(k3*(m1+m2))/m3; my=m1/(m1+m2); FW=10; s=0.4; r=0.4; om=5.73; dv0s=2; dv0=dv0s-my*r*om;
%%starting parameters

      tspan = [0 100];
      tstart = tspan(1);
      tend = tspan(end);
      y0 = [0 dv0 0 0];

      tout = tstart;    % copied from ballode example
      yout = y0;
      teout = [];
      yeout = [];
      ieout = [];

   while tout < tend   % main loop
      [t,y,te,ye,ie] = ode45(@(t,y) eqmi(t,y), [tstart tend], y0, odeset('Events', @ZustandI));
      [t,y,te,ye,ie] = ode45(@(t,y) eqmii(t,y), [tstart tend], y0, odeset('Events', @ZustandII));
      [t,y,te,ye,ie] = ode45(@(t,y) eqmiii(t,y), [tstart tend], y0, odeset('Events', @ZustandIII));

      nt = length(t);        % Accumulate output. Copied from ballode example
      tout = [tout; t(2:nt)];
      yout = [yout; y(2:nt,1:4)];
      teout = [teout; te];
      yeout = [yeout; ye];
      ieout = [ieout; ie];

      y0 = [y(nt,1) y(nt,2) y(nt,3) y(nt,4)];  % new IC and tspan
      tstart = t(nt);
  end
%%Plots
      figure(1);                             % displacement
      plot(teout,yeout(:,1),'b--',teout,yeout(:,3),'r:')
      xlabel('time');
      ylabel('displacement');
      title('mfc');
      legend('v2','v3');

      figure(2);                             % velocity
      plot(teout,yeout(:,2),'r:',teout,yeout(:,4),'b--')
      xlabel('time');
      ylabel('velocity');
      title('mfcvelocity');
      legend('v2dot','v3dot');

      figure(3);                             % v2-v3
      plot(teout,yeout(:,1)-yeout(:,3),'r:')
      xlabel('time');
      ylabel('v2-v3');
      title('mfcunterschied');
      legend('v2-v3');
%%functions
   function dy=eqmi(t,y)
      dy=zeros(4,1);
      dy(1)=y(2);                                             % v2dot
      dy(2)=my*r*om^2*sin(om*t);                              % Gl. 2a
      dy(3)=y(4);                                             % v3dot
      dy(4)=(FW*sin(om*t)-k3*y(3))/m3;                        % Gl. 3a
  end

function dy=eqmii(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                             % v2dot
    dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)-s))/(m1+m2);   % Gl. 2b II
    dy(3)=y(4);                                             % v3dot
    dy(4)=(FW*sin(om*t)-k3*y(3)+k2*(y(1)-y(3)-s))/m3;       % Gl. 3b II
end

function dy=eqmiii(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                             % v2dot
    dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)+s))/(m1+m2);   % Gl. 2b III
    dy(3)=y(4);                                             % v3dot
    dy(4)=(FW*sin(om*t)-k3*y(3)+k2*(y(1)-y(3)+s))/m3;       % Gl. 3b III
end   
%%Events
function [value,isterminal,direction] = ZustandI(t,y)
value = [double(y(1)-y(3)>-s),double(y(1)-y(3)<s)];
isterminal = [1,1];                                     
direction = [-1,1];
end

function [value,isterminal,direction] = ZustandII(t,y)
value = double(y(1)-y(3)>-s);
isterminal = 1;                                     
direction = 0;                                      
end

function [value,isterminal,direction] = ZustandIII(t,y)
value = double(y(1)-y(3)<s);
isterminal = 1;                                     
direction = 0;                                     
end
end

Answer 1

[Jan]

[t,y,te,ye,ie] = ode45(@(t,y) eqmi(t,y), [tstart tend], y0, odeset('Events', @ZustandI));
[t,y,te,ye,ie] = ode45(@(t,y) eqmii(t,y), [tstart tend], y0, odeset('Events', @ZustandII));
[t,y,te,ye,ie] = ode45(@(t,y) eqmiii(t,y), [tstart tend], y0, odeset('Events', @ZustandIII));

Of course this overwrite t,y,te,ye,ie twice. The code computes the complete interval 3 times. But you want to switch between the ODEs.

t   = tstart;
y   = y0;
fcn = @eqmi;
opt = odeset('Events', @ZustandI);
while (t(end) < tend)
  % Run integration until event function stops it:
  % [EDITED, 2018-07-30]: "t" -> "t(end)
  [at, ay, ate, aye, aie] = ode45(fcn, [t(end), tend], y(end, :), opt);

    % Append the new trajectory:
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));

    % Decide for new equation and event function:
    if y(end, :) > ???? 
      fcn = @eqmi;
      opt = odeset('Events', @ZustandI);
    elseif y(end, :) < ???
      fcn = @eqmii;
      opt = odeset('Events', @ZustandII);
    else
      fcn = @eqmiii;
      opt = odeset('Events', @ZustandIII);
    end
  end

I have no idea how to decide for the specific state, but the idea is to integrate only one each time until the event function stops the integration.

Answer 2

[Lennart]

Hi, so I followed your advice and changed my code to the following, but it still doesn't work correctly. The plots I'm getting are nonsensical.

For the if/else conditions I basically used the same statements as for my event functions, just the other way around (all depending on y(1)-y(3) and s). So isn't this basically a double event function because I define when my 3 ODEs should be integrated (if/else) and when they should stop integrating (event functions) at the same time?

my new code, parameters, event functions and equations of motion stay the same as above:

      tspan = [0 100];
      tstart = tspan(1);
      tend = tspan(end);
      y0 = [0 dv0 0 0];

      % set up starting equation
      t = tstart;
      y = y0;
      fcn = @eqmi;
      opt = odeset('Events', @ZustandI);

while t(end) < tend   % main loop
    % Run integration until event function stops it
    [at, ay, ate, aye, aie] = ode45(fcn, [t, tend], y(end, :), opt);

      % Append the new trajectory
      t = cat(1, t, at(2:end));
      y = cat(1, y, ay(2:end,:));

      % Decide for new equation and event function
      if -s<y(end,1)-y(end,3)<s 
        fcn = @eqmi;
        opt = odeset('Events', @ZustandI);
      elseif y(end,1)-y(end,3)<=-s
        fcn = @eqmii;
        opt = odeset('Events', @ZustandII);
      elseif y(end,1)-y(end,3)>=s
        fcn = @eqmiii;
        opt = odeset('Events', @ZustandIII);
      end
  end
  %%Plots
      figure(1);                             % displacement
      plot(t,y(:,1),'b--',t,y(:,3),'r:')
      xlabel('time');
      ylabel('displacement');
      title('mfc');
      legend('v2','v3');

      figure(2);                             % velocity
      plot(t,y(:,2),'r:',t,y(:,4),'b--')
      xlabel('time');
      ylabel('velocity');
      title('mfcvelocity');
      legend('v2dot','v3dot');

      figure(3);                             % v2-v3
      plot(t,y(:,1)-y(:,3),'r:')
      xlabel('time');
      ylabel('v2-v3');
      title('mfcunterschied');
      legend('v2-v3');

I also don't quite understand where I accumulate my new output and new starting parameters y & t in this case represented by the following code (which I copied from the ballode example from matlab) which I removed from my current code:

nt = length(t);        % Accumulate output. Copied from ballode example
    tout = [tout; t(2:nt)];
    yout = [yout; y(2:nt,1:4)];
    teout = [teout; te];
    yeout = [yeout; ye];
    ieout = [ieout; ie];

      y0 = [y(nt,1) y(nt,2) y(nt,3) y(nt,4)];  % new IC and tspan
      tstart = t(nt);

One last question for now. Can I somehowe mark the positions in the plots at which time the equations change from one to another, with a cross or something? This would greatly help understanding the calculations.

[Lennart]

Reviewing it again, it might have something to do with my plot commands. I plotted t & y which I think is incorrect (changing it to at/ay respectivly yields the same result). But if I change it to ate & aye (like in my old code), I get the error message:

Index in position 2 exceeds array bounds.

Error in Testskript2_mehrere_EventFcn (line 66)
    plot(ate,aye(:,1),'b--',ate,aye(:,3),'r:')

[Jan]

In my code at, ay, etc are the data from the current integration. They are accumulated by:

   t = cat(1, t, at(2:end));
   y = cat(1, y, ay(2:end, :));

The first element of at and ay is ignored, because this is the initial value already and contained in t and y before the integration.

By the way: cat(1, x, y] is equivalent to [x; y]

My code had a typo:

[at, ay, ate, aye, aie] = ode45(fcn, [t, tend], y(end, :), opt);
                                      ^

Better:

[at, ay, ate, aye, aie] = ode45(fcn, [t(end), tend], y(end, :), opt);
                                        ^^^^

Now the correct time span is provided.

This will not work:

if -s<y(end,1)-y(end,3)<s

Matlab computes it from left to right. First step:

-s < y(end,1) - y(end,3)

This is either true or false, which are converted to 1 or 0. Then in te second step:

0 < s   or  1 < s

I'm sure you mean:

if (-s < y(end,1) - y(end,3)) && (y(end,1) - y(end,3) < s)

Use the debugger to find out, which index is exceeded. Type this in the command window:

dbstop if error

Run the code again, and when it stops at the error, check the size of e.g. ate.

[Lennart]

I ran the debugger and ate, aye and aie have no values at all, just [ ]. But the ballode example from matlab uses those values to plot the results. And if I'm supposed to plot t,y, the linear plot of the displacements doesn't make any sense.

[Lennart]

Another Thing I found out: The if/else part doesn't do anything. I can delete it or delete parts of it, the final result is the same. The only thing part that makes any difference, is this one:

   fcn = @eqmi;
   opt = odeset('Events', @ZustandI);

I can change fcn and opt so I get slightly different results but it stays on the same equation and event for the rest of the calculation. The changing equations based on events part still doesn't work. eqmi as the starting equation is correct but it doesn't register the corraltion of y(1)-y(3) & s correctly.

[Jan]

@Lennart: Did you read my comment concerning the if condition "-s<y(end,1)-y(end,3)<s"?

[Lennart]

Yeah, it looks like this now:

      % Decide for new equation and event function
      if (-s<y(end,1)-y(end,3) && y(end,1)-y(end,3)<s)
        fcn = @eqmi;
        opt = odeset('Events', @ZustandI);
      elseif y(end,1)-y(end,3)<=-s
        fcn = @eqmii;
        opt = odeset('Events', @ZustandII);
      elseif y(end,1)-y(end,3)>=s
        fcn = @eqmiii;
        opt = odeset('Events', @ZustandIII);
      end

Edit:

I changed the directions of my first event (ZustandI) from [-1,1] to [0,0] which makes a lot more sense because its a logical expression changed to a double and the plot looks much much better. But it still doesn't look right and I'm not sure if it changes the ODEs correctly.

So I'm coming back to my previous question: How can I mark the exact occurances in my plots when it terminates one equation and changes to the next one? It'll help reviewing the code.

[Lennart]

Ok. I take it back. It does seem to work. I finally understood how to detect the event triggers. ate marks the time and aye the solution value of the events at ate. So I put this before my main loop:

    ateout = [];
    ayeout = [];
    aieout = [];

and this inside my main loop:

    ateout = [ateout; ate];
    ayeout = [ayeout; aye];
    aieout = [aieout; aie];

and finally I plot the solution like this:

    figure(1);                             % displacement
    plot(t,y(:,1),t,y(:,3))
    hold on
    plot(ateout,ayeout(:,1),'rx',ateout,ayeout(:,3),'b.')
    xlabel('time');
    ylabel('displacement');
    title('mfc');
    legend('v2','v3');

and I get a nice plot with x and . at the event triggers.

I have to check it manually if it makes sense but it looks promising. Or do you see any obvious mistakes in my code? Thank you very much for your help. It's appreciated. Or do y