中值定理与牛顿法

上一节中，我们根据零点定理，介绍了二分法迭代求解一元方程 $f(x) = 0$ 的方法。只要满足二分法的使用条件，就一定能收敛到一个根上。由于该方法每次都是折半查找，所以每次迭代的误差 $e_{n+1}$ 都不超过上一轮误差 $e_{n}$ 的一半，即 $e_{n+1} ≤ \frac{1}{2}e_{n}$，因此我们常说二分法是线性收敛的。

相比之下，本文将要介绍的牛顿法一般情况下是二阶收敛的，即 $e_{n+1} \approx C e_n^2$，其中 $C$ 是一个常数。虽然它的收敛速度很快，但是迭代初值如果偏离真实的解太远，就很可能不收敛。

1. 牛顿法的基本原理

很多教科书上，都是通过泰勒展开来引出牛顿法的。如果函数 $f(x)$ 在闭区间 $[a, b]$ 上二阶连续，即 $f \in C^2[a,b]$，令 $x_n \in [a, b]$ 是一个接近真实解 $p$ 的点，并且 $f'(x_n) \neq 0$。那么我们在 $x_n$ 处对 $f(x)$ 泰勒展开，有

$$ \begin{equation}\label{f1} f(x) = f(x_n) + (x - x_n)f'(x_n) + \frac{(x - x_n)^2}{2}f^{''}(\xi) \end{equation} $$

最后的余项中 $\xi$ 为 $x_n$ 和 $p$ 之间的数值。因为 $f'(x_n) \neq 0$，上式可以改写为

$$ \begin{equation}\label{f2} \frac{f(x)}{f'(x_n)} = \frac{f(x_n)}{f'(x_n)} + (x - x_n) + (x - x_n)^2\frac{f^{''}(\xi)}{2f'(x_n)} \end{equation} $$

由于 $x_n$ 是一个接近真实解 $p$ 的点，所以 $|x_n - p|$ 是一个很小的数，$|x_n - p|^2$ 就更小了。因此，我们直接抛弃掉最后的余项，又 $f(p) = 0$，有:

$$ \frac{f(p)}{f'(x_n)} = 0 \approx \frac{f(x_n)}{f'(x_n)} + p - x_n \Longrightarrow p \approx x_n - \frac{f(x_n)}{f'(x_n)} $$

由此，我们可以写出迭代关系式：

$$ \begin{equation}\label{f3} x_{n+1} = x_n - \frac{f(x_n)}{f'(x_n)} \end{equation} $$

根据式($\ref{f1}$) 和 ($\ref{f3}$)，可以写出误差的迭代公式 $e_{n+1} = |p - x_{n+1}| = (p - x_n)^2 \left|\frac{f^{''}(\xi)}{2 f'(x_n)} \right| = C \cdot e_{n}^2 $。可见牛顿法是二阶收敛的，相比于线性收敛的二分法来说，它的效率会更高。我们知道式($\ref{f1}$) 中的泰勒展开式只有在 $x_n$ 的邻域内才有意义。如果 $x_n$ 与 $p$ 相差很远，算法就不一定收敛。下面我们先根据迭代式 ($\ref{f3}$) 实现一个牛顿法，再形式化的证明该算法的收敛性 —— 只要初值选的好就一定能收敛。

2. 牛顿法的实现

虽然是牛顿最早想出这个方法的，但我们在各种教科书上看到的迭代公式和使用方式，更接近于拉菲森(Joseph Rephson)的形式，相比之下牛顿的方法更繁琐。据说 Rephson 提出该方法的时候跟牛顿一点关系没有，因为是他先公开这种方法的，只是牛顿名气更大一些而已。所以这里我们用 NewtonRaphson 来命名函数，如下面的代码片段所示。

该函数有 5 个参数，其中 f 是要求解的方程 $f(x) = 0$，df 则是函数 $f(x)$ 的一阶导。x0 是迭代初值。为保证算法一定能够结束，我们用 max_iter 限制最大的迭代次数， tol 判定算法是否收敛。

        template <typename DataType>
        DataType NewtonRaphson(std::function<DataType(DataType)> f,
                               std::function<DataType(DataType)> df,
                               DataType x0, int max_iter = 100, DataType tol = SMALL_VALUE)

我们在一个 for 循环中进行迭代。如下代码片段所示，其中局部变量 re 记录当前的迭代值，初始为 x0。在 for 循环一开始，我们先计算 $f(re)$，若为 0 则 re 就是根直接返回。

        {
            DataType re = x0;
            for (int i = 0; i < max_iter; ++i) {
                DataType y = f(re);
                if (0 == y)
                    return re;

迭代过程中，我们调用函数 df 获得 $f(x)$ 在 re 处的一阶导数。为防止除 0 异常，我们需要一个断言来检查它。

                DataType dydx = df(re);
                assert(0 != dydx);

接下来，代入迭代关系式 ($\ref{f3}$) 计算下一轮的迭代值 $x_{n+1}$。如果 $|\text{re} - \text{x0}|$ 很小，说明算法收敛了，我们将其作为真实解的一个近似返回。

                re = x0 - y / dydx;
                if (std::abs(re - x0) < tol)
                    return re;

最后更新迭代初值 x0。如果 for 循环迭代超过了 max_iter 轮，就将最后一次的迭代值返回。

                x0 = re;
            }
            return re;
        }

实际上，牛顿法会以极快的速度收敛或者发散。所以不需要很多次迭代，我们就可以判定初值是否合理，算法是否收敛。工程上，人们还会加上发散判定提前终止，比如连续出现多次 $f(x_{n+1}) > f(x_n)$ 的情况，很可能迭代过程正在原理真实的解。

3. 牛顿法的收敛原理

在讨论牛顿法的收敛原理之前，我们先来看一下拉格朗日中值定理。如果函数 $f(x)$ 在闭区间 $[a,b]$ 上连续，并且在开区间 $(a,b)$ 上可导，那么必定存在至少一个点 $\xi \in (a, b)$，使得下式成立:

$$ f'(\xi) = \frac{f(b) - f(a)}{b - a} $$

我们的问题是要求解 $f(x) = 0$ 的根 $p$。假设 $b$ 是真实的根，即 $b = p$，以及 $a$ 为当前迭代值 $x_n$。代入上式有:

$$ 0 - f(x_n) = f'(\xi)(p - x_n) \Longrightarrow p = x_n - \frac{f(x_n)}{f'(\xi)} $$

上式与牛顿法的迭代公式($\ref{f3}$)很像，实际就是用 $f'(x_n)$ 近似代替中值定理中未知的 $f'(\xi)$，如此就把一个存在性定理变成了一个构造算法。下面我们来证明:

只要初值选的好就一定能收敛。设函数 $f(x)$ 在闭区间 $[a,b]$ 上二阶连续。如果 $p \in (a, b)$ 满足 $f(p) = 0$ 并且 $f'(p) \neq 0$。那么一定存在 $\delta > 0$，使得任何 $x_0 \in [p - \delta, p+ \delta]$ 内的初值，迭代公式($\ref{f3}$)产生的数列 $\{x_{n}\}_{n = 1}^{\infty}$ 收敛到 $p$。

证:令 $g(x) = x - \frac{f(x)}{f'(x)}$。由于$f(x)$ 是二阶连续的，并且 $f'(p) \neq 0$，所以 $g(x)$ 是连续的，并且其导函数 $g'(x)$ 也是连续的:

$$ g'(x) = 1 - \frac{f'(x)f'(x) - f(x)f''(x)}{[f'(x)]^2} = \frac{f(x)f''(x)}{[f'(x)]^2} $$

因为 $f(p) = 0$，有

$$ g(p) = p - \frac{f(p)}{f'(p)} = p \qquad g'(p) = \frac{f(p)f''(p)}{[f'(p)]^2} = 0 $$

由于$g'(x)$ 是连续的，并且 $g'(p) = 0$。对于 $0 < k < 1$，根据连续函数的性质，一定存在 $\delta > 0$ 使得所有 $x \in [p - \delta, p+ \delta]$，有

$$ |g'(x)| ≤ k $$

根据中值定理，在 $x$ 和 $p$ 之间一定存在一个数 $\xi$，满足 $|g(x) - g(p)| = |g'(\xi)||x - p|$。结合上式有

$$ |g(x) - p| = |g(x) - g(p)| < k |x - p| < |x - p| $$

结合牛顿法的迭代关系式有:

$$ | g(x_n) - p | = | x_{n+1} - p | < | x_{n} - p | $$

所以选择初值 $x_0 \in [p - \delta, p + \delta]$ 时，牛顿法一定收敛。□

上述定理和证明过程，只是证明了存在一个数 $\delta > 0$ 保证算法一定收敛。但并不能给出找到这个 $\delta$ 的方法。

4. 完

牛顿法的逻辑甚至比二分法更简单，只要一直调用迭代关系式就可以了。它是一种二阶收敛的算法，只要初值选的合适就能很快找到解。