Réponse acceptée

Sulaymon Eshkabilov
Sulaymon Eshkabilov le 29 Déc 2022

1 vote

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 commentaires

Sulaymon Eshkabilov
Sulaymon Eshkabilov le 29 Déc 2022
Welcome
Sulaymon Eshkabilov
Sulaymon Eshkabilov le 30 Déc 2022
Déplacé(e) : DGM le 3 Mar 2023
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 le 31 Déc 2022
Déplacé(e) : DGM le 3 Mar 2023
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 le 31 Déc 2022
Déplacé(e) : DGM le 3 Mar 2023
(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 le 31 Déc 2022
Déplacé(e) : DGM le 3 Mar 2023
Most welcome! Happy New Year!
Torsten
Torsten le 31 Déc 2022
Déplacé(e) : DGM le 3 Mar 2023
Thank you for your good wishes. A happy new year for you, too.
Torsten
Torsten le 5 Jan 2023
Déplacé(e) : DGM le 3 Mar 2023
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 le 5 Jan 2023
Déplacé(e) : DGM le 3 Mar 2023
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 le 5 Jan 2023
Déplacé(e) : DGM le 3 Mar 2023
same like here in this file
Torsten
Torsten le 5 Jan 2023
Déplacé(e) : DGM le 3 Mar 2023
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 le 2 Mar 2023
Déplacé(e) : DGM le 3 Mar 2023
May I get the command for graph it every equation alone with t time please?

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Numerical Integration and Differential Equations dans Centre d'aide et File Exchange

Tags

Community Treasure Hunt

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

Start Hunting!

Translated by