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

1 ビュー (過去 30 日間)
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

採用された回答

Torsten
Torsten 2022 年 8 月 23 日
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

その他の回答 (1 件)

Alan Stevens
Alan Stevens 2022 年 8 月 23 日
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 件のコメント
shir hartman
shir hartman 2022 年 8 月 23 日
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 ...)

サインインしてコメントする。

Community Treasure Hunt

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

Start Hunting!

Translated by