SAS: (Qusi-)Newton method Home made code

是 否 +2 论坛币

k人 参与回答

经管之家送您一份

应届毕业生专属福利!

求职就业群

赵安豆老师微信：zhaoandou666

经管之家联合CDA

送您一个全额奖学金名额~ !

立即领取

感谢您参与论坛问题回答

经管之家送您两个论坛币！

+2 论坛币

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. to  Numerical Analysis 5th ed. By Burden Faires           ***/
/* x0--iteration initial value vector                               ***/
/*** N--MAX iteration number  TOL-- 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 value  x0={0.1,0.1,-0.1};

                                          solution is


                                                 X

                                                   0.5
                                             0.0000144
                                             -0.523333


                                          iteration number

                                                     5


2.2.                       initial value  x0={1,1,1};  
                                          solution is


                                                 X

                                                   0.5
                                             0.0000142
                                             -0.523333


                                          iteration number
                                                
                                                    10


2.3.                       initial value  x0={10,7,30};  


                                      solution is


                                                 X

                                                   0.5
                                              0.000015
                                             -0.523333


                                          iteration number
                                                    24

2.4 initial value  x0={10, 7,-1};

                                             solution is


                                                 X

                                             0.5000008
                                              0.000024
                                             -0.523333


                                          iteration number                                                

                                                    14

扫码加我 拉你入群

请注明：姓名-公司-职位

以便审核进群资格，未注明则拒绝

分享0 收藏0 回帖

关键词：Newton Method code Made home Method code Made home