Newton Raphson method for chemical equilibrium system
in [Matlab]
|
Prev: OPC "trend" function
Next: variables in titles ??
From: Gavin on 3 Aug 2010 19:38 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;
From: TideMan on 3 Aug 2010 19:58 On Aug 4, 11:38 am, "Gavin " <wigg...(a)gmail.com> wrote: > 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; Just one comment on your code: Apparently, you are restricting the increment in x to avoid moving into the complex territory of F, yet you are using i as the counter and exponent. This is foolish because next thing you'll be doing something that assumes i is sqrt(-1). Believe me, this happens, and it can be quite tricky to find the bug. Use something else for the counter, like m, for example. Other than that, your code looks good to me. So, if it fails when you add another equation, it points to an error in how you've assembled the Jacobian. Are you sure your algebra is correct? It's easy to make a mistake when doing partial differentiation of a complicated equation. All it takes is one stray negative sign to screw everything up...............
From: Gavin on 3 Aug 2010 23:28 TideMan <mulgor(a)gmail.com> wrote in message <ec2acf45-eb15-4a32-bc4d-4199af5dbbab(a)n19g2000prf.googlegroups.com>... > Just one comment on your code: > Apparently, you are restricting the increment in x to avoid moving > into the complex territory of F, yet you are using i as the counter > and exponent. This is foolish because next thing you'll be doing > something that assumes i is sqrt(-1). Believe me, this happens, and > it can be quite tricky to find the bug. Use something else for the > counter, like m, for example. > > Other than that, your code looks good to me. > So, if it fails when you add another equation, it points to an error > in how you've assembled the Jacobian. Are you sure your algebra is > correct? It's easy to make a mistake when doing partial > differentiation of a complicated equation. All it takes is one stray > negative sign to screw everything up............... TideMan, 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. Gavin
From: TideMan on 3 Aug 2010 23:46 On Aug 4, 3:28 pm, "Gavin " <wigg...(a)gmail.com> wrote: > TideMan <mul...(a)gmail.com> wrote in message <ec2acf45-eb15-4a32-bc4d-4199af5db...(a)n19g2000prf.googlegroups.com>... > > Just one comment on your code: > > Apparently, you are restricting the increment in x to avoid moving > > into the complex territory of F, yet you are using i as the counter > > and exponent. This is foolish because next thing you'll be doing > > something that assumes i is sqrt(-1). Believe me, this happens, and > > it can be quite tricky to find the bug. Use something else for the > > counter, like m, for example. > > > Other than that, your code looks good to me. > > So, if it fails when you add another equation, it points to an error > > in how you've assembled the Jacobian. Are you sure your algebra is > > correct? It's easy to make a mistake when doing partial > > differentiation of a complicated equation. All it takes is one stray > > negative sign to screw everything up............... > > TideMan, > > 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. > > Gavin 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.
From: Roger Stafford on 4 Aug 2010 00:28
TideMan <mulgor(a)gmail.com> wrote in message <126edaf0-e1e5-4760-a737-6d680d5f9dfb(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 |