Find real solution of system of algebraic equations (solved in Maple, but can't do it in Matlab)

Dear Matlab community,

I have solved a system of algebraic equations in Maple and would now like to solve the exact same system in Matlab. I thought this would be quite easy, but turns out I am wrong. Matlab returns an imaginary solution to the system while Maple finds the real one.

To write the script in Matlab I copied the code straight from Maple and then replaced Maple specific commands with Matlab specific commands. The code is attached in a file and you will also find it below. Could anyone urgently help me with finding the real solution in Matlab? Thank you!

The solutions to the system should be

• Pct_Al2O3_i = 55.6042 ...
• Pct_Al_i = 0.1792 ...
• Pct_SiO2 = 0.04636 ...
• Pct_Si_i = 0.09427 ...

clc; clear;

% Initialisation parameters & variables
Pct_Al_i = sym('Pct_Al_i');
Pct_Si_i = sym('Pct_Si_i');
Pct_Al2O3_i = sym('Pct_Al2O3_i');
Pct_SiO2_i = sym('Pct_SiO2_i');

A_int=0.025^2*pi;

T_proc=1600;

Rho_m=7000*10^3; 
Rho_s=2850*10^3;

W_m=500; 
W_s=75;

m_Al=3*10^(-4);
m_Si=3*10^(-4);
m_SiO2=3*10^(-5);
m_Al2O3=3*10^(-5);

Pct_Al_b=0.3;
Pct_Si_b=0;
Pct_SiO2_b=5;
Pct_Al2O3_b=50;

Pct_Al_beq=0.132;
Pct_Si_beq=0.131;
Pct_SiO2_beq=3.13;
Pct_Al2O3_beq=52.12;

AW_Al=26.9815385;
AW_Si=28.085;
AW_O=15.999;
AW_Mg=24.305;
AW_Ca=40.078;
AW_Fe=55.845;

MW_SiO2=AW_Si+2*AW_O;
MW_Al2O3=2*AW_Al+3*AW_O;
MW_MgO=AW_Mg+AW_O;
MW_CaO=AW_Ca+AW_O;

R_cst=8.3144621;

% Mass balance equations
Mass_eq1=0==(Pct_Al_b-Pct_Al_i)+4*AW_Al*m_Si/(3*AW_Si*m_Al)*(Pct_Si_b-Pct_Si_i);

Mass_eq2=0==(Pct_Al_b-Pct_Al_i)-4*Rho_s*m_SiO2*AW_Al/(3*Rho_m*m_Al*MW_SiO2)*(Pct_SiO2_b-Pct_SiO2_i);

Mass_eq3=0==(Pct_Al_b-Pct_Al_i)+2*Rho_s*m_Al2O3*AW_Al/(Rho_m*m_Al*MW_Al2O3)*(Pct_Al2O3_b-Pct_Al2O3_i);

% Equilibrium equation
delta_G0=-720680+133*T_proc;

x_SiO2_i=(Pct_SiO2_i/(100*MW_SiO2))/(Pct_Al2O3_i/(100*MW_Al2O3) + Pct_SiO2_i/(100*MW_SiO2) + (100-Pct_SiO2_i-Pct_Al2O3_i)/(100*MW_CaO));

e_AlonSi=0.058;
e_AlonAl=0.045;
e_SionSi=0.11;
e_SionAl=0.0056;

Gamma_Al_Hry=1;
Gamma_Si_Hry=1;
Gamma_SiO2_Ra=0.01;

ln_h_Al=log(Gamma_Al_Hry)+log(10)*e_SionAl*Pct_Si_i/100+log(Pct_Al_i/100);
ln_h_Si=log(Gamma_Si_Hry)+log(10)*e_AlonSi*Pct_Al_i/100+log(Pct_Si_i/100);
ln_a_Al2O3=log(0.30);
ln_a_SiO2=log(Gamma_SiO2_Ra*x_SiO2_i);

K_cst=-delta_G0/(R_cst*T_proc);
Equil_eq=K_cst==3*ln_h_Si+2*ln_a_Al2O3-4*ln_h_Al-3*ln_a_SiO2;

% Solving
S = solve([Mass_eq1, Mass_eq2, Mass_eq3, Equil_eq], [Pct_Al_i, ...
    Pct_Si_i, Pct_SiO2_i, Pct_Al2O3_i]);

[Geoff Hayes]

Mattias - why not step through the code (using the debugger) and see at which line it deviates from the solution that Maple produces? i.e. if you run both pieces of code side-by-side then at some point the variables will start to differ.

[Walter Roberson]

http://www.mathworks.com/matlabcentral/answers/29922-why-your-question-is-not-urgent-or-an-emergency

Answer 1

[Philip Caplan]

In order to get a real solution to your system of equations, you need to tell MATLAB to assume your variables are real. You can do this by passing the 'real' option to "sym":

% Initialisation parameters & variables
Pct_Al_i = sym('Pct_Al_i','real');
Pct_Si_i = sym('Pct_Si_i','real');
Pct_Al2O3_i = sym('Pct_Al2O3_i','real');
Pct_SiO2_i = sym('Pct_SiO2_i','real');

I get the correct solution you mentioned when using this modification. For more information about assumptions on symbolic variables, please refer to

http://www.mathworks.com/help/symbolic/assumptions-for-symbolic-objects.html