Chapter Contents Previous Next
 General Statistics Examples

## Example 8.2: Newton's Method for Solving Nonlinear Systems of Equations

This example solves a nonlinear system of equations by Newton's method. Let the nonlinear system be represented by

F(x) = 0
where x is a vector and F is a vector-valued, possibly nonlinear, function.

In order to find x such that F goes to 0, an initial estimate x0 is chosen, and Newton's iterative method for converging to the solution is used:

xn+1 = xn - J-1(xn) F(xn)
where J(x) is the Jacobian matrix of partial derivatives of F with respect to x.

For optimization problems, the same method is used, where F(x) is the gradient of the objective function and J(x) becomes the Hessian (Newton-Raphson).

In this example, the system to be solved is

The code is organized into three modules: NEWTON, FUN, and DERIV.
```      /*     Newton's Method to Solve a Nonlinear Function      */
/* The user must supply initial values,                   */
/* and the FUN and DERIV functions.                       */
/* on entry: FUN evaluates the function f in terms of x   */
/* initial values are given to x                          */
/* DERIV evaluates jacobian j                             */
/* tuning variables: CONVERGE, MAXITER.                   */
/* on exit: solution in x, function value in f close to 0 */
/* ITER has number of iterations.                         */

start newton;
run fun;          /* evaluate function at starting values */
do iter=1 to maxiter  /* iterate until maxiter iterations */
while(max(abs(f))>converge);            /* or convergence */
run deriv;                /* evaluate derivatives in j */
delta=-solve(j,f);      /* solve for correction vector */
x=x+delta;                    /* the new approximation */
run fun;                      /* evaluate the function */
end;
finish newton;

maxiter=15;                    /* default maximum iterations */
converge=.000001;           /* default convergence criterion */

/* User-supplied function evaluation */
start fun;
x1=x[1] ;
x2=x[2] ;                           /* extract the values */
f= (x1+x2-x1*x2+2)//
(x1*exp(-x2)-1);                 /* evaluate the function */
finish fun;

/* User-supplied derivatives of the function */
start deriv;
/* evaluate jacobian */
j=((1-x2)||(1-x1) )//(exp(-x2)||(-x1*exp(-x2)));
finish deriv;

do;
print "Solving the system: X1+X2-X1*X2+2=0, X1*EXP(-X2)-1=0" ,;
x={.1, -2};     /* starting values */
run newton;
print x f;
end;
```

The results are shown below.

 The SAS System

 Solving the system: X1+X2-X1*X2+2=0, X1*EXP(-X2)-1=0

 X F 0.0977731 5.3523E-9 -2.325106 6.1501E-8

 Chapter Contents Previous Next Top