去找高斯解。。。哈哈。。

泰勒展开
线性搜索等都可以

@解方程的根@
/*f(x)的泰勒展开:
        f(x)=f(a)+f'(a)/1!*(x-a)+f''(a)/2!*(x-a)^2+ ... +fn(a)/n!*(x-a)^n+...
忽略高阶项得到，f(x)=f(a)+f'(a)*(x-a) .
在x0处有y=f(x0)+f'(x0)*(x-x0).
第一步：
        令0= f(x0)+f'(x0)*(x-x0)
    得到x=x0-f(x0)/f'(x0) .
继续
        x1=x0-f(x0)/f'(x0)
        x2=x1-f(x1)/f'(x1)
        x3=x2-f(x2)/f'(x2)
        ..................
        x[k]=x[k-1]-f(x[k-1])/f'(x[k-1]).
直到|x[k]-x[k-1]|<1e-5 迭代停止*/
        new; cls;
        a=2;
        fn f(x)=x^2-a;

        start=1;
        tol=1e-5; h=1e-6;
        x_1=start;
        i=1;
        do while i<=10000;
              f1x_1=(f(x_1+h)-f(x_1-h))/(2*h);@f(x)的一阶导数@
              x=x_1-f(x_1)/f1x_1;
              if abs(x_1-x)<tol;
                    break;
              endif;
              x_1=x;
              i=i+1;
        endo;
        print x;
        print sqrt(a);

@或者可以把f(x)的一阶导数计算出来@
        new; cls;
        a=2;
        fn f(x)=x^2-a;
        fn f1(x)=2*x; @f(x)的一阶导数@

        start=-2;tol=1e-5;x_1=start;
        i=1;do while i<=10000;
              x=x_1-f(x_1)/f1(x_1);
              if abs(x_1-x)<tol;
                    break;
              endif;
              x_1=x;
           i=i+1;endo;
        print x;

/*lesson2：newton调用命令的使用*/
        new; cls;
        a=2;
        fn f(x)=x^2-a;
        start=-2;
        print newton(&f,start);

proc newton(&f,start);
local tol,x_1,x,i,h,f1x_1,f:proc;
   tol=1e-5; h=1e-6;
   x_1=start;
   i=1;
   do while i<=10000;
      f1x_1=(f(x_1+h)-f(x_1-h))/(2*h);
      x=x_1-f(x_1)/f1x_1;
      if abs(x_1-x)<tol;
         break;
      endif;
      x_1=x;
      i=i+1;
   endo;
   retp(x);
endp;

/*lesson2：Jacobian数值分析*/
        new; cls;
        a=2;
        fn f(x)=x^2-a;
        fn f1(x)=2*x;
        start=-2;
        print newtonJ(&f,&f1,start);

proc newtonJ(&f,&f1,start);
   local tol,x_1,x,i,f:proc,f1:proc;
   tol=1e-5;
   x_1=start;
   i=1;
   do while i<=10000;
      x=x_1-f(x_1)/f1(x_1);
      if abs(x_1-x)<tol;
         break;
      endif;
      x_1=x;
      i=i+1;
   endo;
   retp(x);
endp;

版主辛苦了

谢谢