I want to do stop condition which demand dx=0 in coupled differential equations
2 ビュー (過去 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
0 件のコメント
採用された回答
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
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
その他の回答 (1 件)
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.)
参考
カテゴリ
Help Center および File Exchange で Ordinary Differential Equations についてさらに検索
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!