非线性方程的数值解法

与线性方程组不同，非线性方程 $f(x)=0$ 的求解通常没有直接的解析方法，因此必须依赖迭代数值技术。本章将介绍两种基础且重要的迭代方法：二分法和牛顿迭代法。

不动点迭代与Banach不动点定理

许多求解非线性方程 $f(x)=0$ 的迭代法，都可以归结为一种称为不动点迭代的一般形式。

不动点迭代

不动点迭代法的思想是，首先将方程 $f(x)=0$ 转化为等价的形式：

$$x = g(x) $$

如果 $x^*$ 是这个方程的解，即 $x^* = g(x^*)$，那么我们称 $x^*$ 是函数 $g(x)$ 的一个不动点。

基于这个形式，我们可以构造一个简单的迭代序列：

$$x_{k+1} = g(x_k) $$

从一个初始猜测值 $x_0$ 开始，如果这个序列收敛，即 $\lim_{k \to \infty} x_k = x^*$，那么我们就找到了方程的一个解。

例如，求解 $f(x) = x^2 - 3 = 0$，可以有多种等价形式：

  1. $x = x - (x^2 - 3)$
  2. $x = 3/x$
  3. $x = (x + 3/x) / 2$ （这是求解 $\sqrt{3}$ 的牛顿法形式）

不同的 $g(x)$ 形式会导致迭代序列具有截然不同的收敛性质。

不动点定理 (存在性)

该定理保证了不动点的存在性。

定理：如果函数 $g(x)$ 在闭区间 $[a, b]$ 上连续，并且对于任意 $x \in [a, b]$，都有 $g(x) \in [a, b]$（即 $g(x)$ 将区间 $[a, b]$ 映射到其自身），那么在 $[a, b]$ 内至少存在一个不动点。

这个定理只保证存在性，但没有说明不动点是否唯一，也没有保证迭代序列一定会收敛。

Banach不动点定理 (收敛性)

Banach不动点定理给出了迭代收敛的充分条件，是迭代法收敛性分析的基石。

定义：压缩映射 (Contraction Mapping)
如果存在一个常数 $L$ 满足 $0 \le L < 1$，使得对于定义域内的任意两个不同的点 $x$ 和 $y$，都有：

$$|g(x) - g(y)| \le L|x - y| $$

则称 $g(x)$ 是一个压缩映射。从几何上看，这意味着函数 $g(x)$ 会“拉近”任意两点之间的距离。

如果 $g(x)$ 可导，一个更易于判断的充分条件是，在区间 $[a, b]$ 上始终满足 $|g'(x)| \le L < 1$。

Banach不动点定理：如果 $g(x)$ 是闭区间 $[a, b]$ 上的一个压缩映射，且其值域包含于 $[a, b]$，那么：

  1. 在 $[a, b]$ 内存在唯一的不动点 $x^*$。
  2. 对于任意初始值 $x_0 \in [a, b]$，由 $x_{k+1} = g(x_k)$ 生成的迭代序列都将收敛于 $x^*$。

误差估计：

$$|x_k - x^*| \le \frac{L^k}{1-L} |x_1 - x_0| $$

这表明不动点迭代法至少是线性收敛的，收敛速度取决于压缩常数 $L$ 的大小。$L$ 越小，收敛越快。

这个定理为我们分析和设计迭代算法提供了强大的理论工具。例如，在选择 $f(x)=0$ 的等价形式 $x=g(x)$ 时，我们应尽量选择一个能满足压缩映射条件的 $g(x)$，以保证迭代的收敛性。

唯一性证明 (反证法)：
假设在 $[a, b]$ 内存在两个不同的不动点，记为 $x^*$ 和 $y^*$。根据不动点的定义，我们有：

$$x^* = g(x^*) \quad \text{且} \quad y^* = g(y^*) $$

由于 $g(x)$ 是一个压缩映射，它满足：

$$|g(x^*) - g(y^*)| \le L |x^* - y^*| $$

将 $g(x^*) = x^*$ 和 $g(y^*) = y^*$ 代入上式，得到：

$$|x^* - y^*| \le L |x^* - y^*| $$

因为我们假设 $x^*$ 和 $y^*$ 是不同的，所以 $|x^* - y^*| > 0$。因此，我们可以将不等式两边同时除以 $|x^* - y^*|$，得到：

$$1 \le L $$

但这与压缩映射的条件 $L < 1$ 相矛盾。因此，最初的假设（存在两个不同的不动点）不成立。
这就证明了在满足Banach不动点定理的条件下，不动点是唯一的。

二分法 (Bisection Method)

二分法是最简单、最稳健的求根方法之一。

基本思想

二分法的理论基础是介值定理：如果一个连续函数 $f(x)$ 在区间 $[a, b]$ 的两个端点上的函数值 $f(a)$ 和 $f(b)$ 异号，即 $f(a) \cdot f(b) < 0$，那么在区间 $(a, b)$ 内至少存在一个根。

二分法通过反复将包含根的区间一分为二，并选择保留包含根的那个子区间，从而逐步逼近根。

算法步骤

  1. 初始化：找到一个区间 $[a_0, b_0]$，使得 $f(a_0) \cdot f(b_0) < 0$。设定一个误差容限 $\epsilon$。
  2. 迭代：对于 $k = 0, 1, 2, \ldots$
    a. 计算中点：$x_k = \frac{a_k + b_k}{2}$。
    b. 判断停止条件：如果 $|b_k - a_k| < \epsilon$ 或者 $|f(x_k)| < \epsilon$，则停止迭代，返回 $x_k$ 作为根的近似值。
    c. 更新区间：
    - 如果 $f(a_k) \cdot f(x_k) < 0$，则根在左半部分，令 $[a_{k+1}, b_{k+1}] = [a_k, x_k]$。
    - 否则，根在右半部分，令 $[a_{k+1}, b_{k+1}] = [x_k, b_k]$。

收敛性与误差分析

  • 收敛性：二分法是无条件收敛的。只要初始区间 $[a_0, b_0]$ 包含根，该方法就一定能找到根的近似值。
  • 误差分析：经过 $k$ 次迭代后，区间的长度为：

    $$|b_k - a_k| = \frac{|b_0 - a_0|}{2^k} $$

    设 $x^*$ 是真实解，第 $k$ 次迭代的近似解 $x_k$ 的误差为 $|x_k - x^*|$，其误差界为：

    $$|x_k - x^*| \le \frac{|b_0 - a_0|}{2^k} $$

  • 收敛阶：二分法是线性收敛的，收敛阶为1。误差比率为：

    $$\lim_{k \to \infty} \frac{|x_{k+1} - x^*|}{|x_k - x^*|} = \frac{1}{2} $$

    这意味着每迭代一次，误差大约减半。

优缺点

  • 优点：算法简单，收敛性有保证，对函数要求低（仅需连续）。
  • 缺点：收敛速度慢，无法处理偶数重根（因为根附近函数值不变号），也无法推广到多维情况。

牛顿迭代法 (Newton’s Method)

牛顿法，也称为牛顿-拉夫逊方法，是一种收敛速度极快的迭代方法。

基本思想

牛顿法的思想是利用函数在某一点的泰勒展开。将函数 $f(x)$ 在点 $x_k$ 附近进行一阶泰勒展开：

$$f(x) \approx f(x_k) + f'(x_k)(x - x_k) $$

为了找到 $f(x)=0$ 的根，我们令上式等于0，求解 $x$：

$$f(x_k) + f'(x_k)(x - x_k) = 0 \implies x = x_k - \frac{f(x_k)}{f'(x_k)} $$

我们将这个新解作为下一次迭代的近似值 $x_{k+1}$。几何上，这相当于在点 $(x_k, f(x_k))$ 处作函数的切线，并将切线与x轴的交点作为新的近似根。

迭代公式

$$x_{k+1} = x_k - \frac{f(x_k)}{f'(x_k)} $$

收敛条件

牛顿法的收敛性与初始值的选择密切相关。

  • 局部收敛定理：设 $x^*$ 是 $f(x)=0$ 的单根（即 $f'(x^*) \neq 0$），且 $f(x)$ 在 $x^*$ 附近有二阶连续导数。则存在一个包含 $x^*$ 的邻域 $\delta$，对于任意初始值 $x_0 \in \delta$，牛顿法产生的序列 $\{x_k\}$ 都收敛于 $x^*$。

  • 收敛失败的情况：

    1. $f'(x_k) \approx 0$：切线接近水平，导致 $x_{k+1}$ 离根非常远。
    2. 初始值选择不当：可能导致迭代不收敛或收敛到错误的根。
    3. 函数存在局部极值点。

收敛阶与误差分析

  • 收敛阶：当牛顿法收敛于单根时，它是二次收敛的，收敛阶为2。
  • 误差分析：设 $e_k = x_k - x^*$ 为第 $k$ 次迭代的误差。可以证明：

    $$e_{k+1} \approx \frac{f''(x^*)}{2f'(x^*)} e_k^2 $$

    这意味着每一次迭代后，误差的有效数字位数大约会翻倍，收敛速度非常快。
  • 重根情况：如果 $x^*$ 是 $m$ 重根，牛顿法会退化为线性收敛，收敛比为 $1 - \frac{1}{m}$。

4. 多维情况的推广：牛顿法求解非线性方程组

牛顿法可以自然地推广到求解非线性方程组。考虑方程组：

$$F(X) = \mathbf{0} $$

其中 $X = (x_1, x_2, \ldots, x_n)^T$ 是向量，$F = (f_1, f_2, \ldots, f_n)^T$ 是向量函数。

雅可比矩阵 (Jacobian Matrix)

多维情况下的导数由雅可比矩阵 $J(X)$ 代替：

$$J(X) = \begin{pmatrix} \frac{\partial f_1}{\partial x_1} & \frac{\partial f_1}{\partial x_2} & \cdots & \frac{\partial f_1}{\partial x_n} \\ \frac{\partial f_2}{\partial x_1} & \frac{\partial f_2}{\partial x_2} & \cdots & \frac{\partial f_2}{\partial x_n} \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial f_n}{\partial x_1} & \frac{\partial f_n}{\partial x_2} & \cdots & \frac{\partial f_n}{\partial x_n} \end{pmatrix} $$

迭代公式

将 $F(X)$ 在 $X_k$ 附近进行一阶泰勒展开：

$$F(X) \approx F(X_k) + J(X_k)(X - X_k) $$

令 $F(X) = \mathbf{0}$，得到：

$$J(X_k)(X - X_k) = -F(X_k) $$

这可以看作一个关于变量 $\Delta X_k = X - X_k$ 的线性方程组。求解该方程组得到 $\Delta X_k$，然后更新 $X_k$：

$$X_{k+1} = X_k + \Delta X_k = X_k - [J(X_k)]^{-1} F(X_k) $$

在实际计算中，我们通常不直接求雅可比矩阵的逆，而是通过解线性方程组来获得更新量 $\Delta X_k$：

$$J(X_k) \Delta X_k = -F(X_k) $$

算法步骤

  1. 初始化：选择一个初始猜测值 $X_0$。
  2. 迭代：对于 $k = 0, 1, 2, \ldots$
    a. 计算 $F(X_k)$ 和雅可比矩阵 $J(X_k)$。
    b. 求解线性方程组 $J(X_k) \Delta X_k = -F(X_k)$ 得到 $\Delta X_k$。
    c. 更新解：$X_{k+1} = X_k + \Delta X_k$。
    d. 检查收敛性，例如 $||\Delta X_k||$ 或 $||F(X_{k+1})||$ 是否足够小。

多维牛顿法同样具有二次收敛的特性，但每次迭代都需要计算雅可比矩阵并求解一个线性方程组，计算成本较高。