介值定理与二分法

一元一次方程 $ax + b = 0$ 可以直接写出它的解 $x = -b/a$。一元二次方程 $ax^2 + bx + c = 0$，也可以通过求根公式求解 $x = \frac{-b ± \sqrt{b^2 - 4ac}}{2a}$。一元 $n$ 次的多项式，只要它能写成 $\Pi_{i = 1}^n (x + k_i)$，我们也可以直接得到它的解。但对于那些包含了诸如三角、指数、对数等类型函数的超越方程，比如: $$ \sin(x) + x = 0 \qquad e^x + x = 0 \qquad \log_{10}(x) + x^2 + 1 = 0 $$ 这类方程我们通常写不出解析解，但又普遍存在，还想求解，该怎么办呢? 我们可以通过各类数值计算方法来逼近可能的解，虽然不是完全正确，但可以保证十分精确。

数值求解一般都是采用迭代的方式不断接近解，二分法就是其中一种常用的方法。它有一个重要的特性，就是一定收敛。假设计算机有无限精度，在满足二分法的使用条件的情况下，一直迭代下去，可以无限接近真实的解，也就是误差趋于 0。

1. 二分法原理

如果函数 $f(x)$ 在闭区间 $[a, b]$ 上连续，并且 $f(a) \neq f(b)$，那么对于介于 $f(a), f(b)$ 之间的任意一个实数 $C$，即 $\min(f(a), f(b)) < C < \max(f(a), f(b))$，一定存在一个数 $c \in (a,b)$，使得 $f(c) = C$。这就是著名的介值定理，其几何含义就是，对于一个连续函数，只要起点和终点的值不同，那么它必然经过起点和终点之间的所有值。

介值定理有一个重要的推论。对于在闭区间 $[a, b]$ 上连续的函数 $f(x)$，如果 $f(a) \cdot f(b) < 0$，那么在 $(a,b)$ 中一定存在一点 $c$，使得 $f(c) = 0$。该推论又被称为零点定理。相信很多读者看到这里，就已经能体会到二分法的精髓: 不断的取区间 $[a, b]$ 的中点，得到更小的闭区间，直到闭区间小的不能再小时，就得到了一个数值解。该算法的大致过程如下:

取区间 $[a_i, b_i]$ 的中点 $c$，计算 $f(c)$。若恰好 $f(c) = 0$，那么 $c$ 就是解。
否则，$f(c) \cdot f(a)$ 与 $f(c) \cdot f(b)$ 中一定只有一个小于 0。选择小于 0 的那项，构造新的闭区间 $[a_{i+1}, b_{i+i}]$。
如果闭区间的长度小于一定阈值时，我们就认为 $c$ 就是一个满足精度要求的数值解。否则重复 1 继续迭代。

上述过程中，$a_i, b_i$ 是每次迭代时闭区间的起点和终点，迭代初始时有 $a_0 = a, b_0 = b$。记闭区间 $I_n = [a_n, b_n]$。因为迭代过程中的闭区间具有关系 $I_1 \supset I_2 \supset I_3 \supset \cdots \supset I_n \supset \cdots$。由于每次都取区间的中点构造新的闭区间，所以每次迭代区间长度都会减半，有: $$ \|I_{n+1}\| = |b_{n+1} - a_{n+1}| = \frac{1}{2} |b_n - a_n| = \frac{1}{2^n} | b - a | $$ 因为 $\lim_{n \to \infty} {\frac{1}{2^n}} = 0$，所以 $\lim_{n \to \infty} \|I_{n+1}\| = 0$。根据闭区间套定理，上述迭代过程一定会收敛到一个实数上。记该实数为 $\xi \in (a, b)$。根据连续函数的定义，当 $x \to \xi$ 时，函数的极限 $\lim_{x \to \xi} f(x)$ 也存在，并且 $\lim{x \to \xi} f(x) = f(\xi)$。由于迭代过程中 $a_n ≤ \xi ≤ b_n$，所以 $\lim_{n \to \infty} {f(a_n)}$ 是函数 $f(x)$ 关于 $\xi$ 左极限，$\lim_{n \to \infty} {f(b_n)}$ 是右极限。根据连续函数左右极限相等的性质，有 $\lim_{n \to \infty} {f(a_n)} = \lim_{n \to \infty} {f(b_n)} = f(\xi)$。迭代过程中，我们始终保证 $f(a_n) \cdot f(b_n) \lt 0$，所以只有 $f(\xi) = 0$。

所以，从数学原理上二分法就是绝对稳定的，只要初始区间合适，那么一直迭代下去就一定能找到一个解。理想情况下，一直迭代就一直接近真值。但由于计算机中浮点数的精度有限，所以当误差小到一定程度之后，再继续迭代就没有意义了。

2. 二分法的实现

虽然二分法的数学原理很直观，但实现起来还是有一些细节需要注意的。如下面的代码片段所示，模板函数 Bisection 共有 5 个参数。其中 f, a, b 表示要在闭区间 $[a, b]$ 上求方程 $f(x) = 0$ 的根。为了保证算法一定能够退出，我们通过函数 max_iter 指定了算法的最大迭代次数。最后 tol 是一个很小的数，当迭代闭区间长度小于该值时，我们认为算法收敛到了一个满足精度要求的数值解上了。

        template <typename DataType>
        DataType Bisection(std::function<DataType(DataType)> f, DataType a, DataType b, int max_iter = 100, DataType tol = SMALL_VALUE)

一开始，我们先检查一下闭区间两端的函数值 $f(a), f(b)$，如果恰好为 0，就直接返回。

        {
            DataType A = f(a); DataType B = f(b);
            if (0 == A) return a;
            if (0 == B) return b;

然后检查输入的数据是否满足二分法的要求。为了保持代码的简洁，我们要求用户在调用 Bisection 的时候，保证 $a < b$。判定 $f(a) < f(b)$ 时，我们特意用符号函数 $\text{Sign}(x) = \begin{cases} 1 & x > 0 \\ 0 & x = 0 \\ -1 & x < 0 \end{cases}$ 封装了一下。这样做主要是为了防止浮点数乘法溢出的问题。在迭代过程中，$f(x)$ 接近 0 的时候，浮点数的截止误差可能导致 $f(a) * f(b)$ 输出 0，这会导致迭代逻辑失效。

            assert(a < b);
            assert(Sign(A) * Sign(B) < 0);

接下来就是二分法的迭代循环。我们在 for 循环一开始，先计算了中点 re，这里我们特别采用了 $a + \frac{b - a}{2}$ 而非 $\frac{a + b}{2}$。假如输入参数 $a, b$ 是两个很大的数，那么 $a + b$ 是很有可能溢出的。但 $b - a$ 总是能控制在一个比较的范围内，而且随着迭代的过程在不断接近于 0。

            DataType re = b;
            for (int i = 0; i < max_iter; ++i) {
                DataType half_length = 0.5 * std::abs(b - a);
                re = a + half_length;

然后计算函数值 f(re) 用变量 B 保存，如果碰巧该值为 0，那么 re 就是我们要找的根。否则判定区间长度是否小于阈值。

                B = f(re);
                if (0 == B || half_length < tol)
                    return re;

如果没有满足算法的终止条件，我们就需要更新这些局部变量，缩小闭区间。

                if (Sign(A) * Sign(B) > 0) {
                    A = B; a = re;
                } else {
                    b = re;
                }
            }

最后，超出了最大迭代次数都没有满足收敛条件，我们也将当期的结果作为一个不那么精确的解返回出来。

             return re;
        }

完整的函数实现参见 Bisection。

3. 二分法的误差分析

我们评价一个数值算法的好坏，无外乎两个方面，精度和效率。很多时候还要在两者之间作一些取舍。可以证明，对于在闭区间 $[a,b]$ 上连续的函数 $f$ 满足 $f(a) \cdot f(b) < 0$ 时，通过二分法求解会以 $O(\frac{1}{2^n})$ 的速率收敛到解 $p$ 上。

由于我们在每次迭代时，都取闭区间 $[a_n, b_n]$ 的中点 $p_n = a_n + \frac{b_n - a_n}{2}$。所以有 $b_n - a_n = \frac{1}{2^{n-1}}(b - a)$，并且 $a_n < p < b_n$。所以存在如下不等式关系:

$$ \begin{equation}\label{f1} |p_n - p| = | a_n + \frac{b_n - a_n}{2} - p | < \frac{1}{2}(b_n - a_n) = \frac{b- a}{2^n} \end{equation} $$

所以，当 $n \to \infty$ 时，有 $p_n = p + O(\frac{1}{2^n})$。根据不等式($\ref{f1}$)，我们可以估计出经过多少次迭代就能达到希望的精度。比如，初始闭区间为 $[1, 3]$ 希望精度到达 $10^{-6}$。根据不等式有:

$$ |p_n - p| < \frac{1}{2^n} < 10^{-6} \Rightarrow -n \log_{10}2 < -6 \Rightarrow n > 6 / \log_{10} 2 \approx 19.931568569324174 $$

所以 20 次迭代之后就一定能达到$10^{-6}$ 的精度。这是一个保守的估计，实际可能也不需要 20 次。

4. 完

二分法的逻辑很简单而且绝对可靠，只要初始区间两端的函数值异号，持续迭代下去就一定能找到根。但缺点也很明显，收敛速度比较慢。而且如果区间中有多个根，它只能找到其中一个。下一节我们来讨论一种收敛速度快的牛顿法。