solve function solution issues: spurious solutions

I have the following code to solve for equilibrium concentrations of a chemical reaction. I have 7 unknowns and 7 equations. I try to use the solve function but the answer it provides is a VERY large string of numbers and 'z^6, z^5' etc terms. I am not sure what is wrong. Below is my code followed by the warning followed by the solution of 1 of the variables (all 7 are similar in size and format).

---------------------------------------------------------------------------------------------------------------------------------------------

clear all

%% Initial Conditions and Setup

R = 8.314;

T0 = 298;

P = 10;

Tfuel = 600;

Tair = 300;

x = 75;

m = x./(100-x);

del = 0;

Tguess = 2000;

gn = 65054;

gh = -135643;

gc = -396410 + 285948;

kph = exp(-gh/(R*Tguess));

kpn = exp(-gn/(R*Tguess));

kpc = exp(-gc/(R*Tguess));

a = m;

b = 4-m;

c = (a+b)/2;

d = (13.54*del)/(1-del);

%% Equations

% syms no2 nn2 nh2 nh2o nno nco nco2

o2 = []; n2 = []; h2 = []; h2o = []; no = []; co = []; co2 = [];

eqn = @(no2,nn2,nh2,nh2o,nno,nco,nco2) [a==nco2+nco, 2*b==2*nh2o+2*nh2, a+2*c==2*nco2+nco+nh2o+nno+1*no2, 2*d+7.52*c==nno+2*nn2,...

kph==((nh2o)/(b*sqrt(P)*sqrt(no2/(no2+nh2+nh2o+nno+nn2+nco+nco2)))),...

kpc==((nco2)/(a*sqrt(P)*sqrt(no2/(no2+nh2+nh2o+nno+nn2+nco+nco2)))),...

kpn==((nno)/(sqrt(nn2)*sqrt(no2)))];

S = solve(eqn);

--------------------------------------------------------------------------------------------------

Warning: Possibly spurious solutions.

> In symengine

In mupadengine/evalin (line 123)

In mupadengine/feval (line 182)

In solve (line 293)

In project2 (line 30)

---------------------------------------------------------------------------------------------------

S.nco2

ans =

(632902321880349*10^(1/2)*root(z^6 - (110741925602844662885831111130074544364278775808*10^(1/2)*z^5)/206643248124617443881634579139429536681031199200625 - (82252267417609436537988965206224929949030607884395320604015284576243*z^4)/41124093269650497958202548478968550837856015077606594946973632421875 + (232156115280521206189216435429837602901043707904*10^(1/2)*z^3)/344405413541029073136057631899049227801718665334375 + (5875455492900587090841672373729346777706344561312253391898071269376*z^2)/5874870467092928279743221211281221548265145011086656420996233203125 - (5708990770823839524233143877797980545530986496*10^(1/2)*z)/41328649624923488776326915827885907336206239840125 + 3138550867693340381917894711603833208051177722232017256448/65798549231440796733124077566349681340569624124170551915157811875, z, 1))/274877906944

(632902321880349*10^(1/2)*root(z^6 - (110741925602844662885831111130074544364278775808*10^(1/2)*z^5)/206643248124617443881634579139429536681031199200625 - (82252267417609436537988965206224929949030607884395320604015284576243*z^4)/41124093269650497958202548478968550837856015077606594946973632421875 + (232156115280521206189216435429837602901043707904*10^(1/2)*z^3)/344405413541029073136057631899049227801718665334375 + (5875455492900587090841672373729346777706344561312253391898071269376*z^2)/5874870467092928279743221211281221548265145011086656420996233203125 - (5708990770823839524233143877797980545530986496*10^(1/2)*z)/41328649624923488776326915827885907336206239840125 + 3138550867693340381917894711603833208051177722232017256448/65798549231440796733124077566349681340569624124170551915157811875, z, 2))/274877906944

(632902321880349*10^(1/2)*root(z^6 - (110741925602844662885831111130074544364278775808*10^(1/2)*z^5)/206643248124617443881634579139429536681031199200625 - (82252267417609436537988965206224929949030607884395320604015284576243*z^4)/41124093269650497958202548478968550837856015077606594946973632421875 + (232156115280521206189216435429837602901043707904*10^(1/2)*z^3)/344405413541029073136057631899049227801718665334375 + (5875455492900587090841672373729346777706344561312253391898071269376*z^2)/5874870467092928279743221211281221548265145011086656420996233203125 - (5708990770823839524233143877797980545530986496*10^(1/2)*z)/41328649624923488776326915827885907336206239840125 + 3138550867693340381917894711603833208051177722232017256448/65798549231440796733124077566349681340569624124170551915157811875, z, 3))/274877906944

(632902321880349*10^(1/2)*root(z^6 - (110741925602844662885831111130074544364278775808*10^(1/2)*z^5)/206643248124617443881634579139429536681031199200625 - (82252267417609436537988965206224929949030607884395320604015284576243*z^4)/41124093269650497958202548478968550837856015077606594946973632421875 + (232156115280521206189216435429837602901043707904*10^(1/2)*z^3)/344405413541029073136057631899049227801718665334375 + (5875455492900587090841672373729346777706344561312253391898071269376*z^2)/5874870467092928279743221211281221548265145011086656420996233203125 - (5708990770823839524233143877797980545530986496*10^(1/2)*z)/41328649624923488776326915827885907336206239840125 + 3138550867693340381917894711603833208051177722232017256448/65798549231440796733124077566349681340569624124170551915157811875, z, 4))/274877906944

(632902321880349*10^(1/2)*root(z^6 - (110741925602844662885831111130074544364278775808*10^(1/2)*z^5)/206643248124617443881634579139429536681031199200625 - (82252267417609436537988965206224929949030607884395320604015284576243*z^4)/41124093269650497958202548478968550837856015077606594946973632421875 + (232156115280521206189216435429837602901043707904*10^(1/2)*z^3)/344405413541029073136057631899049227801718665334375 + (5875455492900587090841672373729346777706344561312253391898071269376*z^2)/5874870467092928279743221211281221548265145011086656420996233203125 - (5708990770823839524233143877797980545530986496*10^(1/2)*z)/41328649624923488776326915827885907336206239840125 + 3138550867693340381917894711603833208051177722232017256448/65798549231440796733124077566349681340569624124170551915157811875, z, 5))/274877906944

(632902321880349*10^(1/2)*root(z^6 - (110741925602844662885831111130074544364278775808*10^(1/2)*z^5)/206643248124617443881634579139429536681031199200625 - (82252267417609436537988965206224929949030607884395320604015284576243*z^4)/41124093269650497958202548478968550837856015077606594946973632421875 + (232156115280521206189216435429837602901043707904*10^(1/2)*z^3)/344405413541029073136057631899049227801718665334375 + (5875455492900587090841672373729346777706344561312253391898071269376*z^2)/5874870467092928279743221211281221548265145011086656420996233203125 - (5708990770823839524233143877797980545530986496*10^(1/2)*z)/41328649624923488776326915827885907336206239840125 + 3138550867693340381917894711603833208051177722232017256448/65798549231440796733124077566349681340569624124170551915157811875, z, 6))/274877906944

Answer 1

[Walter Roberson]

you used solve(). solve() attempts to find infinitely precise solutions, formulae for the solution when possible . solve() will not generate numeric approximation in any situation in which it can figure out the formulae .

Your equation turn out to have formulae for the solutions . The formulae turn out to involve the roots of a polynomial of degree 6. polynomial of degree 5 and higher do not always have solutions that can be expressed as "algebraic numbers", which is to say not with addition subtraction multiplication division powers and integral roots such as square roots. polynomial of degree 4 and lower can always have their roots expressed as exact formulae with those operation but degree 5 and higher cannot always be.

But still the roots of polynomials designate definite computable values that can be reasoned about. MATLAB represents the roots of a polynomial as the form root(f(z), z) for some variable zz, and represents one particular one of the roots as root(f(z),z,N) for positive integer N. The order of numbering the roots is undocumented . In the particular case of your equation the solution involves all 6 roots of a polynomial of degree 6.

If you had happened to use 'returnconditions','true' in the equation you probably would get out aa more compact form

[Arjun Anand]

Okay so i should just use the conditions 'returnconditions','true' as part of solve(eqn1, eqn2..., eqn7, 'returnconditions','true') and this would give me a compact numerical answer? Is there another way to do this, maybe by using fsolve?

[Arjun Anand]

I tried using " 'returnconditions', true " but my answers were still very long and in terms of z. Is there any other function i could use to solve this polynomial in z?

[Walter Roberson]

It appears that 'returnconditions', true does not make any difference for this particular situation.

The results you are getting are the solutions. Remember, solve() attempts to find infinitely precise solutions . For example if a solution were sqrt(2) then it would say sqrt(2) not 1.414 .

Question: is it possible that the variables are intended to represent chemical concentrations? Because if so then they would have to be real and non-negative. You would attach assumptions to the variables and attempt to solve:

syms no2 nn2 nh2 nh2o nno nco nco2 positive
S = solve(eqn(no2,nn2,nh2, nh2o,nno,nco,nco2), 'returnconditions', true);
S2 = solve(S1.conditions, 'returnconditions', true);
if isempty(S2.conditions)
    error('No positive real solutions to the original equations')
end

It turns out those equations have no solutions restricted to the non-negative reals. You can get close -- down to just a single value being negative, approximately

nco = 1.4099, nco2 = 1.5901, nh2 = -1.4097, nh2o = 2.4097, nn2 = 7.519979679, nno = 0.4064e-4, no2 = 0.4e-3

If you want to see approximate solutions, then vpa() the results of the solve()

[Arjun Anand]

Yes, these are supposed to represent chemical concentrations so they must be non-negative real. I editied the equations a little bit, and tried to implement your portion of the code however it still gives me a warning:

Undefined variable "S1" or class "S1.conditions".

Error in project2 (line 31)

S2 = solve(S1.conditions, 'returnconditions', true);

Below is my code set up after the initial conditions. I'm not sure how to fix this

syms no2 nn2 nh2 nh2o nno nco nco2 positive
o2 = []; n2 = []; h2 = []; h2o = []; no = []; co = []; co2 = [];
    eqn = @(no2,nn2,nh2,nh2o,nno,nco,nco2) [a==nco2+nco, 2*b==2*nh2o+2*nh2, a+2*c==2*nco2+nco+nh2o+nno+1*no2, 2*d+7.52*c==nno+2*nn2,...
        kph==((nh2o)/(2*b*sqrt(P)*sqrt(2*c/(no2+nh2+nh2o+nno+nn2+nco+nco2)))),...
        kpc==((nco2)/(a*sqrt(P)*sqrt(2*c/(no2+nh2+nh2o+nno+nn2+nco+nco2)))),...
        kpn==((nno)/(sqrt(2*d+7.52*c)*sqrt(2*c)))];
    S = solve(eqn, 'returnconditions', true);
    S2 = solve(S1.conditions, 'returnconditions', true);

[Walter Roberson]

syms no2 nn2 nh2 nh2o nno nco nco2 positive
o2 = []; n2 = []; h2 = []; h2o = []; no = []; co = []; co2 = [];
    eqn = @(no2,nn2,nh2,nh2o,nno,nco,nco2) [a==nco2+nco, 2*b==2*nh2o+2*nh2, a+2*c==2*nco2+nco+nh2o+nno+1*no2, 2*d+7.52*c==nno+2*nn2,...
        kph==((nh2o)/(2*b*sqrt(P)*sqrt(2*c/(no2+nh2+nh2o+nno+nn2+nco+nco2)))),...
        kpc==((nco2)/(a*sqrt(P)*sqrt(2*c/(no2+nh2+nh2o+nno+nn2+nco+nco2)))),...
        kpn==((nno)/(sqrt(2*d+7.52*c)*sqrt(2*c)))];
    S1 = solve(eqn(no2, nn2, nh2, nh2o, nno, nco, nco2), 'returnconditions', true);
    if isempty(S1.conditons)
        error('No positive real solutions to the original equations')
    end
    S2 = solve(S1.conditions, 'returnconditions', true);
    if isempty(S2.conditions)
      error('No positive real solutions to the original equations')
    end

[Arjun Anand]

I'm still unable to get any values, negative or non-negative as S1 has empty fields for all the variables. How would i go about getting values like you did above?

[Walter Roberson]

There are two general solutions to your modified equations. In the general solutions, the only variables that come out real-valued are nn2 and nno, with all of the other variables coming out complex valued.

[Arjun Anand]

If you could show me how you were able to get these solutions that would be helpful, as i could attempt to modify those equations in case they are incorrect and try to find proper solutions.

[Walter Roberson]

I have not implemented the full cross-checking for the parameterized case below:

%% Initial Conditions and Setup
Q = @(v) sym(v, 'r');

R = Q(8.314);
T0 = Q(298);
P = Q(10);
Tfuel = Q(600);
Tair = Q(300);
x = Q(75);
m = x./(100-x);
del = Q(0);
Tguess = Q(2000);
gn = Q(65054);
gh = Q(-135643);
gc = Q(-396410) + Q(285948);
kph = exp(-gh/(R*Tguess));
kpn = exp(-gn/(R*Tguess));
kpc = exp(-gc/(R*Tguess));
a = m;
b = 4-m;
c = (a+b)/2;
d = (Q(13.54)*del)/(1-del);

%% Equations
syms no2 nn2 nh2 nh2o nno nco nco2
eqn = @(no2,nn2,nh2,nh2o,nno,nco,nco2) [a==nco2+nco, 2*b==2*nh2o+2*nh2, a+2*c==2*nco2+nco+nh2o+nno+1*no2, 2*d+7.52*c==nno+2*nn2,...
    kph==((nh2o)/(2*b*sqrt(P)*sqrt(2*c/(no2+nh2+nh2o+nno+nn2+nco+nco2)))),...
    kpc==((nco2)/(a*sqrt(P)*sqrt(2*c/(no2+nh2+nh2o+nno+nn2+nco+nco2)))),...
    kpn==((nno)/(sqrt(2*d+7.52*c)*sqrt(2*c)))];
EQN = eqn(no2, nn2, nh2, nh2o, nno, nco, nco2);
S1 = solve(EQN, 'returnconditions', true);
if isempty(S1.conditions)
    error('No solutions to the original equations :(')
end
if isempty(S1.parameters)
    disp('First equation did not need any parameters to solve!');
    disp(' ')
    S1aug = {S1.nco S1.nco2 S1.nh2 S1.nh2o S1.nn2 S1.nno S1.no2};
    disp('Have a look at the solutions below, which are in the order {nco nco2 nh2 nh2o nn2 nno no2}');
    disp(' ');
    celldisp(S1aug)
    disp(' ');
    disp('Numeric approximations:')
    disp(vpa(S1aug));
    disp(' ' );
    disp('Cross-checking solutions in no parameter case:');
    disp(' ')
    disp(vpa(subs(EQN, {nco nco2 nh2 nh2o nn2 nno no2}, {S1.nco S1.nco2 S1.nh2 S1.nh2o S1.nn2 S1.nno S1.no2}), 5));
    disp(' ')
    disp('Do the signs match on all of the columns? No 3489.5 == -3489.5 for example?')
else
    S2 = solve(S1.conditions, S1.parameters, 'returnconditions', true);
    if isempty(S2.conditions)
        error('First equation needed parameters to solve, but the parameters were unsolvable')
    else
        disp('First equation needed parameters to solve, and the parameters were solvable!');
    end
    S1parms = subs(S1.parameters, S2);
    S1aug = subs({S1.nco S1.nco2 S1.nh2 S1.nh2o S1.nn2 S1.nno S1.no2}, S1.parameters, S1parms);
    disp(' ');
    disp('Have a look at the solutions below, which are in the order {nco nco2 nh2 nh2o nn2 nno no2}')
    disp(' ');
    celldisp(S1aug)
    disp(' ');
    disp('Numeric approximations:');
    disp(' ');
    disp(vpa(S1aug));
end