ASKSAGE: Sage Q&A Forum - Individual question feedhttp://ask.sagemath.org/questions/Q&A Forum for SageenCopyright Sage, 2010. Some rights reserved under creative commons license.Sun, 03 Apr 2016 11:31:12 -0500Solving symbolic equation set gives wrong answers when compared to MATLAB codehttp://ask.sagemath.org/question/32957/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/059085c6-509b-45b7-b9aa-f5f07c7e0194/files/EMEC%20425%20Thermal%20Systems/HW3/Thermal%20Systems%20HW%203.sagews
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/P0)^(vC3+vD3-vA3-vB3)-K3];
end
x=fsolve(@(x) HW3_d(x),[0.3;0.3;0.3])
%% Combustion Reaction: CH8H18 + 37.5*(02 + 3.76*N2) -> 8*CO2 + 9*H2O + 141*N2 + 25*02
% Number of moles
nH2O = 9
nCO2 = 8 - 2*x(1)
nN2 = 141 - x(2) - x(3)
nO2 = 25 + x(1) - x(2) - 2*x(3)
%% 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;
%% Mole Fractions
yH2O=nH2O/ntot
yO2=nO2/ntot
yN2=nN2/ntot
yCO2=nCO2/ntot
yCO=nCO/ntot
yNO=nNO/ntot
yNO2=nNO2/ntot
ytot = yH2O + yO2 + yN2 + yCO2 + yCO + yNO + yNO2
endFri, 01 Apr 2016 19:23:06 -0500http://ask.sagemath.org/question/32957/solving-symbolic-equation-set-gives-wrong-answers-when-compared-to-matlab-code/Comment by schilly for <div class="snippet"><p>Here is the public access link to the sagemath worksheet on the <a href="http://cloud.sagemath.com">cloud.sagemath.com</a> website, which has embedded explanatory text as to what is happening:</p>
<blockquote>
<p><a href="https://cloud.sagemath.com/projects/059085c6-509b-45b7-b9aa-f5f07c7e0194/files/EMEC%20425%20Thermal%20Systems/HW3/Thermal%20Systems%20HW%203.sagews">https://cloud.sagemath.com/projects/0...</a></p>
</blockquote>
<p>also at</p>
<blockquote>
<p><a href="http://preview.tinyurl.com/SageHWhelp">http://preview.tinyurl.com/SageHWhelp</a></p>
</blockquote>
<p>Thank you in advance for the help! This is driving me crazy!</p>
<p>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:</p>
<pre><code>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)
</code></pre>
<p>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.</p>
<p>EDIT again again: I found the MATLAB code. I really need to reorganize my files.</p>
<pre><code>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 ...</code></pre><span class="expander"> <a>(more)</a></span></div>http://ask.sagemath.org/question/32957/solving-symbolic-equation-set-gives-wrong-answers-when-compared-to-matlab-code/?comment=32972#post-id-32972So, 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. Sun, 03 Apr 2016 04:51:53 -0500http://ask.sagemath.org/question/32957/solving-symbolic-equation-set-gives-wrong-answers-when-compared-to-matlab-code/?comment=32972#post-id-32972Answer by nbruin for <div class="snippet"><p>Here is the public access link to the sagemath worksheet on the <a href="http://cloud.sagemath.com">cloud.sagemath.com</a> website, which has embedded explanatory text as to what is happening:</p>
<blockquote>
<p><a href="https://cloud.sagemath.com/projects/059085c6-509b-45b7-b9aa-f5f07c7e0194/files/EMEC%20425%20Thermal%20Systems/HW3/Thermal%20Systems%20HW%203.sagews">https://cloud.sagemath.com/projects/0...</a></p>
</blockquote>
<p>also at</p>
<blockquote>
<p><a href="http://preview.tinyurl.com/SageHWhelp">http://preview.tinyurl.com/SageHWhelp</a></p>
</blockquote>
<p>Thank you in advance for the help! This is driving me crazy!</p>
<p>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:</p>
<pre><code>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)
</code></pre>
<p>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.</p>
<p>EDIT again again: I found the MATLAB code. I really need to reorganize my files.</p>
<pre><code>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 ...</code></pre><span class="expander"> <a>(more)</a></span></div> http://ask.sagemath.org/question/32957/solving-symbolic-equation-set-gives-wrong-answers-when-compared-to-matlab-code/?answer=32975#post-id-32975It 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.Sun, 03 Apr 2016 11:31:12 -0500http://ask.sagemath.org/question/32957/solving-symbolic-equation-set-gives-wrong-answers-when-compared-to-matlab-code/?answer=32975#post-id-32975