多项式方程求根
我们在第一部分中介绍的二分法、 牛顿法、Brent 方法都是实数运算。 很多时候多项式方程是有复数根的，那些算法好像都不能处理。1956 年 Müller 提出了一种算法，能够自然的处理复数， 配合上一节提到的收缩(deflation)约化技术，可以给出多项式方程的所有根。
1. 实数域中的 Müller 算法
Müller 算法与我们熟悉的割线法类似。割线法通过最后两次迭代点构造一条直线， 然后用该直线的根构造新的迭代点。相比之下，Müller 算法使用最后三次迭代点构造一个抛物线，然后用该抛物线的一个根构造新的迭代点。
假设我们现在有三个迭代点 \(x_0, x_1, x_2\)，那么我们可以写出一个二次多项式:
$$ \begin{equation}\label{f1} P(x) = a (x - x_2)^2 + b(x - x_2) + c \end{equation} $$ 如果该多项式经过 \((x_0, f(x_0)), (x_1, f(x_1)), (x_2, f(x_2))\)，那么我么可以写出如下方程组，进而解出多项式的系数 \(a,b,c\): $$ \begin{cases} f(x_0) = a(x_0 - x_2)^2 + b(x_0 - x_2) + c \\ f(x_1) = a(x_1 - x_2)^2 + b(x_1 - x_2) + c \\ f(x_2) = a \cdot 0^2 + b \cdot 0 + c = c \\ \end{cases} \Longrightarrow \begin{cases} c = f(x_2) \\ b = \frac{(x_0 - x_2)^2[f(x_1) - f(x_2)] - (x_1 - x_2)^2[f(x_0) - f(x_2)]}{(x_0 - x_2)(x_1 - x_2)(x_0 - x_1)} \\ a = \frac{(x_1 - x_2)^2[f(x_0) - f(x_2)] - (x_0 - x_2)^2[f(x_1) - f(x_2)]}{(x_0 - x_2)(x_1 - x_2)(x_0 - x_1)} \end{cases} $$我看很多实现，都会定义如下的四个中间变量:
$$ h_0 = x_1 - x_0, \qquad h_1 = x_2 - x_1, \qquad \delta_0 = \frac{f(x_1) - f(x_0)}{h_0}, \qquad \delta_1 = \frac{f(x_2) - f(x_1)}{h_1} $$这样，就可以比较简洁的写出系数 \(a,b,c\)，并且可以省去很多重复的计算:
$$ a = \frac{\delta_1 - \delta_0}{h_1 + h_0} \qquad b = ah_1 + \delta_1 \qquad c = f(x_2) $$有了系数，我们就可以直接通过求根公式写出下一个迭代点。这里出于舍入误差的考虑，我们采用如下的形式:
$$ \begin{equation}\label{f2} x_3 = x_2 - \frac{2c}{b+\text{sign}(b)\sqrt{b^2 - 4ac}} \end{equation} $$参考割线法的实现，我们只要传入三个初始点，并用上式替代割线法的迭代公式，就可以得到一个实数域的 Müller 算法。 参见例程。据说它的收敛阶数是 1.839， 这里我们不再详细分析它的收敛特性了。
2. 复数域的 Müller 算法
递推公式 \(\ref{f2}\) 还是实数域的公式，符号函数 \(\text{sign}\) 只对实数有意义，对于复数而言是不存在正负号的。我看过不少教科书，他们都只提到 Müller 算法天然支持复数根的求解，
但很少有教材配的伪代码说明如何处理复数。后来我悟了，在复数域上除了这个符号函数之外，其它的计算跟实数域的完全一样，只要代码能支持复数运算就行了。
我们用 C++ 的标准库 std::complex
实现了一个复数版本的MullerRoot。
下面左侧是实数域下递推公式的实现，右侧是复数域的。在右侧，我们用 Complex 对标准库 std::complex 简单封装了一下，所有的计算都是在复数域上进行的，
所以目标函数也需要支持复数。
|
|
虽然复数中不存在正负号，但是它可以看做一个二维的矢量，可以计算模长。所以标准库中绝对值函数 std::abs 对复数的计算实际在求其模长。
考虑到舍入误差的问题，我们先计算出 \(b ± \sqrt{b^2 - 4ac}\)，然后选择其中模长较大的那个计算下一个迭代点。这样一来 Müller 算法一次只能输出一个根。
我们还需要配合收缩(deflation)约化技术才能求出多项式方程的所有根。
3. 用 Müller 算法求所有根
要求出所有根，无非是每次用 Müller 算法解出一个根，接着通过综合除法约简一阶，直到多项式降为一次，最后直接写出一次多项式的根。 但我们在一元多项式中就已经提到过，由于每一次求根和约简的过程都是一个近似计算的过程， 存在累积误差的问题，尤其在方程存在重根的情况时更严重。所以人们通常还会把 Müller 算法解出的每一个根，再通过牛顿法细化一下，来削弱累积误差的影响。
下面是我们在复数域上实现的一个求所有根的方法。在左侧，我们定义了一对局部变量。其中 x0, x1, x2 是随便写的三个数，用于每一次调用 Müller 算法的初值。 初值对算法的收敛过程是有一定影响的，但是一般情况下，我们并没有太多的信息，给出一个比较好的初值。工程上很多实现都是用的三个随机数， 也有广泛采样多个数之后，选取多项式函数值最接近 0 的三个数，作为初值。
|
|
上面左侧定义的 P,Q 和指针 p_ptr, q_ptr 用于右侧的综合除法，临时存放约简多项式的。上面右侧是我们的主要循环，这里我们先用 MullerRoot 计算一个粗解， 再用牛顿法 NewtonRaphson 细化削弱累积误差。然后将细化后的解记录到输出数组 roots 中。接着就通过综合除法 SyntheticDivide 约简多项式。 while 循环到一次多项式之后就结束了，在下面的代码片段中，我们可以直接从一次多项式的系数中写出最后一个解。完整的实现 参见例程
if (1 == p_ptr->Degree()) {
Complex r = -p_ptr->At(0) / p_ptr->At(1);
r = NewtonRaphson(r);
roots.push_back(r);
}
}
4. 完
Müller 算法通过二次插值公式提供了求解复数根的方法，在复数域上每次只能求出一个根，但可以通过综合除法约简，逐个计算出所有根。