I want to do stop condition which demand dx=0 in coupled differential equations
2 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
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
ye
%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
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.)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!