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);

Answer 1

[Sam Chak]

Hi @shir hartman

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.

Analysis shows that your system has a stable equilibrium point at 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.

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))

Ts = 1.2273

[Torsten]

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]

I tried now for yo=[1 5] and it still doesnt work :(

[Torsten]

@shir hartman

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]

@TorstenHi , I saw your comment and I tried but it didnt work ... for tfinal =100 it didnt stop .. my code is working but clearly he has some problem because it doesnt work for y0=[1;5]

Thank you anyway :) if you have another idea I would like to hear

[Torsten]

What is the code you are using ?

[shir hartman]

@Torsten

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

[Torsten]

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]

I'm sorry , I tried your correction but what I say is that it doesn't work.. The integration dont stop - I want it will stop when they both stable (you can see its happening on t~20 and not 100).

[Torsten]

Then relax the event from 1e-3 to 1e-2 (see above).

[Sam Chak]

Hi @shir hartman,

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 ?

For example, the trajectory of with is given by

.

It won't reach zero in finite time, and thus you cannot get 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.

[shir hartman]

@Torsten Thank you for all your help!

@Sam Chak Sorry I wasn't clear - I wanted to stop the integration because at some point there is no change in the functions so I wanted to reduce the calculation time. Thanks anyway I think the problem is now solved thanks to @Torsten :)

[Torsten]

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.

[Sam Chak]

Hi @shir hartman

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]

thank you so much! both of you !