采纳的回答

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2022-12-29
Ran in:

1 个投票

Ran in:
There are a few builtin fcn in matlab to solve such ODEs.
(1) Symbolic Solution -see DOC:
% Symbolic
syms x(t) y(t) z(t) p k a1 a2 b2 c d sigma beta
Eqn1 = diff(x(t), t) == p*x*(1-y/k)-a1*x*y;
Eqn2 = diff(y(t), t) == c*a1*x*y-d*y-a2*y*z/(y+a2);
Eqn3 = diff(z(t), t) == sigma*z^2-beta*z^2/(y+b2);
Sols = dsolve([Eqn1, Eqn2, Eqn3])
Warning: Unable to find symbolic solution.
Sols = [ empty sym ]
(2) Numerical Solution - see DOC:
%% Numerical solutions
ICs = [pi; pi/2; 2]; % Initial Conditions
[time, Sol]=ode45(@(t, xyz)FCN(t, xyz), [0, 10], ICs);
plot(time, Sol(:,1), 'r',time, Sol(:,2), 'g', time, Sol(:,3), 'b', 'linewidth', 2)
legend('x(t)', 'y(t)', 'z(t)'); grid on
title('Numerical Solutions')
xlabel('$t$', 'interpreter', 'latex')
ylabel('$x(t), \ y(t), \ z(t)$', 'interpreter', 'latex')
function dxyz = FCN(t, xyz)
p =.1;
k =.2;
a1 =.3;
a2 =.4;
b2 =.5;
c =.6;
d =.7;
sigma =.8;
beta =.9;
dxyz=[p*xyz(1)*(1-xyz(2)/k)-a1*xyz(1).*xyz(2);
c*a1*xyz(1).*xyz(2)-d*xyz(2)-a2*xyz(2).*xyz(3)./(xyz(2)+a2);
sigma*xyz(3).^2-(beta*xyz(3).^2)./(xyz(2)+b2)];
end

11 个评论
显示 9更早的评论 隐藏 9更早的评论

Sulaymon Eshkabilov
Sulaymon Eshkabilov 2022-12-30
移动：DGM 2023-3-3
Note that these two codes must be together in one file - Version 1 :
%% Numerical solutions
ICs = [pi; pi/2; 2]; % Initial Conditions
[time, Sol]=ode45(@(t, xyz)FCN(t, xyz), [0, 10], ICs);
plot(time, Sol(:,1), 'r',time, Sol(:,2), 'g', time, Sol(:,3), 'b', 'linewidth', 2)
legend('x(t)', 'y(t)', 'z(t)'); grid on
title('Numerical Solutions')
xlabel('$t$', 'interpreter', 'latex')
ylabel('$x(t), \ y(t), \ z(t)$', 'interpreter', 'latex')
function dxyz = FCN(t, xyz)
p =.1;
k =.2;
a1 =.3;
a2 =.4;
b2 =.5;
c =.6;
d =.7;
sigma =.8;
beta =.9;
dxyz=[p*xyz(1)*(1-xyz(2)/k)-a1*xyz(1).*xyz(2);
c*a1*xyz(1).*xyz(2)-d*xyz(2)-a2*xyz(2).*xyz(3)./(xyz(2)+a2);
sigma*xyz(3).^2-(beta*xyz(3).^2)./(xyz(2)+b2)];
end
Or it must be in this form in one m-file using function handle - Version 2:
p =.1;
k =.2;
a1 =.3;
a2 =.4;
b2 =.5;
c =.6;
d =.7;
sigma =.8;
beta =.9;
dxyz=@(t, xyz)([p*xyz(1)*(1-xyz(2)/k)-a1*xyz(1).*xyz(2);
c*a1*xyz(1).*xyz(2)-d*xyz(2)-a2*xyz(2).*xyz(3)./(xyz(2)+a2);
sigma*xyz(3).^2-(beta*xyz(3).^2)./(xyz(2)+b2)]);
ICs = [pi; pi/2; 2]; % Initial Conditions
[time, Sol]=ode45(dxyz, [0, 10], ICs);
plot(time, Sol(:,1), 'r',time, Sol(:,2), 'g', time, Sol(:,3), 'b', 'linewidth', 2)
legend('x(t)', 'y(t)', 'z(t)'); grid on
title('Numerical Solutions')
xlabel('$t$', 'interpreter', 'latex')
ylabel('$x(t), \ y(t), \ z(t)$', 'interpreter', 'latex')
Use on these two.
Torsten
Torsten 2022-12-31
移动：DGM 2023-3-3
Ran in:
Ran in:
r1 =.1;
k1 =.2;
a1 =.3;
c1 =.4;
r2 =.5;
k2 =.6;
a2 =.7;
c2 =.8;
dxy=@(t,xy)[r1*xy(1)*(1-xy(1)/k1)*(xy(1)-a1)-r1*c1*xy(1)*xy(2)/k1;...
r2*xy(2)*(1-xy(2)/k2)*(xy(2)-a2)-r2*c2*xy(1)*xy(2)/k2];
ICs = [pi; pi/2]; % Initial Conditions
[time, Sol] = ode45(dxy, [0, 10], ICs);
plot(time, Sol(:,1), 'r',time, Sol(:,2), 'g', 'linewidth', 2)
legend('x(t)', 'y(t)'); grid on
title('Numerical Solutions')
xlabel('$t$', 'interpreter', 'latex')
ylabel('$x(t), \ y(t)$', 'interpreter', 'latex')
Sulaymon Eshkabilov
Sulaymon Eshkabilov 2022-12-31
移动：DGM 2023-3-3
(1) k1 vs. K1 are two different variables.
(2) ERRORs: * missing, x, y are not correct variable names, indices () of xy are missing
dxy=@(t, xy)([r1*x*(1-x/k1)(x-a1)-(r1*c1*x*y/K1);
r2*y*(1-y/k2)(y-a2)-(r2*c2*x*y/K2);
(3) Corrected but not completely
dxy=@(t, xy)([r1*xy*(1-xy/k1)*(xy-a1)-(r1*c1*xy*xy/K1);
r2*xy*(1-xy/k2)*(xy-a2)-(r2*c2*xy*xy/K2);
(4) Now it is your turn to correct indices of xy() and then you will understand what is what for x, y:
dxy=@(t, xy)([r1*xy*(1-xy/k1)*(xy-a1)-(r1*c1*xy*xy/K1); To be Corrected
r2*xy*(1-xy/k2)*(xy-a2)-(r2*c2*xy*xy/K2); To be Corrected
Sulaymon Eshkabilov
Sulaymon Eshkabilov 2022-12-31
移动：DGM 2023-3-3
Most welcome! Happy New Year!
Torsten
Torsten 2022-12-31
移动：DGM 2023-3-3
@Ahmad
Thank you for your good wishes. A happy new year for you, too.
Torsten
Torsten 2023-1-5
移动：DGM 2023-3-3
Ran in:
Ran in:
Since there is no independent variable (time) left, what do you want to plot ? A single point ?
x = 0;
y = 0;
plot(x,y,'o')
Ahmad
Ahmad 2023-1-5
移动：DGM 2023-3-3
I want to graph this nonlinear system with the five steady states (o,o) (a1,0) (k1,o) (0,a1)(0,k1)
everyone alone so i will get 5 graphs
r1 =.1;
k1 =.2;
a1 =.3;
c1 =.4;
r2 =.5;
k2 =.6;
a2 =.7;
c2 =.8;
dxy=@(t,xy)[r1*xy(1)*(1-xy(1)/k1)*(xy(1)-a1)-r1*c1*xy(1)*xy(2)/k1;...
r2*xy(2)*(1-xy(2)/k2)*(xy(2)-a2)-r2*c2*xy(1)*xy(2)/k2];
ICs = [pi; pi/2]; % Initial Conditions
[time, Sol] = ode45(dxy, [0, 10], ICs);
plot(time, Sol(:,1), 'r',time, Sol(:,2), 'g', 'linewidth', 2)
legend('x(t)', 'y(t)'); grid on
title('Numerical Solutions')
xlabel('$t$', 'interpreter', 'latex')
ylabel('$x(t), \ y(t)$', 'interpreter', 'latex')
Ahmad
Ahmad 2023-1-5
移动：DGM 2023-3-3
same like here in this file
Torsten
Torsten 2023-1-5
移动：DGM 2023-3-3
Then call ode45 5 times (maybe in a loop) for the different parameter combinations, save the five results and plot them in figure(1),...,figure(5). What's the problem ?
Ahmad
Ahmad 2023-3-2
移动：DGM 2023-3-3
May I get the command for graph it every equation alone with t time please?

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by