Solving system of ODEs with event functions

Question

Lennart 2018-6-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/405614-solving-system-of-odes-with-event-functions

评论： Torsten 2018-6-27

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

13 个评论
显示 11更早的评论隐藏 11更早的评论

Lennart 2018-6-14

编辑：Lennart 2018-6-14

在 MATLAB Online 中打开

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

Lennart 2018-6-18

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

Lennart 2018-6-19

编辑：Lennart 2018-6-19

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

Lennart 2018-6-27

在 MATLAB Online 中打开

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 2018-6-27

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.

请先登录，再进行评论。

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

Solving system of ODEs with event functions

13 个评论
显示 11更早的评论隐藏 11更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Solving system of ODEs with event functions

13 个评论 显示 11更早的评论隐藏 11更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

13 个评论
显示 11更早的评论隐藏 11更早的评论