How to use Eventfunction for this ODE set?

Question

Angshuman Podder 2021-3-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/776217-how-to-use-eventfunction-for-this-ode-set

编辑： Angshuman Podder 2021-4-17

在 MATLAB Online 中打开

The main variables are R and N.

I want to switch between the functions externally or internally.

I know Event function can help but cannot implement it for this case.

EDIT : (working code based on Jan's template)

clc
clear
close all
y0   = [0,0,0,0]; %random initial nos.
t0   = 0;                  %random no.
tEnd = 10;                 %random no.
opt = odeset('RelTol', 1e-3);
y = y0;
t = t0;
State = 1;
while t(end) < tEnd
    yprime = myFunction(t, y, State);
    
    if yprime(2) > yprime(3)
        State = 1;
        disp('state is changed')
    else
        State = 2;
        disp('state is changed')
    end
    
    
    fcn = @(t,y) myFunction(t, y, State);
    
    opt = odeset('Events', @(t, y) myEvents(t,y,State));
    [at, ay] = ode45(fcn, [t(end), tEnd], y(end, :),opt);
   
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));
end
%solution check
plot(t,y(:,4),'o') %plot results
hold on
z11 =  15*t.^2 - (40/3)*t.^3;
z22 =  50*t.^2 - 1000*t;
plot(t,z11,t, z22) %plot analytical results
figure(2) 
z1 =  30*t - 40*t.^2;
z2 =  100*t - 1000;
plot(t,z1,t, z2) %plot ode set %note the intersection point
function [yprime] = myFunction(t, y, State)
yprime = [5*t;  30*t - 40*t.^2 ; 100*t - 1000 ]; %random funcs.
if State == 1
    yprime(4) = yprime(2);
    disp('state is 1')
else
    yprime(4) = yprime(3);
    disp('state is 2')
end
end
function [val,isterm,dir] = myEvents(t,y,State)
yprime = myFunction(t, y, State);
val      =  yprime(2) - yprime(3);
isterm   =  1;
dir      =  0;
end

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

请先登录，再进行评论。

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

Answer 1

Jan 2021-3-18

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/776217-how-to-use-eventfunction-for-this-ode-set#answer_651482

在 MATLAB Online 中打开

y0   = [1,2,3,4];
t0   = 0;
tEnd = 42;
opt = odeset('Tol', 1e-8);
y = y0;
t = t0;
while t(end) < tEnd
    if -y(2) > y(3)
        State = 1;
    else
        State = 2;
    end
    opt = odeset('Events', @(t, y) myEvents(t, y, State));
    fcn = @(t,y) myFunction(t, y, StaTe);
    
    [at, ay] = ode45(fcn, [t(end), tEnd], y(end, :), opt);
   
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));
end
function myFunction(t, y, State)
if State == 1
    ... implement the first function
else
    ... implement the second function
end
end
function [value,isterminal,direction] = myEvents(t, y, State)
if State == 1
    ... implement the first function
else
    ... implement the second function
end
end

2 个评论
显示无隐藏无

Angshuman Podder 2021-3-20

编辑：Angshuman Podder 2021-3-20

在 MATLAB Online 中打开

@Jan, thanks so much for taking the time to write this up. I understand most of the code except for the Events function. I cannot understand what it is actually doing, so I commented it out. I feel it is telling the solver ode45 to recalculate the step size. I modified the code to test the algorithm is working or not.

Would love your feedback.

P.S. The differential eqns. are dummy eqns.

clc
clear
y0   = [0.1,0.1,0.1,0.21];
t0   = 0;
tEnd = 42;
opt = odeset('RelTol', 1e-3);
y = y0;
t = t0;
State = 1;
while t(end) < tEnd
    yprime = myFunction(t, y, State);
    
    if yprime(2) > yprime(3)
        State = 1;
    else
        State = 2;
    end
    %opt = odeset('Events', @(t, y) myEvents(t, y, State));
    fcn = @(t,y) myFunction(t, y, State);
    
    [at, ay] = ode45(fcn, [t(end), tEnd], y(end, :),opt);
   
    t = cat(1, t, at(2:end));
    y = cat(1, y, ay(2:end, :));
end
plot(t,y)
function [yprime] = myFunction(t, y, State)
yprime = [5*y(1); 4*y(4) + 77*t ; 30*y(4) - 7*t];
if State == 1
    yprime(4) = yprime(2);
    disp('state is 1')
else
    yprime(4) = yprime(3);
    disp('state is 2')
end
end
% 
% function [val,isterm,dir] = myEvents(fun, State)
% if State == 1
%     ... implement the first function
% else
%     ... implement the second function
% val      =  fun(2) - fun(3);
% isterm   =  1;
% dir      =  0;
% end
% end