I want to do stop condition which demand dx=0 in coupled differential equations
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
shir hartman
am 23 Aug. 2022
Kommentiert: shir hartman
am 24 Aug. 2022
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 Kommentare
Akzeptierte Antwort
Torsten
am 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
Weitere Antworten (1)
Alan Stevens
am 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.)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!