Newton's Method for Solving Nonlinear Systems

/****************************************************************/
/*          S A S   S A M P L E   L I B R A R Y                 */
/*                                                              */
/*    NAME: NEWTON                                              */
/*   TITLE: Newton's Method for Solving Nonlinear Systems       */
/* PRODUCT: IML                                                 */
/*  SYSTEM: ALL                                                 */
/*    KEYS: MATRIX  SUGI6                                       */
/*   PROCS: IML                                                 */
/*    DATA:                                                     */
/*                                                              */
/* SUPPORT: Rick Wicklin                UPDATE: SEP2013         */
/*     REF:                                                     */
/*    MISC:                                                     */
/*                                                              */
/****************************************************************/


/**********************************************************/
/*     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 parameters: CONVERGE, MAXITER.                  */
/* On exit: return x such that FUN(x) is close to 0       */
proc iml;
/* User-supplied function evaluation */
start Fun(x);
   x1 = x[1];  x2 = x[2];           /* extract components */
   f1 = x1 + x2 - x1*x2 + 2;
   f2 = x1*exp(-x2) - 1;
   return( f1 // f2 );
finish Fun;

/* User-supplied derivatives of the function */
start Deriv(x);
   x1 = x[1]; x2 = x[2];
   df1dx1 = 1 - x2;
   df1dx2 = 1 - x1;
   df2dx1 = exp(-x2);
   df2dx2 = -x1 * exp(-x2);
   J = (df1dx1 || df1dx2 ) //
       (df2dx1 || df2dx2 );
   return( J );                    /* Jacobian matrix */
finish Deriv;

/* Implementation of Newton's method with default arguments */
/* By default, maximum iterations (maxiter) is 25           */
/*             convergence criterion (converge) is 1e-6     */
start NewtonMethod(x0, maxIter=25, converge=1e-6);
   x = x0;
   f = Fun(x);       /* evaluate function at starting values */
   do iter = 1 to maxiter           /* iterate until maxiter */
      while(max(abs(f))>converge);         /* or convergence */
      J = Deriv(x);                  /* evaluate derivatives */
      delta = -solve(J, f);   /* solve for correction vector */
      x = x + delta;                /* the new approximation */
      f = fun(x);                   /* evaluate the function */
   end;
   /* return missing if no convergence */
   if iter > maxIter then
      x = j(nrow(x0),ncol(x0),.);
   return( x );
finish NewtonMethod;

print "Solve the system: X1+X2-X1*X2+2=0, X1*EXP(-X2)-1=0" ,;
x0 = {.1, -2};     /* starting values */
x = NewtonMethod(x0);
f = Fun(x);
print x f;