ODE45 - want to stop simulation when Steady State occurs
3 次查看(过去 30 天)
显示 更早的评论
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:
Inputvalues:
%%%%%%%%%%%%%%%StartScript.m%%%%%%%%%%%%%%%
global yo j c d g b jm cm gm y y2 tf tf2
j= 0.0031;
c=1.27;
d=0.39;
g=0.8;
b=0.5;
jm=0.0042;
cm=0.9525;
gm=0.7600;
yo =[999 1 0];
%%%%%%%%%%%%%%virttest.m%%%%%%%%%%%%%%%
#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)';
plot(t,Y(:,1),'b',t,Y(:,2),'m',t,Y(:,3),'g',t,Y(:,1)+Y(:,2)+Y(:,3),'k');
tf = t(end);
#second part:
hold on
yt(1)=Y([end],[1]);
yt(2)=Y([end],[2]);
yt(3)=1;
yt(4)=Y([end],[3]);
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');
xlabel('År')
ylabel('Individer')
title('Endemisk SIR model - introduktion af muteret virus med mindsket helbredsrate')
legend(h([1 2 3 4 5]),'S', 'I', 'I*','R','N');
s=[0:70:max(y)];
plot(t(end),s,'.r');
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(1)=(b*sum(y))-j*y(1)*y(2)-d*y(1);
F(2) = j*y(1)*y(2)-((d+c+g)*y(2));
F(3)=g*y(2)-d*y(3);
F = F(:);
%%%%%%%%%%%%Ligeto.m%%%%%%%%%%%%%%
function F = ligeto(y)
global j c d g b jm cm gm
F(1)=(b*sum(y))-d*y(1)-y(1)*(j*y(2)+jm*y(3));
F(2) = j*y(1)*y(2)-((d+c+g)*y(2));
F(3)=jm*y(1)*y(3)-(d+cm+gm)*y(3);
F(4)=(g*y(2))+ gm*y(3)-d*y(4);
F = F(:);
%%%%%%%%%%odeevent1.m%%%%%%%%%%
function [x,isterm,dir] = odeevent1(y)
dy = ligeet(y);
x = norm(dy) - 1e-5;
isterm = 1;
dir = 0;
%%%%%%%%%%odeevent2.m%%%%%%%%%%
function [x,isterm,dir] = odeevent2(y)
dy = ligeto(y);
x = norm(dy) - 1e-5;
isterm = 1;
dir = 0;
0 个评论
回答(1 个)
Venkatachala Sarma
2016-1-12
Hi,
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
0 个评论
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!