用R语言求解非线性方程
从本质上来说,Newtons就是用迭代方式,使近似解(泰勒公式)不断的逼近真实解,当满足精度要求时,即可认为近似解为真实解
下面用R语言实现Newtons法
Newtons<-function(fun,x,ep=1e-5,it_max=100) ##fun为需要求解的方程(组),x为初始解,ep为精度要求,it_max为最大迭代次数
{
index<-0 ##指示是否完成迭代成功,满足精度要求
k<-1 ##迭代次数
while(k<=it_max)
{
x1<-x;obj<-fun(x);
x<-x-solve(obj$J,obj$f) ##obj$J即为方程(组)的jacobj矩阵
norm<-sqrt((x-x1)%*%(x-x1)) ##x(k+1)y与x(k)精度之差
if(norm<ep)
{
index<-1
break
}
k<-k+1
}
obj<-fun(x);
list(root=x,it=k,index=index,FunVal=obj$f,Jacobi=obj$J)
}
现在来构造需要求解的方程(组)和它的jacobj矩阵
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) ##jacobj矩阵,分别对x1和x2求一阶导
list(f=f,J=J)
}
计算结果:
> Newtons(funs,c(0,1))
$root
[1] 1 2
$it
[1] 6
$index
[1] 1
$FunVal
[1] 1.598721e-14 6.217249e-15
$Jacobi
[,1] [,2]
[1,] 2 4
[2,] -1 2
有兴趣的同学,可以用debug(Newtons)来调试计算过程