从本质上来说，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)来调试计算过程

用R语言求解非线性方程