R语言实现牛顿迭代算法

一粒沙

发布于 2019-11-04 17:03:15

13.2K00

代码可运行

文章被收录于专栏：R语言交流中心R语言交流中心

运行总次数：0

代码可运行

我们今天给大家介绍一个用来迭代的算法牛顿迭代法（Newton's method）。单变量下又称为切线法。它是一种在实数域和复数域上近似求解方程的方法。首先我们看下牛顿迭代算法的公式：

其中，Xn+1和Xn主要是指的在n+1和n这个位置的X值。

以上公式的推导可以用泰勒展开公式进行推导。当然上面只是一元一阶的情况的，如果更复杂的那就是多元多阶的情况需要用到更复杂的Hessian矩阵获得以上的通用公式模型。

关于推导的部分，我们就不展开了。接下来我们直接用一个R语言的实例来看下，牛顿迭代是如何工作的。我们看下下面这个例题：

以上就是简单的一元函数求解，当然我们基于我们数学的基础也可以人工展开计算，但是当次幂升到很高，那我们就无从下手了，这时候就可以直接通过牛顿迭代进行获取根。我们就基于上面的例子进行程序设计：

funs=function(x){
  f=x
  J=(2*x^3-4*x^2+3*x-6)/(6*x^2-8*x+3)
  list(f=f,J=J);
                            }
                                                
#Newton迭代法
Newtons=function(fun,x,ep=1e-5,it_max=1000){
 index=0;k=1
 while(k<=it_max){
 norm=abs(fun(x)$J)
 x=x-fun(x)$J
       if(norm<ep){ index=1;break}
       k=k+1
      }
   list(root=x,it=k,index=index)
}                                               
Newtons(funs,1)

通过上面的程序计算就得到我们的根：

上面root就是我们得到的根，it指的迭代的次数，index指的最后的结果1代表找到根；0代表没找到根。

上面是一个简单的例子，那么如果多元的函数我们如何去计算它的根呢，那么就需要我们的Hessian矩阵方法帮我们进行一步转化，然后再生成我们想要的函数。接下来看一下实例：

上面这个例子呢，其实是已经帮大家转化好了对应的矩阵方程，一般我们遇到的那可就不这么简单了，可能是上面的一个方程式组，更可能是只有一个函数方程，需要我们进行下一步的转化，一直到最后转化为我们想要的矩阵方程式。然后构建我们的程序（以下程序引自网友PPT）：

funs=function(x){
  f=c(x[1]^2+x[2]^2-5,(x[1]+1)*x[2]-(3*x[1]+1));
  J=matrix(c(2*x[1],2*x[2],x[2]-3,x[1]+1),nrow=2,byrow=T);
  list(f=f,J=J);
                            }
 
#Newton迭代法
Newtons=function(fun,x,ep=1e-5,it_max=100){
 index=0;k=1
 while(k<=it_max){
 x1=x;obj=fun(x);x=x-solve(obj$J,obj$f);norm=sqrt((x-x1)%*%(x-x1));
       if(norm<ep){ index=1;break}
       k=k+1
      }
   obj=fun(x)
   list(root=x,it=k,index=index,Funval=obj$f)
}
 
Newtons(funs,c(0,1))