§$8.7$ 方程求根与近似计算

方程求根是一个很有实际意义的问题, 在本书中到目前为止也已多次出现. 实际上, 研究迭代生成数列的动机之一就是求方程的近似解. 这方面的一个基本定理就是连续函数的零点存在定理. 在例题 $\mathrm{5.2.1}$ 的注 $2$ 中已经指出, 用闭区间套定理给出的二分法证明同时就提供了一种近似求根方法. 在这一节中我们将介绍 $Newton$ (牛顿) 求根法. 它不仅对方程求根问题是一种有效的计算方法, 同时在整个计算数学领域中都是重要的基本方法. 在数学分析教材中关于方程求根的材料还可以参看 [$14$] 第一卷的第四章第 $5$ 节, [$26$] 的第一卷第八章.

$\mathrm{8.7.1}$ 迭代算法的收敛速度

在第二章的 $\mathrm{2.5.2}$ 小节中已经看到, 在求自然对数的底 $e$ 的近似值时, 不同的计算方法的效果完全不一样. 这是因为在那里的两个数列收敛于 $e$ 的速度有明显的不同. 因此在作迭代计算时必须考虑方法的效率. 在本小节我们将引进迭代算法的阶的概念, 为下一节介绍 $Newton$ 求根法作好准备. 同时还介绍计算圆周率的两种不同的迭代算法.

现在引入迭代算法的阶的定义. 设某个迭代算法的第 $n$ 次计算的误差为 ${\epsilon }_n$ . 若存在一个常数 $\alpha$ , 使得接连两次的误差之间满足递推估计式

$$k{\epsilon }_n^{\alpha }⩽{\epsilon }_{n+1}⩽K{\epsilon }_n^{\alpha },\text{\hspace{1em}(8.19)}$$

其中 $k,K$ 为两个正常数, 则称此算法为 $\alpha$ 阶算法. 其特例之一是存在极限

$$\underset{n\to \infty }{\lim }\frac{{\epsilon }_{n+1}}{{\epsilon }_n^{\alpha }}=A\ne 0,\text{\hspace{1em}(8.20)}$$

则算法的阶为 $\alpha$ 。这时当 $n$ 足够大时成立 ${\epsilon }_{n+1}\approx A{\epsilon }_n^{\alpha }$ 。又若对某个算法只知道公式 $(8.19)$ 右边的不等式成立, 则可以说该算法的阶不低于 $\alpha$ .

对于一般的收敛数列来说, 若有 ${\lim }_{n\to \infty }{x}_n=a$ ，而误差 ${\epsilon }_n=|x_n-a|(n\in {\mathbf{N}}_{+})$ 满足 $(8.19)$ 或 $(8.20)$，则称这个数列的收敛速度为 $\alpha$ 阶.

在 $\alpha =1$ 时的算法称为一阶 (线性) 算法. 不难看出一阶算法的收敛速度是比较慢的. 为简单起见, 不妨设有 ${\epsilon }_{n+1}\approx C{\epsilon }_n$ , $C$ 为某个正常数, 则就可以得到

$${\epsilon }_n\approx C^n{\epsilon }_0.\text{\hspace{1em}(8.21)}$$

可见常数 $C$ 必须小于 $1$, 而且越小越好. 又可以看出, 为了使精度达到 $10^{-n}$ 所需的迭代次数总是 $O(n)$ 的量级.

具有一阶收敛速度的数列在第二章中很多. 从 $(8.20)$ 可见, 只要 $\{x_n\}$ 为无穷小量, 且满足

$$\underset{n\to \infty }{\lim }\left|\frac{x_{n+1}}{x_n}\right|=A\ne 0$$

(参见 $\mathrm{2.7.3}$ 小节的第一组参考题中的题 $4$), 则 $\{x_n\}$ 就是一阶收敛于 $0$ 的数列. (当然在第二章中还有许多收敛速度远低于一阶的数列.) 现在举出一个一阶算法的重要例子. 这就是计算圆周率的 $Archimedes$ - 刘徽算法 (见第二章第一组参考题 $20$).

例$8.7.1$

设 $a_1>b_1>0$ ，并用递推公式

$$a_{n+1}=\frac{2a_n{b}_n}{a_n+b_n},b_{n+1}=\sqrt{a_{n+1}{b}_n},n\in {\mathbf{N}}_{+}\text{\hspace{1em}(8.22)}$$

作迭代, 证明: $\{a_n\}$ 和 $\{b_n\}$ 以一阶速度收敛于同一极限.

证明

在这里只对收敛速度进行分析. 应用与例题 $\mathrm{2.3.5}$ 中类似的方法即可证明 $\{a_n\}$ 和 $\{b_n\}$ 收敛于同一极限. 记此极限为 $A>0$ , 即有

$$\underset{n\to \infty }{\lim }{a}_n=\underset{n\to \infty }{\lim }{b}_n=A,$$

又记

$${\epsilon }_n=a_n-b_n,n\in {\mathbf{N}}_{+},$$

则可以对收敛速度分析如下. 首先进行恒等式运算

$$\begin{array}{rl}& a_{n+1}^2-b_{n+1}^2=a_{n+1}^2-a_{n+1}{b}_n=a_{n+1}(a_{n+1}-b_n)\\ & =a_{n+1}\cdot \frac{2a_n{b}_n-b_n(a_n+b_n)}{a_n+b_n}=\frac{a_{n+1}{b}_n(a_n-b_n)}{a_n+b_n},\end{array}$$

因此得到

$${\epsilon }_{n+1}=\frac{a_{n+1}{b}_n}{(a_{n+1}+b_{n+1})(a_n+b_n)}\cdot {\epsilon }_n.$$

利用 $\{a_n\}$ 和 $\{b_n\}$ 收敛于同一极限 $A>0$ ，就有近似估计

$${\epsilon }_{n+1}\approx \frac{1}{4}{\epsilon }_n,$$

因此收敛速度为一阶.

注$1$

若取 $a_1=2\sqrt{3},b_1=3,$ 即单位圆的外切和内接正六边形的半周长, 可以证明极限 $A=\pi$ . 这就是计算圆周率的 $Archimedes$ - 刘徽算法. 从以上分析有

$${\left(\frac{1}{4}\right)}^5=\frac{1}{2^{10}}\approx 0.001,$$

可见这种算法每迭代 $5$ 次大致可以增加 $3$ 位新的有效数字。这个估计与用这个算法求 $\pi$ 的大量实际计算完全符合.

注$2$

$Archimedes$ - 刘徽算法有个变形：令 $A_n=1/a_n,B_n=1/b_n$ ，有递推式：

$$A_{n+1}=\frac{A_n+B_n}{2},B_{n+1}=\sqrt{A_{n+1}{B}_n}.\text{\hspace{1em}(8.23)}$$

与算法 $(8.22)$ 比较, 每次迭代的计算量少得多. 但是这个改进并没有提高收敛速度. 若仍令 ${\epsilon }_n=B_n-A_n$ , 则还是得到 ${\epsilon }_{n+1}\approx \frac{1}{4}{\epsilon }_n$ .

当 $\alpha =2$ 时算法为二阶收敛 (也称为平方收敛). 这时情况大不相同. 为简便起见, 只考虑误差递推估计的单侧不等式

$${\epsilon }_{n+1}⩽K{\epsilon }^2.\text{\hspace{1em}(8.24)}$$

上式两边乘以常数 $K$ , 就有 $K{\epsilon }_{n+1}⩽(K{\epsilon }_n{)}^2$ . 因此可以继续做下去, 得到

$$K{\epsilon }_{n+1}⩽(K{\epsilon }_n{)}^2⩽(K{\epsilon }_{n-1}{)}^4⩽\cdots ⩽(K{\epsilon }_0{)}^{2^{n+1}}.$$

这样就得到

$${\epsilon }_n⩽\frac{1}{K}(K{\epsilon }_0{)}^{2^n}.\text{\hspace{1em}(8.25)}$$

这里的常数 $K$ 不一定要小于 $1$, 只要取初始值足够好, 使得 $K{\epsilon }_0<1$ 即可. 这时在公式 $(8.25)$ 右边的表达式收敛于 $0$ 的速度是非常快的.

下一个例题是对在第二章中的例题 $\mathrm{2.3.5}$ (其中的极限是 $Gauss$ 的算术几何平均值) 进行收敛速度的分析.

例$8.7.2$

从 $0<b<a$ 出发, 设 $a_0=a,b_0=b$ , 并用递推公式

$$a_n=\frac{a_{n-1}+b_{n-1}}{2},b_n=\sqrt{a_{n-1}{b}_{n-1}},n\in {\mathbf{N}}_{+}\text{\hspace{1em}(8.26)}$$

作迭代, 证明: $\{a_n\}$ 和 $\{b_n\}$ 以二阶速度收敛于同一极限.

证明

记极限为 $AG(a,b)$ , ${\epsilon }_n=a_n-b_n$ , $n\in {\mathbf{N}}_{+}$ . 利用恒等式

$$a_{n+1}^2-b_{n+1}^2=\frac{1}{4}(a_n-b_n{)}^2,$$

就可以得到

$${\epsilon }_{n+1}=\frac{1}{4(a_{n+1}+b_{n+1})}{\epsilon }_n^2\approx \frac{1}{8AG(a,b)}{\epsilon }_n^2,$$

因此是二阶算法.

□

注

迭代公式 $(8.26)$ 与 $(8.23)$ 很相似, 但实际上收敛速度完全不同.

$Salamin$-$Brent$算法

在 $1976$ 年出现了 $Salamin$-$Brent$ 算法. 它就是以 $(8.26)$ 和上述分析为基础的. 由于收敛速度快, 因此成为一类重要的算法, 具有广泛的应用. 自此以后计算圆周率 $\pi$ 的新纪录很多都是用这类新算法得到的. 目前已经有计算 $\pi$ 的任意高阶的算法. 下面列出计算 $\pi$ 的二阶算法中的一个算法, 以及它的计算效果. 关于它以及其他有关材料可以从 [$4$, $1$] 中找到, 还可以参看 [$60$] 第六章的圆周率及其计算.

$Salamin$-$Brent$算法：令 $a_0=1$, $b_0=s_0=\frac{1}{\sqrt{2}}$ ，并用递推公式

$$\begin{array}{l}{a}_n=\frac{a_{n-1}+b_{n-1}}{2},\\ b_n=\sqrt{a_{n-1}{b}_{n-1}},\\ s_n=s_{n-1}-2^n(a_n^2-b_n^2),\\ p_n=\frac{2a_n^2}{s_n}\end{array}$$

作迭代, 则 $\{p_n\}$ 二阶收敛于 $\pi$ .

它的计算结果为

迭代次数	1	2	3	4	5	6	7	8	9	$\cdots$	25
有效位数	1	4	9	20	42	85	173	347	697	$\cdots$	$⩾4.5\times 10^6$

从上面的计算结果可以看到, 每迭代一次, 有效位数几乎增加一倍. 实际上, 这是二阶算法的共同特征. 从公式 $(8.24)$ 可以知道, 若常数 $K<1$ , 则一定如此. 这在下面的简单例题中可以看得很清楚.

例$8.7.3$

从初始值 $x_0=1$ 开始, 用迭代算法

$$x_{n+1}=\frac{1}{2}{x}_n+\frac{1}{x_n}$$

求无理数 $\sqrt{2}$ 的近似值, 观察有效位数的增长情况.

解

前几个值很容易计算:

$$x_2=1.5,x_3\approx 1.417,x_4\approx \mathrm{1.414216.}$$

与 $\sqrt{2}=1.4142135623\cdots$ 比较, 可见有效位数分别为 $1$, $3$, $6$. 为了继续计算下去, 我们需要使用如 $Mathematica$ 那样的软件. 具体地说, 即每次将 $x_n$ 与 $\sqrt{2}$ 都计算到足够多的位数, 然后进行比较, 从而确定第 $n$ 次近似值 $x_n$ 的有效位数. 这里只列出实际计算结果:

迭代次数	1	2	3	4	5	6	7	8
有效位数	1	3	6	12	25	49	98	196

实际上, 不难直接验证在这个例题中的迭代算法确实为二阶算法. 用 §$2.5$ 的方法知数列 $\{x_n\}$ 从 $x_2$ 起严格单调减少收敛于 $\sqrt{2}$ . 然后估计迭代误差

$$0<{\epsilon }_{n+1}=x_{n+1}-\sqrt{2}=\frac{1}{2}{x}_n+\frac{1}{x_n}-\sqrt{2}=\frac{(x_n-\sqrt{2}{)}^2}{2x_n}\approx \frac{1}{2\sqrt{2}}{\epsilon }_n^2,$$

可见恰好为二阶算法.

我们即将看到, 这个开平方根的算法就是 $Newton$ 求根法的一个例子.

$\mathrm{8.7.2}$ $Newton$ 求根法

可以从不同的角度来导出 $Newton$ 求根法 (也称为 $Newton$-$Raphson$ (拉弗森) 求根法). 第一种推导方法完全来自于 $Newton$ 求根法的几何意义. 从图 $8.11$ 上可以看出已有 $f(a)<0$ 和 $f(b)>0$ , 因此从连续函数的零点存在定理知道在区间 $(a,b)$ 内方程 $f(x)=0$ 一定有根. 在图上记这个根为 $\xi$ . 问题是如何计算出根 $\xi$ 的近似值.

以 $b=x_0$ 为初值, 在点 $B(b,f(b))$ 作曲线 $y=f(x)$ 的切线, 与 $x$ 轴的交点记为 $x_1$ . 然后再在点 $(x_1,f(x_1))$ 作曲线 $y=f(x)$ 的切线, 与 $x$ 轴交于点 $x_2$ . 如此进行下去, 不难得到一般的递推公式是

图8.11

$$x_{n+1}=x_n-\frac{f(x_n)}{f^{\prime }(x_n)}.\text{\hspace{1em}(8.27)}$$

从图$8.11$可以看出, 所得到的数列 $\{x_n\}$ 有可能会很快收敛到方程 $f(x)=0$ 的根 $\xi$ . 这就是 $Newton$ 求根法, 也称为 $Newton$ 切线法.

但这里实际上有许多问题需要研究. 首先, 若设方程 $f(x)=0$ 在区间 $(a,b)$ 中只有唯一的一个根 $\xi$ , 是否一定成立

$$\underset{n\to \infty }{\lim }{x}_n=\xi .$$

其次, 这个算法的收敛速度如何? 此外还有一个如何取初始值的问题. 在图 $8.11$ 上从点 $A(a,f(a))$ 也作出了曲线 $y=f(x)$ 的切线, 但它与 $x$ 轴的交点 (在图上未画出) 很可能会越出函数 $f(x)$ 的定义域.

可以从完全不同的角度导出 $Newton$ 求根法. 例如, 假定已得到根 $\xi$ 的第 $n$ 个近似值 $x_n$ , 又设 $f$ 二阶可微, 则可以写出带 $Lagrange$ 余项的 $Taylor$ 展开式:

$$0=f(\xi )=f(x_n)+f^{\prime }(x_n)(\xi -x_n)+\frac{f^{\prime \prime }({\theta }_n)}{2}(\xi -x_n{)}^2,\text{\hspace{1em}(8.28)}$$

其中的中值 ${\theta }_n$ 在 $x_n$ 和 $\xi$ 之间. 在公式 $(8.28)$ 中弃去右边的最后一项, 并将由此得到的根 $\xi$ 的近似值记为 $x_{n+1}$ , 就得到了前面已有的公式 $(8.27)$. 在将 $Newton$ 求根法推广到高维空间时用这种局部线性化的方法没有本质上的困难.

现在给出保证 $Newton$ 求根法能够成功的一组充分条件, 并证明它的收敛速度恰好是二阶. 为直观起见, 下面所述的条件均与图 $8.11$ 一致. 读者可自行写出与这组充分条件平行的其他充分条件.

命题$8.7.1$

设函数 $f(x)$ 在区间 $[a,b]$ 上二阶连续可微，且满足以下条件

(1) $f(a)<0,f(b)>0$ ;

(2) $f^{\prime }(x)>0,\forall x\in [a,b];$

(3) $f^{\prime \prime }(x)>0,\forall x\in [a,b];$

(4) 取 $x_0=b$ 为初始值;

那么就有

(1) 方程 $f(x)=0$ 在区间 $(a,b)$ 内存在唯一的根 $\xi$ ;

(2) 由递推公式 $(8.27)$ 得到的数列 $\{x_n\}$ 是在区间 $[a,b]$ 中的严格单调减少数列, 且以 $\xi$ 为极限;

(3) 数列 $\{x_n\}$ 的收敛速度为二阶.

证明

从连续函数的零点存在定理知道方程 $f(x)=0$ 在 $(a,b)$ 中有根. 从 $f^{\prime }(x)$ 在区间上处处大于 $0$ 可知函数 $f(x)$ 严格单调增加, 因此方程的根唯一.

从公式 $(8.28)$ 和 $f^{\prime \prime }(x)>0,\forall x\in [a,b]$ 可见, 成立

$$f(x_n)+f^{\prime }(x_n)(\xi -x_n)<0,\forall n\in {\mathbf{N}}_{+}.$$

由此即可得到

$$\xi <x_n-\frac{f(x_n)}{f^{\prime }(x_n)}=x_{n+1},\forall n\in {\mathbf{N}}_{+}.$$

从递推公式 $(8.27)$ 又可看出只要 $f(x_n)>0$ ，就有

$$x_{n+1}<x_n$$

成立. 由于初始值 $x_0=b$ 处有 $f(b)>0$ , 而当 $x_n>\xi$ 时就有 $f(x_n)>0$ , 因此就保证了数列 $\{x_n\}$ 是严格单调减少数列, 且以 $\xi$ 为下界. 记该数列的极限为 $\eta$ , 在公式 $(8.27)$ 两边令 $n\to \infty$ , 得到

$$\eta =\eta -\frac{f(\eta )}{f^{\prime }(\eta )},$$

可见 $f(\eta )=0$ 。由于方程 $f(x)=0$ 在 $[a,b]$ 中的根唯一，因此 $\eta =\xi$ 。这样就证明了 $Newton$ 求根法所得到的迭代数列 $\{x_n\}$ 收敛于方程在区间 $[a,b]$ 中的唯一根。

现在用 $Taylor$ 公式 $(8.28)$ 估计迭代误差

$$\begin{array}{rl}0<x_{n+1}-\xi & =x_n-\frac{f(x_n)}{f^{\prime }(x_n)}-\xi \\ & =\frac{(x_n-\xi )f^{\prime }(x_n)-f(x_n)}{f^{\prime }(x_n)}\\ & =\frac{f^{\prime \prime }({\theta }_n)(x_n-\xi {)}^2}{2f^{\prime }(x_n)},\end{array}$$

其中 ${\theta }_n\in (\xi ,x_n)$ .利用 $f^{\prime \prime }$ 连续且处处大于 $0$，就有

$$\underset{n\to \infty }{\lim }\frac{x_{n+1}-\xi }{(x_n-\xi {)}^2}=\frac{f^{\prime \prime }(\xi )}{2f^{\prime }(\xi )}\ne 0,$$

因此 $\{x_n\}$ 二阶收敛于 $\xi$ .

注$1$

容易看到在命题中的某些条件可以放宽, 而仍保证迭代数列为二阶收敛或不低于二阶收敛. 但另一方面也可以举出各种例子, 说明在 $4$ 个条件中的某些条件不成立时, $Newton$ 求根法可能失败. 一般而言, $Newton$ 求根法在有拐点、重根和多根时会有困难.

不难得到 $Newton$ 求根法的先验估计和事后估计 (可以与命题 $\mathrm{3.4.4}$, 即压缩映射原理中的结果作比较). 由于 $f$ 在 $[a,b]$ 上二阶连续可微, 因此在命题的条件下 $f^{\prime \prime }$ 在 $[a,b]$ 上有最大值 $M>0$ , 而 $f^{\prime }(x)$ 在 $[a,b]$ 上有最小值 $m=f^{\prime }(a)>0$ . 从证明中关于迭代误差的估计有

$$\left|x_{n+1}-\xi \right|⩽\frac{M}{2m}{\left|x_n-\xi \right|}^2,$$

代入公式 $(8.25)$ 就得到

$$\left|x_n-\xi \right|⩽\frac{2m}{M}{\left(\frac{M}{2m}\left|x_0-\xi \right|\right)}^{2^n}.$$

在计算前, 用它可以估计为达到指定精度所需的迭代次数, 即先验估计.

注意：为了满足 $\frac{M}{2m}|x_0-\xi |<1$ 的要求，区间 $[a,b]$ 的长度应当满足不等式

$$b-a<\frac{2m}{M}.$$

但从命题的证明知道, 即使这个条件不成立, 迭代数列仍然二阶收敛.

再利用 $Lagrange$ 中值定理得到

$$|f(x_n)-f(\xi )|=|f(x_n)|⩾m|x_n-\xi |,$$

又用 $Taylor$ 公式有

$$\begin{array}{rl}f(x_n)& =f(x_{n-1})+f^{\prime }(x_{n-1})(x_n-x_{n-1})+\frac{f^{\prime \prime }(\eta )}{2}(x_n-x_{n-1}{)}^2\\ & =\frac{f^{\prime \prime }(\eta )}{2}(x_n-x_{n-1}{)}^2,\end{array}$$

其中 $\eta$ 在 $x_n$ 与 $x_{n-1}$ 之间，就可以得到

$$|x_n-\xi |⩽\frac{1}{m}|f(x_n)|⩽\frac{M}{2m}|x_n-x_{n-1}{|}^2.$$

这可用于从相继两次计算的结果去估计当时的误差大小, 即事后估计.

例$8.7.4$

用 $Newton$ 求根法导出开平方根的迭代算法.

解

设要求正数 $A$ 的平方根. 取函数 $f(x)=x^2-A$ , 求 $\sqrt{A}$ 的问题就成为方程求根问题了. 将 $f$ 代入迭代公式 $(8.27)$ 中, 就得到

$$x_{n+1}=x_n-\frac{x_n^2-A}{2x_n}=\frac{x_n}{2}+\frac{A}{2x_n}.$$

注$1$

在 $A=2$ 时就得到例题 $\mathrm{8.7.3}$ 中的算法. 实际上本题的算法已在 $\mathrm{2.6.3}$ 小节的练习题 $8$ 中出现. 其中练习题 $9$ 还给出了求平方根的一个三阶算法.

注$2$

这个算法的历史可以上溯到古代巴比伦文明 (见 [$35$]). 其思路可能是: 如果 $x$ 是 $\sqrt{A}$ 的一个近似值, 那么 $A/x$ 也是一个近似值, 这两个近似值的乘积等于 $A$ , 而且分别在 $\sqrt{A}$ 的两侧, 因此取算术平均值可能会得到更好的结果.

$\mathrm{8.7.3}$ 练习题

本节的计算题应当根据在学习中所能使用的计算工具来安排. 下面的题中只有第 $1$ 题是计算题, 且只需用计算器.

$1$. 用 $Newton$ 求根法计算 (要求精确到 $0.0001$):

(1) $x^3-2x^2-4x-7=0$ 在 $[3,4]$ 之间的根的近似值；

(2) $\sin x=1-x$ 的根的近似值.

$2$. 在命题 $\mathrm{8.7.1}$ 中给出了保证 $Newton$ 求根法成功的充分条件, 其中共有 $4$ 项要求. 试举出例子, 说明不满足其中的某些要求时, 用 $Newton$ 求根法有可能失败.

$3$. 证明: 在例题 $\mathrm{5.2.1}$ 中提供的二分法是方程求根的一阶算法.

$4$. 证明: $\mathrm{2.6.3}$ 小节的题 $9$ 中的算法是求平方根的三阶算法.

$5$. 设 $A>0$ , 为计算立方根 $\sqrt[3]{A}$ , 从 $x_0>0$ 出发用递推公式

$$x_{n+1}=\frac{x_n(x_n^3+2A)}{2x_n^3+A}$$

作迭代, 证明数列 $\{x_n\}$ 收敛, 求出其极限, 并确定其收敛的阶.

$6$. 用 $Newton$ 求根法设计一个求 $\sqrt[k]{A}(A>0,k>0)$ 的迭代算法, 并对其收敛速度作出分析.

$7$. 用 $Newton$ 求根法设计一个求 $\frac{1}{A}$ 的迭代算法, 其中只用加法和乘法运算, 并对其收敛速度作出分析.