ODE45 - want to stop simulation when Steady State occurs

조회 수: 7 (최근 30일)
Marc Birkedal Jensen
Marc Birkedal Jensen 2016년 1월 7일
답변: Venkatachala Sarma 2016년 1월 12일
So with inspiration from the post "steady state criterion" I have tried to get my simulation to stop when steady state occurs, using an event with ODE45. However it doesn't stop when steady state occurs, it just moves on to the end time defined, which should change to the time of steady state when the script runs. I have two almost similar parts in the script doing almost the same, the second part works great if I fiddle in the numbers missing from the first part, but the first part keeps on getting a t value of 1000 - I know this value should be 30.8 with the current input. Can't seem to figure out why it doesn't work.
These are my scripts:
global yo j c d g b jm cm gm y y2 tf tf2
j= 0.0031;
yo =[999 1 0];
#First part of the script:
global yo y y2 tf tf2
to = 0;
odefun1 = @(t,y) ligeet(y);
efun1 = @(t,y) odeevent1(y);
options = odeset('Events',efun1);
sol = ode45(odefun1,[to 1000],yo,options);
t = (sol.x)';
Y = (sol.y)';
tf = t(end);
#second part:
hold on
odefun2 = @(t2,y2) ligeto(y2);
efun2 = @(t2,y2) odeevent2(y2);
options = odeset('Events',efun2);
sol2 = ode45(odefun2,[t(end) 10000],yt,options);
t2 = sol2.x;
t2 = t2';
Y2 = (sol2.y)';
h =plot(t2,Y2(:,1),'b',t2,Y2(:,2),'m',t2,Y2(:,3),'c',t2,Y2(:,4),'g',t2,Y2(:,1)+Y2(:,2)+Y2(:,3)+Y2(:,4),'k');
title('Endemisk SIR model - introduktion af muteret virus med mindsket helbredsrate')
legend(h([1 2 3 4 5]),'S', 'I', 'I*','R','N');
axis([0 t2(end) 0 max(max(y))+min(max(y))]);
tf2 = t2(end);
%%%%%%%%%%%%Ligeet.m%%%%%%%%%%%%%% function F = ligeet(y)
global j c d g b
F(2) = j*y(1)*y(2)-((d+c+g)*y(2));
F = F(:);
function F = ligeto(y)
global j c d g b jm cm gm
F(2) = j*y(1)*y(2)-((d+c+g)*y(2));
F(4)=(g*y(2))+ gm*y(3)-d*y(4);
F = F(:);
function [x,isterm,dir] = odeevent1(y)
dy = ligeet(y);
x = norm(dy) - 1e-5;
isterm = 1;
dir = 0;
function [x,isterm,dir] = odeevent2(y)
dy = ligeto(y);
x = norm(dy) - 1e-5;
isterm = 1;
dir = 0;

답변 (1개)

Venkatachala Sarma
Venkatachala Sarma 2016년 1월 12일
I initially got an error while running your code. I changed the 'y' to 'sol2.y' in the last but 1 and last but 3 lines. Then it started working. I am able to observe the plot only for 1000 and not for 10000 as should have been the case for the second part of the code.
The first plot gives a response for 4 graphs and the second plot gives a few red dots at the maximum value of the solution 2's output.
Are you getting any error? Because in the above code, y is not defined before using it. MATLAB will treat it as a null ([]). So the line in which s is defined, will give an error.
-Venkatachala Sarma


Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by