Newton Raphson method for chemical equilibrium system
in [Matlab]
|
Prev: OPC "trend" function
Next: variables in titles ??
From: TideMan on 4 Aug 2010 00:50 On Aug 4, 4:28 pm, "Roger Stafford" <ellieandrogerxy...(a)mindspring.com.invalid> wrote: > TideMan <mul...(a)gmail.com> wrote in message <126edaf0-e1e5-4760-a737-6d680d5f9...(a)m17g2000prl.googlegroups.com>... > > No, I'm not interested in looking at your code. > > All I can tell you is that in 40 years since graduation, I've done a > > lot of work using Newton Raphson and every time I've had a problem it > > has been because I got the algebra wrong in assembling the Jacobian. > > First, you need to check your partial differentiation. > > Next, you need to check that you have correctly transcribed the > > partial differentials into Matlab code. > > - - - - - - - - - - > It is conceivable to me that the Newton-Raphson method for multiple dimensions could run into convergence problems with real-life 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. > > 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 Newton-Raphson could arise. It would depend on the nature of the reaction rates involved. > > Roger Stafford Well, in my experience with water waves (which admittedly are generally pretty smooth in time and space), if something as robust as Newton Raphson blows up it is pointing to a problem - the waves have broken or the numerical scheme has gone unstable - but usually the problem is not with Newton Raphson itself, which is what the OP is assuming.
From: Steve Amphlett on 4 Aug 2010 03:58 "Gavin " <wigging(a)gmail.com> wrote in message <i3a98v$lck$1(a)fred.mathworks.com>... > 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 A-matrix 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 > Raphson helps, but I am unsure of how to apply this to my problem. > > function solution = Newton(Function,Jacobian,x_i,tol) > x = x_i; > error = 2*tol; > while error > tol > F = feval(Function,x); > error1 = max(abs(F)); > error2 = max(abs(F)); > J = feval(Jacobian,x); > dx = J\(-F); > i=1; > while error2 >= error1||~isreal(F) > xnew = x+0.5^i*dx; > F = feval(Function,xnew); > error2 = max(abs(F)); > i = i+1; > end > x = xnew; > F = feval(Function,x); > error = max(abs(F)); > disp(['error = ' num2str(error)]) > end > solution = x; Gavin, 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: http://elib.uni-stuttgart.de/opus/volltexte/2006/2725/pdf/Diss_Grill.pdf - Steve
From: Gavin on 4 Aug 2010 10:00 Thanks Steve, but I can't read German. I'm trying to solve the following system of equations: where: CH4 = N(1), H2O = N(2), CO = N(3), H2 = N(4), CO2 = N(5), O2 = N(6) CH4+CO+CO2-C = 0 4*CH4+2*H2O+2*H2-H = 0 H2O+CO+2*CO2+2*O2-O = 0 -LnKp1+2*log(P/sum(N))+log(N(3))+3*log(N(4))-log(N(1))-log(N(2)) = 0 -LnKp2+log(N(5))+log(N(4))-log(N(3))-log(N(2)) = 0 -LnKp3-0.5*log(P/sum(N))+log(N(2))-log(N(4))-0.5*log(N(6)) = 0 Known inputs are LnKp, C, H, O. Where LnKp is the log of the equilibrium constant for that particular reaction at the specified temperature.
From: Torsten Hennig on 4 Aug 2010 06:55 > > Thanks Steve, but I can't read German. I'm trying to > solve the following system of equations: > where: CH4 = N(1), H2O = N(2), CO = N(3), H2 = N(4), > CO2 = N(5), O2 = N(6) > > CH4+CO+CO2-C = 0 > 4*CH4+2*H2O+2*H2-H = 0 > H2O+CO+2*CO2+2*O2-O = 0 > -LnKp1+2*log(P/sum(N))+log(N(3))+3*log(N(4))-log(N(1)) > -log(N(2)) = 0 > -LnKp2+log(N(5))+log(N(4))-log(N(3))-log(N(2)) = 0 > -LnKp3-0.5*log(P/sum(N))+log(N(2))-log(N(4))-0.5*log(N > (6)) = 0 > > Known inputs are LnKp, C, H, O. Where LnKp is the > log of the equilibrium constant for that particular > reaction at the specified temperature. The first thing I'd do is to substitute N(i) = N~(i)^2 for i=1,...,6 to avoid negative arguments within the logarithm. Best wishes Torsten.
From: Maxx Chatsko on 4 Aug 2010 11:05
"Gavin " <wigging(a)gmail.com> I'm pretty handy with pchem and differentials. Would you mind sending me your code? Maxx |