# A system of nonlinear equations with three variables

3 views (last 30 days)
Mei Cheng on 29 Jan 2023
Commented: Walter Roberson on 29 Jan 2023
Hi,
I am solving a system of nonlinear equations with two variable x, y, and a changing parameter T. I want to solve this problem using loop.
The description of nonlinear equations is desribed as below code:
% A changing parameters 273<T<700;
T=linespace(273,700,1000);
% two variables 0.02<x<0.2, 0.02<y<0.2
x=linespace(0.02,0.2,1000);
y=linespace(0.02,0.2,1000);
% All equations is below
% SO equation
SO = 5.8e30.*exp(-2.562./T).*((0.2-x)./(x-0.02))^(-2);
% also S has another expression
SO = 0.577.*(0.4-2.*x);
%Mo
MO = 0.6;
% yy
yy = (MO+SO)./(1-MO.*SO);
% SS equations is silimar to SO
SS = 5.8e29.*exp(-2.562./T).*((0.2-y)./(y-0.02))^(-2);
% also SS has another expression
SS = 0.577.*(0.4-2.*y);
% M equation
MM = 0.6+4.7e-5.*T;
% YY
YY= (MM+SS)./(1-MM.*SS);
% Plot a figure
plot (T-273, -2.3.*(YY-yy));

Torsten on 29 Jan 2023
The equations for x and y (eqnx and eqny) are polynomials of degree 3 in x resp. y.
They might have several real solutions x and y.
I arbitrarily took the last one in the list MATLAB returns.
You might want to examine sol_x and sol_y in the below loop in more detail.
syms x y T
eqnx = 5.8e30*exp(-2.562/T)*((0.2-x)/(x-0.02))^(-2) == 0.577*(0.4-2*x);
solx = solve(eqnx,x,'MaxDegree',3);
eqny = 5.8e29*exp(-2.562/T)*((0.2-y)/(y-0.02))^(-2) == 0.577*(0.4-2*y);
soly = solve(eqny,y,'MaxDegree',3);
format long
T_array=linspace(273,700,1000);
for i = 1:numel(T_array)
sol_x = double(subs(solx,T,T_array(i)));
sol_x_real = sol_x(abs(imag(sol_x))<1e-15);
x_array(i) = sol_x_real(end);
SO1(i) = 0.577*(0.4-2*x_array(i));
sol_y = double(subs(soly,T,T_array(i)));
sol_y_real = sol_y(abs(imag(sol_y))<1e-15);
y_array(i) = sol_y_real(end);
SS1(i) = 0.577*(0.4-2*y_array(i));
end
%MO
MO = 0.6;
% yy
yy = (MO+SO1)./(1-MO*SO1);
% M equation
MM = 0.6+4.7e-5*T_array;
% YY
YY= (MM + SS1)./(1-MM.*SS1);
% Plot a figure
plot(T_array-273, -2.3*(YY-yy))
Warning: Imaginary parts of complex X and/or Y arguments ignored.
Walter Roberson on 29 Jan 2023
abs(imag(sol_x))<1e-15 is a test to be sure that sol_x is real-valued "to within round-off error" . That is, sometimes calculations that mathematically "should" be real-valued, involve intermediate calculations that are complex-valued, and sometimes due to round-off errors and the fact you are using finite precision, the imaginary parts do not exactly cancel out when mathematically they should, and you are left with a small imaginary part in practice.
abs(imag(sol_x))<1e-15 is not a test that attempts to isolate x and y to be within th range 0.1 to 0.2 . If you need that then you can test
sol_x >= 0.1 & sol_x <= 0.2

Sulaymon Eshkabilov on 29 Jan 2023
In your code, there are a few errors with matrix operations. It is more computationally to use vectorized approach instead of for .. loop operation.
% A changing parameters 273<T<700;
T=linspace(273,700,1000);
% two variables 0.02<x<0.2, 0.02<y<0.2
x=linspace(0.02,0.2,1000);
y=linspace(0.02,0.2,1000);
% SO equation (1)
SO1 = 5.8e30*exp(-2.562./T).*((0.2-x)./(x-0.02)).^(-2);
% Another expression (2): SO
SO2 = 0.577*(0.4-2.*x);
%MO
MO = 0.6;
% yy
yy = (MO+SO1)./(1-MO*SO1);
% SS equation (1)
SS1 = 5.8e29*exp(-2.562./T).*((0.2-y)./(y-0.02)).^(-2);
% Another expression (2): SS
SS2 = 0.577*(0.4-2*y);
% M equation
MM = 0.6+4.7e-5*T;
% YY
YY= (MM + SS1)./(1-MM.*SS1);
% Plot a figure
plot(T-273, -2.3.*(YY-yy));
Mei Cheng on 29 Jan 2023
@Sulaymon Eshkabilov thanks for your reply. in this code, SO1 is the same with SO2, SS1 is the same with SS2. It means the variables x, T, or y, T are correlated.
In your revised code, i just see the
yy = (MO+SO1)./(1-MO*SO1);
YY= (MM + SS1)./(1-MM.*SS1);
Here, do SO1 and SS1 mean the SO and SS? Or SO1 and SS1 is constrained by SO2 and SS2, respectively?
From the simulation results, it seems a bit strange.

### Categories

Find more on Systems of Nonlinear Equations in Help Center and File Exchange

R2021a

### Community Treasure Hunt

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

Start Hunting!

Translated by