用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)來調試計算過程


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM