A system of nonlinear equations with three variables

Hi,

I thank you for your kind help in advance!

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));

Answer 1

[Torsten]

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.

[Mei Cheng]

@Torsten thanks for your reply. your code looks good.

here, we need to notice that, 0.02<x<0.2, 0.02<y<0.2, it means x and y cannot be imaginary number, and they must be real numbers.

However, in the code, x and y can be imaganary number sometimes. Could you constrain the range of x and y? or avoid the generation of imaginary numbers? Thanks!

[Torsten]

Could you constrain the range of x and y? or avoid the generation of imaginary numbers? Thanks!

If you try to understand the code, you will notice that I did this by choosing the solutions for x and y for which abs(imag(x)) resp. abs(imag(y)) < 1e-15.

Still in some cases, there are multiple real-valued solutions for x and y (at most 3 for a polynomial of degree 3 in x resp. y). That's why I told you to examine sol_x resp. sol_y to extract the physical senseful ones for your application.

[Mei Cheng]

@Torsten thanks for your reply. Yes, I notice that you try to choose the solution for x and y.

In actual, x and y are in the range of about 0.1-0.2. And you define abs(imag(sol_x))<1e-15; I don't know how to change the range in code.

Also, only one solution has a physical meaning. The extraction of x and y is difficult for me.

Could you please help define the range of x and y, or extract the useful number (x and y) automatically?

Thanks again!

[Walter Roberson]

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

Answer 2

[Sulaymon Eshkabilov]

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]

@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.