SAS: (Qusi-)Newton method Home made code-经管之家官网！

Hi, Guys:
Here I show how to solve equation system/equation by quasi-newton method. Other method like
newton, secant method are similar.This is inspired by a friend's help posting. Hope this is helpful.
Note: Newton algorithm is not always convergent. You'd better draw the picture to guess solution and
let inital value close to real solution.
Just my 2 cents.
Happy holiday and 2010.
1. Code
/*********************************************************************/
/**Quai-Newton method (Broyden) *********/
/**Ref. toNumerical Analysis 5th ed. By Burden Faires ***/
/* x0--iteration initial value vector ***/
/*** N--MAX iteration numberTOL-- Tolerence */
/*We(gzjb) show quasi newton method by slove equation systems: */
/*****f1(x)= 3x1-cos(x2x3)-1/1=0 */
/*****f2(x)= x1^2-81(x2+0.1)^2+sinx3+1.04=0 */
/**** f3(x)= exp(-x1x2)+20x3+9.466667=0 */
/* i.e. F(x)=0, where F(x)=(f1(x),f2(x),f3(x))', x=(x1,x2,x3)' a column vector**/
/** You can solve other equation/equ sys using same algorithm */
/** Advantage of Broydon: Only calculate inverse of Jacobian matrix ONCE */
/** SAS version: SAS 9.1.3 Service Pack 4 *****/
/************************************************************************/
proc iml;
reset noprint;
start quasi_newton(x0,N,TOL);
A0=j(3,3,0); v=j(3,1,0);
* calculate Jacobian matrix at x0;
A0[1,1]=3; A0[1,2]=x0[3]*sin(x0[2]*x0[3]); A0[1,3]=x0[2]*sin(x0[2]*x0[3]);
A0[2,1]=2*x0[1]; A0[2,2]= -162*(x0[2]+0.1); A0[2,3]= cos(x0[3]);
A0[3,1]=-x0[2]*exp(-x0[1]*x0[2]); A0[3,2] = -x0[1]*exp(-x0[1]*x0[2]); A0[3,3]= 20;
* calculate F(x0);
v[1]=3*x0[1]-cos(x0[2]*x0[3])-1/2;
v[2]= x0[1]**2-81*(x0[2]+0.1)**2+sin(x0[3])+1.06;
v[3]= exp(-x0[1]*x0[2])+20*x0[3]+9.466667;

A=ginv(A0);
s=-A*v;x=x0+s;
* iteration algorithm;
do k=2 to N;
w=v;
v[1]= 3*x[1]-cos(x[2]*x[3])-1/2;
v[2]= x[1]**2-81*(x[2]+0.1)**2+sin(x[3])+1.06;
v[3]= exp(-x[1]*x[2])+20*x[3]+9.466667;
y=v-w;
z=-A*y;
p=-t(s)*z;
ut=t(s)*A;
if abs(p)< 0.0001*TOL then do;
print x; abort;
end;
else A=A+(s+z)*ut/p;
s=-A*v;
x=x+s;
if sqrt(ssq(s)) < TOL then do;
print 'solution is'; print x;
print 'iteration number'; print k;
abort;
end;
end;
print 'Maximum number of iteration exceed';
finish quasi_newton;

*x0={0.1,0.1,-0.1};
*x0={1,1,1};
x0={10,7,30};
N=500; TOL=0.0001;
run quasi_newton(x0,N,TOL);
quit;

2. Running results
2.1. initial valuex0={0.1,0.1,-0.1};
solution is
X
0.5
0.0000144
-0.523333
iteration number
5
2.2. initial valuex0={1,1,1};
solution is
X
0.5
0.0000142
-0.523333
iteration number

10
2.3. initial valuex0={10,7,30};
solution is
X
0.5
0.000015
-0.523333
iteration number
24
2.4 initial valuex0={10, 7,-1};
solution is
X
0.5000008
0.000024
-0.523333
iteration number
14