I dont understand the error
显示 更早的评论
I would like it to stop calculating the ode when both functions meet the requirement that dx=0 , and I did:
alpha=0.5;
beta=0.5;
r1=2;
r2=3;
s1=1;
s2=1;
t0 = 0;
tfinal = 100;
y0 = [1;1];
AnonFun = @(t,y)diag([2+0.5*y(2)-1*y(1),3+0.5*y(1)-1*y(2)])*y;
if (alpha>0)&&(beta>0)
Opt=odeset('Events',@(t,y)myEvent1(t,y,AnonFun));
Other
Opt=odeset('Events',@(t,y)myEvent2(t,y,AnonFun));
end
[t,y,te,ye,ie] = ode23(AnonFun,[t0 tfinal],y0,Opt);
plot (t,y)
function [value, isterminal, direction] = myEvent1(t,y,AnonFun)
value = AnonFun(t,y)-1.0e-3;
isterminal = 1; % Stop the integration
direction = -1;
end
function [value, isterminal, direction] = myEvent2(t,y,AnonFun)
value=abs(AnonFun(t,y))-0.001;
isterminal = 1; % Stop the integration
direction = -1;
end
But when I change the vector y0 to [1;5] for example, I got this message:
Index exceeds the number of array elements. Index must not exceed 1.
Error in odezero (line 142) if any(isterminal(indzc))
Error in ode23 (line 335) odezero(@ntrp23,eventFcn,eventArgs,valt,t,y,tnew,ynew,t0,h,f,idxNonNegative);
Error in LogisticGrowthForTwoSpecies (line 17)
[t,y,te,ye,ie] = ode23(AnonFun,[t0 tfinal],y0,Opt);
采纳的回答
Edit: After understanding what you really want in your latest clarification. There is a simpler and intuitive way to code the program such that the simulation only run for approximately 3 or 4 times the Settling Time, 
.
. This allows 
or 
of the plot window to show the transient trajectories of the states from the initial values to the steady-state values. This method is only meaningful for systems that have stable equilibrium points.
or of the plot window to show the transient trajectories of the states from the initial values to the steady-state values. This method is only meaningful for systems that have stable equilibrium points.Analysis shows that your system has a stable equilibrium point at 
and 
and three other unstable equilibrium points at 
, 
, and 
.
and and three other unstable equilibrium points at , , and . Note: I changed the initial values for 
and 
because 
and 
, and to show you the difference between the settling time approach and the event function approach.
and because and , and to show you the difference between the settling time approach and the event function approach.alpha = 0.5;
beta = 0.5;
r1 = 2;
r2 = 3;
s1 = 1;
s2 = 1;
t0 = 0;
tfinal = 5; % Adjust this parameter roughly 4 times the Settling Time, Ts
y0 = [6 6];
AnonFun = @(t,y) diag([2 + 0.5*y(2) - 1*y(1), 3 + 0.5*y(1) - 1*y(2)])*y;
% AnonFun = @(t,y) [(2 + 0.5*y(2) - 1*y(1))*y(1);
% (3 + 0.5*y(1) - 1*y(2))*y(2)];
[t, y] = ode23(AnonFun, [t0 tfinal], y0);
plot(t, y), grid on
y1 = y(:,1);
y1e = 14/3;
idx = find(t > 1 & y1/y1e > 0.98 & y1/y1e < 1.02); % applying 2% criterion after 1 sec
Ts = t(idx(1))
14 个评论
From previous questions, the OP wants to stop integration if dy1/dt and dy2/dt are 0.
Thus stopping if dy1/dt or dy2/dt are 0 as in the above code is wrong.
Use
value = norm(AnonFun(t,y))-1e-3
instead of
value = AnonFun(t,y)-1.0e-3;
value = abs(AnonFun(t,y))-0.001;
shir hartman
2022-8-24
Read my comment.
Use
value = abs(AnonFun(t,y))-1.0e-3;
isterminal = [1;1];
direction = [-1;-1];
if you want to stop integration if dy1/dt = 0 or dy2/dt = 0.
Use
value = norm(AnonFun(t,y))-1.0e-3;
isterminal = 1;
direction = -1;
if you want to stop integration if dy1/dt = 0 and dy2/dt = 0.
shir hartman
2022-8-24
Torsten
2022-8-24
What is the code you are using ?
shir hartman
2022-8-25
Why don't you change the code if I gave you the correction two times already ?
And reread my comment for that you know why you get different results for the two event functions below.
alpha=0.5;
beta=0.5;
r1=2;
r2=3;
s1=1;
s2=1;
t0 = 0;
tfinal = 100;
y0 = [1;5];
AnonFun = @(t,y)diag([2+0.5*y(2)-1*y(1),3+0.5*y(1)-1*y(2)])*y;
if (alpha>0)&&(beta>0)
Opt=odeset('Events',@(t,y)myEvent1(t,y,AnonFun));
else
Opt=odeset('Events',@(t,y)myEvent2(t,y,AnonFun));
end
[t,y,te,ye,ie] = ode45(AnonFun,[t0 tfinal],y0,Opt);
plot (t,y)
function [value, isterminal, direction] = myEvent1(t,y,AnonFun)
value = norm(AnonFun(t,y))-1.0e-2;
isterminal = 1; % Stop the integration
direction = -1;
end
function [value, isterminal, direction] = myEvent2(t,y,AnonFun)
value=abs(AnonFun(t,y))-1.0e-2;
isterminal = [1;1]; % Stop the integration
direction = [-1;-1];
end
shir hartman
2022-8-25
Torsten
2022-8-25
Then relax the event from 1e-3 to 1e-2 (see above).
Sam Chak
2022-8-26
You probably need some mathematics to explain the behavior of dynamics (you are expert in your dynamics). The question is, do the states actually reach the stable equilibrium point 
in finite time or asymptotically converge to the point at 
?
in finite time or asymptotically converge to the point at ? For example, the trajectory of 
with 
is given by
with is given by 
.It won't reach zero in finite time, and thus you cannot get 
in finite time.
in finite time.If you can tell us why you want to stop ode45 for stable equilibrium point in finite time, perhaps we can help.
Note that it is possible to stop ode45 when state achieves a certain condition, say when 
reaches and remains within a given error band (
). The time at which the state attains this is called the Settling Time, 
. But it is not mathematically true that the system reaches 
in finite time.
reaches and remains within a given error band (). The time at which the state attains this is called the Settling Time, . But it is not mathematically true that the system reaches in finite time.
shir hartman
2022-8-26
Torsten
2022-8-26
What Sam Chak means is:
If the solution is e.g. f(t) = exp(-t), then it seems that the function doesn't change for t big enough. But y=0 is never reached. So you must define a threshold thr (e.g.1e-2) to stop when f(t) (or its derivative) is < thr.
I have updated my Answer to show you an alternative approach. But the following uses the Event function approach to force stop the ode45, so that you can see which one suits your needs.
alpha = 0.5;
beta = 0.5;
r1 = 2;
r2 = 3;
s1 = 1;
s2 = 1;
t0 = 0;
tfinal = 10;
y0 = [6 6];
AnonFun = @(t,y) diag([2 + 0.5*y(2) - 1*y(1), 3 + 0.5*y(1) - 1*y(2)])*y;
% AnonFun = @(t,y) [(2 + 0.5*y(2) - 1*y(1))*y(1);
% (3 + 0.5*y(1) - 1*y(2))*y(2)];
if (alpha > 0) && (beta > 0)
Opt = odeset('Events', @(t, y) myEvent1(t, y, AnonFun));
else
Opt = odeset('Events', @(t, y) myEvent2(t, y, AnonFun));
end
[t, y, te, ye, ie] = ode23(AnonFun, [t0 tfinal], y0, Opt);
plot(t, y), grid on
function [value, isterminal, direction] = myEvent1(t, y, AnonFun)
value = norm(AnonFun(t,y)) - 1e-2;
isterminal = 1; % Stop the integration
direction = -1;
end
function [value, isterminal, direction] = myEvent2(t, y, AnonFun)
value = norm(AnonFun(t,y)) - 1e-2;
isterminal = 1; % Stop the integration
direction = -1;
end
shir hartman
2022-8-29
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
