A system of nonlinear equations with three variables (part II)

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. But i met some code problems, particularly for the solutions of x and y.

The description of nonlinear equations is desribed as below code:

% a-b_CNS model at the steady state_Cheng Mei_2022

clear all

% effective normal stress (Pa)
p.sigma=5e7;p.To=273.15+20;p.R=8.314;p.K=6e10;p.Lt=8e-4;p.lambda=0.0625;
p.phio=0.32;p.dsb=2e-6;p.dbulk=2e-5;p.qsb=0.4;p.phic=0.2;p.phio=0.02;
p.qbulk=0.7;p.H=0.577;p.n=1.7;p.m=3;p.An=0.0759;p.At=4.4*p.An;
p.Ea=213e3;p.omega=3.69e-7;p.muo=0.6;p.ao=0.006;p.xo = 0.1813;
p.V0=1e-6;p.V1=0.1e-6;p.gammao=0.02;p.M = 2;p.N = 1;


% T=linspace(273.15,600,1000);
% x=linspace(0.1,0.2,1000);
% y=linspace(0.1,0.2,1000);
syms x y T positive

% At reference (prestep) velocity of V0
fn =((p.phic-x)./(x-p.phio)).^(-p.M);
eqnx = p.Lt.*p.lambda./p.V0*p.An.*exp(-p.Ea./p.R./T).*(p.sigma).^p.n./(p.dsb).^p.m.*fn==p.H.*(p.qsb-2.*x).^p.N;
solx = solve(eqnx,x,'MaxDegree', 3);

% At a velocity (poststep) of V1
fnn=((p.phic-y)./(y-p.phio)).^(-p.M);
eqny= p.Lt.*p.lambda./p.V1*p.An.*exp(-p.Ea./p.R./T).*(p.sigma).^p.n./(p.dsb).^p.m.*fnn==p.H.*(p.qsb-2.*y).^p.N;
soly = solve(eqny,y,'MaxDegree', 3);

format long
T_array=linspace(273,600,100);
for i = 1:numel(T_array)
  sol_x = subs(solx,T,T_array(i));
  tf = isAlways(sol_x >= 0.1 & sol_x <= 0.2);
  x_array(i) = sol_x;
  tanso(i) = 0.577*(0.4-2*x_array(i));
  sol_y = subs(soly,T,T_array(i));
  tff =isAlways(sol_y >= 0.1 & sol_y <= 0.2);
  y_array(i) = sol_y;     
  tansi(i) = 0.577*(0.4-2*y_array(i));
end

musi=p.muo+p.ao*T_array./p.To.*log(p.V1./p.V0);
muso=p.muo;
bb=-tansi.*(1+musi.^2)./(1-musi.*tansi)./(1-musi.*tansi).*(1-(p.V1./p.V0).^(1/(p.M+p.N)))./log(p.V1./p.V0);

plot(T_array-273.15,bb,'k-',LineWidth=2);
%xlabel(''Temperature (°C)')

[John D'Errico]

Is this not the third or 4th time you have asked this question?

[Mei Cheng]

@John D'Errico haha, they are similar! but there are still problems about my questions. thanks!

[Torsten]

@Mei Cheng

sol_x and sol_y are vectors of length 3 (the roots of a polynomial of degree 3 in x resp. y).

Thus setting

x_array(i) = sol_x;
y_array(i) = sol_y;  

will throw an error because you try to assign a vector (sol_x, sol_y) to a scalar (x_array(i), y_array(i)).

That's why I repeatedly told you to examine sol_x and sol_y to see which of the three values makes physical sense and thus which of the three solution you should save to x_array(i) resp. y_array(i).

[Mei Cheng]

@Torsten thanks for your reply.

I understand what you said as below. I am sorry that i dont know how to examine sol_x and sol_y one by one in three values. Actually, there must be one solution is physcially reasonable. 0.1< x <0.2, 0.1<y<0.2. Do you want how to examine each solution? or how to constrain the range of x, y.

If possible, you can edit the below code directly. Thanks again!

T_array=linspace(273,600,100);
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);
  tf = isAlways(sol_x >= 0.1 & sol_x <= 0.2);
  x_array(i) = sol_x_real(end);
  tanso(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);
  tff =isAlways(sol_y >= 0.1 & sol_y <= 0.2);
  y_array(i) = sol_y_real(end);     
  tansi(i) = 0.577*(0.4-2*y_array(i));
end

[Torsten]

Look at what is printed to console with the below code after the solutions with small imaginary parts and those which don't satisfy your bound constraints have been deleted by the commands

sol_x_real = sol_x(abs(imag(sol_x))<1e-15 & real(sol_x) >= 0.1 & real(sol_x) <=0.2)
sol_y_real = sol_y(abs(imag(sol_y))<1e-15 & real(sol_y) >= 0.1 & real(sol_y) <= 0.2)

Maybe it will only be the empty array - maybe there will still be more than one solution left so that you must specify the index of the solution which should be saved in x_array(i) resp. y_array(i).

T_array=linspace(273,600,100);
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 & real(sol_x) >= 0.1 & real(sol_x) <=0.2)
    x_array(i) = sol_x_real(end);
    tanso(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 & real(sol_y) >= 0.1 & real(sol_y) <= 0.2)
    y_array(i) = sol_y_real(end);     
    tansi(i) = 0.577*(0.4-2*y_array(i));
end

[Mei Cheng]

@Torsten thank you so much! you are a very kind man! Though i cannot fix my problem, I really appreaciate your help. I guess i will try to use "fsolve" to find a better solution. thanks again!

[Torsten]

What do you expect from using "fsolve" ?

The three solutions for x and y you get from the above symbolic approach are the only ones your system has.

If they don't satisfy your constraints, then there are no solutions that satisfy your constraints.