| RMVL 2.4.0 Robotic Manipulation and Vision Library

载入中...
搜索中...
未找到
非线性方程（组）数值解与迭代法

涉及 不动点迭代 与 牛顿迭代 两种非线性方程组数值解的求法

作者
  赵曦

日期
  2024/01/12

版本
  1.0

上一篇教程：非线性最小二乘
下一篇教程：常微分方程（组）数值解与 Runge-Kutta 算法

## 1 不动点迭代

对于一个函数 \(y=f(x)\)，希望求其零点，即 \(f(x)=0\)，可以将其改写成等价的

\[x=\varphi(x)\tag{1-1}\]

因此，选定一个初值 \(x_0\)，对 \(\text{(1-1)}\) 构建迭代公式

\[x_{k+1}=\varphi(x_k)\quad(k=0,1,\cdots)\tag{1-2}\]

如果该迭代最终收敛于精确解，则这种方法被称为不动点迭代法，需要特别注意的是，选择的迭代函数 \(\varphi(x)\)不能发散，也就是需要满足一定的收敛条件，这里给出局部收敛的判断方法。令 \((x)=0\)的精确解为 \(x^*\)，若满足

\[\left|\varphi'(x^*)\right|<1\tag{1-3}\]

则迭代公式 \(x=\varphi(x)\) 在 \(x^*\) 处局部收敛。

示例

用不动点迭代法求解

\[x^3+4x^2-10=0\]

解答

上式可以改写为

\[x=\frac12\sqrt{10-x^3}\]

，因此可以构建如下迭代方程

\[x_{k+1}=\frac12\sqrt{10-x_k^3}\]

## 2 Newton 迭代

对于方程 \(f(x)=0\)，其精确解为 \(x^*\)，对于某一次迭代结果 \(x_k\)，将 \(f(x^*)\)在 \(x_k\)处展开，有

\[0=f(x^*)=f(x_k)+f'(x_k)(x^*-x_k)+\frac12f''(x_k)(x^*-x_k)^2+\cdots\tag{2-1}\]

即

\[\begin{align} x^*&=x_k-\frac{f(x_k)}{f'(x_k)}-\frac{f''(x_k)}{2f'(x_k)}(x^*-x_k)^2-\cdots\\ &\approx x_k-\frac{f(x_k)}{f'(x_k)} \end{align}\tag{2-2}\]

因此可以建立迭代方程

\[f(x_{k+1})=x_k-\frac{f(x_k)}{f'(x_k)}\tag{2-3}\]

上式就是 Newton 迭代的迭代公式，对于存在连续二阶导数的函数 \(f(x)\)，Newton 迭代在用于求解单根（一重零点）的时候，至少具有 2 阶收敛性。在用于求解重根（ \(m\)重零点）的时候，Newton 迭代仅具备 1 阶收敛性，需要改写为

\[f(x_{k+1})=x_k-m\frac{f(x_k)}{f'(x_k)}\tag{2-3}\]

方可实现 2 阶收敛。

## 3 Newton 迭代的简化

对于计算机算法而言，普通的 Newton 迭代需要求导 \(f'(x_k)\)的操作，这一步可使用向后差商来代替

\[f'(x_k)\approx\frac{f(x_k-x_{k-1})}{x_k-x_{k-1}}=f\left[x_k,x_{k-1}\right]\tag{3-1}\]

此时的 Newton 迭代被称为离散 Newton 迭代，也称为弦截法，表示为如下形式，可参考 rm::NonlinearSolver

\[f(x_{k+1})=x_k-\frac{f(x_k)}{f(x_k-x_{k-1})}(x_k-x_{k-1})\tag{3-2}\]

## 4 使用示例

- C++

  #include <rmvl/algorithm.hpp>
  #include <cstdio>
  int main() {
  auto foo = rm::NonlinearSolver([](double x) { return x * x - 4; });
  printf("x1 = %.2f, x2 = %.2f\n", foo(-5), foo(5));
  // x1 = -2.00, x2 = 2.00
  }
  rm::NonlinearSolver
  非线性方程求解器
  定义 numcal.hpp:152

- Python

  import rm
  foo = rm.NonlinearSolver(lambda x: x * x - 4)
  print("x1 = {:.2f}, x2 = {:.2f}".format(foo(-5), foo(5)))
  # x1 = -2.00, x2 = 2.00