I want to do stop condition which demand dx=0 in coupled differential equations

1 view (last 30 days)
I have two equations and I want to stop the integration when those two reach to dx=0.
the equations :
dx/dt=2*x-x^2+0.5*x*y
dy/dt=3*y-y^2+0.5*x*y
my code:
alpha=0.5;
beta=0.5;
r1=2;
r2=3;
s1=1;
s2=1;
t0 = 0;
tfinal = 10;
y0 = [1; 1];
[t,y] = ode23(@Two,[t0 tfinal],y0);
yp(1) = (r1+alpha*y(2))*y(1)-s1*(transpose(y(1))*y(1));
yp(2) = (r2 + beta*y(1))*y(2)-s2*(transpose(y(2))*y(2));
plot(t,y)
grid on
title('Two species')
xlabel('time')
ylabel('Population')
legend('X(t)','Y(t)')
function yp = Two(t,y)
yp = diag([2+0.5*y(2)-1*y(1),3+0.5*y(1)-1*y(2)])*y;
end

Accepted Answer

Torsten
Torsten on 23 Aug 2022
alpha=0.5;
beta=0.5;
r1=2;
r2=3;
s1=1;
s2=1;
t0 = 0;
tfinal = 10;
y0 = [1; 1];
AnonFun = @(t,y)diag([2+0.5*y(2)-1*y(1),3+0.5*y(1)-1*y(2)])*y;
Opt=odeset('Events',@(t,x)myEvent(t,x,AnonFun));
[t,y,te,ye,ie] = ode23(AnonFun,[t0 tfinal],y0,Opt);
te
te = 9.8627
ye
ye = 1×2
4.6658 5.3343
%yp(1) = (r1+alpha*y(2))*y(1)-s1*(transpose(y(1))*y(1));
%yp(2) = (r2 + beta*y(1))*y(2)-s2*(transpose(y(2))*y(2));
plot(t,y)
grid on
title('Two species')
xlabel('time')
ylabel('Population')
legend('X(t)','Y(t)')
function [value, isterminal, direction] = myEvent(t,x,AnonFun)
value = norm(AnonFun(t,x))-1.0e-2;
isterminal = 1; % Stop the integration
direction = -1;
end

More Answers (1)

Alan Stevens
Alan Stevens on 23 Aug 2022
You could do something like the following (doesn't actually stop the integration, but only plots the relevant bit):
alpha=0.5;
beta=0.5;
r1=2;
r2=3;
s1=1;
s2=1;
t0 = 0;
tfinal = 10;
y0 = [1; 1];
[t,y] = ode23(@Two,[t0 tfinal],y0);
yp(1) = (r1+alpha*y(2))*y(1)-s1*(transpose(y(1))*y(1));
yp(2) = (r2 + beta*y(1))*y(2)-s2*(transpose(y(2))*y(2));
plot(t,y)
grid on
title('Two species')
xlabel('time')
ylabel('Population')
legend('X(t)','Y(t)')
function yp = Two(~,y)
if y(1)>14/3
y = [NaN; NaN];
end
yp = diag([2+0.5*y(2)-1*y(1),3+0.5*y(1)-1*y(2)])*y;
end
(The value of 14/3 comes from setting your two ode's to zero and solving for the resulting value of x.)
  1 Comment
shir hartman
shir hartman on 23 Aug 2022
Thank you , but if I may - I want it to be more specific.
Is there a way to use odeset ? like that :
Opt=odeset("Events",@(t,x)myEvent(t,x,AnonFun));
[t,x]=ode45(AnonFun,[0,5],1,Opt);
function [value, isterminal, direction] = myEvent(t,x,AnonFun)
value = AnonFun(t,x)-1.0e-4;
isterminal = 1; % Stop the integration
direction = -1;
end
(this is from example of one function ...)

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!