# How to use Eventfunction for this ODE set?

3 Ansichten (letzte 30 Tage)
Angshuman Podder am 18 Mär. 2021
Bearbeitet: Angshuman Podder am 17 Apr. 2021
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 Kommentare-2 ältere Kommentare anzeigen-2 ältere Kommentare ausblenden

Melden Sie sich an, um zu kommentieren.

### Antworten (1)

Jan am 18 Mär. 2021
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 KommentareKeine anzeigenKeine ausblenden
Angshuman Podder am 20 Mär. 2021
Bearbeitet: Angshuman Podder am 20 Mär. 2021
@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.
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
Jan am 20 Mär. 2021
The event function does not influence the step size, but detects events. This is used to stop the integration and to change the state.
If e.g. the function to be integrated changes, when a component of the trajectory gets 0, this is detected in the event function.

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by