【机器学习】数值分析02——任意方程求根

任意方程求根

全文目录

（博客园）机器学习

（Github）MachineLearning Math

1.简介

方程和函数是代数数学中最为重要的内容之一，从初中直到大学，我们都在研究着方程与函数，甚至我们将图形代数化，从而发展出了代数几何、解析几何的内容。而在方程与函数中，我们研究其性质最多的，往往就是方程的根（零点），即使是研究方程的极值点、鞍点等，我们无非也只是研究其微商的零点。
我们在初等数学中已经学过许多简单初等函数、线性方程的求解方法，在本文中，我们重点讨论任意方程，尤其是计算困难的非线性方程的求根方法。

2.方程

2.1分类和介绍

方程就是指含有未知数的等式。是表示两个数学式（如两个数、函数、量、运算）之间相等关系的一种等式，使等式成立的未知数的值称为“解”或“根”。在这里，根据一些性质的不同，我们将方程分成以下几类：

• 单个方程
  • 线性方程：本质是等式两边乘以任何相同的非零数，方程的本质都不受影响。通常认为只含有一次项的方程。
  • 非线性方程：是因变量与自变量之间的关系不是线性的关系的方程。
    • 多项式方程
    • 超越方程：指含有未知量的超越式（指数、对数、三角函数、反三角函数等）的方程。换言之，超越方程中都有无法用含有未知数的多项式、分式或开方表示的式子。
• 多个方程
  • 线性方程组
  • 非线性方程组

2.2方程的零点（根、解）

若有一个值或一些值能够使得方程 \(f(x)=0\) 成立，那么这个值就被成为方程的解，也常常被叫做零点和根。

若方程有且只有一个解\(x^*\)，那么我们称方程有单根\(x^*\)。

若对于方程\(f(x)=0\)，有\(f(x^*) = 0,f^{'}(x^*)=f^{''}(x^*)=\cdots=f^{(k)}(x^*)=0,f^{(k+1)}(x^*)\neq0\),那么称\(x^*\)为方程的k+1重根

PS：若方程是简单幂函数多项式组成，那么方程的解的数量应和最高此项的数值一致，因为存在虚根。

3.求根方法

求根的方法基本上大同小异，都是通过区间去逼近方程的根的点。

首先我们说一个定理1：对于实函数方程\(f(x)=0\)，当\(x\in(a,b)\)，且\(f(x)\)在\(x\in(a,b)\)时单调且连续，若\(f(a)\cdot f(b)<0\)，则方程在\(x\in(a,b)\)有且只有一个根。

3.1二分法

3.1.1普通二分法的原理及操作

二分法和算法中的二分搜索法非常的类似，取定有根区间\([a,b]\)进行对分，求得\(mid = \frac{a+b}{2}\)进行判断含根区间，如果\(f(a)\cdot f(mid)<0\),则令\(b=mid\)；反之若\(f(b)\cdot f(mid)<0\),则令\(a=mid\)。当\(|b_n-a_n|<\epsilon\)停止计算返回结果。

产生的截断误差为\(|e_{n-1}| = x_{n+1} - x^*\leq[b_n - a_n] = \frac{b_0 - a_0}{2^n}\)。

可以计算出最小迭代次数为\(n = \frac{lg(c_0-a_0)-lg\epsilon}{lg2}\)

代码实现（更多语言的代码见仓库中Code文件夹）：

private static double epsilon = 0.001;
// func为函数，写法如x=>x*x+2*x-1，a，b必须为有效的含根开区间
public static double Binary(Func<double, double> func, double a, double b)
{
    var f1 = func.Invoke(a);
    var f2 = func.Invoke(b);
    if (f1 * f2 > 0)
        throw new ArgumentException("此区间无根或根不唯一");
    double mid = (a + b) / (double)2;
    var fm = func.Invoke(mid);
    if (fm == 0)
        return fm;
    if (f1 * fm < 0)
        b = mid;
    else if (f2 * fm < 0)
        a = mid;
    if (Math.Abs(b - a) <= epsilon)
        return (a + b) / (double)2;
    return Binary(func, a, b);
}

3.1.2普通二分法准确度及速度分析

假设\(\left[a,b\right]\)是二分法的初始区间，在进行n此二分之后，得到的最终根的分布区间\([a_n,b_n]\)的长度为\((b-a)/2^n\)，我们将其中点作为根的最优估计值，与真实值之间的误差不超过区间长度的一半。

\[e = \left|x_c-r\right|<\frac{b-a}{2^{n+1}} \]

若\(e<\frac{1}{2}\times10^{-p}\),则精确到小数点后p位，

事实上我们也可以简单的计算出函数执行的次数为n+2次。

3.2浅谈迭代

3.2.1 迭代是什么

迭代是重复反馈过程的活动，其目的通常是为了逼近所需目标或结果。每一次对过程的重复称为一次“迭代”，而每一次迭代得到的结果会作为下一次迭代的初始值。
重复执行一系列运算步骤，从前面的量依次求出后面的量的过程。此过程的每一次结果，都是由对前一次所得结果施行相同的运算步骤得到的。这篇文章中对于迭代有一个非常不错的概述究竟什么是迭代？

3.2.2 不动点的定义

各位可以尝试以下操作，

1. 随意输入一个数字\(\lambda\)
2. 然后对其进行\(cos\lambda\)运算，
3. 将运算结果作为新的值传回\(cosx\)函数之中

当你一直重复以上三个操作，你会发现数字最后会定格在0.73908513左右不动了。这是一个非常有趣的现象，事实上我们的这个操作就是\(x_n = cosx_{n-1}\)，而最后不变化的数值实际上就是\(x = cosx\)这个方程的解。

这就是我们不动点的一种实际情况，不动点原理是数学上一个重要的原理，也叫压缩映像原理或巴拿赫（Banach）不动点定理，完整的表达：完备的度量空间上，到自身的一个压缩映射存在唯一的不动点。用初等数学可以这么理解：连续映射\(f\)的定义域包含值域，则存在一个\(x\)使得\(f(x)=x\)。

若某函数满足\(f(\lambda)=\lambda,\lambda \in R\)，我们就称\(\lambda\)为函数的一阶不动点。同样的，我们推广一下，若\(f(f(\lambda))=\lambda,\lambda \in R\)，则称为二阶不动点，一阶不动点必定是二阶不动点。

不动点的存在定理：若某函数\(y=f(x)\)与\(y=x\)存在至少一个交点，那么函数必然存在不动点。

3.2.3 函数的相似性

我们常常说图形之间的相似，事实上函数也有其相似的定义。若函数\(f(x)\)对其换元，令\(t = \varphi(x),x=\varphi^{-1}(t),f(x)=>f(\varphi^{-1}(t))\)。我们的不动点是一个基于迭代的过程，迭代就是让\(f(x_0) = x_1 = t,f(x_1)=f(f(x_0))=f(t)=x_2...\)的过程，如果我们换个角度去思考，迭代不就是和我们刚刚提到的换元是一个道理吗？那么我们思考一下，如果对于二次迭代\(f(f(x))\)，如果说我们把函数完全看成是一个复合函数，那么我们现在需要做的就是试图将函数写成\(f(\varphi^{-1}(t))\)的形式。

现在做出如下变换，将外层的\(f(x)\)换元为\(f(\varphi^{-1}(t))\)，那么二次迭代式就变成了\(g(t) = \varphi(f(\varphi^{-1}(t)))\)。随后对我们的\(g(x)\)再进行迭代就会得到\(g(g(t)) = \varphi(f(f(\varphi^{-1}(t)))\)，利用数学归纳法计算n次迭代式，那么就会得到\(g_n(t) = \varphi(f_n(\varphi^{-1}(t)))\)。

我们在这给出总结和定义：若函数\(g(x),f(x),\varphi(x)\)，若有\(g(x)=\varphi(f(\varphi^{-1}(x)))\)，那么称函数\(g(x)\)和\(f(x)\)通过\(\varphi(x)\)相似，记作\(f\sim g\)，\(\varphi(x)\)称为相似的桥函数，且之前我们的过程得出了一个重要的结论就是\(g\sim f=>g^n\sim f^n\)，同时，相似函数的不动点是完全一致的，证明如下：

\[若函数g(x)的不动点为x_0，则有 \]

\[g(x_0)=x_0=\varphi(f(\varphi^{-1}(x_0))) \]

\[根据反函数性质有\varphi(\varphi^{-1}(x))=x \]

\[故f(\varphi^{-1}(x_0)) = \varphi^{-1}(x_0)=x_0,\varphi^{-1}(x_0)是f(x)的一个不动点 \]

通过上述得出的不动点，很容易通过刚才讲的内容得出\(\varphi^{-1}(x_0)\)实际上也是\(g(x)\)换元后得到的不动点的位置，因此相似函数的不动点是一一对应的。

3.2.4 迭代收敛速度

迭代是一种逼近、利用数列、函数收敛的性质进行求解的方法，那么对于收敛的速度，我们很难直接从数值看出速度，因为收敛过程中，误差的变化是越来越小的，因此我们再次进行一个比阶。

\[\exists p\in N \And C>0,\displaystyle \lim_{k \to \infty}\frac{e_{k+1}}{e_k^p} = C \]

其中：

\[e_k = |x_k-x^*| \]

\[p= \begin{cases} p=1,线性收敛\\ p=2,平方收敛\\ p>2,高阶收敛 \end{cases} \]

并且有一个定理：若\(\varphi(x)\)在所求的根\(x^*\)邻域有连续的二阶导数，并且有\(|\varphi^{'}(x)|<1\)，有以下特点：

1. 当\(\varphi^{'}(x)\neq0\)，那么迭代过程线性收敛

2. 当\(\varphi^{'}(x)=0,\varphi^{''}(x)\neq0\)，那么迭代是平方收敛的。

证明过程留给读者，只需要利用泰勒展开是便可以证明该定理。

3.3不动点迭代法

3.3.1不动点迭代法的原理及操作

迭代法是将方程求根问题转换成了函数求交点问题再转换成数列求极限的问题，这个过程中，对方程进行一个巧妙的变换之后，方程就可以在若干次迭代之后解出一个近似解。

操作方法如下：首先将方程\(f(x)=0\)改写成\(x = g(x)\)的形式，这样就可以将方程的解看成是函数\(y=x\)和\(y=g(x)\)的交点。给定初始值\(x_0\)后，则\(x_1 = g(x_0)\)，不断重复这个过程，若\(\displaystyle \lim_{k \to \infty}x_k\)存在，则迭代便可以达到使得\(x_k\)趋于交点，而这个交点就是我们刚才讲了那么久的不动点。

我们即使构造出了\(x=g(x)\)的迭代式，往往由于我们的变换只是一些简单的移项、通分等四则运算获得，使得迭代式也并不是一个很容易被计算的函数，这个时候，我们之前讲到的那么多函数相似就派上了用场，我们可以利用函数相似性去寻找一个更为简单的函数\(\varphi(x)\sim g(x)\)，由于相似的性质，不动点并不会发生变化。不过，我们需要保证最终的迭代函数在指定的含根区间内存在一阶导数且导数值\(|\varphi^{'}(x)|<1\)，否则迭代函数将会不收敛。

迭代过程收敛的情况如下图所示：

如果说我们的导数不满足条件，那么我们的迭代将会发散：

代码样例

// func是迭代函数而不是实际函数
public static double Iterative(Func<double, double> func, double x)
{
    double temp = func.Invoke(x);
    if (Math.Abs(temp - x) <= epsilon)
    {
        return temp;
    }
    x = temp;
    return Iterative(func, x);
}

3.3.2迭代法的收敛性证明

在这里，我们将证明迭代法求根的合理性和可行性。

前提条件：设函数\(\varphi(x), x\in[a,b]\)有连续的一阶导数，并且满足以下条件：

• \(\forall x\in [a,b],\varphi(x)\in[a,b]\)

• \(\exists L \in (0,1),\forall x\in[a,b],|\varphi^{'}(x)|\leq L\)

要证明和解决以下命题和问题：

• 在\(x\in[a,b],\exists x^*,\varphi(x^*) = 0\)

• \(\forall x_0\in[a,b]\)，迭代过程\(x_{k+1} = \varphi(x_k)\)均收敛与\(x^*\)

• 求解误差分析式

现在开始证明

1：证明在区间内有且只有一个根存在：

证：在\(x \in [a,b]\)时，\(\varphi^{'}(x)\)存在，所以有\(\varphi(x)\)连续，于是可以作\(g(x) = x - \varphi(x)\)，易知\(g(x)\)连续。

又因为\(\varphi(x) \in [a,b]\)，且\(g(a)*g(b)<0\),故存在实根，使得\(x=\varphi(x)\)。

利用反证法：若在[a,b]上还有一实根\(\bar{x}\)，那么通过拉格朗日中值定理必定有：

\[x^*-\bar{x} = \varphi(x^*)-\varphi(\bar{x}) = \varphi^{'}(\xi)(x^*-\bar{x})\Longrightarrow \varphi^{'}(\xi) = 1 \]

显然与假设不符合。

2：证明这个根收敛于\(x^*\)

根据拉格朗日中值定理，容易知道

\[x^*-x_{k+1} = \varphi(x^*)-\varphi(x_k) = \varphi^{'}(\xi)(x^*-x_k) \]

又因\(|\varphi^{'}(x)|\leq L\)，故易得：

\[|x^*-x_{k+1}|\leq|L(x^*-x_k)|=|L^2(x^*-x_{k-1})=\cdots=|L^k(x^*-x_0)| \]

因为\(\displaystyle \lim_{k \to \infty}L^k=0\)，故\(\displaystyle \lim_{k \to \infty}|x^*-x_{k+1}|=0\)(绝对值必定非负)，得\(x^*=x_{k+1}(k\to\infty)\)

3:迭代法的误差式：

设某一次迭代后误差为\(\epsilon\)，则可以推出：

\[|x_{k+1}-x_{k}| = |(x^*-x_k)-(x^*-x_{k+1})\geq|x^*-x_k|-|x^*-x_{k+1}|\geq|x^*-x_k|-L|x^*-x_k| \]

其中从\(|x^*-x_{k+1}| \leq L|x^*-x_k|\)则是因为对于一个可以收敛的函数而言，每一次迭代后，误差总是减少的。

故可以轻松的计算出误差估计式为：

\[|x^*-x_k|\leq\frac{1}{1-L}|x_{k+1}-x_k| \]

通过中值定理，又可以推出：

\[|x_{k+1}-x_k| = |\varphi'(\xi)(x_k-x_{k-1})| \]

因为\(|\varphi'(\xi)|\leq L\)，将上式递推后放缩成第二种误差估计式：

\[|x^*-x_k|\leq\frac{L^k}{1-L}|x_1-x_0| \]

其中L可以近似的看作是一个足够接近方程解的函数导数的值。

3.3.3不动点迭代法的收敛速度

首先我们先直接给出一个结论：在不动点迭代过程中，第i步迭代和第i+1步迭代的表达式分别为\(e_i,e_{i+1}\)，且有\(\displaystyle \lim_{ i \to \infty}\frac{e_{i+1}}{e_i}=S<1\neq0\)，因此不动点迭代法是一个线性收敛的算法。

如何能证明我们的算法是一个线性收敛的函数呢？实际上非常简单。证明过程如下：

\[若\varphi(x)连续可微,不动点的一个足够近似的值为x_0，x_{i+1}和x_i是两次迭代的估计值 \]

\[根据中值定理，\exists r\in (x_i,r),\frac{\varphi(x_i)-x_0}{x_i-x_0}=\varphi'(r) \]

\[\varphi(x_i) = x_{i+1},故x_{i+1}-r = \varphi'(r)(x_i-r) \]

\[e_{i+1}=\left|\varphi'(r)\right|e_i，即S=\left|\varphi'(r)\right|为常数 \]

我们可以在这里举几个例子，例如我们希望求函数\(f(x) = x^3+x-1=0\)的根，倘若我们取\(\varphi(x) = 1-x^3\)，根据我们误差分析的式子，包括我们收敛分析的公式来说，后一步的误差将会是\(e_{i+1}=-3r^2e_i\)，我们预估的根大概是0.6823左右，可以发现我们的收敛速度的绝对值已经大于1，意味着每一步的误差都在扩大，因此函数不会收敛。如果我们取\(\varphi(x) = \sqrt{1-x}\)，我们继续计算出的收敛速度大约是0.7左右，很明显误差在以线性缩小。如果再次修改我们的递推式\(\varphi(x) = \frac{1+3x^3}{1+3x^2}\)，我们的S几乎是等于0的，因此也解释了这种方法为何是最快收敛到我们的根。

//todo 图片

再举一个例子，函数\(f(x) = cosx-sinx=0\)的迭代求根过程，我先给定我们的近似值大约是0.7853982。这里似乎很难构造我们的迭代递推式，事实上只需要方程两边同时加x就可以凑出\(x=cosx-sinx+x=\varphi(x)\)。我们利用计算机编写好程序并计算，最后得到了这样一组迭代过程中的一些数据：

不难看出我们的\(\frac{e_i}{e_{i-1}}\)在一个恒定值0.414保持了相当的次数，这也是因为\(S=\left|\varphi'(r)\right|=|1-sinr-cosr|\approx0.414\)，印证了我们的推论是完全正确的，当我们不动点迭代在一个区间内，往往会以线性模式进行误差的缩小。

迭代的终止

3.4牛顿法

牛顿在整个微积分、数值分析中都有着举足轻重的地位，这里阐述的牛顿法，全名牛顿-拉夫逊法，就是Taylor展开式的一部分。牛顿迭代法的核心思想就是：设法将一个非线性方程\(f(x)=0\)转化为某种线性方程去求解。在数学意义上，我们知道泰勒公式可以将任意函数以简单幂函数的形式表示出来，而在几何意义上，我们是利用切线与X轴的交点去进行迭代，在根处，切线必过零点。若设\(f(x)=0\)的近似解为\(x_k\)，则方程可以用一阶Taylor展开进行近似表示，牛顿迭代法的核心公式如下：

\[p_1(x) = f(x_k)+f^{'}(x_k)(x-x_k) \]

3.4.1普通牛顿法

从上述公式中我们知道，其中\(p_1(x)\)就是泰勒多项式的表达，将其看成是方程\(f(x)=0\)，通过迭代的思想，从而得到了一个线性方程：

\[x_{k+1} = x_k - \frac{f(x_k)}{f^{'}(x_k)} \]

这就是普通牛顿法的迭代公式。在几何意义上，他就是一个切线与x轴的交点去逼近零点，如图：

我们不断的通过迭代这个过程，就能让我们的值越来越精确。

3.4.2牛顿下山法

所有的迭代法都有一个无法逃过的缺点，如果选取的初始值离根太远，则会导致迭代过程次数过多，并且有可能导致迭代发散，因为迭代法只具有局部收敛性。

为了避免迭代失败或时间过长，我们加上这样一个条件用于保证数值稳定下降

\[|f(x_{k+1})|<|f(x_k)| \]

将这个条件和牛顿法结合再一起，即再保证数值稳定下降的过程中，利用牛顿法加快迭代，这就是牛顿下山法。具体的操作如下：

将牛顿法的结果\(\bar{x}_{k+1}\)与前一步的迭代值\(x_k\)进行加权平均作为新的迭代值\(x_{k+1}\)，

则有：

\[x_{k+1} = \lambda\bar{x}_{k+1} + (1-\lambda)x_k \]

或

\[x_{k+1} = x_k - \lambda\frac{f(x_k)}{f^{'}(x_k)} \]

其中\(\lambda(0\leq\lambda<1)\)称为下山因子，它的值是一个逐步试探的过程，可以从1开始取值，一旦满足\(|f(x_{k+1})|<|f(x_k)|\)则称为下山成功，否则需要另选初始值\(x_0\)进行试算。

3.4.3简单牛顿法

牛顿法可以说是一个非常有效的计算方法，准确度和迭代次数上都要比普通的迭代法要好得多，但是牛顿法最大的问题是我们需要求方程的导数，对于某些极其复杂的函数而言，导数是无法通过人工的方式计算，假如我们使用微积分的方式去求解导数，这对整个程序的性能会有比较大的影响，因此我们可以利用一个常数值\(\lambda\)来代替方程中的导数项。此时迭代公式为：

\[x_{k+1} =x_k - \frac{f(x_k)}{\lambda} \]

不过对于常数\(\lambda\)的取值是有限制的，因为我们需要保证迭代函数的收敛性，如果函数不收敛于\(x^*\)，那么一切都没有意义。于是有：

\[\varphi(x) = x - \frac{f(x)}{\lambda}\Longrightarrow\varphi^{'}(x) = 1-\frac{f^{'}(x)}{\lambda} \]

牛顿迭代法的收敛性遵循前文中普通迭代法的收敛性，于是可以得到：

\[|\varphi^{'}(x)| = |1-\frac{f^{'}(x)}{\lambda}|\Longrightarrow0<\frac{f^{'}(x)}{c}<2 \]

这样我们就可以很轻松的确定下c的值了。

简单牛顿法的几何意义就简单许多了，和我们之前讨论的普通迭代法一致，只不过普通迭代法是将函数值和\(y=x\)进行处理变换，而简单牛顿法则是和\(y= \lambda(x-x_k)+f(x_k)\)进行变换，本质是一致的。

///这里只写普通牛顿法，另外的函数由读者自己思考
// 其中f1x为导数
public static double Newtown(Func<double, double> fx, Func<double, double> f1x, double x)
{
    var temp = f1x.Invoke(x);
    if (temp == 0)
    {
        throw new ArgumentException();
    }
    x = x - fx.Invoke(x) / temp;
    if (Math.Abs(fx.Invoke(x)) <= epsilon)
    {
        return x;
    }
    return Newtown(fx, f1x, x);
}

3.4.4 改进牛顿法——弦截法

弦截法是牛顿法的一种改进，保留了牛顿法中收敛速度快的优点，还克服了再牛顿法中需要计算函数导数值\(f^{'}\)的缺点。

弦截法中，我们用差商

\[\frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}} \]

去代替牛顿法中的导数值，于是可以得到以下离散化的迭代

\[x_{k+1} = x_k - \frac{f(x_k)}{f(x_k)-f(x_{k-1})}(x_k-x_{k-1}) \]

这种方法叫做双点弦截法，如图所示：

从图中可以知道，弦截法一直利用两点之间的连线作为迭代的内容，那么，他的合理性在哪里呢？

整个迭代法都离不开中值定理，这里也是这样，事实上，这个差商之所以可以对导数值进行替代，是因为中值定理中说过，连续函数中两函数上的点的连线的斜率必为两点之间某一点的导数值，并且迭代过程中，这两点的中值会越来越逼近函数零点，事实上这已经说明了弦截法是牛顿法的改进方法了。

不过，如果将函数中\(x_{k-1}\)用一个定点\(x_0\)代替，这种方法叫做单点弦截法，几何意义如图所示：

//这里就对双点弦截法进行编程
public static double StringCut(Func<double, double> func, double x1, double x2)
{
    var temp = x1 - (func.Invoke(x1) / (func.Invoke(x1) - func.Invoke(x2))) * (x1 - x2);
    x2 = x1;
    x1 = temp;
    if (Math.Abs(func.Invoke(x1)) <= epsilon)
    {
        return x1;
    }
    return StringCut(func, x1, x2);
}

3.4.5 牛顿法的收敛性

牛顿法的收敛性有可能是一阶收敛也有可能是二阶收敛，具体要看我们的迭代函数。我们首先讲讲牛顿法二次收敛的情形。

假定我们牛顿法的迭代函数为\(g(x) = x-\frac{f(x)}{f'(x)}\)，那么我们迭代函数的导数为\(g'(x) = \frac{f(x)f''(x)}{f'(x)^2}\)，设根为r，那么易知\(g'(r)=0\)，迭代函数是收敛的。

我们再将\(f(r)\)进行泰勒展开，取\(x_i\)是一个足够接近根的一次迭代结果，并且代入\(f(r)=0\)则有

\[f(r) = f(x_i)+(r-x_i)f'(x_i)+\frac{(r-x_i)^2}{2}f''(\xi),\xi \in (x_i,r) \]

\[f(r) = 0 =>x_i-\frac{f(x_i)}{f'(x_i)} -r = \frac{(r-x_i)^2}{2}\frac{f''(\xi)}{f'(x_i)} \]

若\(f'(x_i)\neq0\)对此函数迭代式在进行一次迭代可以得到以下式子

\[x_{i+1}-r=e_i^2\frac{f''(\xi)}{2f'(x_i)} \]

\[e_{i+1} = e_i^2\left|\frac{f''(\xi)}{2f'(x_i)}\right| \]

这样我们就得到了牛顿法的二阶收敛的定义，对于二阶收敛，我们要求是需要在该处的一阶导数值不为0。

不过我们的牛顿法并不是总为二次收敛的，在面临多重根的时候，牛顿法往往会倒退到线性收敛的速度。，例如函数\(f(x)=x^m\)利用牛顿法求根的时候，牛顿公式如下：

\[x_{i+1} = x_i-\frac{x_i^m}{mx_i^{m-1}}=\frac{m-1}{m}x_i \]

唯一的根是0，但是0是一个m-1重根，因此得到

\[e_{i+1} = Se_i,S=\frac{m-1}{m} \]

我们这里可以给出一个定义，若m+1阶连续可微函数有一个m重根，那么利用牛顿法会线性局部收敛，误差式如上。

收敛迭代的加速与精度

利用迭代法进行收敛的计算时，往往面临一个很麻烦的问题，就是当运行足够多次迭代的时候，我们的精度往往会不降反升，这是因为一种类似摇摆的问题出现，例如我规定此刻你站在坐标1上，我希望你走到坐标0上，你每次可以向左或者向右走三步，那么第一次你向左走三步之后，得到的迭代结果时-2，反而离精确值更远了。这是一个令人糟心的问题。这里我们对这种现象进行一次解释。

收敛的精度

y还是x？

事实上我们的误差分为了前向误差和后向误差，什么是前向误差呢？实际上就是我们要求解的根与实际的差值，假设\(f(r)=0\)这个方程求根，前向误差定义为\(\left|r-x_i\right|\)，而后向误差指代的就是我们在函数值上的接近程度，也就是y轴的接近程度，我们定义为\(\left|f(x_i)\right|\)。

对于计算机而言，存储的精度是有极限的，假定我们现在有一台古老的计算机，他的存储精度最多就到小数点后4位，那么假设我们的迭代函数是\(f(x)=x^2\)，我们知道在精度范围内，0.001是他能取到最小具有有效数字的精度解，如果我们现在通过迭代取到的数字是0.001的话，我们会发现一个惊人的问题，再进行一次迭代之后，得到\(f(0.001)=0.000001\)，由于我们的计算机无法存储足够的精度，也就是说这个值会被计算机看成是0，但我们知道这并不是我们的最优解。

这其实就是我们后向误差足够小，但是前向误差并不够小导致的。我们在实际的代码编写的时候，一定要注意的一个地方就是我们的误差在什么时候应当利用前向误差进行终止，什么时候又应该利用后向误差进行求解。

敏感性

读者可以尝试将迭代初始值设置成我们方程实际的根进行迭代，按着我们正常人的理解，他应该立即停止迭代直接返回我们的初始值即可。但是算法往往是不尽人意的，在绝大多数工具或语言的内置求根函数中，即使我们传入的初始值是根，往往得到的结果反而是一个不精确的解。例如著名的威尔金森多项式：

\[W(x) = (x-1)(x-2)...(x-20) \]

如果我们利用这个已经因式分解好了的式子进行求根运算，这似乎没有必要，我们可以得到一个非常工整的方程根，他就是1-20之间所有的自然数组成，一旦我们将其展开后投入运算，由于存在许多大数字，很多大数字会因为存储、计算等过程中产生损失和误差，这些误差实际上是很小的数字，但是却对我们的计算结果产生了很大的误差，最终我们难以得到一个精确解。这就是对于求根过程中对数字扰动的敏感性问题。

我们进一步去探索一下究竟是什么导致了误差放大，并且能粗略的计算出这些误差究竟有多少，我们假设是寻找一个函数\(f(x)=0\)的根r，同时对我们的函数做出一个很小的扰动\(\epsilon g(x)\)，其中我们的系数\(\epsilon\)很小，用\(\Delta r\)表示根的变化。那么会有以下式子。

\[f(r+\Delta r)+\epsilon g(r+\Delta r) = 0 \]

将上式统统泰勒展开

\[f(r)+(\Delta r)f'(r)+\epsilon g(r)+\epsilon(\Delta r)g'(r)+O((\Delta r)^2) = 0,其中\displaystyle O((\Delta r)^2) \to 0,\Delta r \to 0 \]

忽略后面的高阶无穷小，同时r又是根，于是可以得到一个关于r的变化值的近似

\[\Delta r (f'(r)+\epsilon g'(r)) \approx -f(r)-\epsilon g(r) =-\epsilon g(r) \]

又因为\(\epsilon\)也是一个很小的值（需要远小于\(f'(r)\)），我们直接认为他就是0，得到了\(\Delta r\)的最终估计:

\[\Delta r \approx -\epsilon\frac{g(r)}{f'(r)} \]

这就是我们根的敏感性公式，利用这个公式，我们可以得出，如果出现了某些数字的扰动，根会如何变化。

这里举一个黄皮书上的例子，如果函数\(P(x) = (x-1)(x-2)...(x-6)-10^{-6}x^7\)，求出它的最大根。我们肉眼可见最大根应该是在6附近，如果没有最后那一项，我们的最大根就是6，那么问题就在于如果我去除最后一项，对根会有多大的影响呢？

那么我们直接设近似函数为\(f(x)=(x-1)(x-2)...(x-6)\)，令我们的\(\epsilon g(x)=-10^{-6}x^7\)。通过我们的根敏感性公式计算一番之后，得到了

\[\Delta r \approx -\epsilon\frac{g(r)}{f'(r)} = 2332.8\times 10^{-6} \]

因此我们估计的最大根应当是6.0023328左右，读者可以自行使用迭代法去求解最大根，最后得出的结果与这个估计值是极度接近的。这里我们定义一个新东西误差放大因子，也就是我们的\(\left|\frac{g(r)}{rf'(r)}\right|\)。

普通迭代加速

对于迭代过程，利用中值定理有

\[x^*-x_{k+1} = \varphi(x^*)-\varphi(x_k) = \varphi^{'}(\xi)(x^*-x_k) \]

当\(\Delta x \rightarrow0\)，我们将\(\varphi^{'}(\xi)\)看成定值\(\lambda(\lambda<1)\)

于是有

\[\lambda(x^*-x_k) = x^*-x_{k+1}\longrightarrow x^* = \frac{1}{1-\lambda}\bar{x}_{k+1}-\frac{\lambda}{1-\lambda}x_k \]

故可以推出最终的迭代公式为

\[\begin{cases} \bar{x}_{k+1} = \varphi(x_k)\\ x_{k+1} = \bar{x}_{k+1} + \frac{\lambda}{1-\lambda}(\bar{x}_{k+1}-x_k) \end{cases} \]

利用这种方式的好处就是再求得\(x_k\)时，利用加权的方式，使得迭代法变得有点类似像牛顿迭代法一样变成切线性质的线而不是x=c或y=c的线，你经过画图和简单的代数运算后，你会发现\(\bar{x}_{k+1}<x_{k+1}\)，也就是说达到了我们的加速目的。

Atiken加速法

Atiken加速法其实就是再普通加速迭代公式上在进行一次迭代，这里直接写出公式：

\[\begin{cases} \bar{x}_{k+1} = \varphi(x_k)\\ \bar{\bar{x}}_{k+1} = \varphi(\bar{x}_{k+1})\\ x_{k+1} = \bar{\bar{x}}_{k+1} - \frac{(\bar{\bar{x}}_{k+1}-\bar{x}_{k+1})}{\bar{\bar{x}}_{k+1} - 2\bar{x}_{k+1}+x_k} \end{cases} \]

练习题

1. 试着不借助计算器计算\(\sqrt{2}\)的值
2. 试着证明\(f(x)=ax+b\)在牛顿法中会一步收敛
3. 编程题：由一个高度为10m的圆柱构成的发射井顶部是一个半球，体积\(400m^3\)，确定发射井底部班级
4. 设函数为\(f(x)=x^n-ax^{n-1},g(x)=x^n\)，利用敏感性公式手动估计\(f_\epsilon(x) = x^n-ax^{n-1}+\epsilon x^n\)的非零根

Reference

《数值分析》（原书第二版） —— Timothy Sauer

《数值计算方法》（清华大学出版社） —— 吕同富等

About Me

作　　者：WarrenRyan
出　　处：https://www.cnblogs.com/WarrenRyan/
本文对应视频：BiliBili
关于作者：热爱数学、热爱机器学习，喜欢弹钢琴的不知名小菜鸡。
版权声明：本文版权归作者所有，欢迎转载，但未经作者同意必须保留此段声明，且在文章页面明显位置给出原文链接。若需商用，则必须联系作者获得授权。
特此声明：所有评论和私信都会在第一时间回复。也欢迎园子的大大们指正错误，共同进步。或者直接私信我
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是作者坚持原创和持续写作的最大动力！


博主一些其他平台：
微信公众号：寤言不寐

QQ群：236996657
BiBili——小陈的学习记录
Github——StevenEco
BiBili——记录学习的小陈（计算机考研纪实）
掘金——小陈的学习记录
知乎——小陈的学习记录

联系方式

电子邮件：cxtionch@live.com

社交媒体联系二维码：