http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456
MATLAB Central Newsreader  Newton Raphson method for chemical equilibrium system
Feed for thread: Newton Raphson method for chemical equilibrium system
enus
©19942014 by MathWorks, Inc.
webmaster@mathworks.com
MATLAB Central Newsreader
http://blogs.law.harvard.edu/tech/rss
60
MathWorks
http://www.mathworks.com/images/membrane_icon.gif

Tue, 03 Aug 2010 23:38:07 +0000
Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768234
Gavin
I am trying to solve a system of equations relating to a chemical equilibrium problem in the form or A*x = B, five equations and five unknowns. The Amatrix is 5x5 as is the Jacobian. I have no problem solving the current system; however, when I try to solve for a larger system I start running into problems. I was wondering if anyone could provide some suggestions on how to use this method for a larger system of equations, thus solving for more unknowns. My function for the Newton Raphson method is shown below. Please let me know how I could improve the code. It currently works well even with bad initial guesses for the 5x5 system, but fails when the system is larger than 5x5. When I run the problem for a 6x6 system, it converges but gives a warning that the matrix is close to singular or badly scaled. How can I fix this? I've read that incorporating a damping method with Newton <br>
Raphson helps, but I am unsure of how to apply this to my problem. <br>
<br>
function solution = Newton(Function,Jacobian,x_i,tol)<br>
x = x_i;<br>
error = 2*tol;<br>
while error > tol<br>
F = feval(Function,x);<br>
error1 = max(abs(F));<br>
error2 = max(abs(F));<br>
J = feval(Jacobian,x);<br>
dx = J\(F);<br>
i=1;<br>
while error2 >= error1~isreal(F)<br>
xnew = x+0.5^i*dx;<br>
F = feval(Function,xnew);<br>
error2 = max(abs(F));<br>
i = i+1;<br>
end<br>
x = xnew;<br>
F = feval(Function,x);<br>
error = max(abs(F));<br>
disp(['error = ' num2str(error)])<br>
end<br>
solution = x;

Tue, 03 Aug 2010 23:58:19 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768239
TideMan
On Aug 4, 11:38 am, "Gavin " <wigg...@gmail.com> wrote:<br>
> I am trying to solve a system of equations relating to a chemical equilibrium problem in the form or A*x = B, five equations and five unknowns. The Amatrix is 5x5 as is the Jacobian. I have no problem solving the current system; however, when I try to solve for a larger system I start running into problems. I was wondering if anyone could provide some suggestions on how to use this method for a larger system of equations, thus solving for more unknowns. My function for the Newton Raphson method is shown below. Please let me know how I could improve the code. It currently works well even with bad initial guesses for the 5x5 system, but fails when the system is larger than 5x5. When I run the problem for a 6x6 system, it converges but gives a warning that the matrix is close to singular or badly scaled. How can I fix this? I've read that incorporating a damping method with Newton<br>
> Raphson helps, but I am unsure of how to apply this to my problem.<br>
><br>
> function solution = Newton(Function,Jacobian,x_i,tol)<br>
> x = x_i;<br>
> error = 2*tol;<br>
> while error > tol<br>
> F = feval(Function,x);<br>
> error1 = max(abs(F));<br>
> error2 = max(abs(F));<br>
> J = feval(Jacobian,x);<br>
> dx = J\(F);<br>
> i=1;<br>
> while error2 >= error1~isreal(F)<br>
> xnew = x+0.5^i*dx;<br>
> F = feval(Function,xnew);<br>
> error2 = max(abs(F));<br>
> i = i+1;<br>
> end<br>
> x = xnew;<br>
> F = feval(Function,x);<br>
> error = max(abs(F));<br>
> disp(['error = ' num2str(error)])<br>
> end<br>
> solution = x;<br>
<br>
Just one comment on your code:<br>
Apparently, you are restricting the increment in x to avoid moving<br>
into the complex territory of F, yet you are using i as the counter<br>
and exponent. This is foolish because next thing you'll be doing<br>
something that assumes i is sqrt(1). Believe me, this happens, and<br>
it can be quite tricky to find the bug. Use something else for the<br>
counter, like m, for example.<br>
<br>
Other than that, your code looks good to me.<br>
So, if it fails when you add another equation, it points to an error<br>
in how you've assembled the Jacobian. Are you sure your algebra is<br>
correct? It's easy to make a mistake when doing partial<br>
differentiation of a complicated equation. All it takes is one stray<br>
negative sign to screw everything up...............

Wed, 04 Aug 2010 03:28:04 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768280
Gavin
TideMan <mulgor@gmail.com> wrote in message <ec2acf45eb154a32bc4d4199af5dbbab@n19g2000prf.googlegroups.com>...<br>
<br>
> Just one comment on your code:<br>
> Apparently, you are restricting the increment in x to avoid moving<br>
> into the complex territory of F, yet you are using i as the counter<br>
> and exponent. This is foolish because next thing you'll be doing<br>
> something that assumes i is sqrt(1). Believe me, this happens, and<br>
> it can be quite tricky to find the bug. Use something else for the<br>
> counter, like m, for example.<br>
> <br>
> Other than that, your code looks good to me.<br>
> So, if it fails when you add another equation, it points to an error<br>
> in how you've assembled the Jacobian. Are you sure your algebra is<br>
> correct? It's easy to make a mistake when doing partial<br>
> differentiation of a complicated equation. All it takes is one stray<br>
> negative sign to screw everything up...............<br>
<br>
TideMan,<br>
<br>
I changed the " i " in the code to " m " to avoid the assumption of i is sqrt(1). This does not fix my problem. I also checked my Jacobian and everything seems to be calculated correctly. My problem seems to arise when certain values in the Jacobian become very large or very small compared to the other values in the matrix. I can provide my entire code if you are interested.<br>
<br>
Gavin

Wed, 04 Aug 2010 03:46:39 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768284
TideMan
On Aug 4, 3:28 pm, "Gavin " <wigg...@gmail.com> wrote:<br>
> TideMan <mul...@gmail.com> wrote in message <ec2acf45eb154a32bc4d4199af5db...@n19g2000prf.googlegroups.com>...<br>
> > Just one comment on your code:<br>
> > Apparently, you are restricting the increment in x to avoid moving<br>
> > into the complex territory of F, yet you are using i as the counter<br>
> > and exponent. This is foolish because next thing you'll be doing<br>
> > something that assumes i is sqrt(1). Believe me, this happens, and<br>
> > it can be quite tricky to find the bug. Use something else for the<br>
> > counter, like m, for example.<br>
><br>
> > Other than that, your code looks good to me.<br>
> > So, if it fails when you add another equation, it points to an error<br>
> > in how you've assembled the Jacobian. Are you sure your algebra is<br>
> > correct? It's easy to make a mistake when doing partial<br>
> > differentiation of a complicated equation. All it takes is one stray<br>
> > negative sign to screw everything up...............<br>
><br>
> TideMan,<br>
><br>
> I changed the " i " in the code to " m " to avoid the assumption of i is sqrt(1). This does not fix my problem. I also checked my Jacobian and everything seems to be calculated correctly. My problem seems to arise when certain values in the Jacobian become very large or very small compared to the other values in the matrix. I can provide my entire code if you are interested.<br>
><br>
> Gavin<br>
<br>
No, I'm not interested in looking at your code.<br>
All I can tell you is that in 40 years since graduation, I've done a<br>
lot of work using Newton Raphson and every time I've had a problem it<br>
has been because I got the algebra wrong in assembling the Jacobian.<br>
First, you need to check your partial differentiation.<br>
Next, you need to check that you have correctly transcribed the<br>
partial differentials into Matlab code.

Wed, 04 Aug 2010 04:28:04 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768297
Roger Stafford
TideMan <mulgor@gmail.com> wrote in message <126edaf0e1e54760a7376d680d5f9dfb@m17g2000prl.googlegroups.com>...<br>
> No, I'm not interested in looking at your code.<br>
> All I can tell you is that in 40 years since graduation, I've done a<br>
> lot of work using Newton Raphson and every time I've had a problem it<br>
> has been because I got the algebra wrong in assembling the Jacobian.<br>
> First, you need to check your partial differentiation.<br>
> Next, you need to check that you have correctly transcribed the<br>
> partial differentials into Matlab code.<br>
         <br>
It is conceivable to me that the NewtonRaphson method for multiple dimensions could run into convergence problems with reallife situations. To take an example from geometry suppose you are seeking the point of tangency between a plane and a sphere. Suppose further that you wish to start the iteration with a point (x,y,z) near the tangency point but not necessarily constrained to either the plane or the sphere. (Admittedly this is not a wise approach for this problem.) The closer you get to the tangency point, the more nearly singular the Jacobian gets. Though I haven't tried this problem I suspect there would be real trouble finding the tangency point accurately using this approach.<br>
<br>
My training in chemical equilibria is too ancient to be certain, but I speculate that there may be some situations involving a nearly unstable point of equilibrium where an analogous difficulty with NewtonRaphson could arise. It would depend on the nature of the reaction rates involved.<br>
<br>
Roger Stafford

Wed, 04 Aug 2010 04:50:02 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768304
TideMan
On Aug 4, 4:28 pm, "Roger Stafford"<br>
<ellieandrogerxy...@mindspring.com.invalid> wrote:<br>
> TideMan <mul...@gmail.com> wrote in message <126edaf0e1e54760a7376d680d5f9...@m17g2000prl.googlegroups.com>...<br>
> > No, I'm not interested in looking at your code.<br>
> > All I can tell you is that in 40 years since graduation, I've done a<br>
> > lot of work using Newton Raphson and every time I've had a problem it<br>
> > has been because I got the algebra wrong in assembling the Jacobian.<br>
> > First, you need to check your partial differentiation.<br>
> > Next, you need to check that you have correctly transcribed the<br>
> > partial differentials into Matlab code.<br>
><br>
>          <br>
> It is conceivable to me that the NewtonRaphson method for multiple dimensions could run into convergence problems with reallife situations. To take an example from geometry suppose you are seeking the point of tangency between a plane and a sphere. Suppose further that you wish to start the iteration with a point (x,y,z) near the tangency point but not necessarily constrained to either the plane or the sphere. (Admittedly this is not a wise approach for this problem.) The closer you get to the tangency point, the more nearly singular the Jacobian gets. Though I haven't tried this problem I suspect there would be real trouble finding the tangency point accurately using this approach.<br>
><br>
> My training in chemical equilibria is too ancient to be certain, but I speculate that there may be some situations involving a nearly unstable point of equilibrium where an analogous difficulty with NewtonRaphson could arise. It would depend on the nature of the reaction rates involved.<br>
><br>
> Roger Stafford<br>
<br>
Well, in my experience with water waves (which admittedly are<br>
generally pretty smooth in time and space), if something as robust as<br>
Newton Raphson blows up it is pointing to a problem  the waves have<br>
broken or the numerical scheme has gone unstable  but usually the<br>
problem is not with Newton Raphson itself, which is what the OP is<br>
assuming.

Wed, 04 Aug 2010 07:58:06 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768334
Steve Amphlett
"Gavin " <wigging@gmail.com> wrote in message <i3a98v$lck$1@fred.mathworks.com>...<br>
> I am trying to solve a system of equations relating to a chemical equilibrium problem in the form or A*x = B, five equations and five unknowns. The Amatrix is 5x5 as is the Jacobian. I have no problem solving the current system; however, when I try to solve for a larger system I start running into problems. I was wondering if anyone could provide some suggestions on how to use this method for a larger system of equations, thus solving for more unknowns. My function for the Newton Raphson method is shown below. Please let me know how I could improve the code. It currently works well even with bad initial guesses for the 5x5 system, but fails when the system is larger than 5x5. When I run the problem for a 6x6 system, it converges but gives a warning that the matrix is close to singular or badly scaled. How can I fix this? I've read that incorporating a damping method with <br>
Newton <br>
> Raphson helps, but I am unsure of how to apply this to my problem. <br>
> <br>
> function solution = Newton(Function,Jacobian,x_i,tol)<br>
> x = x_i;<br>
> error = 2*tol;<br>
> while error > tol<br>
> F = feval(Function,x);<br>
> error1 = max(abs(F));<br>
> error2 = max(abs(F));<br>
> J = feval(Jacobian,x);<br>
> dx = J\(F);<br>
> i=1;<br>
> while error2 >= error1~isreal(F)<br>
> xnew = x+0.5^i*dx;<br>
> F = feval(Function,xnew);<br>
> error2 = max(abs(F));<br>
> i = i+1;<br>
> end<br>
> x = xnew;<br>
> F = feval(Function,x);<br>
> error = max(abs(F));<br>
> disp(['error = ' num2str(error)])<br>
> end<br>
> solution = x;<br>
<br>
Gavin,<br>
<br>
I assume you are still looking to solve your thermochemical equilibrium equations. If your technical German is up to it, you may find this link interesting:<br>
<br>
<a href="http://elib.unistuttgart.de/opus/volltexte/2006/2725/pdf/Diss_Grill.pdf">http://elib.unistuttgart.de/opus/volltexte/2006/2725/pdf/Diss_Grill.pdf</a><br>
<br>
 Steve

Wed, 04 Aug 2010 14:00:11 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768406
Gavin
<br>
Thanks Steve, but I can't read German. I'm trying to solve the following system of equations:<br>
where: CH4 = N(1), H2O = N(2), CO = N(3), H2 = N(4), CO2 = N(5), O2 = N(6)<br>
<br>
CH4+CO+CO2C = 0<br>
4*CH4+2*H2O+2*H2H = 0<br>
H2O+CO+2*CO2+2*O2O = 0<br>
LnKp1+2*log(P/sum(N))+log(N(3))+3*log(N(4))log(N(1))log(N(2)) = 0<br>
LnKp2+log(N(5))+log(N(4))log(N(3))log(N(2)) = 0<br>
LnKp30.5*log(P/sum(N))+log(N(2))log(N(4))0.5*log(N(6)) = 0 <br>
<br>
Known inputs are LnKp, C, H, O. Where LnKp is the log of the equilibrium constant for that particular reaction at the specified temperature.

Wed, 04 Aug 2010 14:55:44 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768431
Torsten Hennig
> <br>
> Thanks Steve, but I can't read German. I'm trying to<br>
> solve the following system of equations:<br>
> where: CH4 = N(1), H2O = N(2), CO = N(3), H2 = N(4),<br>
> CO2 = N(5), O2 = N(6)<br>
> <br>
> CH4+CO+CO2C = 0<br>
> 4*CH4+2*H2O+2*H2H = 0<br>
> H2O+CO+2*CO2+2*O2O = 0<br>
> LnKp1+2*log(P/sum(N))+log(N(3))+3*log(N(4))log(N(1))<br>
> log(N(2)) = 0<br>
> LnKp2+log(N(5))+log(N(4))log(N(3))log(N(2)) = 0<br>
> LnKp30.5*log(P/sum(N))+log(N(2))log(N(4))0.5*log(N<br>
> (6)) = 0 <br>
> <br>
> Known inputs are LnKp, C, H, O. Where LnKp is the<br>
> log of the equilibrium constant for that particular<br>
> reaction at the specified temperature.<br>
<br>
The first thing I'd do is to substitute N(i) = N~(i)^2<br>
for i=1,...,6 to avoid negative arguments within<br>
the logarithm.<br>
<br>
Best wishes<br>
Torsten.

Wed, 04 Aug 2010 15:05:26 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768438
Maxx Chatsko
"Gavin " <wigging@gmail.com><br>
I'm pretty handy with pchem and differentials. Would you mind sending me your code? <br>
Maxx

Wed, 04 Aug 2010 15:19:05 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768441
Maxx Chatsko
"Maxx Chatsko" <chatskom@chemimage.com>

Wed, 04 Aug 2010 23:40:23 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768641
Gavin
"Maxx Chatsko" <chatskom@chemimage.com> wrote in message <i3c0d8$1o4$1@fred.mathworks.com>...<br>
> "Maxx Chatsko" <chatskom@chemimage.com><br>
<br>
Maxx, <br>
I emailed the code to you. Let me know what you think. Thanks.

Thu, 05 Aug 2010 03:03:05 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768676
Gavin
> The first thing I'd do is to substitute N(i) = N~(i)^2<br>
> for i=1,...,6 to avoid negative arguments within<br>
> the logarithm.<br>
> <br>
> Best wishes<br>
> Torsten.<br>
<br>
Torsten,<br>
<br>
Substituting N(i) = N~(i)^2 for i=1...6 does not solve my problem. I still receive warnings for ill conditioned matrix. It also increases the number of iterations needed to solve the system.<br>
<br>
Gavin

Thu, 05 Aug 2010 06:26:14 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768703
Torsten Hennig
> > The first thing I'd do is to substitute N(i) =<br>
> N~(i)^2<br>
> > for i=1,...,6 to avoid negative arguments within<br>
> > the logarithm.<br>
> > <br>
> > Best wishes<br>
> > Torsten.<br>
> <br>
> Torsten,<br>
> <br>
> Substituting N(i) = N~(i)^2 for i=1...6 does not<br>
> solve my problem. I still receive warnings for ill<br>
> conditioned matrix. It also increases the number of<br>
> iterations needed to solve the system.<br>
> <br>
> Gavin<br>
<br>
Did you adjust the Jacobian matrix according to the substitution ?<br>
<br>
Best wishes<br>
Torsten.

Thu, 05 Aug 2010 07:39:22 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#768714
Darren Rowland
Torsten makes a good point about the need for updating the Jacobian.<br>
I would also recommend writing a simple routine which calculates the Jacobian numerically, then compare your analytical version with the numerical.<br>
<br>
Another common approach when working with chemical equilibria is to work only with the logarithms of the concentrations, i.e. let M = log(N) then run the NR scheme on M. You will need to update your Jacobian for this but the changes are only minor.<br>
<br>
Hth<br>
Darren<br>
<br>
@ Maxx, <br>
Trying to access www.chemimage.com I got the following error<br>
<br>
Content Encoding Error<br>
The page you are trying to view cannot be shown because it uses an invalid or unsupported form of compression.<br>
<br>
This is using Firefox 3.6.8

Fri, 06 Aug 2010 04:06:04 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#769047
Gavin
> > Gavin<br>
> <br>
> Did you adjust the Jacobian matrix according to the substitution ?<br>
> <br>
> Best wishes<br>
> Torsten.<br>
<br>
Torsten,<br>
I have now adjusted the jacobian accordingly. However, I'm having issues with making the rest of the code consistent with changing to the N(i)^2 terms. I've posted links to my files below. Let me know how you would fix this. I am running the model with inputs of nC=2, nH=2, nO=4, T=500..1500, P=1. It runs fine at temperatures greater than 1200K but not so well below 1200K. The answers at 1300K should be yCH4=4.36e8, yH2O=0.24882, yCO=0.24882, yH2=0.08451, yCO2=0.41785.<br>
Gavin<br>
<br>
I've posted my original mfile here... https://files.me.com/wigging/fh1src<br>
and my modified (i.e. N(i)^2 ) file here... https://files.me.com/wigging/47n7xa

Fri, 06 Aug 2010 09:12:44 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#769086
Torsten Hennig
> > > Gavin<br>
> > <br>
> > Did you adjust the Jacobian matrix according to the<br>
> substitution ?<br>
> > <br>
> > Best wishes<br>
> > Torsten.<br>
> <br>
> Torsten,<br>
> I have now adjusted the jacobian accordingly.<br>
> However, I'm having issues with making the rest of<br>
> f the code consistent with changing to the N(i)^2<br>
> terms. I've posted links to my files below. Let me<br>
> know how you would fix this. I am running the model<br>
> with inputs of nC=2, nH=2, nO=4, T=500..1500, P=1.<br>
> It runs fine at temperatures greater than 1200K but<br>
> t not so well below 1200K. The answers at 1300K<br>
> should be yCH4=4.36e8, yH2O=0.24882, yCO=0.24882,<br>
> yH2=0.08451, yCO2=0.41785.<br>
> Gavin<br>
> <br>
> I've posted my original mfile here...<br>
> https://files.me.com/wigging/fh1src<br>
> and my modified (i.e. N(i)^2 ) file here...<br>
> https://files.me.com/wigging/47n7xa<br>
<br>
In every expression in the function F, you will have <br>
to replace N(i) by N(i)^2.<br>
This gives new functions F1,...,F6 which you will<br>
have to partially differentiate with respect to <br>
the N(i) to put the result in function Jac.<br>
<br>
Best wishes<br>
Torsten.

Fri, 06 Aug 2010 23:42:05 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#769364
Gavin
> In every expression in the function F, you will have <br>
> to replace N(i) by N(i)^2.<br>
> This gives new functions F1,...,F6 which you will<br>
> have to partially differentiate with respect to <br>
> the N(i) to put the result in function Jac.<br>
> <br>
> Best wishes<br>
> Torsten.<br>
<br>
Torsten,<br>
I've made the N(i)^2 changes to the F function but now I'm getting imaginary numbers. Something is definitely not right. I've updated the jacobian accordingly too. I've posted the latest file here:<br>
https://files.me.com/wigging/6f6j7j<br>
<br>
Gavin

Sat, 07 Aug 2010 10:15:30 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#769412
Torsten Hennig
I did not check your equations in detail, just<br>
found two errors that I marked.<br>
<br>
Best wishes<br>
Torsten.<br>
<br>
<br>
%%%%%%%%%%%%%% MAIN FUNCTION using Newton Raphson method %%%%%%%%%%%%%%%%<br>
function gasf_Newton_6sJ<br>
global T<br>
global E<br>
global Stoichiometry<br>
global P<br>
<br>
% input the flows of elements entering the system<br>
nC = input('enter > nC, mol = ');<br>
nH = input('enter > nH, mol = ');<br>
nO = input('enter > nO, mol = ');<br>
E = [nC; nH; nO];<br>
% set the temperature (K)<br>
T = input('enter > T, K = ');<br>
% set the pressure (Po = 1 bar)<br>
P = input('enter > P, bar = ');<br>
%set the stiociometry matrix<br>
Stoichiometry = [1 0 1 0 1 0;...<br>
4 2 0 2 0 0;...<br>
0 1 1 0 2 2];<br>
<br>
% solve the system of equations<br>
[N] = Newton(@f,@jac,[1 1 1 1 1 1]',1e6);<br>
<br>
% outputs of system, moles and total moles<br>
CH4 = sqrt(N(1));<br>
H2O = sqrt(N(2));<br>
CO = sqrt(N(3));<br>
H2 = sqrt(N(4));<br>
CO2 = sqrt(N(5));<br>
O2 = sqrt(N(6));<br>
nt = sqrt(sum(N));<br>
<br>
>> nt must be the sum of the sqrt, not the <br>
>> sqrt of the sum.<br>
<br>
% mole fractions<br>
yCH4 = CH4/nt<br>
yH2O = H2O/nt<br>
yCO = CO/nt<br>
yH2 = H2/nt<br>
yCO2 = CO2/nt<br>
yO2 = O2/nt<br>
ntt = yCH4 + yH2O + yCO + yH2 + yCO2 + yO2<br>
<br>
%%%%%%%%%%%%%%%%%%%%%%%%%% Jacobian Matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
function [Jac] = jac(N)<br>
<br>
Jac(1:3,:) = [2*N(1), 0, 2*N(3), 0, 2*N(5), 0 ;...<br>
8*N(1), 4*N(2), 0, 4*N(4), 0, 0 ;...<br>
0, 2*N(2), 2*N(3), 0, 4*N(5), 4*N(6)]; <br>
<br>
% CH4 + H2O > CO + 3*H2<br>
% lnKp1 = 2*ln(P/nt) + ln(nCO) + 3*ln(nH2)  ln(nCH4)  ln(nH2O)<br>
Jac(4,1) = (4*N(1))/sum(N)2/N(1);<br>
Jac(4,2) = (4*N(2))/sum(N)2/N(2);<br>
Jac(4,3) = (4*N(3))/sum(N)+2/N(3);<br>
Jac(4,4) = (4*N(4))/sum(N)+6/N(4);<br>
Jac(4,5) = (4*N(5))/sum(N);<br>
Jac(4,6) = (4*N(6))/sum(N);<br>
% CO + H2O > CO2 + H2<br>
% lnKp2 = ln(nCO2) + ln(nH2)  ln(nCO)  ln(nH2O)<br>
Jac(5,1) = 0;<br>
Jac(5,2) = 2/N(2);<br>
Jac(5,3) = 2/N(3);<br>
Jac(5,4) = 2/N(4);<br>
Jac(5,5) = 2/N(5);<br>
Jac(5,6) = 0;<br>
% H2 + 0.5*O2 > H2O<br>
% lnKp4 = 0.5*ln(P/nt) + ln(nH2O)  ln(nH2)  0.5*ln(nO2)<br>
Jac(6,1) = N(1)/sum(N);<br>
Jac(6,2) = N(2)/sum(N)+2/N(2);<br>
Jac(6,3) = N(3)/sum(N);<br>
Jac(6,4) = N(4)/sum(N)2/N(4);<br>
Jac(6,5) = N(5)/sum(N);<br>
Jac(6,6) = N(6)/sum(N)1/N(6);<br>
<br>
>> You still use sum(N) in the Jacroutine although <br>
>> you changed N to N.^2.<br>
<br>
%%%%%%%%%%%%%%%%%%%%% Equilibrium Constants, lnKp %%%%%%%%%%%%%%%%%%%%%%%<br>
function [LnKp1 LnKp2 LnKp3] = LnKp(T)<br>
% database of temperatures, K<br>
Tdb = [300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 ...<br>
1800 1900 2000];<br>
<br>
% database of lnKp values from GasEq<br>
lnKp1 = [56.864 35.9835 23.2035 14.5385 8.2618 3.5002 0.2359 3.2434...<br>
5.716 7.7828 9.536 11.04 12.3441 13.4842 14.4895 15.3812 16.1773 16.8919];<br>
lnKp2 = [11.4575 7.3537 4.9302 3.3489 2.2479 1.4448 0.8381 0.3667 0.0081...<br>
0.3117 0.5614 0.7696 0.9449 1.0942 1.2225 1.3333 1.4296 1.5141];<br>
lnKp3 = [91.6186 67.3317 52.7002 42.9072 35.8863 30.6019 26.479 23.1708...<br>
20.4575 18.191 16.2691 14.6193 13.187 11.9323 10.8234 9.8368 8.9531...<br>
8.1571];<br>
<br>
% interpolate equilibrium constant<br>
LnKp1 = interp1(Tdb,lnKp1,T,'spline','extrap');<br>
LnKp2 = interp1(Tdb,lnKp2,T,'spline','extrap');<br>
LnKp3 = interp1(Tdb,lnKp3,T,'spline','extrap');<br>
<br>
%%%%%%%%%%%%%%%%%%%%% Set of Equations to Solve %%%%%%%%%%%%%%%%%%%%%%%%%<br>
function [eq] = f(N)<br>
global Stoichiometry<br>
global E<br>
global T<br>
global P<br>
<br>
% the element balance equations<br>
eq(1:3) = Stoichiometry*N.^2  E.^2;<br>
% the equilibrium constants<br>
[LnKp1 LnKp2 LnKp3] = LnKp(T);<br>
<br>
% equilibrium equations<br>
eq(4) = LnKp1+2*log(P/sum(N.^2))+log(N(3)^2)+3*log(N(4)^2)log(N(1)^2)log(N(2)^2);<br>
eq(5) = LnKp2+log(N(5)^2)+log(N(4)^2)log(N(3)^2)log(N(2)^2);<br>
eq(6) = LnKp30.5*log(P/sum(N.^2))+log(N(2)^2)log(N(4)^2)0.5*log(N(6)^2);<br>
eq = eq';<br>
<br>
%%%%%%%%%%%%%%%%%%%%%%% Newton Raphson Method %%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
% modified from example in Numercial Methods in Chemical Engineering<br>
function solution = Newton(Function,Jacobian,x_i,tol)<br>
x = x_i;<br>
error = 2*tol;<br>
iter = 0;<br>
while error > tol<br>
F = feval(Function,x);<br>
error1 = max(abs(F));<br>
error2 = max(abs(F));<br>
J = feval(Jacobian,x);<br>
dx = J\(F);<br>
m = 1;<br>
while error2 >= error1~isreal(F)<br>
xnew = x+0.5^m*dx;<br>
F = feval(Function,xnew);<br>
error2 = max(abs(F));<br>
m = m+1;<br>
end<br>
x = xnew;<br>
F = feval(Function,x);<br>
error = max(abs(F));<br>
iter = iter+1;<br>
disp(['error = ' num2str(error)])<br>
disp(['iter = ' num2str(iter)])<br>
end<br>
solution = x;

Mon, 09 Aug 2010 01:44:06 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#769625
Gavin
> % outputs of system, moles and total moles<br>
> CH4 = sqrt(N(1));<br>
> H2O = sqrt(N(2));<br>
> CO = sqrt(N(3));<br>
> H2 = sqrt(N(4));<br>
> CO2 = sqrt(N(5));<br>
> O2 = sqrt(N(6));<br>
> nt = sqrt(sum(N));<br>
> >> nt must be the sum of the sqrt, not the <br>
> >> sqrt of the sum.<br>
> <br>
> Jac(6,4) = N(4)/sum(N)2/N(4);<br>
> Jac(6,5) = N(5)/sum(N);<br>
> Jac(6,6) = N(6)/sum(N)1/N(6);<br>
> >> You still use sum(N) in the Jacroutine although <br>
> >> you changed N to N.^2.<br>
<br>
Torsten,<br>
I've made the changes as you stated above and still get the wrong solutions. Imaginary numbers are still showing up somehow. Uploaded the latest mfile and the jacobian mfile that I'm using.<br>
Gavin<br>
<br>
mfile for model... https://files.me.com/wigging/10xkqn<br>
jacobian mfile... https://files.me.com/wigging/cpg7p8

Mon, 09 Aug 2010 06:44:55 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#769649
Torsten Hennig
%%%%%%%%%%%%%% MAIN FUNCTION using Newton Raphson method %%%%%%%%%%%%%%%%<br>
function gasf_Newton_6sJ<br>
global T<br>
global E<br>
global Stoichiometry<br>
global P<br>
<br>
% input the flows of elements entering the system<br>
nC = input('enter > nC, mol = ');<br>
nH = input('enter > nH, mol = ');<br>
nO = input('enter > nO, mol = ');<br>
E = [nC; nH; nO];<br>
% set the temperature (K)<br>
T = input('enter > T, K = ');<br>
% set the pressure (Po = 1 bar)<br>
P = input('enter > P, bar = ');<br>
%set the stiociometry matrix<br>
Stoichiometry = [1 0 1 0 1 0;...<br>
4 2 0 2 0 0;...<br>
0 1 1 0 2 2];<br>
<br>
% solve the system of equations<br>
[N] = Newton(@f,@jac,[1 1 1 1 1 1]',1e6);<br>
<br>
% outputs of system, moles and total moles<br>
CH4 = sqrt(N(1));<br>
H2O = sqrt(N(2));<br>
CO = sqrt(N(3));<br>
H2 = sqrt(N(4));<br>
CO2 = sqrt(N(5));<br>
O2 = sqrt(N(6));<br>
nt = CH4+H2O+CO+H2+CO2+O2;<br>
<br>
>> I was wrong in my last reply.<br>
>> The last seven lines should read<br>
>> CH4 = N(1)^2<br>
>> H2O = N(2)^2<br>
>> CO = N(3)^2<br>
>> H2 = N(4)^2<br>
>> CO2 = N(5)^2<br>
>> O2 = N(6)^2<br>
>> nt = CH4+H2O+CO+H2+CO2+O2<br>
<br>
% mole fractions<br>
yCH4 = CH4/nt<br>
yH2O = H2O/nt<br>
yCO = CO/nt<br>
yH2 = H2/nt<br>
yCO2 = CO2/nt<br>
yO2 = O2/nt<br>
ntt = yCH4 + yH2O + yCO + yH2 + yCO2 + yO2<br>
<br>
%%%%%%%%%%%%%%%%%%%%%%%%%% Jacobian Matrix %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
function [Jac] = jac(N)<br>
<br>
sum = N(1)^2+N(2)^2+N(3)^2+N(4)^2+N(5)^2+N(6)^2;<br>
<br>
Jac(1:3,:) = [2*N(1), 0, 2*N(3), 0, 2*N(5), 0 ;...<br>
8*N(1), 4*N(2), 0, 4*N(4), 0, 0 ;...<br>
0, 2*N(2), 2*N(3), 0, 4*N(5), 4*N(6)]; <br>
<br>
% CH4 + H2O > CO + 3*H2<br>
% lnKp1 = 2*ln(P/nt) + ln(nCO) + 3*ln(nH2)  ln(nCH4)  ln(nH2O)<br>
Jac(4,1) = (4*N(1))/sum2/N(1);<br>
Jac(4,2) = (4*N(2))/sum2/N(2);<br>
Jac(4,3) = (4*N(3))/sum+2/N(3);<br>
Jac(4,4) = (4*N(4))/sum+6/N(4);<br>
Jac(4,5) = (4*N(5))/sum;<br>
Jac(4,6) = (4*N(6))/sum;<br>
% CO + H2O > CO2 + H2<br>
% lnKp2 = ln(nCO2) + ln(nH2)  ln(nCO)  ln(nH2O)<br>
Jac(5,1) = 0;<br>
Jac(5,2) = 2/N(2);<br>
Jac(5,3) = 2/N(3);<br>
Jac(5,4) = 2/N(4);<br>
Jac(5,5) = 2/N(5);<br>
Jac(5,6) = 0;<br>
% H2 + 0.5*O2 > H2O<br>
% lnKp4 = 0.5*ln(P/nt) + ln(nH2O)  ln(nH2)  0.5*ln(nO2)<br>
Jac(6,1) = N(1)/sum;<br>
Jac(6,2) = N(2)/sum+2/N(2);<br>
Jac(6,3) = N(3)/sum;<br>
Jac(6,4) = N(4)/sum2/N(4);<br>
Jac(6,5) = N(5)/sum;<br>
Jac(6,6) = N(6)/sum1/N(6);<br>
<br>
%%%%%%%%%%%%%%%%%%%%% Equilibrium Constants, lnKp %%%%%%%%%%%%%%%%%%%%%%%<br>
function [LnKp1 LnKp2 LnKp3] = LnKp(T)<br>
% database of temperatures, K<br>
Tdb = [300 400 500 600 700 800 900 1000 1100 1200 1300 1400 1500 1600 1700 ...<br>
1800 1900 2000];<br>
<br>
% database of lnKp values from GasEq<br>
lnKp1 = [56.864 35.9835 23.2035 14.5385 8.2618 3.5002 0.2359 3.2434...<br>
5.716 7.7828 9.536 11.04 12.3441 13.4842 14.4895 15.3812 16.1773 16.8919];<br>
lnKp2 = [11.4575 7.3537 4.9302 3.3489 2.2479 1.4448 0.8381 0.3667 0.0081...<br>
0.3117 0.5614 0.7696 0.9449 1.0942 1.2225 1.3333 1.4296 1.5141];<br>
lnKp3 = [91.6186 67.3317 52.7002 42.9072 35.8863 30.6019 26.479 23.1708...<br>
20.4575 18.191 16.2691 14.6193 13.187 11.9323 10.8234 9.8368 8.9531...<br>
8.1571];<br>
<br>
% interpolate equilibrium constant<br>
LnKp1 = interp1(Tdb,lnKp1,T,'spline','extrap');<br>
LnKp2 = interp1(Tdb,lnKp2,T,'spline','extrap');<br>
LnKp3 = interp1(Tdb,lnKp3,T,'spline','extrap');<br>
<br>
%%%%%%%%%%%%%%%%%%%%% Set of Equations to Solve %%%%%%%%%%%%%%%%%%%%%%%%%<br>
function [eq] = f(N)<br>
global Stoichiometry<br>
global E<br>
global T<br>
global P<br>
<br>
% the element balance equations<br>
eq(1:3) = Stoichiometry*N.^2  E.^2;<br>
<br>
>> If your first version (before transforming<br>
>> the variables) was correct, this should read<br>
>> eq(1:3) = Stoichiometry*N.^2  E;<br>
<br>
% the equilibrium constants<br>
[LnKp1 LnKp2 LnKp3] = LnKp(T);<br>
<br>
% equilibrium equations<br>
eq(4) = LnKp1+2*log(P/sum(N.^2))+log(N(3)^2)+3*log(N(4)^2)log(N(1)^2)log(N(2)^2);<br>
eq(5) = LnKp2+log(N(5)^2)+log(N(4)^2)log(N(3)^2)log(N(2)^2);<br>
eq(6) = LnKp30.5*log(P/sum(N.^2))+log(N(2)^2)log(N(4)^2)0.5*log(N(6)^2);<br>
eq = eq';<br>
<br>
%%%%%%%%%%%%%%%%%%%%%%% Newton Raphson Method %%%%%%%%%%%%%%%%%%%%%%%%%%%<br>
% modified from example in Numercial Methods in Chemical Engineering<br>
function solution = Newton(Function,Jacobian,x_i,tol)<br>
x = x_i;<br>
error = 2*tol;<br>
iter = 0;<br>
while error > tol<br>
F = feval(Function,x);<br>
error1 = max(abs(F));<br>
error2 = max(abs(F));<br>
J = feval(Jacobian,x);<br>
dx = J\(F);<br>
m = 1;<br>
while error2 >= error1~isreal(F)<br>
xnew = x+0.5^m*dx;<br>
F = feval(Function,xnew);<br>
error2 = max(abs(F));<br>
m = m+1;<br>
end<br>
x = xnew;<br>
F = feval(Function,x);<br>
error = max(abs(F));<br>
iter = iter+1;<br>
disp(['error = ' num2str(error)])<br>
disp(['iter = ' num2str(iter)])<br>
end<br>
solution = x;<br>
<br>
Best wishes<br>
Torsten.

Mon, 09 Aug 2010 09:28:06 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#769681
Gavin
> % outputs of system, moles and total moles<br>
> CH4 = sqrt(N(1));<br>
> H2O = sqrt(N(2));<br>
> CO = sqrt(N(3));<br>
> H2 = sqrt(N(4));<br>
> CO2 = sqrt(N(5));<br>
> O2 = sqrt(N(6));<br>
> nt = CH4+H2O+CO+H2+CO2+O2;<br>
> <br>
> >> I was wrong in my last reply.<br>
> >> The last seven lines should read<br>
> >> CH4 = N(1)^2<br>
> >> H2O = N(2)^2<br>
> >> CO = N(3)^2<br>
> >> H2 = N(4)^2<br>
> >> CO2 = N(5)^2<br>
> >> O2 = N(6)^2<br>
> >> nt = CH4+H2O+CO+H2+CO2+O2<br>
> <br>
> % the element balance equations<br>
> eq(1:3) = Stoichiometry*N.^2  E.^2;<br>
> <br>
> >> If your first version (before transforming<br>
> >> the variables) was correct, this should read<br>
> >> eq(1:3) = Stoichiometry*N.^2  E;<br>
<br>
Torsten,<br>
I made the changes which seems to have taken care of the imaginary number problem. The model also runs from T = 700K1500K with no warnings (original model only ran from T = 1200K1500K with no warnings). However, the model fails to converge at T = 600K, and at T = 300K500K matlab starts to give warnings again. The warning given at T = 300K is as follows:<br>
<br>
> In gasf_Newton_6sJ>Newton at 126<br>
In gasf_Newton_6sJ at 23<br>
error = 1.4951e06<br>
iter = 162<br>
Warning: Matrix is close to singular or badly scaled.<br>
Results may be inaccurate. RCOND = 1.137846e35.<br>
<br>
It would be great if I could avoid these warnings. My goal is to fix this so I can avoid problems when I use the model for a larger system. I've posted the current mfile if you want to run it and see if you get the same warning. My inputs are:<br>
enter > nC, mol = 2<br>
enter > nH, mol = 2<br>
enter > nO, mol = 4<br>
enter > T, K = 300<br>
enter > P, bar = 1<br>
<br>
latest mfile is here... https://files.me.com/wigging/s2d25v<br>
<br>
Gavin

Mon, 09 Aug 2010 09:55:23 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#769687
Torsten Hennig
> > % outputs of system, moles and total moles<br>
> > CH4 = sqrt(N(1));<br>
> > H2O = sqrt(N(2));<br>
> > CO = sqrt(N(3));<br>
> > H2 = sqrt(N(4));<br>
> > CO2 = sqrt(N(5));<br>
> > O2 = sqrt(N(6));<br>
> > nt = CH4+H2O+CO+H2+CO2+O2;<br>
> > <br>
> > >> I was wrong in my last reply.<br>
> > >> The last seven lines should read<br>
> > >> CH4 = N(1)^2<br>
> > >> H2O = N(2)^2<br>
> > >> CO = N(3)^2<br>
> > >> H2 = N(4)^2<br>
> > >> CO2 = N(5)^2<br>
> > >> O2 = N(6)^2<br>
> > >> nt = CH4+H2O+CO+H2+CO2+O2<br>
> > <br>
> > % the element balance equations<br>
> > eq(1:3) = Stoichiometry*N.^2  E.^2;<br>
> > <br>
> > >> If your first version (before transforming<br>
> > >> the variables) was correct, this should read<br>
> > >> eq(1:3) = Stoichiometry*N.^2  E;<br>
> <br>
> Torsten,<br>
> I made the changes which seems to have taken care of<br>
> the imaginary number problem. The model also runs<br>
> from T = 700K1500K with no warnings (original model<br>
> only ran from T = 1200K1500K with no warnings).<br>
> However, the model fails to converge at T = 600K,<br>
> , and at T = 300K500K matlab starts to give warnings<br>
> again. The warning given at T = 300K is as follows:<br>
> <br>
> > In gasf_Newton_6sJ>Newton at 126<br>
> In gasf_Newton_6sJ at 23<br>
> error = 1.4951e06<br>
> iter = 162<br>
> Warning: Matrix is close to singular or badly<br>
> dly scaled.<br>
> Results may be inaccurate. RCOND =<br>
> ccurate. RCOND = 1.137846e35.<br>
> <br>
> It would be great if I could avoid these warnings.<br>
> My goal is to fix this so I can avoid problems when<br>
> n I use the model for a larger system. I've posted<br>
> the current mfile if you want to run it and see if<br>
> you get the same warning. My inputs are:<br>
> enter > nC, mol = 2<br>
> enter > nH, mol = 2<br>
> enter > nO, mol = 4<br>
> enter > T, K = 300<br>
> enter > P, bar = 1<br>
> <br>
> latest mfile is here...<br>
> https://files.me.com/wigging/s2d25v<br>
> <br>
> Gavin<br>
<br>
Just out of curiosity:<br>
Why don't you use MATLAB's fsolve for your <br>
nonlinear system of equations ?<br>
<br>
Best wishes<br>
Torsten.

Mon, 09 Aug 2010 12:47:04 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#769739
Gavin
I have tried using fsolve but ended up with incorrect results. Feel free to try with my system of equations and let me know if you get it to work correctly. I found that the newton raphson method gave more consistent answers and was easier to fix if problems arose. <br>
<br>
Gavin

Tue, 10 Aug 2010 13:57:04 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#770087
Gavin
> > Gavin<br>
> <br>
> Just out of curiosity:<br>
> Why don't you use MATLAB's fsolve for your <br>
> nonlinear system of equations ?<br>
> <br>
> Best wishes<br>
> Torsten.<br>
<br>
Torsten,<br>
<br>
I created a model using fsolve but still have problems at lower temperature ranges. At higher temperatures it will run but often gives incorrect answers. The warning that it gives at lower temperatures is the following...<br>
>>Optimizer appears to be converging to a point which is not a root.<br>
Norm of relative change in X is less than max(options.TolX^2,eps) but<br>
sumofsquares of function values is greater than or equal to sqrt(options.TolFun)<br>
Try again with a new starting guess.<br>
I receive much better results when implementing the newton raphson method that we have been discussing.<br>
<br>
Does anyone know how I can make the newton raphson model stable at lower temperature ranges?<br>
<br>
Gavin

Wed, 23 May 2012 10:05:09 +0000
Re: Newton Raphson method for chemical equilibrium system
http://www.mathworks.com/matlabcentral/newsreader/view_thread/288456#877683
Zaki
Dear Gavin,<br>
<br>
Has this problem been solve? <br>
If yes, which approach did you choose? The Fsolve or the newton Rahpson?<br>
<br>
Thanks.<br>
<br>
"Gavin" wrote in message <i3rlrg$lsm$1@fred.mathworks.com>...<br>
> > > Gavin<br>
> > <br>
> > Just out of curiosity:<br>
> > Why don't you use MATLAB's fsolve for your <br>
> > nonlinear system of equations ?<br>
> > <br>
> > Best wishes<br>
> > Torsten.<br>
> <br>
> Torsten,<br>
> <br>
> I created a model using fsolve but still have problems at lower temperature ranges. At higher temperatures it will run but often gives incorrect answers. The warning that it gives at lower temperatures is the following...<br>
> >>Optimizer appears to be converging to a point which is not a root.<br>
> Norm of relative change in X is less than max(options.TolX^2,eps) but<br>
> sumofsquares of function values is greater than or equal to sqrt(options.TolFun)<br>
> Try again with a new starting guess.<br>
> I receive much better results when implementing the newton raphson method that we have been discussing.<br>
> <br>
> Does anyone know how I can make the newton raphson model stable at lower temperature ranges?<br>
> <br>
> Gavin