7.1 非线性方程求根的基本问题#

类似于线性代数方程组 Ax=b 求解的问题, 非线性方程求解的问题也可以提为 $f(x)=y$ , 其中 $f: R^{n} \to R^{m}$ 是一个非线性函数, $y \in R^{m}$ 是一个给定的向量, $x \in R^{n}$ 是一未知向量. 这一问题也可以描述如下: 在映射 f 下, y 的原像 x 是什么? 一般地, 我们可以令 $g(x)=f(x)-y$ , 原问题也就等价于求解 $g(x)=0$ , 后者也可以认为是非线性求根问题的一般形式. 方程 $f(x)=0$ 的根, 称为函数 f 的零点, 函数的零点可以是实零点, 也可以是复零点. 若函数 f 有零点 $x^{*}$ , 并且 $f(x)=(x-x^{*})^{m}g(x)$ , 其中 $g(x^{*}) \neq 0$ , m 是正整数, 则称 $x^{*}$ 是函数 f 的 m 重零点.

非线性方程求根问题可能有解, 也可能无解, 不论是求实根还是求复根. 对于多项式函数, 我们有如下熟知的代数基本定理:

讨论其他非线性函数的零点的个数, 甚至零点的存在性, 都是一个复杂的问题, 但在某些情况下我们可以保证零点的存在, 例如中值定理.

下面的例子表明, 就算非线性函数的形式是固定的, 讨论它的零点个数的问题也会因为其中参数值的变化而变得异常复杂.

例 7.1.3 考察非线性方程 $f(x)=\frac{x\sin x}{x^{2}+1}=\alpha$ 的根的个数, 其中 $\alpha$ 是某个给定常数.

函数 $f$ 的图像如图7-1所示. 易知给定的非线性方程的根的个数与 $\alpha$ 的值有着非常复杂的关系.

图7-1 函数 $f(x) = \frac{x\sin x}{x^2 + 1}$ 的图像

通常情况下非线性方程只能求近似根. 求近似根经常采用下面的步骤: 首先确定非线性方程的有根区间, 即包含根的区间, 或者给出根的一个粗略的近似值; 然后, 根据精度要求把区间缩小到一定程度, 或者利用这个粗略的近似值把足够精度的根计算出来. 如果需要的精度很高, 则不同的方法, 其计算过程快慢可能相差非常大. 概括起来, 非线性方程求根需要考虑如下的三个问题.

(1) 根的存在性问题, 即能不能在理论上确定有没有根 (或实根), 有几个根.

(2) 有根区间的确定, 或者是给出一个粗略的近似值; 可以采用作图、分析函数性态, 或者大范围算法.

(3) 求出足够精度的近似根, 通过某些条件判别根的精度.

非线性方程求根的精度可能是很不相同的,例如下面的例子.

例 7.1.4 求解方程 $f(x)=\mathrm{e}^{x}-\alpha=0,\alpha$ 是一给定参数.#

对于不同的 $\alpha$ 值, 想要达到相同的近似根的精度, 困难程度是不一样的. 例如图7-2, 对应 $\alpha = 2$ 和 $\alpha = 0.1$ .

图 7-2 非线性方程求根的不同困难程度

在实际计算中, $e^{x}$ 的计算会有误差, 所以 $f(x)$ 的值可能在图示的虚线带内变化, 其近似根则在虚线带与实轴相交的区间内变化. 容易推出, 该区间的长度应正比于 $u/|f'(x^{*})|$ , 其中 u 为机器精度, 这个量可以用来衡量求根的难易程度.

当 $f'(x^*)$ 数值很小时, 称方程的根 $x^*$ 是坏条件的. 对于坏条件的根, 仅依靠 $|f(x)| \leqslant \varepsilon$ 来判定根的精度是不够的, 一般还应加上条件 $|x - x^*| \leqslant \varepsilon$ . 在重根处, 总有 $f'(x^*) = 0$ , 因此重根总是坏条件的.

若设定求根的精度 $\varepsilon = 0$ ，即在理论上探讨近似算法能否收敛到精确根，一般情况下我们就得到了一个无限的序列。该序列收敛的快慢是算法好坏的一个标准。记 $e_k = x_k - x^*$ 为第 $k$ 步迭代的误差，其中 $x^*, x_k$ 分别为真解和第 $k$ 步的近似解。一个迭代序列称为 $r$ 阶收敛的，或者收敛阶为 $r, r \geqslant 0$ ，如果 $e_k \to 0$ ，并且

这里, $C$ 是大于零的常数. 当 $r = 1,2,3$ 时, 分别称序列是线性收敛、平方收敛和立方收敛的; 如果 $2 > r > 1$ , 或者

称序列是超线性收敛的. 一个 $r$ 阶 $(r > 1)$ 收敛序列, 其近似解的有效位数基本按照 $r$ 倍的速度扩大. 能够产生 $r$ 阶收敛序列的方法也称为 $r$ 阶算法.

例 7.1.5 考察如下几个序列的收敛速度:

其中, $F_{k}$ 是 Fibonacci 序列, 即满足递推关系 $F_{k + 2} = F_{k + 1} + F_{k}, F_{1} = F_{2} = 1$ 的序列.

可以证明四个序列都收敛到 0, 且前两个是线性收敛的, 后两个分别是超线性收敛和平方收敛的, 如图 7-3 所示. 事实上, 利用通项公式 $F_{n} = \frac{1}{\sqrt{5}} \left( \left( \frac{1 + \sqrt{5}}{2} \right)^{n} - \left( \frac{1 - \sqrt{5}}{2} \right)^{n} \right)$ 可以证明, 第三个序列的收敛阶为 $r = \frac{1 + \sqrt{5}}{2} \approx 1.618$ .

图 7-3 不同的序列收敛速度

由于在求解方程之前，精确的根并不清楚，一个合理的收敛估计可以设定为 $|x_{k + 1} - x_k| / |x_k|^r$ 。若有某个 $r$ ，使得该项变化不大，可以认为它就是算法的收敛速度。在计算过程中，若 $|x_{k + 1} - x_k| / |x_k|$ 非常小，迭代没有显著变化，继续迭代变得没有必要，因此可以用它来控制迭代的继续与否。若对求解的方程的非线性性态并不十分了解，可以同时加上条件 $|f(x_k)|$ 足够小作为收敛条件。

7.2 二 分 法#

设有单变量连续函数 $f(x)$ ，若存在某一连续区间 $[a, b]$ 使得 $f(a)f(b) < 0$ ，即函数值异号，则 $f$ 在区间 $[a, b]$ 中至少存在一个零点.

二分法的基本思想就是每次把区间二等分, 给出两个等分区间中有根的那个区间, 达到把区间缩小的目的. 二分法的具体做法如下.

记 $a_1 = a, b_1 = b$ . 令 $x_1 = (a_1 + b_1) / 2$ , 考虑下面各种情形. 若 $f(x_1) = 0$ , 则 $x_1$ 是函数 $f$ 的零点; 不然, $f(x_1)$ 必定与函数值 $f(a)$ 和 $f(b)$ 中的一个异号, 这样就给出一个新的有根区间, $[x_1, b_1]$ 或者 $[a_1, x_1]$ , 不管哪一种情形, 我们都把它记为 $[a_2, b_2]$ . 反复施行这个做法, 我们可以得到一个区间套

该区间套中每一个区间长度都只有它的前一个的一半. 由操作的方式, 如果该过程无限, 则每个区间有如下的几个性质:

(1) $f(a_{k})f(b_{k}) < 0$ ，至少存在一个 $x^{*}$ ，对于所有 $k, x^{*} \in [a_{k}, b_{k}]$

(2) $b_{k} - a_{k} = (b - a) / 2^{k - 1};$

(3) $x_{k} = (a_{k} + b_{k}) / 2,$ 并且 $|x_{k} - x^{*}|\leqslant (b - a) / 2^{k}$

由最后一个性质, 二分法是线性收敛的. 如果指定精度 $\varepsilon$ , 则最多需要迭代步数为

其中, 记号 $\lceil x \rceil$ 代表 x 往上取整.

算法 7.2.1 (二分法)

(1) 给定初始区间 $[a, b]$ ，满足 $f(a)f(b) < 0$ ，以及计算精度 $\varepsilon$ ;

(2) 令 $c = a / 2 + b / 2$ ;

(3) 若 $b - a \leqslant \varepsilon$ 或者 $|f(c)| \leqslant \varepsilon$ , 停止算法;

(4) 若 $\operatorname{sign}(f(a)) = \operatorname{sign}(f(c))$ , 令 $a = c$ , 否则 $b = c$ ; 转步骤 (2).

这里不写 $c = (a + b)/2$ 或者 $f(a)f(c) > 0$ ，而是用算法中的式子，这里考虑了在计算中尽可能避免溢出的情形，即 a, b 都很大，或者 $f(a), f(c)$ 都很大的情形.

例 7.2.2 用二分法求解非线性方程 $f(x)=x^{3}-x-1=0$ .

易知, 函数 $f(x)$ 满足 $f(1) = -1 < 0$ , $f(2) = 5 > 0$ , 所以 $f$ 在区间 [1,2] 内有一根. 用二分法迭代, 结果见表 7-1. 因此, 方程的近似根为 $x^{*} = \frac{1.32422 + 1.32520}{2} = 1.3247$ .

表 7-1 二分法计算结果

二分法对函数的要求低, 仅需要函数的连续性, 编程简单. 但是二分法只用到了函数值的符号, 而非函数值的大小, 因此收敛速度并不快. 另一方面, 二分法不能求复根, 一般不能求偶数重根, 也不能直接应用到多变量的情形.

7.3 不动点迭代方法#

线性方程组 Ax = b 可以写成某个等价的 $x = Bx + g$ ，从而产生对应的迭代方法 $\boldsymbol{x}^{(k+1)} = \boldsymbol{B}\boldsymbol{x}^{(k)} + \boldsymbol{g}$ ，或者是雅可比方法、高斯-赛德尔方法，或者其他的方法。类似于线性的情形，一个非线性方程 $f(x) = 0$ 也可以写成等价的方程 $x = \varphi(x)$ ，从而产生迭代算法$x_{k+1} = \varphi(x_k)$ . 不同的是, 这种产生函数 $\varphi$ 的方式可以有很多种, 甚至可以说是无穷无尽的. 若把函数 $\varphi$ 看成是一个映射, 求解 $x^{*}$ 满足 $f(x^{*}) = 0$ 相当于求解 $x^{*} = \varphi(x^{*})$ , 即求在映射 $\varphi$ 下不动的点. 因此, 方程 $x = \varphi(x)$ 称为不动点方程, $x^{*}$ 称为函数 $\varphi$ 的不动点, 该问题也称为不动点问题.

例 7.3.1 非线性方程 $f(x)=x^{3}-x-1=0$ 可以有如下四个不同的改造方式, 从而得到不同的迭代方法.

(1) $x = x^{3} - 1$ , 对应迭代方法 $x_{k + 1} = \varphi_1(x_k) = x_k^3 - 1$ ;

(2) $x = \sqrt[3]{x + 1}$ , 对应迭代方法 $x_{k + 1} = \varphi_2(x_k) = \sqrt[3]{x_k + 1}$ ;

(3) $x = \frac{1}{x^2 - 1}$ , 对应迭代方法 $x_{k + 1} = \varphi_3(x_k) = \frac{1}{x_k^2 - 1}$ ;

(4) $x = x - \frac{x^3 - x - 1}{3x^2 - 1}$ , 对应迭代方法 $x_{k + 1} = \varphi_4(x_k) = x_k - \frac{x_k^3 - x_k - 1}{3x_k^2 - 1}$ .

同前面的例子, 函数 $f(x)$ 满足 $f(1) = -1 < 0$ , $f(2) = 5 > 0$ , 所以 $f$ 在区间 $[1,2]$ 内有一根. 取 $x_0 = 1.5$ , 用上述的四个迭代方法计算的结果见表7-2. 可以看到, 不动点函数 $\varphi_1$ 对应的迭代发散到 $+\infty$ ; $\varphi_2$ 对应的迭代是收敛的, 但是速度不是很快, 迭代6步才有5个有效数字; $\varphi_3$ 对应的迭代不收敛, 实际上再计算下去, 到第12步就会碰到除零的情况; $\varphi_4$ 对应的迭代收敛最快, 到第6步显示的数字已全为有效数字.

表 7-2 $x^{3} - x - 1 = 0$ 的四个不同的迭代方法

下面讨论不动点函数 $\varphi$ 具有什么性质时, 可以使得不动点迭代收敛.

定理 7.3.2　设一元函数 $\varphi(x)$ 在区间 $[a,b]$ 上一阶连续可导, 且

(1) $a \leqslant \varphi(x) \leqslant b$ 对一切 $x \in [a, b]$ 成立;

(2) 存在常数 $L, 0 \leqslant L < 1$ , 使得 $|\varphi'(x)| \leqslant L$ 对一切 $x \in [a, b]$ 成立.

则成立如下结论:

(a) 对任何 $x_0 \in [a, b]$ , 由 $x_{k+1} = \varphi(x_k)$ 产生的迭代序列 $\{x_k\}$ 必定收敛于函数 $\varphi(x)$ 在区间 $[a, b]$ 上的唯一不动点, 即 $x^* = \varphi(x^*)$ ;

(b) 序列 $\{x_{k}\}$ 的收敛速度有估计

证：(a) 本定理的条件保证了 $\varphi$ 不动点的存在性和唯一性. 事实上, 若令 $g(x) = x - \varphi(x)$ , 则函数 $g$ 在区间 $[a, b]$ 上连续, 且有 $g(a) \leqslant 0$ , $g(b) \geqslant 0$ . 因此, 要么 $a, b$ 中有一个是 $g$ 的零点, 要么存在 $x^{*} \in (a, b)$ , 使得 $g(x^{*}) = 0$ . 函数 $g$ 的零点实际上就是 $\varphi$ 的不动点.

下面证明不动点的唯一性. 若还有 $\bar{x}^{*}$ 满足 $\bar{x}^{*} = \varphi(\bar{x}^{*})$ ，则

其中, $\xi$ 位于 $x^{*}$ 和 $\bar{x}^{*}$ 之间. 但我们有 $L < 1$ , 上式仅当 $x^{*} = \bar{x}^{*}$ 时成立.

由条件, 若 $x_{0} \in [a, b]$ , 则有 $x_{1} \in [a, b]$ , 依此类推则有, 对于任意 $k, x_{k} \in [a, b]$ . 因为

其中, $\xi_{k - 1}$ 位于 $x_{k - 1}$ 与 $x^{*}$ 之间. 递推下去就有

因此, $x_{k}\to x^{*}$

(b) 由上面的推导, 我们有 $|x_{k+1} - x^*| \leqslant L|x_k - x^*$ . 因此,

同样, 由中值定理,

所以，

由收敛速度的估计可知, 当 $L$ 较小时, 收敛较快, 反之, 当 $L$ 很靠近 1 时收敛很慢. 若 $L \geqslant 1$ , 则迭代不收敛. 从收敛的估计式, 我们可以得到算法需要达到某个指定精度所需的最大迭代步数, 当然实际所需的迭代步数可能比理论估计出来的少些. 该定理的运用有时是比较困难的, 我们需要给定区间 $[a, b]$ , 并估计给定的不动点函数在整个区间上的导数的最大值 $\max_{x \in [a, b]} |\varphi'(x)|$ . 倘若我们有 $|\varphi'(x^*)| < 1$ , 其中 $x^*$ 是不动点, 则存在 $x^*$ 的某个邻域 $N(x^*) = [x^* - \delta, x^* + \delta]$ , 使得对于任意 $x \in N(x^*)$ , 有 $|\varphi'(x)| \leqslant L < 1$ . 这样, 取 $x_0 \in N(x^*)$ , 不动点迭代即可收敛. 由于 $\delta$ 预先并不清楚, 因此无法说明 $x_0$ 是否足够好, 能够在区间 $N(x^*)$ 内, 这种在解的某个小邻域内收敛的性质称为局部收敛性, 若没有这种要求, 或者可以清楚地给出邻域的表达, 则称为全局收敛性.

例 7.3.1 的几个函数的图像如图 7-4 所示, 其局部的函数图像说明了定理的结论. 图中阴影区域的边界直线方程为 y = x 和 $x + y = 2x^{*}$ , 即通过点 $(x^{*}, x^{*})$ 且斜率为 $\pm 1$ 的两条直线.

图 7-4 不动点迭代的局部收敛性

例 7.3.3 求解方程 $x \ln x = 1$ 的根.

方程 $x\ln x = 1$ 等价于下面的几个不动点方程

其中最后一个方程是在原方程两边加上 x，对左边分解因式后得到的.

首先, $1.5 \ln 1.5 - 1 = -0.3918 < 0$ , $2 \ln 2 - 1 = \ln 4 - 1 > 0$ , 因此原方程在区间 [1.5,2] 中有一根. 记四个迭代函数分别为 $\phi_1, \phi_2, \phi_3, \phi_4$ . 由上面的判别方法, 我们有

因此, 对于所有 $x \in [1.5, 2]$ 有

采用第四个迭代函数, 取初值 $x_0 = 1.75$ 计算, 则在第4步就达到机器精度了, 计算结果列在表7-3中. 对其他三个函数的迭代过程也一并列在表7-3中. 计算的结果说明了先前的估计是有效的, 可以用来预测迭代的进程.

表 7-3 方程 $x \ln x = 1$ 的不同迭代法

7.4 迭代加速#

假设有不动点迭代 $x_{k + 1} = \varphi (x_k)$ ，且设不动点为 $x^{*}$ 。在前面不动点迭代的收敛分析中，我们有

其中, $\xi_{k}$ 是位于 $x_{k}$ 和 $x^{*}$ 之间的某个点. 倘若 $\varphi'(x)$ 变化不大, 能够得到估计 $\varphi'(\xi_k) \approx L$ , 那么从上式可以解出

上式可以从两方面来解释．一方面，它可以理解为利用两步迭代以及一些收敛分析，可以得到这两步迭代近似值的某种平均，它将是更好的近似；另一方面，该式可以看成是一个新的迭代方法，对应的新不动点函数为

在上述假定下, $\bar{\varphi}'(x)$ 几乎为零, 因此迭代收敛速度非常快. 当然, 这有赖于原来不动点函数 $\varphi(x)$ 导数估计的准确性.

例 7.4.1 加速例 7.3.3 中的函数 $\phi_{3}(x)$ 的迭代.#

和例 7.3.3 的计算过程类似, 我们有 $0.53 \approx \phi_{3}^{\prime}(1.5) \geqslant \phi_{3}^{\prime}(x) \geqslant \phi_{3}^{\prime}(2) \approx 0.44$ . 因此, 导数变化不大. 令 $L = \frac{0.53 + 0.44}{2} = 0.485$ , 加速迭代方法为

同样, 令 $x_0 = 1.75$ , 可以得到如表 7-4 的计算结果. 原来需要 41 迭代步, 经过加速后仅需 8 步.

表 7-4 加速方法计算效果

下面介绍的艾特肯 (Aitken) 加速方法可以在不估计导数的情形下达到加速的效果, 所需的假设仍是 $\varphi'(x)$ 变化不大.

根据前面的方法, 我们有

其中, $\xi_{k}$ 位于 $x_{k}$ 与 $x^{*}$ 之间, $\xi_{k - 1}$ 位于 $x_{k - 1}$ 与 $x^{*}$ 之间. 假若 $\varphi^{\prime}(x)$ 变化不大, 可以认为 $\varphi^{\prime}(\xi_{k})\approx \varphi^{\prime}(\xi_{k - 1})$ ，那么就有

从上式可以解得

假如原来有一个迭代序列 $\{x_{k}\}$ , 则该序列任意相邻的三项 $x_{k-1}, x_{k}, x_{k+1}$ 都可以依照上述公式得到一个值, 记为 $\bar{x}_{k+1}$ . 一般地, 序列 $\{\bar{x}_k\}$ 的收敛速度要比原序列 $\{x_k\}$ 快得多, 特别是原来的序列只有线性收敛时.

艾特肯加速方法虽然避开了导数 $\varphi'(x)$ 的估计, 但是它也有一定的缺陷. 当迭代序列靠近解时, $x_{k-1}, x_k, x_{k+1}$ 都是靠近 $x^*$ 的数, 公式 (7.22) 中就出现了分子分母都很靠近 0 的情况. 一般在这里, 我们可以采用双精度的计算方式.

艾特肯加速也可以采用下面的方式来进行:

关于艾特肯加速方法, 我们有如下的定理, 证明略.

例 7.4.3 对例 7.3.3 中的迭代函数 $\phi_{2}(x)$ 做艾特肯加速.

已知 $\phi_2(x) = \exp \left(\frac{1}{x}\right)$ , 因此它的艾特肯加速为

取初值 $x_{0} = 1.75$ ，可以得到如表 7-5 的计算结果。原来需要 51 迭代步，经过加速后仅需 3 步，但第四步就出现了除零的情形。艾特肯加速方法每一步的计算量约为原迭代每一步计算量的两倍。

表 7-5 艾特肯加速效果

7.5 牛顿法#

牛顿法是求解非线性方程的一个重要方法, 有时也称为牛顿-拉弗森 (Newton-Raphson) 方法, 可以由多种途径导出.

假设需要求解非线性方程 $f(x) = 0$ ，其中 $f$ 是一个二阶连续可微函数。若已有方程的某一个近似根 $x_{k}$ ，函数 $f$ 在近似根 $x_{k}$ 处的泰勒展开为

其中, $\xi_{k}$ 是位于 $x_{k}$ 与 $x$ 之间的某个数. 取泰勒展开的前两项, 称函数 $h$ 是函数 $f$ 的线性化:

当 $f'(x_k) \neq 0$ , 求解 $h(x) = 0$ , 把它的根作为 $f$ 的更好的近似根, 则有

这实际上就是迭代函数为 $\varphi(x) = x - \frac{f(x)}{f'(x)}$ 的不动点迭代方法. 例7.3.1和例7.3.3中的最后一个迭代方法都是牛顿法.

牛顿法的几何意义如图 7-5 所示. 求解非线性方程 $f(x) = 0$ 等同于求曲线 $y = f(x)$ 和 $x$ 坐标轴的交点. 如果有近似解 $x_{k}$ , 在曲线上就存在一点 $(x_{k}, f(x_{k}))$ , 曲线在该点的切线方程为 $y - f(x_{k}) = f'(x_{k})(x - x_{k})$ . 该切线和 $x$ 坐标轴的交点就是 $(x_{k+1}, 0)$ . 因此, 牛顿法也称为切线法.

图 7-5 牛顿法的几何意义

牛顿法也可以看成一种特殊的加速方法. 方程 $f(x) = 0$ 等价于 $x = x + f(x) = \phi(x)$ , 而这是一种不动点迭代的形式. 若令 $L = \phi'(x_k) = 1 + f'(x_k)$ , 则加速方法为

此即牛顿迭代方法.

关于牛顿法的局部收敛性, 我们有如下的定理.

证：牛顿法是对应于函数 $\varphi (x) = x - \frac{f(x)}{f'(x)}$ 的不动点迭代.我们有

若 $f'(x) \neq 0$ ，则有 $\varphi'(x^{*}) = 0$ 。因此，牛顿法是局部收敛的。

和不动点收敛性类似, 对于牛顿法迭代, 我们有

其中, $\xi_{k}$ 位于 $x_{k}$ 和 $x^{*}$ 之间. 因此

牛顿法有下面的全局收敛性理论.

证：由 $f'(x)$ 连续且非零，知 $f'(x)$ 定号，所以函数 f 在区间 $[a, b]$ 中严格单增或者严格单减。因此 $f(x) = 0$ 在该区间中的根是唯一的。考虑函数 f 的一阶和二阶导数的符号，由于一、二阶导数皆定号，可得如下四种情形：

(1) $f'(x)>0$ 且 $f''(x)>0;$

(2) $f'(x)>0$ 且 $f''(x)<0;$

(3) $f'(x) < 0$ 且 $f''(x) > 0;$

(4) $f'(x) < 0$ 且 $f''(x) < 0$ .

下面只考虑第一种情形, 其他三种情形类似可证.

在情形 (1) 下, $f(a) < 0$ 且 $f(b) > 0$ , 函数单调递增. 由于 $f''(x_0) > 0$ , $f(x_0) > 0 = f(x^*)$ , 因此 $x_0 > x^*$ . 下面用归纳法证明, 对于所有 $k, x_k > x_{k+1} > x^*$ . 首先, 由牛顿法迭代公式, 以及假设 $f(x_k) > 0$ ,

利用函数 f 的泰勒展开以及 $f''(x) > 0$ , 有

此即

上式左边也就是 $x_{k + 1}$ , 即 $x_{k + 1} > x^{*}$ . 综合上面两式, 对所有 $k$ , 成立 $x_{k} > x_{k + 1} > x^{*}$ . 所以, $\{x_{k}\}$ 是一个单调下降且有下界 $x^{*}$ 的序列, 必有极限, 记该极限为 $\tilde{x}$ . 在牛顿法迭代式两边取极限, 有

所以, $\tilde{x}$ 是 $f(x)=0$ 的根, 再由根的唯一性, $\tilde{x}=x^{*}$ .

上面的定理实际上要求初始点要有一定的性质, 下面的定理仅需要初始点落在某一区间内, 证明略.

例 7.5.4 设计一个算法计算 $\sqrt{a}$ ，其中 a > 0.

计算 $\sqrt{a}$ 相当于求解方程 $x^{2}-a=0$ 的正根. 牛顿法应用在该问题上, 有

其中 $x_{0}>0$ 可以预先给定. 容易证明,

因此

所以上述迭代法是二阶收敛的. 以 a = 2 为例, 有如表 7-6 的计算结果

表 7-6 $\sqrt{2}$ 的不同近似值

如果函数 $f$ 有重根, 则牛顿法一般不是二阶收敛的. 设 $f(x) = (x - x^{*})^{m}g(x)$ , 其中, $m \geqslant 2$ 且 $g(x^{*}) \neq 0$ . 我们有 $f'(x) = m(x - x^{*})^{m-1}g(x) + (x - x^{*})^{m}g'(x)$ , 所以

$x^{*}$ 是函数 $\mu (x) = \frac{f(x)}{f'(x)}$ 的单根, 运用牛顿法于函数 $\mu (x)$ , 我们得到

此时迭代仍有二阶收敛性质, 只是每一迭代步运算量较大, 还需要计算 f 的二阶导数.

由前面的定理, 牛顿法的全局收敛性不管是对函数还是对初始点, 要求都是比较高的, 一般情形不易验证. 牛顿下山法是有效降低这些要求的一种技巧. 牛顿下山法的基本思想如下.

在一定的条件下, 求解方程 $f(x) = 0$ 可以等价地看成求函数 $|f(x)|$ 的最小点. 若把函数 $|f(x)|$ 的图像想象为许多山峰的话, 求极小点就相当于找到山谷谷底. 牛顿法如果不收敛, 通常是在两面 (或更多) 山坡之间跳跃, 每次都跳过谷底. 从函数的角度来讲, 这个现象出现是因为每次修正迭代点 $x_{k}$ 时, 修正的幅度太大了.

牛顿下山法的迭代公式为

其中，

即牛顿法对应的参数 $\lambda$ 恒为 1, 而下山法的参数 $\lambda$ 是每一步变化的, $\lambda$ 取满足函数值下降条件 $|f(x_{k+1})| < |f(x_k)|$ 的 $1, \frac{1}{2}, \frac{1}{4}, \cdots$ 中最大的那个值.

(1) 给定初始值 $x_{0}$ ，精度 $\varepsilon, k = 0;$

(2) 若 $|f(x_k)| \leqslant \varepsilon$ , 近似解为 $x_k$ , 停止迭代;

(3) 令 $d_{k} = -\frac{f(x_{k})}{f'(x_{k})}, \lambda = 1;$

(4) 若 $|f(x_k + \lambda d_k)| < |f(x_k)|$ ，则 $x_{k+1} = x_k + \lambda d_k$ ，转 (5); 否则， $\lambda = \frac{1}{2}\lambda$ ，重复步骤 (4);

(5) $k = k + 1$ , 转 (2).

牛顿下山法对初始点没有特别的要求, 因此整个算法对初始点的依赖大大减小, 原来用牛顿法不收敛的问题用下山法就可能收敛了. 下山法是一种技巧, 不仅在牛顿法中可以应用, 在其他各种求根方法中或求极值的方法中都可以应用. 鉴于下山法的广泛使用, 我们把牛顿下山法的程序附在下面.

function [x,it,convg] = newton(x0,f,g,maxit,tol)
% find the zero of function f, with gradient g provided
% Usage: [x,it,convg] = newton(x0,f,g,maxit,tol)
    if nargin<5, tol = 1e-10;
    if nargin<4, maxit = 100;
    end; end
    x = x0;
    fx = feval(f,x);
    convg = 0;
    it = 1;
    while ~convg,
    it = it + 1;
    if norm(fx)<=tol,
    fprintf('Newton Iteration successes!!\n');
    convg = 1;
    return;
    end
    d = - feval(g,x) \ fx;
    lambda = 1;

function v = f(x)
    v = x.^2 + sin(10*x) - 1;
文件 g.m

lsdone = 0;
while ~lsdone,
    xn = x + lambda * d;
    fn = feval(f,xn);
    if abs(fn) < abs(fx),
    lsdone = 1;
    else
    lambda = 1/2 * lambda;
    if lambda <= eps,
    convg = -1;
    error('line search fails!!');
    end
    end
end
x = xn;
fx = fn;
if it > maxit,
    convg = 0;
    error('Newton method needs more iterations.!!');
end

例如, 若我们需要计算方程 $f(x) = x^{2} + \sin 10x - 1 = 0$ 的根, 该函数的图像如图7-6所示.

图7-6 函数 $f(x) = x^{2} + \sin 10x - 1$ 的图像

我们需要编写两个 Matlab 文件 (函数文件及其导数文件): 文件 f.m

function v = g(x)
    v = 2*x + 10*cos(10*x);
调用时, 需要在命令行上输入
>> [x, it, convg] = newton(30,'f','g')
Newton Iteration successes!!
x =
-0.412101013664971
it =
12
convg =
1
或者
>> [x, it, convg] = newton(-30,'f','g',100,1e-16)
Newton Iteration successes!!
x =
-1.41420927316939
it =
15
convg =
1

7.6 割线法#

牛顿法是一个二阶收敛的方法, 但是它要求计算导数, 在推广到多变量情形时会有很大的计算量. 割线法使用近似计算的方法代替牛顿法中的导数, 从而使算法保持较快的收敛速度, 同时计算量较小且不需要计算导数. 下面考虑单变量的情形.

假设计算非线性方程 $f(x) = 0$ 的根, 并且已有近似值 $x_{k-1}$ 和 $x_k$ . 则通过曲线 $y = f(x)$ 上的两点 $(x_{k-1}, f(x_{k-1}))$ , $(x_k, f(x_k))$ 可作一条直线, 其方程为

令该直线和 x 轴的交点横坐标为 $x_{k+1}$ ，有

再用 $x_{k}, x_{k + 1}$ 两点做同样的事情, 进行迭代, 此即为割线法. 该方法相当于在牛顿法中的切线被割线所替代. 图7-7所示为割线法的几何意义.

例 7.6.1 用割线法求方程 $x^{3} + x^{2} - 1 = 0$ 的解.

易知, 方程在区间 $[0,1]$ 中有一根, 取 $x_{0}=0, x_{1}=1$ . 割线法迭代公式为

图 7-7 割线法的几何意义

计算结果如表 7-7, 我们把牛顿法的计算结果也列出来. 可以看到, 割线法比牛顿法稍慢.

表 7-7 牛顿法和割线法的计算效果

割线法每一个计算步骤都牵涉到前面两个点, 如果没有两个初始值, 则割线法无法开始计算. 这时, 若固定其中一个点, 例如

称为单点割线法, 原来的方法相应地称为双点割线法. 单点割线法也可以认为是一种不动点迭代. 此外, 双点割线法的两个初始点不需要函数值异号, 这和二分法是不一样的.

关于双点割线法有如下的收敛性定理, 证明略.

7.7 非线性方程组简介#

考虑如下的方程组

其中, $\boldsymbol{x}=(x_{1},x_{2},\cdots,x_{n})^{\mathrm{T}}\in\mathbb{R}^{n}$ 为未知变量, $f_{i}$ 是定义在某区域 $D\subset R^{n}$ 上的 n 元实函数. 若至少有一个 $f_{i}$ 不是线性函数, 则称上述方程组为非线性方程组. 这里, 我们只考虑方程个数和自变量个数相同的情形.

和线性代数方程组类似, 非线性方程组有如下的高斯-赛德尔迭代方法:

算法 7.7.1 (非线性方程组的 GS 方法)

(1) 给定 $\boldsymbol{x}^{(0)} = (x_{1}^{(0)}, x_{2}^{(0)}, \cdots, x_{n}^{(0)})^{\mathrm{T}}, k = 0,$ 控制精度 $\varepsilon;$

(2) 对 $j=1,2,\cdots,n,$ 以 $t=x_{j}^{(k)}$ 为初始值求解如下问题:

并把其解记为 $x_{j}^{(k+1)}$ ;

(3) 若 $\| \pmb{x}^{(k + 1)} - \pmb{x}^{(k)}\| \leqslant \varepsilon$ ，则迭代收敛，方程组的近似解为 $\pmb{x}^{(k + 1)}$ ；否则， $k = k + 1$ ，转步骤 (2).

其中, 步骤 (2) 中求解的是一个单变量的非线性方程, 可以用前面的各种方法求解. 鉴于它处在内层迭代, 因此该问题求解不需要很精确, 通常用牛顿法迭代一两步就可以了. 读者也可以类似地给出雅可比方法.

非线性方程组的不动点方法也有同样的算法框架．首先引入下面的记号．记向量函数 $\boldsymbol{F}(\boldsymbol{x})=(f_{1}(\boldsymbol{x}),f_{2}(\boldsymbol{x}),\cdots,f_{n}(\boldsymbol{x}))^{\mathrm{T}}:\mathbb{R}^{n}\to\mathbb{R}^{n}$ ，其中 $f_{i}(\boldsymbol{x})=f_{i}(x_{1},x_{2},\cdots,x_{n})$ 是x的函数.若对于某个 $x\in R^{n}$ ，存在x的一个邻域 $\{z||z-x|\leqslant\varepsilon,\varepsilon>0\}\subset D$ ，则称x是D的内点.对于n元向量函数F，若存在矩阵 $\boldsymbol{A}(\boldsymbol{x})\in\mathbb{R}^{n\times n}$ ，使得

则称 F 在 x 点处可微, $A(x)$ 称为函数 F 在 x 的雅可比矩阵, 并与导数相类似, 记 $F'(x) = A(x)$ . 可以证明, F 在 x 点可微的充分必要条件是每个分量函数 $f_i(x)$ 关于每个变量 $x_j$ 可微, 此时, 我们有

若有函数 $\Phi(x): \mathbb{R}^n \to \mathbb{R}^n$ , 使得 $x = \Phi(x)$ 等价于 $F(x) = 0$ , 我们就可以构造一个向量函数的不动点迭代: $x^{(k+1)} = \Phi(x^{(k)})$ . 向量函数的不动点迭代有如下的收敛性定理 7.7.2, 也称为压缩映像定理.

非线性方程组的牛顿法可以用类似的方法导出. 设非线性函数 $F(\boldsymbol{x}) : \mathbb{R}^{n} \to \mathbb{R}^{n}$ 在区域 D 内连续可微, 且它的分量函数为 $f_{i}(\boldsymbol{x}), i = 1, 2, \cdots, n$ . 假如在某一迭代步有近似值 $\boldsymbol{x}^{(k)} = (x_{1}^{(k)}, x_{2}^{(k)}, \cdots, x_{n}^{(k)})^{\mathrm{T}} \in D$ , 做函数 $f_{i}(\boldsymbol{x})$ 在 $\boldsymbol{x}^{(k)}$ 点的泰勒展开, 并只取线性部分, 有

写成向量形式, 即为

令上式左边为零, 若 $\boldsymbol{F}^{\prime}(\boldsymbol{x}^{(k)})$ 可逆, 可从上式解出 x. 把求得的解作为新的迭代点 $\boldsymbol{x}^{(k+1)}$ :

此即为非线性方程组求根的牛顿法. 若向量函数 $F$ 在解 $\pmb{x}^*$ 的某个邻域内二阶连续可微, 则方程组的牛顿法有着相同的局部二阶收敛性质, 在此不再列出.

例 7.7.3 求解非线性方程组#

采用不动点迭代方法, 我们有迭代格式

改造成高斯-赛德尔方式, 可得格式

而牛顿法的迭代格式为

以 $\boldsymbol{x}^{(0)}=(-1,1)^{\mathrm{T}}$ 为初始点, 精度为 Matlab 中的 $eps=2.2204\times10^{-16}$ , 用上面三种不同方式迭代, 计算结果如下表 7-8 三种不同迭代方法的迭代步数分别为 47 步、23 步和 6 步. 虽然计算精度有 16 位小数, 由于篇幅关系, 我们只列出了前面的 8 位. 对于非线性方程求根问题, 不同算法或者不同的初始点都可能得到方程不同的解.

表 7-8 方程组的不同迭代效果比较

牛顿法有很好的收敛性质, 但是它在理论上要求较高, 且在计算上每一步要求计算 $n^2$ 个偏导数, 这在 $n$ 比较大时是不可接受的. 拟牛顿法就类似于单变量情形的割线法, 采用某些近似方式逼近导数, 在不精确计算导数的情形下尽量保持较快的收敛速度.

设在空间 $\mathbb{R}^n\times \mathbb{R}^n$ 中有点 $(\pmb{x}^{\prime},\pmb{F}^{\prime})$ 和点 $(\pmb{x}^{\prime \prime},\pmb{F}^{\prime \prime})$ ，其中 $F^{\prime} = F(x^{\prime}),F^{\prime \prime} = F(x^{\prime \prime}),F$ 是某个多变量的向量函数．通过 $(\pmb{x}^{\prime},\pmb{F}^{\prime})$ 的任一直线方程为 $z - F^{\prime} = B(x - x^{\prime}),(x,z)$ 是直线上点的坐标，其中 $B\in \mathbb{R}^{n\times n}$ ．若另一点 $(\pmb{x}^{\prime \prime},\pmb{F}^{\prime \prime})$ 也在该直线上，则矩阵 $B$ 满足 $F'' - F' = B(x'' - x')$ .可以看到，该方程在退化到单变量时，就相当于用割线斜率代替切线斜率.多变量情形下，这个方程就称为拟牛顿方程，我们用满足拟牛顿方程的矩阵 $B$ 来代替向量函数 $F(x)$ 的雅可比矩阵.可以说，拟牛顿方程就是 $F(x)$ 的雅可比矩阵在这两点连线方向上的近似.

设 $(\pmb{x}^{(k)},\pmb{F}^{(k)})$ ， $(\pmb{x}^{(k + 1)},\pmb{F}^{(k + 1)})$ 分别是第 $k$ 步和 $k + 1$ 步的点．当第 $k$ 步迭代完成后，就应该考虑矩阵 $\pmb{B}$ 的更新，以保证拟牛顿方程 $\pmb{F}^{(k + 1)} - \pmb{F}^{(k)} = \pmb{B}(\pmb{x}^{(k + 1)} - \pmb{x}^{(k)})$ 成立．记第 $k$ 步和第 $k + 1$ 步的拟牛顿矩阵分别为 $\pmb{B}_k$ 和 $\pmb{B}_{k + 1}$ .从 $\pmb{B}_k$ 得到新矩阵 $\pmb{B}_{k + 1}$ 的修正方式称为拟牛顿修正公式.一般地，为了记号的方便，在拟牛顿方法中，我们记 $\pmb{s}^{(k)} = \pmb{x}^{(k + 1)} - \pmb{x}^{(k)}$ 记 $\pmb{y}^{(k)} = \pmb{F}^{(k + 1)} - \pmb{F}^{(k)}$ . $\pmb{s}^{(k)}$ 也称为位移.

可以想象, 满足拟牛顿方程的矩阵可以有无穷多个. 下面我们列举出几个常用的拟牛顿修正公式. 下面的公式称为秩一修正 (Rank 1 Update) 公式:

又如, BFGS(Broyden-Fletcher-Goldfarb-Shanno) 修正公式

这样, 一个拟牛顿算法通常有以下的算法框架:

算法 7.7.4 (拟牛顿法)

(1) 给定初始值 $\boldsymbol{x}^{(0)}$ ，初始矩阵 $B_{0}$ ，控制精度 $\varepsilon, k = 0;$

(2) 若 $\| F(\pmb{x}^{(k)})\| \leqslant \varepsilon$ 或者 $k\geqslant 1$ 且 $\| \pmb{x}^{(k)} - \pmb{x}^{(k - 1)}\| \leqslant \varepsilon ,\pmb{x}^{(k)}$ 为近似解，停止迭代；

(3) 求解 $\boldsymbol{B}_{k}\boldsymbol{s}^{(k)} = -\boldsymbol{F}(\boldsymbol{x}^{(k)})$ ;

(4) 令 $\pmb{x}^{(k + 1)} = \pmb{x}^{(k)} + \lambda_k\pmb{s}^{(k)}$ ;

(5) $\boldsymbol{y}^{(k)} = \boldsymbol{F}(\boldsymbol{x}^{(k+1)}) - \boldsymbol{F}(\boldsymbol{x}^{(k)})$ ;

(6) 利用秩一修正公式或 BFGS 修正公式计算 $B_{k+1}$ ;

(7) $k = k + 1$ , 转 (2).

在上面的算法中, 第 (4) 步 $\lambda_{k}$ 可以恒取 1 或者可以每步调节, 相当于牛顿方法和牛顿下山方法, 调节的方式也和牛顿下山法类似. 而在初始步中, $B_{0}$ 一般取为单位矩阵, 这样不仅可以节省存储, 而且由于下面叙述的一些原因, 甚至可以不计算 $B_{k}$ 的逆.

下面定理给出的公式也称为低秩修正公式.

该公式可以直接用逆矩阵的定义证明. 这个公式的一个重要应用是当 $m$ 远小于 $n$ 时, 本来需要计算一个 $n$ 阶矩阵的逆, 利用这个公式则仅需要计算一个 $m$ 阶矩阵的逆.

若记拟牛顿算法中的矩阵 $B_{k}$ 的逆为 $H_{k}, B_{k+1}$ 的逆为 $H_{k+1}$ ，则利用 Sherman-Morrison 公式，可以导出从 $H_{k}$ 到 $H_{k+1}$ 的修正公式。对应于上面的秩一修正，逆矩阵的修正公式为

对应于上面的 BFGS 修正, 逆矩阵的修正公式为

这样, 上述的拟牛顿方法就可以改造成如下的方法, 从而在每一迭代步中不需去计算一个线性方程组的解.

算法 7.7.6 (改进拟牛顿法)

(1) 给定初始值 $\boldsymbol{x}^{(0)}$ ，初始矩阵 $H_{0}=I$ （单位矩阵），控制精度 $\varepsilon, k=0;$

(2) 若 $\| F(x^{(k)})\| \leqslant \varepsilon$ 或者 $k\geqslant 1$ 且 $\| x^{(k)} - x^{(k - 1)}\| \leqslant \varepsilon ,x^{(k)}$ 为近似解，停止迭代；

(3) 求解 $\boldsymbol{s}^{(k)} = -\boldsymbol{H}_{k}\boldsymbol{F}(\boldsymbol{x}^{(k)})$ ;

(4) 令 $\pmb{x}^{(k + 1)} = \pmb{x}^{(k)} + \lambda_k\pmb{s}^{(k)}$ ;

(5) $\boldsymbol{y}^{(k)} = \boldsymbol{F}(\boldsymbol{x}^{(k+1)}) - \boldsymbol{F}(\boldsymbol{x}^{(k)})$ ;

(6) 利用逆矩阵的秩一修正公式或 BFGS 修正公式计算 $H_{k+1}$ ;

(7) $k = k + 1$ , 转步骤 (2).

7.8 非线性最小二乘问题#

在 §4.4 节中, 我们给出了一个较为简单的非线性最小二乘问题. 下面我们描述一般的非线性最小二乘问题及其常用的算法. 假设有数据 $(t_{i}, y_{i}), i = 1, 2, \cdots, N$ , 我们需要找到某种特定形式的函数 $y = f(t)$ ，使之能够最好地反映数据 $(t_i, y_i)$ ，即使下面的误差达到最小：

若存在基函数 $\phi_{k}(t)$ , $k = 1, 2, \cdots, n$ ，使得 $f(t)$ 具有形式 $f(t) = x_{1}\phi_{1}(t) + x_{2}\phi_{2}(t) + \cdots + x_{n}\phi_{n}(t)$ ，其中 $x_{1}, x_{2}, \cdots, x_{n}$ 是待定参数，那么我们就得到了一个线性最小二乘问题。否则，我们通常得到一个非线性最小二乘问题。例如 $f(t) = x_{1} + tx_{2} + x_{3}e^{-x_{4}t}$ 。

假设有某一非线性最小二乘问题, 其拟合函数为 $y = f(t, x_{1}, x_{2}, \cdots, x_{n})$ , 需要拟合的数据为 $(t_{i}, y_{i}), i = 1, 2, \cdots, N$ . 我们需要求解如下的最小化问题

类似于线性最小二乘拟合, 我们可以得到如下的稳定点的条件:

记 $r_{i}(\boldsymbol{x}) = y_{i} - f(t_{i}, x_{1}, x_{2}, \cdots, x_{n})$ ，且 $\boldsymbol{r}(\boldsymbol{x}) = (r_{1}(\boldsymbol{x}), r_{2}(\boldsymbol{x}), \cdots, r_{n}(\boldsymbol{x}))^{\mathrm{T}}$ 是一个向量函数，则 $F(x) = \frac{1}{2} r(\boldsymbol{x})^{\mathrm{T}} r(\boldsymbol{x})$ 。经过计算，上面的稳定点条件可以改写成

其中， $J(x) = r^{\prime}(\pmb{x}) = \left(\frac{\partial r_{i}(\pmb{x})}{\partial x_{j}}\right)$ 是函数 $r(\pmb{x})$ 的雅可比矩阵.我们还有

如果我们计算非线性方程 $F'(x)=0$ 的根, 而已有近似值 $x^{(k)}$ , 采用牛顿法, 我们将需要计算

由于矩阵 $J(x)^{\mathrm{T}}J(x)$ 至少是半正定的, 如果忽略 $F''(x)$ 后面的一项, 例如 $r(\boldsymbol{x})$ 接近线性的情形, 我们就得到了称为高斯-牛顿算法的方法.

(1) 给定初始点 $\boldsymbol{x}^{(0)}$ ，控制精度 $\varepsilon, k = 0;$

(2) 如果 $\| \pmb{r}(\pmb{x}^{(k)})\| \leqslant \varepsilon, \| \pmb{J}(\pmb{x})^{\mathrm{T}}\pmb{r}(\pmb{x})\| \leqslant \varepsilon$ 或者 $k > 1$ 且 $\| \pmb{x}^{(k)} - \pmb{x}^{(k - 1)}\| \leqslant \varepsilon,$ 则 $\pmb{x}^{(k)}$ 为近似解，停止迭代；

(3) 求解 $\boldsymbol{J}(\boldsymbol{x}^{(k)})^{\mathrm{T}}\boldsymbol{J}(\boldsymbol{x}^{(k)})\boldsymbol{s}^{(k)} = -\boldsymbol{J}(\boldsymbol{x}^{(k)})^{\mathrm{T}}\boldsymbol{r}(\boldsymbol{x}^{(k)})$ ;

(4) $\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\lambda_{k}\boldsymbol{s}^{(k)},k=k+1$ 转(2).

$\lambda_{k}$ 可以恒取1或者每步调节, 当 $J(x^{*})$ 是一个满秩矩阵, 其中 $x^{*}$ 是解时, 可以证明, 在适当条件下, 若 $x^{(0)}$ 足够靠近 $x^{*}$ , 则高斯-牛顿方法局部收敛. 实际上, 该算法在每一步都要求矩阵 $J(x^{(k)})^{\mathrm{T}} J(x^{(k)})$ 可逆, 或者等价地, 矩阵 $J(x^{(k)})$ 列满秩.

如果该条件不成立, 或者问题本身使得 $F''(x)$ 的第二项不能够忽略, 我们通常可以采用 Levenberg-Marquardt 算法, 或者简称 LM 算法.

算法 7.8.2 (LM 算法)

(1) 给定初始点 $\boldsymbol{x}^{(0)}$ ，控制精度 $\varepsilon, k = 0;$

(2) 如果 $\| r(\pmb{x}^{(k)})\| \leqslant \varepsilon, \| \pmb{J}(\pmb{x})^{\mathrm{T}}\pmb{r}(\pmb{x})\| \leqslant \varepsilon$ 或者 $k > 1$ 且 $\| \pmb{x}^{(k)} - \pmb{x}^{(k - 1)}\| \leqslant \varepsilon,$ 则 $\pmb{x}^{(k)}$ 为近似解，停止迭代；

(3) 求解 $(\boldsymbol{J}(\boldsymbol{x}^{(k)})^{\mathrm{T}}\boldsymbol{J}(\boldsymbol{x}^{(k)}) + \mu_{k}D_{k})s^{(k)} = -\boldsymbol{J}(\boldsymbol{x}^{(k)})^{\mathrm{T}}\boldsymbol{r}(\boldsymbol{x}^{(k)})$ ;

(4) $\boldsymbol{x}^{(k+1)}=\boldsymbol{x}^{(k)}+\lambda_{k}\boldsymbol{s}^{(k)},k=k+1,$ 转 (2).

在该算法中, $D_{k}$ 是一个对角调比矩阵, $\mu_{k}$ 是一个每一步都可以调节的参数. 如果没有实行调比, 可取 $D_{k}$ 恒为单位矩阵; $\mu_{k}$ 可以取为 $\| r(x)\|$ 或者 $\| r(x)\|^2$ 的常数倍, 或者由其他方式得到, 例如信赖域方法. 关于这两者其他合适的选取, 我们不再在此展开, 请读者参考最优化方面的书籍.

7.9 大范围求解方法#

通常来讲, 一个算法的局部收敛性总是比全局收敛性容易验证, 但是我们没有很好的方法确定一个初始点怎样才算是足够靠近解, 使得我们可以确保迭代收敛. 大范围求解方法就是一种把近似解引导向精确解或者足够精确的近似解的方法. 下面介绍的同伦算法是其中的一种.

假设我们要求解非线性方程 $f(x) = 0$ ，而我们知道另一方程 $g(x) = 0$ 的根，例如，若设 $f(x) = 0$ 有近似根 $a$ ，可令 $g(x) = x - a$ 。现令函数 $h(t,x) = tf(x) + (1 - t)g(x)$ ，则有 $h(0,x) = g(x), h(1,x) = f(x)$ 。若对于任意参数 $t \in [0,1]$ ，我们把方程 $h(t,x) = 0$ 的解记为 $x(t)$ ，则在某些条件下， $x(t)$ 是 $t$ 的连续函数而 $x(0)$ 就是 $g(x) = 0$ 的根。这样，我们可以采取顺藤摸瓜的方法，令 $t$ 从 0 到 1 逐渐增大，从而求得 $x(1)$ ，也就是方程 $f(x) = 0$ 的根。

利用微分方程的方法, 可以实现这一想法, 请参见第 9 章. 设 $x = x(t)$ 是 $t$ 的连续函数, 满足 $h(t, x(t)) = 0$ . 方程两边对 $t$ 求导, 有

因此, $x'(t) = -[h_x(t,x(t))]^{-1}h_t(t,x(t))$ .若 $g(x) = 0$ 的根为 $a$ ，则我们得到一个常微分方程初值问题：

令 $x'(t) = \frac{x(t + \delta t) - x(t)}{\delta t}$ , 则有

若在 $t + \delta t$ 点处用牛顿法迭代一步, 有

但是, 该方法也有不足的地方. 在某些情况下, $x(t)$ 可能不是 $t$ 的连续函数, $x$ 和 $t$ 的对应关系可能不存在, 可能有多值, 甚至其他更复杂的情况.

例 7.9.1 考虑下面三个方程的同伦算法.

(1) 求 $f(x) = x^3 + x - 6$ 的根, 取初值 $x_0 = 1$ ;

(2) 求 $f(x) = x^2 - 6x + 8$ 的根, 取初值 $x_0 = 1$ ;

(3) 求 $f(x) = x^3 + 5x^2 + x + 5$ 的根, 取初值 $x_0 = 1$ ;

三个函数的解曲线 $x = x(t)$ 如图 7-8 所示.

图 7-8 三个方程同伦算法的解曲线

(1) 对于第一种情形, $h(t,x) = t(x^3 + x - 6) + (1 - t)(x - 1)$ . 当 $t = 0$ 时有唯一实根 $x(0) = 1$ ; 当 $0 < t \leqslant 1$ 时, $h_x'(t,x) = 3tx^2 + 1 > 0$ , 函数关于 $x$ 严格单增, 也有唯一实根. 求根曲线是一条连续曲线, 并且 $t$ 和 $x$ 单值对应.

(2) 对于第二种情形, $h(t,x) = t(x^2 - 6x + 8) + (1 - t)(x - 1)$ . 易知, $h(0,x) = 0$ 和 $h(1,x) = 0$ 有实根; 但 $h\left(\frac{1}{2}, x\right) = \frac{1}{2}\left(x^2 - 5x + 7\right) = \frac{1}{2}\left(x - \frac{5}{2}\right)^2 + \frac{3}{8} = 0$ 无实根. 求根曲线 $x = x(t)$ 在中间断开.

(3) 对于第三种情形, $h(t,x) = t(x^3 + 5x^2 + x + 5) + (1 - t)(x - 1)$ . $h(0,x) = x - 1 = 0$ 和 $h(1,x) = x^3 + 5x^2 + x + 5 = (x^2 + 1)(x - 5) = 0$ 都有唯一实根. 但当 $t = \frac{1}{5}$ 时, $h\left(\frac{1}{5}, x\right) = \frac{1}{5}(x^3 + 5x^2 + 5x + 1) = \frac{1}{5}(x + 1)(x^2 + 4x - 1)$ 有三个实根. 实际上, 求实根的曲线 $x = x(t)$ 有“迂回”的现象.

关于同伦算法如何处理这些更复杂的情形, 我们这里不再展开, 请读者参看其他书籍. 例 7.9.2 利用同伦算法计算非线性方程 (组) 的根:

(1) 求 $f(x)=x^{3}+x-6$ 的根, 取初值 $x_{0}=1;$

(2) 求

$\left\{\begin{array}{l}x_{1} + \cos x_{2} - 1 = 0,\-sin x_{1} + x_{2} - 1 = 0, \end{array} \right.$ 的根，取初值 $(-1,1)^{\mathrm{T}}$ 编写同伦算法的程序如下：function $\mathbf{x} = \text{homo}(\mathbf{f},\mathbf{g},\mathbf{x}\mathbf{0},\mathbf{dt})$ if nargin $< 4$ dt $= 0.01$ end n = length(x0); I = eye(n); told $= 0$

xold = x0;
xp = x0;
for t = dt:dt:1,
    x = xold - dt * ([told*feval(g,xold)+(1-told)*I] \ [feval(f,xold)-(xold-xp)]);
    x = x - [t*feval(g,x)+(1-t)*I] \ [t*feval(f,x)+(1-t)*(x-xp)];
    xold = x;
    told = t;
end

function v = f(x)
    v = x.^3 + x - 6;

function v = g(x)
    v = 3 * x.^2 + 1;

function v = f1(x)
    v = [x(1)+cos(x(2))-1
    -sin(x(1))+x(2)-1];

function v = g1(x)
    v = [1 -sin(x(2))
    cos(x(1)) 1];

其中, 程序后面包含了需求根的两个方程的函数及导数, 分别为 f,g 和 f1,g1. 调用过程如下:
>> x=homo('f','g',1)
x =
    1.63436529303759
>> x=homo('f1','g1',[-1 1])
x =
    0.262119524902768
    0.740871739227832

在算法中, 我们每一步都有一步牛顿法修正, 真正实现算法时, 可以求得一个较好的初始点, 然后再用牛顿方法, 而不必每一步都有牛顿法修正.

若要计算方程 $f(x) = 0$ 的根, 且初始值为 $x_0$ . 令 $h(t,x) = f(x) - \mathrm{e}^{-t}f(x_0)$ , 则 $h(0,x) = f(x) - f(x_0)$ 有根 $x_0$ ; 而当 $t \to +\infty$ 时, $\lim_{t\to +\infty}h(t,x) = f(x)$ . 因此, 同伦算法需要我们计算曲线 $x(t)$ 在 $t\to +\infty$ 的情形. 用上面微分方程的方法, 我们有

两边对 t 求导, 可得

但是, $0 = h(t,x(t)) = f(x(t)) - \mathrm{e}^{-t}f(x_0)$ , 所以

即 $x'(t) = -[f'(x(t))]^{-1} f(x(t))$ . 用近似计算的方法, 令 $x'(t) = \frac{x(t + \delta t) - x(t)}{\delta t}$ , $\delta t = 1$ , 则

上式即为牛顿方法.

评 注#

非线性方程求根和非线性最小二乘问题、非线性最优化问题等有着非常密切的联系, 牛顿法和拟牛顿法在这些问题上都是非常快速有效的算法, 尤其是在多变量的情形. 其他常用的方法在这些问题上也是可以相互借鉴的, 例如, 求解非线性最优化问题的有些被称为现代计算方法的算法——遗传算法或者模拟退火算法, 它们的基本思想可以用到非线性方程求根上.

作为非线性问题的算法, 我们总是期望在问题退化为线性时, 算法能有极其快速的收敛, 甚至有限步或一步就收敛, 这样, 当问题的非线性不是很严重时, 算法也能有非常满意的效果.

习题七#

1. 设方程 $10 - 2x - \cos x = 0$ 的迭代法为 $x = \frac{1}{2}(10 - \cos x)$ , 说明对于任意初值此迭代收敛, 并估计要求近似根具有 10 位有效数字时大约要迭代多少步.

2. 试证明, 对于任意初值 $x_0$ , 迭代格式 $x_{k+1} = \cos x_k$ 都收敛于方程 $x = \cos x$ 的同一实根.

3. 方程 $x^{3} - 2x - 2 = 0$ 在2附近有一实根, 把方程写成下面的等价的形式, 并建立相应的迭代格式,

(a) $x = \sqrt[3]{2x + 2};$

(b) $x = \frac{1}{2}(x^{3} - 2);$

(c) $x = \frac{2}{x^{2} - 2};$

试判别它们的收敛性. 选取一个最有效的格式进行计算.

4. 将 $x = \tan x$ 化为合适的迭代格式, 求解 x = 4.5 附近的根.

5. 求解下面的不动点方程, 若其迭代不收敛, 请加以改造, 给出一个收敛的迭代格式:

6. 求解下面非线性方程在区间 [2,3] 中的根, 精确到 4 位小数:

7. 若用牛顿法求解方程 $f(x)=\sin(x^{3})=0$ ，效果如何？你有什么加速的方法.

8. 求解方程 $f(x)=0$ 的 Halley 方法如下:

其中, $f_{n} = f(x_{n})$ . 说明这个公式是把牛顿法应用在 $f(x) / \sqrt{f'(x)} = 0$ 得到的. 编程实现该方法.

9. 求出多项式 $p(x)=63x^{5}-70x^{3}+15x-1$ 在区间 $[-1,1]$ 中的所有实根.

10. 用牛顿法或弦截法计算方程 $f(x) = 3x^{3} - 8x^{2} - 8x - 11 = 0$ 的某个近似根, 使误差具有精度 $10^{-4}$ .

11. 求解下面的非线性方程组, 取初值 $\boldsymbol{x}_{0} = (0.8, 0.4)^{\mathrm{T}}$ :

12. 求解下面的非线性方程组

13. 用二分法、不动点方法、牛顿法、割线法等求解下面各个问题, 列表比较各算法的性能.

(a) $x^{5}-3x-10=0;$

(b) $\sin 10x + 2\cos x - x - 3 = 0;$

(c) $x + \arctan x = 3;$

(d) $(x + 2)\ln(x^{2} + x + 1) + 1 = 0.$

数值实验七#

1. 编写牛顿方法和拟牛顿方法的程序来求解下面方程组:

2. 尝试用各种方法计算下面方程在区间 $[-10,10]$ 中的所有根:

第8章 矩阵特征值与特征向量的计算#