From: TideMan <>
Newsgroups: comp.soft-sys.matlab
Subject: Re: Newton Raphson method for chemical equilibrium system
Date: Tue, 3 Aug 2010 16:58:19 -0700 (PDT)
Lines: 53
Message-ID: <>
References: <i3a98v$lck$>
Mime-Version: 1.0
Content-Type: text/plain; charset=ISO-8859-1
Content-Transfer-Encoding: quoted-printable
X-Trace: 1280879900 16272 (3 Aug 2010 23:58:20 GMT)
NNTP-Posting-Date: Tue, 3 Aug 2010 23:58:20 +0000 (UTC)
Injection-Info:; posting-host=; 
User-Agent: G2/1.0
X-HTTP-UserAgent: Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US; rv: 
	Gecko/20100722 Firefox/3.6.8,gzip(gfe)
Xref: comp.soft-sys.matlab:658806

On Aug 4, 11:38 am, "Gavin " <> 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...............