采纳的回答

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 中查找有关 Programming 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by