[科学计算]非线性方程的数值求解

实际问题中，线性方程甚至线性方程组都不是最常见的形式，世界的规律是错综复杂的，涉及更多的是非线性方程，而求解非线性方程的问题就显得颇为重要。那么这一节中，我们的出发点就是如何求解非线性方程$f(x)=0$ 的根.

1. 对分区间法

相必大家记得这么一个定理：$若f(x)\in C[a,b]，且f(a)\cdot f(b)<0，则f(x)在(a,b)上必有一根。$这一句话便是“对分区间法”的依据。

操作步骤

如下图所示

1. 选取当前区间的对分点 $x_1=\frac{a+b}{2}$ ；
2. 判断$f(x_1)$的正负，若 $f(x_1)<0$则 $a_1=x_1,b_2=b$;若$f(x_1)>0$则 $a_1=a,b_2=x_1.$
3. 重复步骤1、2，直到满足条件 $|x_k-x^{\ast }|<\epsilon$ 或者 $|f(x_k)|<\epsilon$ .

需要注意的是对于不同的情况，需要合理选择不同的停机条件。

误差分析

那么通过区间对分的方法求方程的根所带来的误差究竟有多大呢？在第1步中，$|x_1-x^{\ast }|\le \frac{b-a}{2}$, 在第2步中，$|x_2-x^{\ast }|\le \frac{b-a}{2^2}$, 以此类推，可知第$k$步二分区间得到的$x_n$ 的误差为 $|x_k-x^{\ast }|\le \frac{b-a}{2^k}$.

因此，在给定误差范围 $\epsilon$ 时，我们可以计算出需要二分区间的次数 $k$ .
$$\begin{array}{r}|x_k-x^{\ast }|\le \frac{b-a}{2^k}\le \epsilon \\ k\ge \frac{\ln (b-a)-\ln \epsilon }{\ln 2}\end{array}$$

该方法的不足之处

1. 不能解决有重根的情况

2. 迭代法

$$\begin{array}{r}\text{等价变换：}f(x)=0=>x=g(x)\end{array}$$

将 $f(x)=0$ 求根问题转化为 $x=g(x)$ 的不动点问题。从一个初值 $x_0$ 出发，依次计算 $x_1=g(x_0),x_2=g(x_1),\cdots$ 得出一个序列 { x i } \{x_i\} {xi​} ,如果这个序列是收敛于 $x^{\ast }$, 即 $\underset{k\to +\infty }{\lim }{x}_k=\underset{k\to +\infty }{\lim }{x}_{k+1}=\underset{k\to +\infty }{\lim }g(x_k)=x^{\ast }$,如果此时的 $g(x)$ 是连续函数，根据连续的定义可得出 $g(x^{\ast })=x^{\ast }$,也就是说 $x^{\ast }$ 是 $g(x)$的不动点，也就是 $f(x)=0$ 的根。

那么如何保证由 $g(x)$ 构造出来的序列 { x i } \{x_i\} {xi​} 是收敛的呢？这里有一个定理来保证这件事情。

全局收敛定理： 考 虑 方 程 x = g ( x ) , g ( x ) ∈ C [ a , b ] . 若 : ( 1 ) x ∈ [ a , b ] 时 , g ( x ) ∈ [ a , b ] . ( 2 ) ∃ L , 0 ≤ L < 1 , 使 得 ∣ g ′ ( x ) ∣ ≤ L < 1. ∀ x ∈ [ a , b ] 成 立 . 则 任 取 x 0 ∈ [ a , b ] , 由 x k + 1 = g ( x k ) 产 生 的 序 列 { x i } 一 定 收 敛 于 g ( x ) 在 [ a , b ] 上 唯 一 的 不 动 点 , 并 且 有 误 差 估 计 式 ∣ x ∗ − x k ∣ ≤ 1 1 − L ∣ x k + 1 − x k ∣ ∣ x ∗ − x k ∣ ≤ L k 1 − L ∣ x 1 − x 0 ∣ 考虑方程 x=g(x),g(x) \in C[a,b].若: \\ (1) x \in [a,b]时,g(x)\in [a,b]. \\ (2)\exist L,0\le L<1,使得|g'(x)| \le L < 1.\forall x \in [a,b]成立. \\则任取x_0 \in [a,b],由x_{k+1}=g(x_k)产生的序列\{x_i\}一定收敛于g(x)在[a,b]上唯一的不动点,并且有误差估计式\\ |x^*-x_k|\le \frac{1}{1-L}|x_{k+1}-x_k| \\ |x^* - x_k| \le \frac{L^k}{1-L}|x_1-x_0| 考虑方程x=g(x),g(x)∈C[a,b].若:(1)x∈[a,b]时,g(x)∈[a,b].(2)∃L,0≤L<1,使得∣g′(x)∣≤L<1.∀x∈[a,b]成立.则任取x0​∈[a,b],由xk+1​=g(xk​)产生的序列{xi​}一定收敛于g(x)在[a,b]上唯一的不动点,并且有误差估计式∣x∗−xk​∣≤1−L1​∣xk+1​−xk​∣∣x∗−xk​∣≤1−LLk​∣x1​−x0​∣

上面的定理给出的初值 $x_0$ 的范围相对较大，我们还可以利用局部收敛性进行简化定理： 设 x ∗ 为 g ( x ) 的 不 动 点 , 并 且 g ′ ( x ) 在 x ∗ 的 某 个 邻 域 R 内 连 续 , 还 有 ∣ g ′ ( x ∗ ) ∣ < 1 , 则 任 意 x 0 ∈ R ， 迭 代 产 生 的 序 列 { x i } ∈ R 收 敛 到 x ∗ . 设x^*为g(x)的不动点,并且g'(x)在x^*的某个邻域R内连续,还有|g'(x^*)|<1, 则任意x_0\in R，迭代产生的序列\{x_i\}\in R\\收敛到x^*. 设x∗为g(x)的不动点,并且g′(x)在x∗的某个邻域R内连续,还有∣g′(x∗)∣<1,则任意x0​∈R，迭代产生的序列{xi​}∈R收敛到x∗.

迭代法的收敛阶

对于迭代过程的收敛快慢，我们可以用收敛阶来度量。什么是收敛阶呢？我们这么定义：设迭代公式 $x_{k+1}=g(x_k)$ 构造的序列 { x i } \{x_i\} {xi​} 收敛到 $g(x)$ 的不动点 $x^{\ast }$，令 $e_k=x_k-x^{\ast }$ ，若 $\underset{k\to +\infty }{\lim }\frac{|e_{k+1}|}{|e_k{|}^p}=C>0$，则称该迭代格式为 $p$阶收敛，$C$ 为渐进误差常数。

一方面，我们可以根据定义来确定收敛阶、渐进误差常数。但实际上是很难操作的，那么在平常我们该如何确定收敛阶和渐进误差常数呢？我们有下面这个定理：

$$\begin{gathered} 设x^{\ast }为g(x)的不动点,若g(x)\in C^p[B_{\delta }(x^{\ast })],p\ge 2,g^{\prime }(x^{\ast })=\cdots =g^{p-1}(x^{\ast })=0,且g^p(x^{\ast })\ne 0. \\ 则x_{k+1}=g(x_k)在B_{\delta }(x^{\ast })内p阶收敛. \end{gathered}$$
$$\begin{array}{rl}证明：令g(x_k)在x^{\ast }处作Taylor展开& \\ g(x_k)& =g(x^{\ast })+g^{\prime }(x^{\ast })(x_k-x^{\ast })+\frac{g^{(2)}(x^{\ast })}{2!}(x_k-x^{\ast }{)}^2+\cdots +\frac{g^{(p-1)}(x^{\ast })}{(p-1)!}(x_k-x^{\ast }{)}^{p-1}+\frac{g^{(p)}(x^{\ast })}{p!}(x_k-x^{\ast }{)}^p\\ x_{k+1}& =x^{\ast }+\frac{g^{(p)}(x^{\ast })}{p!}(x_k-x^{\ast }{)}^p\\ \frac{|x_{k+1}-x^{\ast }|}{|x_k-x^{\ast }{|}^p}& =\frac{g^{(p)}(x^{\ast })}{p!}=C\end{array}$$

几种常见的迭代法

可以看出，迭代法的关键在于选取合适的迭代函数g(x)。对于 $f(x)=0$ ，最简单的一种构造方法便是 $x=x+f(x)$，迭代格式为 $x_{k+1}=x_k+f(x_k)$ 。但是一般情况下，这种格式不会收敛或者收敛的很慢。

为此，我们可以对它进行迭代-加速的改进处理。令 ${\hat{x}}_{k+1}=x_k+f(x_k)$
$$\begin{array}{rl}{x}_{k+1}& =\frac{1}{1-q}{\hat{x}}_{k+1}-\frac{q}{1-q}{x}_k\\ x_{k+1}& =\frac{1}{1-q}[x_k+f(x_k)]-\frac{q}{1-q}{x}_k\\ x_{k+1}& =x_k+\frac{1}{1-q}f(x_k)\end{array}$$
更一般的形式为 $x_{k+1}=x_k+\frac{1}{M}f(x_k)$.

牛顿迭代法

简单来说，该方法就是利用Taylor公式将 $f(x)$ 在 $x^{\ast }$ 附近的某一点 $x_0$ 处进行一阶展开表示整理而来的。

$$f(x)=f(x_0)+f^{\prime }(x_0)(x-x^{\ast })+\frac{f^{\prime \prime }(\xi )}{2!}(x-x^{\ast }{)}^2,\xi 在x_0与x之间$$

忽略余项，将 $x^{\ast }$ 代入上式得 $f(x^{\ast })=f(x_0)+f^{\prime }(x_0)(x^{\ast }-x_0)=0$ ，整理可得 $x^{\ast }=x_0-\frac{f(x_0)}{f^{\prime }(x_0)}$ ， 从而牛顿迭代公式可以这样定义
$$\begin{array}{r}{x}_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_k)}\end{array}$$
我们可以从几何意义上去理解牛顿迭代法，如图

假设此时并不知道 $x^{\ast }$ ，我们有一个初始值 $x_0$ ,根据 $f(x)$ 我们可以求出 $f(x_0),f^{\prime }(x_0)$ ，图中蓝线的斜率为 $f^{\prime }(x_0)=\frac{f(x_0)}{x_0-x_1}$ ，从而可以确定 $x_1$的位置，可以这样迭代下去，直到 $|f(x_k)|<\epsilon$ .

这时，我们要问了：由牛顿迭代构造的序列 { x i } \{x_i\} {xi​} 一定收敛吗？初值是不是好选择？它收敛的值是 $g(x)$ 的不动点吗？

首先，在序列 { x i } \{x_i\} {xi​} 收敛的条件下，最后一个问题可以解决，只要 $g(x)$ 满足连续就可以了。由于 $g(x)=x-\frac{f(x)}{f^{\prime }(x)}$ ，所以只需要 $f(x)\in C^1[a,b],且f^{\prime }(x)\ne 0$。

第一个问题该怎么解决呢？这就回到了全局收敛定理，它是从不动点的角度出发讨论条件，这个条件有点难，我们换一种从根的角度讨论的思路：

收 敛 的 充 分 条 件 ： 设 f ( x ) ∈ C 2 [ a , b ] . 若 满 足 ： ( 1 ) f ( a ) f ( b ) < 0 ( 2 ) 在 [ a , b ] 上 f ′ ′ ( x ) 不 变 号 且 f ′ ( x ) ≠ 0 ( 3 ) 选 取 x 0 ∈ [ a , b ] , 需 要 满 足 f ( x 0 ) f ′ ′ ( x 0 ) > 0 则 由 牛 顿 迭 代 法 构 造 的 序 列 { x i } 收 敛 于 f ( x ) 在 [ a , b ] 上 的 唯 一 根 收敛的充分条件：设f(x)\in C^2[a,b].若满足：\\(1)f(a)f(b)<0\\(2) 在[a,b]上f''(x)不变号且f'(x) \neq 0 \\ (3) 选取x_0 \in [a,b],需要满足f(x_0)f''(x_0)>0\\ 则由牛顿迭代法构造的序列\{x_i\}收敛于f(x)在[a,b]上的唯一根 收敛的充分条件：设f(x)∈C2[a,b].若满足：(1)f(a)f(b)<0(2)在[a,b]上f′′(x)不变号且f′(x)​=0(3)选取x0​∈[a,b],需要满足f(x0​)f′′(x0​)>0则由牛顿迭代法构造的序列{xi​}收敛于f(x)在[a,b]上的唯一根

上面定理中前两个条件还好说，是一个硬要求，但是第三个条件对我们提出的第二个问题就不友好了。因此利用迭代法局部收敛的定理，做下面的改进： 局 部 收 敛 定 理 ： 设 f ( x ) ∈ C 2 [ a , b ] , 若 x ∗ 为 f ( x ) 在 [ a , b ] 上 的 根 ， 且 f ′ ( x ∗ ) ≠ 0. 则 存 在 B δ ( x ∗ ) , 使 得 ∀ x 0 ∈ B δ ( x ∗ ) , 由 牛 顿 迭 代 法 产 生 的 序 列 { x i } 收 敛 到 x ∗ , 且 满 足 lim ⁡ k → + ∞ x ∗ − x k + 1 ( x ∗ − x k ) 2 = − f ′ ′ ( x ∗ ) 2 f ′ ( x ∗ ) 局部收敛定理：设f(x)\in C^2[a,b],若x^*为f(x)在[a,b]上的根，且f'(x^*)\neq 0.则存在B_{\delta}(x^*),使得\forall x_0 \in B_{\delta}(x^*),\\由牛顿迭代法产生的序列\{x_i\}收敛到x^*,且满足\displaystyle \lim_{k \rightarrow +\infty} \frac{x^*-x_{k+1}}{(x^*-x_k)^2}=-\frac{f''(x^*)}{2f'(x^*)} 局部收敛定理：设f(x)∈C2[a,b],若x∗为f(x)在[a,b]上的根，且f′(x∗)​=0.则存在Bδ​(x∗),使得∀x0​∈Bδ​(x∗),由牛顿迭代法产生的序列{xi​}收敛到x∗,且满足k→+∞lim​(x∗−xk​)2x∗−xk+1​​=−2f′(x∗)f′′(x∗)​

在这里 $x_0$ 的选择是任意的，但是是在邻域内，这个邻域的寻找还是有点问题的。因此牛顿迭代法的收敛性速度还是依赖于 $x_0$的选择。

简化牛顿迭代法

这个的改进是针对 $f‘(x_k)$ 比较难计算而提出的。将 $f^{\prime }(x_k)$ 都替换为 $f^{\prime }(x_0)$ ，不用每次迭代都计算导数，节省大量计算难度
$$\begin{array}{r}{x}_{k+1}=x_k-\frac{f(x_k)}{f^{\prime }(x_0)}\end{array}$$
更一般的形式是不计算导数，而直接用一个常数代替
$$\begin{array}{r}{x}_{k+1}=x_k-Cf(x_k)\end{array}$$

牛顿下山法

在前面提到的方法中，并不能保证每次的迭代 $|f(x_{k+1})|<|f(x_k)|$ ，为了更好地确保迭代的收敛，我们添加一条件，使之具有单调性，即若 $|f(x_{k+1})|\ge |f(x_k)|$ ,需要在 $x_{k+1}$ 和 $x_k$ 之间重新选择 ${\hat{x}}_{k+1}$ 点 。

$$\begin{array}{rl}{\hat{x}}_{k+1}& =\lambda x_{k+1}+(1-\lambda )x_k\\ {\hat{x}}_{k+1}& =\lambda (x_k-\frac{f(x_k)}{f^{\prime }(x_k)})+(1-\lambda )x_k\\ {\hat{x}}_{k+1}& =x_k-\lambda \frac{f(x_k)}{f^{\prime }(x_k)}\end{array}$$
默认情况下 $\lambda =1$ ，如果在迭代中由 $x_k$ 得到的 $x_{k+1}$ 不能使 $|f|$ 减小 ，则将 $\lambda$ 减半 计算直到满足条件。

正割法、抛物线法

正割法是将迭代格式中的 $f^{\prime }(x_k)$ 使用差商来近似，避免求导的过程。此时 $f^{\prime }(x_k)\approx \frac{f(x_k)-f(x_{k-1})}{x_k-x_{k-1}}$ .
$$\begin{array}{rl}{x}_{k+1}& =x_k-\frac{f(x_k)}{f^{\prime }(x_k)}\\ x_{k+1}& \approx x_k-\frac{f(x_k)(x_k-x_{k-1})}{f(x_k)-f(x_{k-1})}\end{array}$$
这个时候则需要两个初始值。

抛物线法在这里暂不介绍。

问题与思考

1. 似乎这些方法对于多根情况没有考虑，迭代法在多根的情况下是否适用呢？我的想法可能是需要将多根转化为一个个局域单根的情况来解决，那样的话就涉及到另一个问题，如何确定根的个数和每一个根的局部范围？
2. 在迭代法中，重根情况下牛顿迭代法还收敛吗？收敛速度怎么样？如何加速啊？