R语言用newton法解非线性方程组

Newton法是一种常用的求解非线性方程组的方法，也适用于R语言。假设我们要求解如下非线性方程组：

$$
\begin{cases}
f_1(x_1, x_2, \cdots, x_n) = 0 \\
f_2(x_1, x_2, \cdots, x_n) = 0 \\
\vdots \\
f_n(x_1, x_2, \cdots, x_n) = 0
\end{cases}
$$

其中，$x_1, x_2, \cdots, x_n$ 是未知数，$f_1, f_2, \cdots, f_n$ 是已知函数。Newton法的基本思路是：从一个初始点 $(x^{(0)})$ 开始迭代，使用一阶泰勒展开式来逼近非线性方程组的解，得到迭代公式：

$$
x^{(k+1)} = x^{(k)} - J^{-1}(x^{(k)})F(x^{(k)})
$$

其中，$J(x^{(k)})$ 是 $F(x^{(k)})$ 的雅可比矩阵，$F(x^{(k)}) = (f_1(x^{(k)}), f_2(x^{(k)}), \cdots, f_n(x^{(k)}))^T$ 是方程组在点 $x^{(k)}$ 处的取值。

在 R 语言中，我们可以先定义一个函数来计算雅可比矩阵 $J(x)$ 和方程组 $F(x)$ 的值，再使用 while 循环进行迭代求解。具体代码如下：

# 定义方程组
f <- function(x) {
  res <- c(x[1]^2 + x[2]^2 - 1, x[1]*x[2] - 1)
  return(res)
}

# 定义雅可比矩阵
J <- function(x) {
  res <- matrix(c(2*x[1], 2*x[2], x[2], x[1]), nrow = 2, ncol = 2, byrow = TRUE)
  return(res)
}

# 初始解
x0 <- c(1, 1)

# 迭代求解
tol <- 1e-8  # 精度要求
maxiter <- 1000  # 最大迭代次数
k <- 0
while (k <= maxiter) {
  Fk <- f(x0)
  Jinv <- solve(J(x0))
  x1 <- x0 - Jinv %*% Fk
  if (max(abs(x1 - x0)) < tol) break
  x0 <- x1
  k <- k + 1
}
cat("解为：", x1, "\n")