Solving system of ODEs with event functions

Hello guys, I'm trying to solve a system of two different equation of motions dv and dw which change depending on the displacements v(1) and w(1). I got the tip to solve it with event functions which I think is the right way to go. I got three different conditions on which my eqm will change, leaving me with 6 different functions which you can see at the end of my code below. The problem is, it doesn't work. I don't get any error messages but the calculation will never get to an end.

I have a few ideas why it won't work but I didn't got a solution yet. So here are my questions, maybe you can help me. I tried to comment everything in my code for you:

1) I need to set a starting equation before the event functions begin to get the first values for v(1) and w(1) on which the events depend. Otherwise I get an error message (not enough input parameters). You can see it in my code just after the beginning of the main loop. I tried to use different stiff and nonstiff ODEs but everytime it just calculates for ages on [t,w] = ode15s(@(t,w) eqm3ai(t,w), tspan, w0);. Where is the problem with that and how can I correct that?

2) I want my event-functions to stop integrating if the conditions (also seen in my code down below) are met. I know that normally it stops integrating if the desired value crosses zero from the desired direction but this is a little problematic in my case because in one case it should include v(1)-w(1)=s and in another case it shouldn't include v(1)-w(1)=s. But maybe I oversee something important. How can I include my conditions seen down below?

3) Which ODE-solver is the best for me? I tried a few (stiff and nonstiff) and none worked.

4) A little redundant but is my event setup more or less correct for my case?

Thank you in advance and sorry for my bad english.

-------------------------code begins here---------------------------------------

    function mfc

    %%Parameters

    m1=4; m2=4; m3=34.6; k3=350; k2=(k3*(m1+m2))/m3; my=m1/(m1+m2); s=0.25 r=0.25; om=5.37; dv0s=5; dv0=dv0s-my*r*om; FW=20;

    tspan = [0 10];     % Initial conditions
    v0 = [0 dv0];
    w0 = [0 0];

    tout = tspan(1);    % copied from ballode example
    vout = v0.';
    wout = w0.';
    teout = [];
    veout = [];
    weout = [];
    ieout = [];
% main loop
while tout(end) <= tspan(end)
while tout(end) == 0            % starting equations, for first v(1), w(1)
    [t,v] = ode15s(@(t,v) eqm2ai(t,v), tspan, v0);
    [t,w] = ode15s(@(t,w) eqm3ai(t,w), tspan, w0);
end
    [t,v,te,ve,ie] = ode15s(@(t,v) eqm2ai(t,v), tspan, v0, odeset('Events', @ZustandI));
    [t,w,te,we,ie] = ode15s(@(t,w) eqm3ai(t,w), tspan, w0, odeset('Events', @ZustandI));
    [t,v,te,ve,ie] = ode15s(@(t,v) eqm2bii(t,v), tspan, v0, odeset('Events', @ZustandII));
    [t,w,te,we,ie] = ode15s(@(t,w) eqm3bii(t,w), tspan, w0, odeset('Events', @ZustandII));
    [t,v,te,ve,ie] = ode15s(@(t,v) eqm2biii(t,v), tspan, v0,odeset('Events', @ZustandIII));
    [t,w,te,we,ie] = ode15s(@(t,w) eqm3biii(t,w), tspan, w0, odeset('Events', @ZustandIII));

    nt = length(t);        % Accumulate output. Copied from ballode example
    tout = [tout; t(2:nt)];
    vout = [vout; v(2:nt,:)];
    wout = [wout; w(2:nt,:)];
    teout = [teout; te];
    veout = [veout; ve];
    weout = [weout; we];
    ieout = [ieout; ie];

    v0 = [y(nt,1) y(nt,2)];     % new IC and tspan. Copied from ballode
    w0 = [y(nt,1) y(nt,2)];
    tspan(1) = t(nt);  
end

    figure;
    plot(teout,veout(:,1),'ro')
    hold on
    plot(teout,weout(:,1),'ro')
    xlabel('time');
    ylabel('displacement');
    title('placeholder');
    legend('y_1','y_2');
    hold off
    box on
% functions------------------------------------------------------------
function dv = eqm2ai(t,v)
    dv = zeros(2,1);
    dv(1) = v(2);
    dv(2) = my*r*om^2*sin(om*t);
end
function dw = eqm3ai(t,w)
    dw = zeros(2,1);
    dw(1) = w(2);
    dw(2) = (FW*sin(om*t)*t-k3*w(1))/m3;
end  
function dv = eqm2bii(t,v)
    dv = zeros(2,1);
    dv(1) = v(2);
    dv(2) = my*r*om^2*sin(om*t)-(k2*(v(1)-w(1)-s))/(m1+m2);
end
function dw = eqm3bii(t,v)
    dw = zeros(2,1);
    dw(1) = w(2);
    dw(2) = (FW*sin(om*t)-k3*w(1)+k2*(v(1)-w(1)-s))/m3;
end   
function dv = eqm2biii(t,v)
    dv = zeros(2,1);
    dv(1) = v(2);
    dv(2) = my*r*om^2*sin(om*t)-(k2*(v(1)-w(1)+s))/(m1+m2);
end
function dw = eqm3biii(t,v)
    dw = zeros(2,1);
    dw(1) = w(2);
    dw(2) = (FW*sin(om*t)-k3*w(1)+k2*(v(1)-w(1)+s))/m3;
end 
% events --------------------------------------------------------------
function [value,isterminal,direction] = ZustandI(t,v,w)
value = all(v(1)-w(1)<=-s && v(1)-w(1)>=s);
isterminal = 1;                        
direction = 0;                        
end
function [value,isterminal,direction] = ZustandII(t,v,w)
value = all(v(1)-w(1)>-s);
isterminal = 1;                        
direction = 0;                        
end
function [value,isterminal,direction] = ZustandIII(t,v,w)
value = all(v(1)-w(1)<s);                         
isterminal = 1;                        
direction = 0;                       
end
end

[Torsten]

The first thing you'll have to change is to put the equations for v and w together in one function - thus solving for v(1), v(2) and w(1), w(2) simultaneously, not one after the other.

Come back with the modified code after this change.

Best wishes

Torsten.

[Lennart]

I don't quite understand how to do that. I guess you mean this part of my code:

function dv = eqm2ai(t,v)
    dv = zeros(2,1);
    dv(1) = v(2);
    dv(2) = my*r*om^2*sin(om*t);
end
function dw = eqm3ai(t,w)
    dw = zeros(2,1);
    dw(1) = w(2);
    dw(2) = (FW*sin(om*t)*t-k3*w(1))/m3;
end

I didn't find a way to properly set them up in one function without error messages. Could you ellaborate?

[Torsten]

function dy=eqmai(t,y)
    dy=zeros(4,1)
    dy(1)=y(2);
    dy(2)=my*r*om^2*sin(om*t);
    dy(3)=y(4);
    dy(4)=(FW*sin(om*t)*t-k3*y(3))/m3;
end

[Lennart]

Thank you. Seems much less bloated that way. Still calculates forever on my "starting equation" though. But here is my new code. The parameters stay the same from above.

    tspan = [0 10];
    y0 = [0 dv0 0 0];

    tout = tspan(1);    % copied from ballode example
    yout = y0.';
    teout = [];
    yeout = [];
    ieout = [];

   % main loop
   while tout(end) <= tspan(end)

   % starting equations, for first y(1) & y(3)
   while tout(end) == 0            
         [t,y] = ode45(@(t,y) eqmai(t,y), tspan, y0);
   end

    [t,y,te,ye,ie] = ode45(@(t,y) eqmai(t,y), tspan, y0, odeset('Events', @ZustandI));
    [t,y,te,ye,ie] = ode45(@(t,y) eqmbii(t,y), tspan, y0, odeset('Events', @ZustandII));
    [t,y,te,ye,ie] = ode45(@(t,y) eqmbiii(t,y), tspan, y0, odeset('Events', @ZustandIII));

    nt = length(t);        % Accumulate output. Copied from ballode example
    tout = [tout; t(2:nt)];
    yout = [yout; y(2:nt,:,2:nt,:)];       % unsure about this
    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, also unsure about this
    tspan(1) = t(nt);   
end

    figure;
    plot(teout,veout(:,1),'b--',teout,yeout(:,3),'r:')
    %plot(teout,veout(:,2),'b--',teout,yeout(:,4),'r:')
    hold on
    xlabel('time');
    ylabel('displacement');
    title('placeholder');
    legend('v_1','w_1');
    hold off
    box on %functions----------------------------------------------------
function dy=eqmai(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                             % v(2)
    dy(2)=my*r*om^2*sin(om*t);                              % Gl. 2a
    dy(3)=y(4);                                             % w(2)
    dy(4)=(FW*sin(om*t)*t-k3*y(3))/m3;                      % Gl. 3a
end
function dy=eqmbii(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                             % v(2)
    dy(2)=my*r*om^2*sin(om*t)-(k2*(v(1)-w(1)-s))/(m1+m2);   % Gl. 2b II
    dy(3)=y(4);                                             % w(2)
    dy(4)=(FW*sin(om*t)-k3*w(1)+k2*(v(1)-w(1)-s))/m3;       % Gl. 3b II
end
function dy=eqmbiii(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                             % v(2)
    dy(2)=my*r*om^2*sin(om*t)-(k2*(v(1)-w(1)+s))/(m1+m2);   % Gl. 2b III
    dy(3)=y(4);                                             % w(2)
    dy(4)=(FW*sin(om*t)-k3*w(1)+k2*(v(1)-w(1)+s))/m3;       % Gl. 3b III
end  
% events --------------------------------------------------------
function [value,isterminal,direction] = ZustandI(t,y)
value = all(y(1)-y(3)<=-s || y(1)-y(3)>=s);
isterminal = 1;             
direction = 0;               
end
function [value,isterminal,direction] = ZustandII(t,y)
value = all(y(1)-y(3)>-s);
isterminal = 1;                            
direction = 0;                          
end
function [value,isterminal,direction] = ZustandIII(t,y)
value = all(y(1)-y(3)<s);
isterminal = 1;              
direction = 0;              
end
end

[Torsten]

You don't reset tout(end) from 0 after solving your starting equation - so it will be solved over and over again.

Further I don't understand why you call eqmai, eqmbii and eqmbiii one after the other - I thought it depends on y(1)-y(3) which one is to be called next.

Best wishes

Torsten.

[Lennart]

Ok it seems to me that I can just delete my starting equation or just insert a break command inside the "while t==0 loop" without drawbacks (for now) but than I get this error: Undefined function 'sign' for input arguments of type 'logical'.

Error in odezero (line 46) indzc = find((sign(vL) ~= sign(vR)) & (direction .* (vR - vL) >= 0));

Error in ode45 (line 406) odezero(@ntrp45,eventFcn,eventArgs,valt,t,y,tnew,ynew,t0,h,f,idxNonNegative);

Error in mechanical_frequency_converter (line 38) [t,y,te,ye,ie] = ode45(@(t,y) eqmi(t,y), tspan, y0, odeset('Events', @ZustandI));

which seem to idicate that I can't have the value-statements from the eventfunctions in my current form. So I tried to rewrite them so that it detects zero-crossings in the desired direction

e.g.

function [value,isterminal,direction] = ZustandIa(t,y)

value = y(1)-y(3)+s;

isterminal = 1;

direction = -1;

end

for the first event. But then I get the next error:

Index in position 3 exceeds array bounds (must not exceed 1).

Error in Testskript (line 45) yout = [yout; y(2:nt,:,2:nt,:)]; % unsure about this

so I'm getting further but still no end is near.

Also I don't know what you mean with calling my 3 equations one after another. I thought inside my main loop, the desired equation will be called automatically only if the event-statement is true? Or in other words how could I rewrite my code to do exatly that?

thanks for your help

[Torsten]

yout = y0;
...
yout = [yout; y(2:nt,1:4)];

And calling "ode45" with "eqmbii" after "eqmai" means that you know that for your initial conditions you have -s<v(1)-w(1)<s and that you will switch to the condition v(1)-w(1)<=-s. Are you sure this will happen ?

[Lennart]

No I don't know that will happen. I know eqmi is the starting function but I don't know if eqmii or eqmiii will follow it. But I want my system to detect that automatically everytime another event condition is fulfilled. How can I do that?

Also I think something is fishy with my main loop because after changing my event values from not supported logical impressions (y(1)-y(3)>-s etc.) to zero crossings like described in my last answer will yield me the next error:

Error using odearguments (line 83) The last entry in tspan must be different from the first entry.

Error in ode45 (line 115) odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in Testskript (line 34) [t,y,te,ye,ie] = ode45(@(t,y) eqmi(t,y), tspan, y0, odeset('Events', @ZustandIa));

I also changed thee first eventfunction into two eventfunctions (Ia and Ib) for zero crossings in each direction respectively.

If I delete the while loop completely, the calculation will finish but I get an empty plot.

I could add the altered code if desired.

Thanks

Lennart

[Steven Lord]

I've skimmed the code you've posted here and in previous questions, and I'm not sure reading through the code will be sufficient to understand what you're trying to do. Can you post the text description and/or equations (no code) of the problem you're trying to solve? Starting from the original problem description we may be able to offer a more efficient way to solve that problem than calling ode45 multiple times.

[Torsten]

@Steven Lord:

The original code can serve as a problem description:

https://de.mathworks.com/matlabcentral/answers/402546-how-to-solve-differential-equations-of-motion-of-second-order-which-are-changing-depending-on-if-el

@Lennart Hartmann: You set tspan(1)=t(end) which makes the first and last entry of tspan identical.

You should only use one event function for all cases with

value = (y(1)-y(3)-s)*(y(1)-y(3)+s)

This covers zero crossings at s and -s.

[Lennart]

Hey @Torsten, I implemented everything you said and a bit more and it seems to work. Thank you for your help. Much appreciated. Now I only need to check if the end results are correct and if it properly changes back and forth to the three different functions but it looks kinda ok for the first time. I might get back to this this thread in the future but I hope it's not necessary.

I leave my new code in the next answer down below, maybe you find a new error I don't see otherwise it might help another person in the future with a similar problem.

[Lennart]

new working code:

function mfc

%% Parameters

    m1 = 8; m2 =  8; m3 =  100; k3 =  200; k2 = (k3*(m1+m2))/m3; my = m1/(m1+m2);
    s = 0.2; r = 0.2; om = 5.37; dv0s = 5; dv0 = dv0s-my*r*om; FW = 10;

%% calculation

    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', @ZustandI));
    [t,y,te,ye,ie] = ode45(@(t,y) eqmiii(t,y), [tstart tend], y0, odeset('Events', @ZustandI));

    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
% plot
    figure;                             % displacement
    box on
    hold on
    plot(teout,yeout(:,3),'r:')     %plot(teout,yeout(:,1),'b--')
    xlabel('time');
    ylabel('displacement');
    title('mfc');
    legend('v2','v3');

    figure;                             % velocity
    box on
    hold on
    plot(teout,yeout(:,2),'r:',teout,yeout(:,4),'b--')
    xlabel('time');
    ylabel('velocity');
    title('mfctime');
    legend('v2dot','v3dot');  
%%functions
function dy=eqmi(t,y)
    dy=zeros(4,1);
    dy(1)=y(2);                                             % v(2)
    dy(2)=my*r*om^2*sin(om*t);                              % Gl. 2a
    dy(3)=y(4);                                             % w(2)
    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);                                             % v(2)
    dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)-s))/(m1+m2);   % Gl. 2b II
    dy(3)=y(4);                                             % w(2)
    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);                                             % v(2)
    dy(2)=my*r*om^2*sin(om*t)-(k2*(y(1)-y(3)+s))/(m1+m2);   % Gl. 2b III
    dy(3)=y(4);                                             % w(2)
    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 = (y(1)-y(3)+s)*(y(1)-y(3)-s);
isterminal = 1;                                     % 
direction = 0;                                      %
end
end

[Torsten]

The results you collect in your solution vectors tout,yout,... all stem from the set of equations defined in eqmiii(t,y), regardless of the relationship between y1-y3 and s. So your solution can't be correct - sorry.

Best wishes

Torsten.