# Solving symbolic equation set gives wrong answers when compared to MATLAB code

Here is the public access link to the sagemath worksheet on the cloud.sagemath.com website, which has embedded explanatory text as to what is happening:

https://cloud.sagemath.com/projects/0...

also at

http://preview.tinyurl.com/SageHWhelp

Thank you in advance for the help! This is driving me crazy!

EDIT: I guess I should add some more explanatory text here about the sage worksheet. The answers in that are all correct up until the very last cell that contains the following:

f1(x1,x2,x3)= (yC1^vC1*yD1^vD1) / (yA1^vA1)*(P/P0)^(vC1+vD1-vA1)-K1
f2(x1,x2,x3)= (yC2^vC2) / (yA2^vA2*yB2^vB2)*(P/P0)^(vC2-vB2-vA2)-K2
f3(x1,x2,x3)= (yC3^vC3) / (yA3^vA3*yB3^vB3)*(P/P0)^(vC3-vB3-vA3)-K3
solutions = solve([f1, f2,f3],x1,x2,x3,solution_dict=True)


EDIT again: Wow that is embarrassing, I put up the wrong MATLAB code for comparison. It also appears that I have lost the MATLAB code to compare with. I guess the question is still there, though: Is sage not capable of solving these equations? I can tell you that I compared each one of the final equations f1 through f3 and they are all correct going into the solver as far as values go. It's just that the solver spits back different values than I expect. I had no less than 4 others check the numbers going in, including my professor, although I will say that none of us are very familiar with sage or python.

EDIT again again: I found the MATLAB code. I really need to reorganize my files.

function HW3_d()
%% Embedded Function

function F = HW3_d(x)

%% Given Conditions
P0=100;     % Ambient pressure (kPa)
T0=298.15;  % Ambient temperature (K)
P=500;      % Products pressure (kPa)
T=1600;     % Products temperature (K)

%% Combustion Reaction: CH8H18 + 37.5*(02 + 3.76*N2) -> 8*CO2 + 9*H2O + 141*N2 + 25*02
nO2=25+x(1)-x(2)-2*x(3);
nN2=141-x(2)-x(3);
nH2O=9;
nCO2=8-2*x(1);
% constants are initial concentrations

%% Dissociation Reactions

% Reaction #1: 2CO2 -> 2CO + O2
% Reaction #2: N2 + 02 -> 2NO
% Reaction #3: N2 + 202 -> 2NO2

% Number of Moles
nCO = 2*x(1);
nNO = 2*x(2);
nNO2 = 2*x(3);
ntot = nCO2 + nH2O + nN2 + nO2 + nCO + nNO + nNO2;

%% Equilibrium constants at 1600 K from Table A.11
K1=exp(-21.656); K2=exp(-10.547); K3=exp(-20.126);

%% Mole Fractions
yA1=nCO2/ntot; yB1=0; yC1=nCO/ntot; yD1=nO2/ntot;
yA2=nN2/ntot; yB2=nO2/ntot; yC2=nNO/ntot; yD2=0;
yA3=nN2/ntot; yB3=nO2/ntot; yC3=nNO2/ntot; yD3=0;

%% Stoichiometric Coefficients
vA1=2; vB1=0; vC1=2; vD1=1;
vA2=1; vB2=1; vC2=2; vD2=0;
vA3=1; vB3=2; vC3=2; vD3=0;

%% Equilibrium Constants
F=[(yC1^vC1*yD1^vD1)/(yA1^vA1*yB1^vB1)*(P/P0)^(vC1+vD1-vA1-vB1)-K1;
(yC2^vC2*yD2^vD2)/(yA2^vA2*yB2^vB2)*(P/P0)^(vC2+vD2-vA2-vB2)-K2;
(yC3^vC3*yD3^vD3)/(yA3^vA3*yB3^vB3)*(P ...
edit retag close merge delete

So, I have no idea how to help, especially what exactly is wrong? I only see a lot of matlab code I can't make use of and a worksheet that shows a bunch of solutions.

( 2016-04-03 04:51:53 -0600 )edit

Sort by » oldest newest most voted

It seems to me sage isn't doing anything particularly wrong. However, I think you're using the wrong tool for the job: solve is for exact expressions (and has its limitations). When you're putting floats in, the underlying package that does the real work, maxima, will replace them with rational numbers that are close by and do exact arithmetic. That gets rid of rounding errors during the computations, but it does mean that maxima ends up solving a system of equations with slightly different coefficients than you had in mind. Unfortunately, maxima also find a solution that corresponds to the vanishing of the denominator, probably because it cleared denominators before trying to solve (it has negative coefficients, however). Your homework probably expects you to use a numerical method.

more