Example 3: Solve with Newton-Raphson Algorithm

* SOLVE TWO NONLINEAR EQUATIONS WITH TWO UNKNOWNS;

* solve 2 non-linear equations for x and y;
* F1:  x*x + y*y = 2;
* F2:  x*x - y*y = 1;

PROC IML;
niter = 9;  * enter: niter = number of interations;

* set initial estimates of the solution vector : xx = | x y |;
xx=1:2; xx[1]=1; xx[2]=-2;

* Initialize two computation vectors ;
incr=J(1,2,.); tmp=J(1,2,.);

fn=J(2,1,.);  * initialize function matrix (column vector);
DR=J(2,2,.);  * initialize 2x2 derivative matrix;
print fn dr;
*quit;

* initialize output matrix to collect
  the solution estimates and function values at each iteration ;

out=J(niter,4,.);

DO i=1 TO niter;

* enter the two nonlinear functions in the
  following notation to solve for x and y;

* x*x + y*y - 2 = 0;
* x*x - y*y - 1 = 0;

FN[1,]=xx[1]**2 + xx[2]**2 - 2;  * evaluate functions at x and y;
FN[2,]=xx[1]**2 - xx[2]**2 - 1;  * for iteration i;
print i niter fn;

out[i,] = xx||fn`; * collect values of estimates and functions;
print out;

* compute matrix of partial derivatives
  at values of xx for iteration i;
*                                          Partial Derivaties;
DR[1,1]= 2*xx[1]; DR[1,2]= 2*xx[2];     * | dF1/dx  dF1/dy  | ;
DR[2,1]= 2*xx[2]; DR[2,2]=-2*xx[1];     * | dF2/dx  dF2/dy  | ;

print dr, fn, xx;

incr = INV(DR)*fn;  * compute the increment to add to solution;
tmp = xx`-incr ;
xx = tmp`;          * update parameter estimates;
print incr, xx;

END;

* place results in a SAS dataset;

CREATE results FROM out[colname={"x","y","F1","F2"}];
APPEND FROM out;

RUN;
quit;

PROC PRINT DATA=results NOobs n; RUN;