R语言 牛顿-拉夫森迭代法求方程组

牛顿-拉夫森迭代法：

xk+1=xk−[f′(x)]−1f(x)

其中，f′(x)为Jacobi矩阵。
例，设非线性方程组为：
x2+y2=1,
y2=2x
求方程组的解。
Jocabi行列式：

R代码如下：

---

fun <- function(x){  f <- c(x[1]^2+x[2]^2-1, x[2]^2-2*x[1])  joc <- matrix(c(2*x[1],-2,2*x[2],2*x[2]),nr=2)  list(f=f,jac=jac)}Newton <-function(fun, x,eps = 1e-5){  k <- 0  repeat{  k <- k+1    x1 <- x    obj <- fun(x)    x <- x - solve(obj$jac,obj$f)    if((x-x1)%*%(x-x1)<eps){      return(list(root=x,iter=k))      break    }  }}

---

最后设初始值为c(1,1.2). 注: 选择初始值必须式Jacobi行列式不为零。

Newton(fun,c(1,1.2)

$root
[1] 0.4142136 0.9101797

$iter
[1] 4
而利用rootSolve包解方程组multiroot(model,c(1,1.2)解出结果与上述结果一致，而迭代次数为5.