Solving Systems of Nonlinear Algebraic Equations using Newton's Method

Solving Systems of Nonlinear Algebraic Equations using Newton's Method

The standard form for these problems is that we are trying to find an X that satisfies

F(X) = 0

where both F and X are vectors.

Newton's method is based on the Taylor Expansion:

F_i(X) = F_i(X_o) + Sum { (dF_i/dX_k) * DX_k } + higher order terms

where DX=X-X_o , and which in matrix form is

F(X) = F(X_o) + J(X_o) * DX + higher order terms

where the elements of the Jacobian matrix J evaluated at X_o are given by

J_ik(X_o) = (dF_i/dX_k)|_Xo

so if we want F(X) to vanish, and the higher order terms are negligible, we should solve

0 = F(X_o) + J(X_o) * (X-X_o)

i.e.

J(X_o) * DX = -F(X_o)

This is a linear system of equations, so we can solve it for DX=(X-X_o). So here is Newton's algorithm:

0) Make an initial guess at the solution X

1) Evaluate J(X) and F(X)

2) Solve J*DX = -F

3) X = X + DX

4) if for all k, |DX_k| < RTOL * |X_k | + ATOL we have converged.

5) if not, go back to step 1

It can be proved that this algorithm converges quadratically (i.e. very rapidly!) if the initial guess is good enough. However, this method can be very poor if the initial guess is not good enough. It is advisable to build in a limit on the maximum number of iterations; if the program is taking a lot of iterations it is unlikely to ever converge.

The user must supply a function which evaluates F(X) and J(X) for any input X. Any error in the function which computes F or J will generally make the method fail completely, in fact sometimes even minor numerical errors (e.g. from using numerical derivatives rather than exact analytical derivatives to compute J) can cause serious problems. So this function must be checked very carefully!

For a more detailed discussion, see the chapter on solving systems of equations in Numerical Recipes in C. For a careful discussion of Newton's method in one dimension, see the Course Notes.