Solve the following simultaneous equations:

2 views (last 30 days)
Mr.DDWW on 27 Apr 2022
Commented: Torsten on 7 May 2022
I am trying to make a code and plot and Draw the diagram of X and Y vs. t for 0<t<5
here is my code. I think there is somthoing worng
theta=0.10895;
YF=0.0667;
X=0;
alpha= 0.29;
beta= 0.68;
y1=450;
y2=11.25;
t=0; Y=0.0667; Z=0;
f1=@(t,y)[-X/theta+(1+alpha)*y1(1-X)*Y^2+beta*y1(1-X)*Z^2;(YF-Y/theta)+(1+alpha)*y1(1-X)*Y^2-y2*Y;-Z/theta+beta*y1(1-X)*Z^2+2*alpha*y1(1-X)*Y^2-((y2*Z/beta));];
[T,u]=ode45(f1,[100 120],[0 0 0.0667]);
plot(T,u(:,1),'-',T,u(:,3),'-.',T,u(:,3),'.');
Mr.DDWW on 28 Apr 2022
% I tried to run the code but it is not working. here is the code
clc;clf;clear all;close all;
syms X(t) Y(t) Z(t) YF alpha beta gamma1 gamma2 T theta
% The given data
YF=0.0667;
alpha= 0.29;
beta= 0.68;
gamma1=450;
gamma2=11.25;
theta=0.10895;
Ode1 = diff(X) == -X/theta + (1+alpha)*gamma1*(1-X)*Y^2 + beta*gamma1*(1-X)*Z^2;
Ode2 = diff(Y) == (YF-Y)/theta+(1-alpha)*gamma1*(1-X)*Y^2-gamma2*Y;
Ode3 = diff(Z) == -Z/theta+beta*gamma1*(1-X)*Z^2+2*alpha*gamma1*(1-X)*Y^2-(gamma2*Z)/beta;
Odes = [Ode1;Ode2; Ode3];
[VF, Subs] = OdeToVectorField ([Ode1, Ode2,Ode3]);
M = matlabFunction (VF, 'Vars',{'T','Y'});
% The given condtions
[t,y]=ode45(M, [0 5],[0.0667 0 0]);
figure
plot(t,y(:,1:2))
gride on
title('X and Y')
legend (string(subs))
[t,y]=ode45(M, [100 120],[0.0667 0 0]);
figure
plot(t,y(:,1:2))
gride on
title('X and Y')
legend (string(subs))
[t,y]=ode45(M, [100 120],[0.0667 0 0]);
figure
plot(t,y(:,1),t,y(:,2))
gride on
title('X vs Y')
legend (string(subs))

Jan on 27 Apr 2022
Edited: Jan on 27 Apr 2022
f1 = @(t,y) [-X/theta+(1+alpha)*y1(1-X)*Y^2+beta*y1(1-X)*Z^2; ...
... % ^ ^ * are missing
(YF-Y/theta)+(1+alpha)*y1(1-X)*Y^2-y2*Y; ...
... % ^^^^^^^^^^^^ ^ and a missing * again
-Z/theta+beta*y1(1-X)*Z^2+2*alpha*y1(1-X)*Y^2-((y2*Z/beta));];
... % ^ ^
y1(1-X) would be indexing in the vector y1. You mean y1 * (1 - X) at all 5 locations.
At the 2nd marked location the closing parenthesis is at the wrong position:
(YF - Y) / theta
Using some spaces would increase the readability and support the debugging:
f1 = @(t,y) ...
[-X / theta + (1+alpha) * y1 * (1-X) * Y^2 + beta * y1 * (1-X) * Z^2; ...
(YF-Y) / theta + (1+alpha) * y1 * (1-X) * Y^2 - y2 * Y; ...
-Z / theta + beta * y1 * (1-X) * Z^2 + 2 * alpha * y1 * (1-X) * Y^2 - y2 * Z / beta];

Torsten on 27 Apr 2022
Edited: Torsten on 28 Apr 2022
theta=0.10895;
YF=0.0667;
alpha= 0.29;
beta= 0.68;
gamma1=450;
gamma2=11.25;
X0=0; Y0=0.0667; Z0=0;
f=@(t,y)[-y(1)/theta+(1+alpha)*gamma1*(1-y(1))*y(2)^2+beta*gamma1*(1-y(1))*y(3)^2;...
(YF-y(2))/theta+(1-alpha)*gamma1*(1-y(1))*y(2)^2-gamma2*y(2);...
-y(3)/theta+beta*gamma1*(1-y(1))*y(3)^2+2*alpha*gamma1*(1-y(1))*y(2)^2-gamma2*y(3)/beta];
[T,Y]=ode45(f,[100 120],[X0,Y0,Z0]);
plot(T,Y(:,1),'-',T,Y(:,3),'-.',T,Y(:,3),'.');
Torsten on 7 May 2022
The MATLAB solvers expect the unknowns be defined as a vector (here: y), not as a collection of scalar variables (here: X, Y and Z). So one has to decide which of your solution variables (X, Y and Z) to placed at which position of this vector. I decided to take X = y(1), Y = y(2), Z = y(3).