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;