【学习记录-R】牛顿迭代解非线性方程组及单个方程

newtons非线性方程组迭代过程:

 

思路分解：

1. index=1,max=100，n=0,每迭代一次n++,在100次以内算出index=1标记；100次以内无法算出符合精度的解，index=0标记
2. Function(x)自建函数f,进行J()*f()计算，输出list A，含多个结果，A$J,A$f
3. 迭代，输入x,y原始值,精度ep,函数f,自建函数y=x,x=y-J[(y)]-1 * f(y),用x-y与ep进行比较
4. 输出结果

 参考：

用R语言求解非线性方程 - R爱好者 - 博客园 (cnblogs.com)

R语言中%*%运算符 - 小鲨鱼2018 - 博客园 (cnblogs.com)

例题：

  

 （我）难理解的地方：

1.精度的表示，开始看到%*%符号下意识以为是矩阵乘积，但是其结果无法用于做判断条件，后来发现代码中的%*%运算实际上是两个向量在进行，精度判断条件为:

[x(k+1)-x(k)]%*%[x(k+1)-x(k)]<ep,即每个解的第k+1次迭代结果与上一次差值的和的开方。

 

2.雅可比矩阵的计算

调用numDeriv包计算：

library(numDeriv)
fw = function(x){
y1 = x[1]^2+x[2]^2-7*x[1]-9.4
y2 = x[2]^2*x[1]-54.648
c(y1,y2)
}

f<-function(x)
{
  J<-jacobian(func=fw,x)#雅各比矩阵
  fx<-c(x[1]^2+x[2]^2-7*x[1]-9.4,x[2]^2*x[1]-54.648) #方程组
  list(fx=fx,J=J)
}

参考：R中计算jacobian矩阵和hessian矩阵 | Keep Learning (biocuckoo.org)

手动计算：

代码：
 

f<-function(x)
{
  fx<-c(x[1]^2+x[2]^2-7*x[1]-9.4,x[2]^2*x[1]-54.648) ##方程组
  J<-matrix(c(2*x[1]-7,2*x[2],x[2]^2,x[2]*x[1]-1),nrow=2,byrow=T) ##jacobj矩阵，分别对x(x[1]),y(x[2])求一阶导,按行
  list(fx=fx,J=J)
}

完整代码：

library(numDeriv)
fw = function(x){
  y1 = x[1]^2+x[2]^2-7*x[1]-9.4
  y2 = x[2]^2*x[1]-x[2]-54.648
  c(y1,y2)
}

f<-function(x)
{
  J<-jacobian(func=fw,x)#雅各比矩阵
  fx<-c(x[1]^2+x[2]^2-7*x[1]-9.4,x[2]^2*x[1]-x[2]-54.648) #方程组
  list(fx=fx,J=J)
}

Newtons<-function(f,x,ep=1e-15,max=50) ##f需要求解的方程(组)和它的jacobj矩阵,x为初始解，ep为精度要求，max为最大迭代次数
{
  index<-1 #在n>max之前跳出，则为成功找到
  n<-1 #迭代次数
  while(norm>ep)
  {
    x1<-x
    obj<-f(x1)
    x<-x1-solve(obj$J,obj$fx) ##obj$J即为方程(组)的jacobj矩阵
    norm<-sqrt((x-x1)%*%(x-x1)) ##x(k+1)y与x(k)精度之差
    n=n+1
    if(n>max)
    {
      index<-0
      break
    }
  }
  obj<-f(x)
  list(root=x,it=n,index=index,FunVal=obj$fx)
}

结果：

  

单个方程组迭代过程：

解方程例题 （方程组同理）：

已知总体的分布密度为：，x的样本值为：

0.1，0.2，0.9，0.8，0.7，0.7，求α的估计量

 

 

 

 牛顿迭代解方程代码：

f<-function(x,mean){
  x1=x-(x+1)*(x+2)+mean*(x+2)^2
  list(x1=x1)
}
newton<-function(f,mean,x1,x2=-999,ep=10^(-5),max=50){
  index=1;n=0
  while(abs(x1-x2)>ep){
    n=n+1
    x2<-x1
    x<-f(x2,mean)
    x1<-x$x1#如果直接x1<-f(),x1不能参与x1-x2
    if(n>=max){
      index=0;list(x1=x1,index=index,n=n)
      break
    }
  }
  list(x1=x1,index=index,n=n)
}

 结果：