【NA】代数方程求根算法

• 前面说过，对于非线性方程 $$f(x)=0$$ 若 $f(x)$ 是多项式函数，那么上述方程被称为代数方程，已经完成的二分法、Aitken迭代法以及牛顿法都适用于代数方程的求解。但针对多项式函数的特殊性，可以寻求一些更加特殊、有效的算法。

代数方程の牛顿法.

• 考虑代数方程 $$f(x)=a_0\cdot x^0+a_1\cdot x^1+\cdots +a_n\cdot x^n=\sum _{i=0}^n{a}_i\cdot x^i=0$$应用牛顿迭代公式$$x_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)}\text{\hspace{1em}(Newton.1)}$$于是问题转化为计算函数 $f(x)$ 在给定点 $x_k$ 的函数值 $f(x_k)$ 和导数值 $f^{\prime }(x_k).$

秦九韶算法.

• 考虑多项式$$f(x)=a_0\cdot x^0+a_1\cdot x^1+\cdots +a_n\cdot x^n$$其中系数均为实数，以一次式 $(x-x_0)$ 除 $f(x)$ 得到商 $p(x)$，显然余数为 $f(x_0)$，即下面的等式成立$$f(x)=f(x_0)+(x-x_0)\cdot p(x)\text{\hspace{1em}(1)}$$因为等式两边次数需要一致，所以 $p(x)$ 是 $n-1$ 次多项式，设为$$p(x)=b_0\cdot x^0+b_1\cdot x^1+\cdots +b_{n-1}\cdot x^{n-1}\text{\hspace{1em}(2)}$$将其代入 $(1)$ 式，可得
  $$\begin{array}{rl}{a}_0\cdot x^0+a_1\cdot x^1+\cdots +a_n\cdot x^n& =f(x_0)+\sum _{i=0}^{n-1}{b}_i\cdot x^{i+1}-x_0\cdot \sum _{i=0}^{n-1}{b}_i\cdot x^i\\ & =f(x_0)-x_0\cdot b_0+\sum _{i=1}^{n-1}(b_{i-1}-x_0\cdot b_i)\cdot x^i+b_{n-1}\cdot x^n\end{array}$$对应次数的系数应该相等，所以有$$a_0=f(x_0)-x_0\cdot b_0$$$$a_i=b_{i-1}-x_0\cdot b_i,i\in [1,n-1]$$$$a_n=b_{n-1}$$将上述三公式变形，得到$$b_{n-1}=a_n$$$$b_i=a_{i+1}+x_0\cdot b_{i+1},i\in [0,n-2]$$$$f(x_0)=a_0+x_0\cdot b_0\text{\hspace{1em}(3)}$$$(3)$ 式提供了快速计算 $f(x_0)$ 的方法，称为秦九韶算法，也叫Horner算法。
• 对 $(1)$ 式求导，得$$f^{\prime }(x)=p(x)$$因此有 $f^{\prime }(x_0)=p(x_0)$，也就是说一阶导数值是多项式 $p(x)$ 在 $x_0$ 的值，再次应用上述秦九韶算法，将 $p(x)$ 表示如下$$p(x)=p(x_0)+(x-x_0)\cdot q(x)$$$$q(x)=c_0\cdot x^0+c_1\cdot x^1+\cdots +c_{n-2}\cdot x^{n-2}\text{\hspace{1em}(4)}$$和上面推导出 $f(x_0)$ 的过程一样，我们可以得到$$c_{n-2}=b_{n-1}$$$$c_i=b_{i+1}+x_0\cdot c_{i+1},i\in [0,n-3]$$$$f^{\prime }(x_0)=b_0+x_0\cdot c_0\text{\hspace{1em}(5)}$$由此得到了一阶导数值 $f^{\prime }(x_0)$，从而可以使用牛顿迭代公式求解方程的近似解。
• 直观来看，秦九韶算法所做的事情就是将形如$$p(x)=\sum _{i=0}^n{a}_i\cdot x^i$$这类多项式函数重写为$$p(x)=(\cdots (a_n\cdot x+a_{n-1})x+\cdots )x+a_0$$最终所需要的乘法此时仅为 $n$ 次，而原本计算单个 $x^n$ 就需要 $n$ 次乘法。秦九韶算法在很多国外教材中被称为霍纳规则Horner’s rule，下图是Anany 182-183中的描述：
  
  
  
• 至此，我们完成了对于代数方程的牛顿迭代公式 $(Newton.1)$ 的具体化。

贝尔斯多夫方法.

• 给定多项式函数 $$f(x)=\sum _{i=0}^n{a}_i\cdot x^i$$若能从中分离出一个二次因式$${\omega }^{\ast }(x)=x^2+u^{\ast }x+v^{\ast }$$此二次因式的零点即确定了一对共轭复根。Bairstov方法的策略就是从这对共轭复根出发，通过适当的迭代收敛到精确解，通常情况下，二次因式也可以是某个近似的二次因式$$\omega (x)=x^2+ux+v\text{\hspace{1em}(6)}$$
• 用二次因式 $\omega (x)$ 除 $f(x)$，记商为 $p(x)$，余项为 $r(x)$，显然 $p(x)$ 是 $n-2$ 次多项式函数，$r(x)$ 是一次函数，不妨设为$$r(x)=r_0+r_1\cdot x$$于是可以将 $f(x)$ 表示如下：$$f(x)=\omega (x)\cdot p(x)+r(x)=(x^2+ux+v)\cdot p(x)+r_1\cdot x+r_0\text{\hspace{1em}(7)}$$通过比较等号两边相同次数项系数，可以发现 $r_0,r_1$ 是 $u,v$ 的函数，分别记为$$r_0=g(u,v)$$$$r_1=h(u,v)$$贝尔斯多夫方法的策略就是逐步修改 $u,v$ 的值，使得余项 $r_0,r_1$ 变得很小。
• 考察方程组$$\{\begin{array}{rl}& r_0=g(u,v)=0\\ & r_1=h(u,v)=0\end{array}$$假设上述方程组有解 $(u^{\ast },v^{\ast })$，将上述方程组左边使用多元函数Taylor展开在 $(u,v)$ 展开到一次项(注意这里的 $u,v$ 代表的是本次迭代时的值)，则有$$\{\begin{array}{rl}& g(u,v)+\frac{\partial g}{\partial u}\cdot \Delta u+\frac{\partial g}{\partial v}\cdot \Delta v=0\\ & h(u,v)+\frac{\partial h}{\partial u}\cdot \Delta u+\frac{\partial h}{\partial v}\cdot \Delta v=0\end{array}\text{\hspace{1em}(8)}$$对于方程组 $(8)$ 计算出 $(\Delta u,\Delta v)$即可得到改进的二次因式$$\omega (x)=x^2+(u+\Delta u)x+(v+\Delta v)$$于是目前的关键问题就是如何计算出 $(8)$ 式中的 $6$ 个系数。

$r_0=g(u,v);r_1=h(u,v)$

• 设 $$p(x)=\sum _{i=0}^{n-2}{b}_i\cdot x^i$$$$f(x)=\sum _{i=0}^n{a}_i\cdot x^i$$代入 $(7)$ 式可得$$\sum _{i=0}^n{a}_i\cdot x^i=r_0+r_1\cdot x+(x^2+ux+v)\cdot \sum _{i=0}^{n-2}{b}_i\cdot x^i$$比较等式两边相等次数项系数，可得$$a_n=b_{n-2}$$$$a_{n-1}=u\cdot b_{n-2}+b_{n-3}$$$$a_i=v\cdot b_i+u\cdot b_{i-1}+b_{i-2},i\in [2,n-2]$$$$a_1=v\cdot b_1+u\cdot b_0+r_1$$$$a_0=v\cdot b_0+r_0$$将其变形则得到$$b_{n-2}=a_n$$$$b_{n-3}=a_{n-1}-u\cdot b_{n-2}$$$$b_i=a_{i+2}-v\cdot b_{i+2}-u\cdot b_{i+1},i\in [0,n-4]$$$$r_1=a_1-v\cdot b_1-u\cdot b_0$$$$r_0=a_0-v\cdot b_0$$

$\frac{\partial g}{\partial v},\frac{\partial h}{\partial v}$

• 将 $(7)$ 对 $v$ 求导，得到$$0=p(x)+(x^2+ux+v)\frac{\partial p}{\partial v}+s_{1}x+s_0$$$$s_1=\frac{\partial h}{\partial v}$$$$s_0=\frac{\partial g}{\partial v}$$
• 因此可以将 $p(x)$ 表示为$$p(x)=-(x^2+ux+v)\frac{\partial p}{\partial v}-s_1\cdot x-s_0\text{\hspace{1em}(9)}$$$(9)$ 的意义在于，用 $\omega (x)=x^2+ux+v$ 除 $p(x)$，可以得到商为 $n-4$ 次多项式 $\frac{\partial p}{\partial v}$，余数为 $s_1\cdot x+s_0$，上面的表述中忽略了负号。
• 设$$\frac{\partial p}{\partial v}=\sum _{i=0}^{n-4}{c}_i\cdot x^i$$采用和上面一致的比较相等次数项系数的方法，可以得到 $s_0,s_1.$

$\frac{\partial g}{\partial u},\frac{\partial h}{\partial u}$

• 将 $(7)$ 对 $u$ 求导，得到$$0=x\cdot p(x)+(x^2+ux+v)\frac{\partial p}{\partial u}+\frac{\partial h}{\partial u}\cdot x+\frac{\partial g}{\partial u}$$即$$x\cdot p(x)=-(x^2+ux+v)\frac{\partial p}{\partial u}-\frac{\partial h}{\partial u}\cdot x-\frac{\partial g}{\partial u}$$另外根据 $(9)$ 式可知$$x\cdot p(x)=-(x^2+ux+v)\cdot x\cdot \frac{\partial p}{\partial v}-(s_1\cdot x-s_0)\cdot x=-(x^2+ux+v)\cdot (x\cdot \frac{\partial p}{\partial v}+s_1)+u{s}_0\cdot x+s_1\cdot x+v\cdot s_1$$对比可知$$\frac{\partial h}{\partial u}=-(u{s}_0+s_1)$$$$\frac{\partial g}{\partial u}=-v{s}_1$$

---

• 关于 $$\frac{\partial f}{\partial u}=\frac{\partial f}{\partial v}=0$$因为 $f(x)$ 的取值是和 $u,v$ 无关的，在 $(7)$ 式中改变 $u,v$ 会影响到 $r_1=g(u,v),r_0=h(u,v)$ 的取值，但不会影响整个函数值。

---