【笔记】多种方法求解方程及方程组
断清秋
·
·
算法·理论
不动点法解非线性方程
对于非线性方程 f(x)=0，可以通过变形将其变为等价的 g(x)=x 形式，其中满足 g(x^*)=x^* 的 x^* 称为 g 的不动点。显然，由于两个方程等价，x^* 即为原方程 f(x)=0 的解。
那么问题转化为求解 x^*。
我们通过以下方式求解：
- 猜测一个近似值 x_0；
- 将 x_0 代入得到 x_1=g(x_0)；
- 再将 x_1 代入原式，以此类推，最终得到 x_{n+1}=g(x_n)。
- 如果数列 \left\{x_n\right\} 收敛于某个值 x^*，那么两边取极限可知 x^* 即为我们所求的答案。
当然，这种做法依赖于迭代函数 g(x) 的收敛性。具体地，迭代函数收敛的充分条件为：在不动点 x^* 的邻域内，迭代函数 g(x) 的导数的绝对值必须小于 1，即 |g^{\prime}(x)|<1。（显然，绝对值越小收敛越快）
总体来说，这种方法的步骤依次为：
- 构造迭代函数 g(x)；
- 判断迭代函数 g(x) 在不动点邻域内的收敛性；
- 迭代求方程解的近似值。
例题
证明：当 x_0 = 1.5 时，迭代法 x_{k+1}=\sqrt{\frac{10}{4+x_k}} 收敛于方程 f(x) = x^3 + 4x^2 - 10 = 0 在区间 [1,2] 内的唯一实根 x^*。
用上述不动点迭代法求满足精度要求 |x_k - x_{k-1}| \leq 10^{-5} 的近似根。运算过程中注意有效数字位数。
解答
先证明迭代法收敛。
考虑 x=g(x)=\sqrt{\frac{10}{4+x}}，化简得 x^3+4x^2-10=0。
考虑对于 x \in [1,2]，显然有 \sqrt{\frac{5}{3}} \le g(x) \le \sqrt{2}，满足 g(x) \in [1,2]。
对 g(x) 求导，g'(x)=-\frac{\sqrt{10}}{2(x+4)^{\frac{3}{2}}}，在 [1,2] 上满足 \max_{x \in [1,2]} |g'(x)|=|g'(1)|=\frac{\sqrt{2}}{10}<1，综上可知题中迭代法收敛于 f(x) 在 [1,2] 上唯一实根。
接下来按迭代法求解近似根，取 x_0=1.5。
x_1=\sqrt{\frac{10}{4+x_0}}=\sqrt{\frac{10}{5.5}} \approx 1.348400
x_2=\sqrt{\frac{10}{4+x_1}}=\sqrt{\frac{10}{5.348400}} \approx 1.367376
x_3=\sqrt{\frac{10}{4+x_2}}=\sqrt{\frac{10}{5.367376}} \approx 1.364957
x_4=\sqrt{\frac{10}{4+x_3}}=\sqrt{\frac{10}{5.364957}} \approx 1.365265
x_5=\sqrt{\frac{10}{4+x_4}}=\sqrt{\frac{10}{5.365265}} \approx 1.365226
x_6=\sqrt{\frac{10}{4+x_5}}=\sqrt{\frac{10}{5.365226}} \approx 1.365231
代入 $x_6=1.365231$，有 $x_6^3+4x_6^2-10=0.00001629 \approx 0$.
因此，原方程的近似根是 $1.365231$。
## Newton 切线迭代法解非线性方程
Newton 切线迭代法本质就是选取特定构造函数 $g(x)=x-\dfrac{f(x)}{f^{\prime}(x)}$ 的不动点法。这是一个十分巧妙的构造，在初值 $x_0$ 选取接近 $x^*$ 时保证了 $g(x)$ 的收敛性，且速度更快（平方收敛）。
### 例题
当 $x_0 = 1.5$ 时，用 Newton 切线迭代法求方程 $f(x) = x^3 + 4x^2 - 10 = 0$ 方程根的近似值，满足精度要求 $|x_k - x_{k-1}| \leq 10^{-5}$。
### 解答
设 $f(x)=x^3+4x^2-10$，则 $f'(x)=3x^2+8x$。
牛顿迭代公式：$x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}$。
$x_1=x_0-\frac{f(x_0)}{f'(x_0)}=1.373333
x_2=x_1-\frac{f(x_1)}{f'(x_1)}=1.365262
x_3=x_2-\frac{f(x_2)}{f'(x_2)}=1.365230
x_4=x_3-\frac{f(x_3)}{f'(x_3)}=1.36522998... \approx 1.365230
## Doolittle LU 分解法解线性方程组
Doolittle LU 分解法是在线性代数中求解线性方程组 $A\mathbf{x}=\mathbf{b}$ 的一种常用方法。具体地，它将系数矩阵 $A$ 分解成两个矩阵的乘积：单位下三角矩阵 $L$ 和普通上三角矩阵 $U$。
通过这种分解方式，可以将原方程转化为 $LU\mathbf{x}=\mathbf{b}$。如果令 $U\mathbf{x}=\mathbf{y}$，那么只需依次解 $L\mathbf{y}=\mathbf{b}$ 和 $U\mathbf{x}=\mathbf{y}$ 两个方程组即可。利用 $U,L$ 两个矩阵的性质，这是容易的。
可以看出，对于系数矩阵 $A$ 相同但 $\mathbf{b}$ 不同的多个线性方程组，先进行 LU 分解可以大大降低时间复杂度。
那么难点在于如何求解 $L,U$ 矩阵。
先写出一般形式，假设 $A$ 是 $n$ 阶方阵，那么有
$$
A = LU =
\begin{bmatrix}
1 & 0 & \cdots & 0 \\
l_{21} & 1 & \cdots & 0 \\
\vdots & \vdots & \ddots & \vdots \\
l_{n1} & l_{n2} & \cdots & 1
\end{bmatrix}
\begin{bmatrix}
u_{11} & u_{12} & \cdots & u_{1n} \\
0 & u_{22} & \cdots & u_{2n} \\
\vdots & \vdots & \ddots & \vdots \\
0 & 0 & \cdots & u_{nn}
\end{bmatrix}
$$
一般地，我们通过以下步骤计算：
对于 $k = 1, 2, \dots, n$ 依次执行以下两步：
1. 计算 $U$ 的第 $k$ 行元素 $(j = k, k+1, \dots, n)
u_{kj} = a_{kj} - \sum\limits_{t=1}^{k-1} l_{kt} u_{tj}
- 计算 L 的第 k 列元素 (i = k+1,k+2, \dots, n)
l_{ik} = \frac{a_{ik} - \sum_{t=1}^{k-1} l_{it} u_{tk}}{u_{kk}}
简单来说，就是先算 U 的第 1 行，再算 L 的第 1 列，再算 U 的第 2 行……以此类推。
显然地，要满足 u_{kk} \neq 0。
例题
使用 Doolittle LU 分解法解方程组
4x_1 + 2x_2 - 2x_3 = 10 \\
2x_1 + 2x_2 - 3x_3 = 5 \\
-2x_1 -3x_2 + 14x_3 = 4
\end{cases}
解答
A = \begin{bmatrix} 4 & 2 & -2 \\ 2 & 2 & -3 \\ -2 & -3 & 14 \end{bmatrix},
\mathbf{b} = \begin{bmatrix} 10 \\ 5 \\ 4 \end{bmatrix}
求解 U 的第 1 行：u_{11}=a_{11}=4,u_{12}=a_{12}=2,u_{13}=a_{13}=-2；
求解 L 的第 1 列：l_{11}=1,l_{21} = \dfrac{a_{21}}{u_{11}} = \dfrac{1}{2},l_{31} = \dfrac{a_{31}}{u_{11}} = -\dfrac{1}{2}；
求解 U 的第 2 行：u_{21}=0,u_{22}=a_{22}-l_{21}u_{12}=1,u_{23}=a_{23}-l_{21}u_{13}=-2；
求解 L 的第 2 列：l_{12}=0,l_{22} = 1,l_{32} = \dfrac{a_{32}-l_{31}u_{12}}{u_{22}} = -2；
求解 U 的第 3 行：u_{31}=0,u_{32}=0,u_{33}=a_{33}-(l_{31}u_{13} + l_{32}u_{23})=9；
求解 L 的第 3 列：l_{13}=0,l_{23} = 0,l_{33} = 1；
可以求出 L,U 矩阵如下：
L = \begin{bmatrix}
1 & 0 & 0 \\
\frac{1}{2} & 1 & 0 \\
-\frac{1}{2} & -2 & 1
\end{bmatrix},
\quad
U = \begin{bmatrix}
4 & 2 & -2 \\
0 & 1 & -2 \\
0 & 0 & 9
\end{bmatrix}
解 L\mathbf{y} = \mathbf{b}，即求$
\begin{bmatrix}
1 & 0 & 0 \
\frac{1}{2} & 1 & 0 \
-\frac{1}{2} & -2 & 1
\end{bmatrix}
\begin{bmatrix} y_1 \ y_2 \ y_3 \end{bmatrix}=\begin{bmatrix} 10 \ 5 \ 4 \end{bmatrix}
解得 $y_1=10,y_2=0,y_3=9$，故 $\mathbf{y}=\begin{bmatrix} 10 \\ 0 \\ 9 \end{bmatrix}$。
解 $U\mathbf{x} = \mathbf{y}$，即求
\begin{bmatrix}
4 & 2 & -2 \
0 & 1 & -2 \
0 & 0 & 9
\end{bmatrix}
\begin{bmatrix} x_1 \ x_2 \ x_3 \end{bmatrix}=\begin{bmatrix} 10 \ 0 \ 9 \end{bmatrix}
解得 $x_1=2,x_2=2,x_3=1$，写成向量形式即 $\mathbf{x}=\begin{bmatrix}2 \\ 2 \\ 1 \end{bmatrix}$。
## Gauss-Seidel 迭代法解线性方程组
与 Doolittle 分解法不同，Gauss-Seidel 迭代法不是通过有限的步骤算出精确解，而是从一个猜测的初始解出发，通过不断迭代更新，让解一步步逼近真实值，直到误差满足要求。
对于 $n$ 元线性方程组 $A\mathbf{x}=\mathbf{b}$， 考虑将第 $i$ 个方程变形，将 $x_i$ 孤立在等号左边，把其他的 $x_j$ 移到等号右边，可以得到如下式子：$ x_i = \dfrac{1}{a_{ii}} \left( b_i - \sum\limits_{j \neq i} a_{ij} x_j \right)$。
这样就可以得到一个迭代公式，能够通过旧的 $x$ 值算出新的 $x$ 值。
为了保证收敛，通常要求系数矩阵 $A$ 是**严格对角占优**的。
即对于每一行，对角线元素的绝对值要大于该行其他所有元素绝对值之和： $ |a_{ii}| > \sum_{j \neq i} |a_{ij}| $。(如果原矩阵不满足条件，可以尝试通过交换行来达成目的)
可以发现，对于第 $x_i$ 项的计算，我们应当使用本轮迭代已经算出的 $x_1,\ldots,x_{i-1}$ 的值和上一轮迭代中算出的 $x_{i+1},\ldots,x_n$ 的值。这样我们可以更快地应用本轮迭代得到的新信息，较传统的 Jacobi 迭代法收敛更快。
由上述结论，对于 $i = 1, 2, \dots, n$，可以将迭代方程改写为：
$$
x_i^{(k+1)} = \frac{1}{a_{ii}} \left( b_i - \sum_{j=1}^{i-1} a_{ij} x_j^{(k+1)} - \sum_{j=i+1}^{n} a_{ij} x_j^{(k)} \right)
$$
一般地，我们需要设定初始猜测向量 $x^{(0)}$（第 0 轮迭代向量），通常设为全 0 向量 $[0, 0, \dots, 0]^T$。
除此以外，还要设定允许误差 $\epsilon$。我们只需要不断迭代直到满足停止条件 $ ||x^{(k+1)} - x^{(k)}|| < \epsilon $，即可得到要求误差范围内的近似解。
### 例题
用Gauss-Seidel迭代法解方程组
$\begin{cases}
5x_1 + 2x_2 + x_3 = -12 \\
-x_1 + 4x_2 + 2x_3 = 20 \\
2x_1 - 3x_2 + 10x_3 = 3
\end{cases}
初值从 (0, 0, 0)^T 开始迭代，要求当 \|\mathbf{x}^{(k)} - \mathbf{x}^{(k-1)}\|_\infty < 10^{-3} 时终止计算。(\|\mathbf{v}\|_\infty 表示向量 \mathbf{v} 中绝对值最大的元素)
解答
一、写出迭代函数
将原方程组改写为 Gauss-Seidel 迭代格式：
x_1^{(k)} = \frac{-12 - 2x_2^{(k-1)} - x_3^{(k-1)}}{5}
x_2^{(k)} = \frac{20 + x_1^{(k)} - 2x_3^{(k-1)}}{4}
x_3^{(k)} = \frac{3 - 2x_1^{(k)} + 3x_2^{(k)}}{10}
其中 k 为迭代次数，初值 \mathbf{x}^{(0)} = (0, 0, 0)^T。
二、判断迭代法是否收敛
系数矩阵为：A = \begin{bmatrix}
5 & 2 & 1 \\
-1 & 4 & 2 \\
2 & -3 & 10
\end{bmatrix}
判断是否严格对角占优：
第一行：|5| > |2| + |1|
第二行：|4| > |-1| + |2|
第三行：|10| > |2| + |-3|
因此，系数矩阵严格对角占优，Gauss-Seidel 迭代法收敛。
三、迭代求解过程
从初值 \mathbf{x}^{(0)} = (0, 0, 0) 开始迭代：
第 1 次迭代：
x_1^{(1)} = \frac{-12 - 2 \times 0 - 0}{5} = -2.4
x_2^{(1)} = \frac{20 + (-2.4) - 2 \times 0}{4} = 4.4
x_3^{(1)} = \frac{3 - 2 \times (-2.4) + 3 \times 4.4}{10} = 2.1
\mathbf{x}^{(1)} = (-2.4, 4.4, 2.1)$，$\|\mathbf{x}^{(1)} - \mathbf{x}^{(0)}\|_{\infty} = 4.4 > 10^{-3}
第 2 次迭代：
x_1^{(2)} = \frac{-12 - 2 \times 4.4 - 2.1}{5} = -4.58
x_2^{(2)} = \frac{20 + (-4.58) - 2 \times 2.1}{4} = 2.805
x_3^{(2)} = \frac{3 - 2 \times (-4.58) + 3 \times 2.805}{10} = 2.0575
\mathbf{x}^{(2)} = (-4.58, 2.805, 2.0575)$，$\|\mathbf{x}^{(2)} - \mathbf{x}^{(1)}\|_{\infty} = 2.18 > 10^{-3}
第 3 次迭代：
x_1^{(3)} = \frac{-12 - 2 \times 2.805 - 2.0575}{5} = -3.9335
x_2^{(3)} = \frac{20 + (-3.9335) - 2 \times 2.0575}{4} = 2.987875
x_3^{(3)} = \frac{3 - 2 \times (-3.9335) + 3 \times 2.987875}{10} = 1.9830625
\mathbf{x}^{(3)} = (-3.9335, 2.987875, 1.9830625)$，$\|\mathbf{x}^{(3)} - \mathbf{x}^{(2)}\|_{\infty} = 0.6465 > 10^{-3}
第 4 次迭代：
x_1^{(4)} = \frac{-12 - 2 \times 2.987875 - 1.9830625}{5} = -3.9917625
x_2^{(4)} = \frac{20 + (-3.9917625) - 2 \times 1.9830625}{4} = 3.010528125
x_3^{(4)} = \frac{3 - 2 \times (-3.9917625) + 3 \times 3.010528125}{10} = 2.0015109375
\mathbf{x}^{(4)} = (-3.9917625, 3.010528125, 2.0015109375)$，$\|\mathbf{x}^{(4)} - \mathbf{x}^{(3)}\|_{\infty} = 0.0582625 > 10^{-3}
第 5 次迭代：
x_1^{(5)} = \frac{-12 - 2 \times 3.010528125 - 2.0015109375}{5} = -4.0045134375
x_2^{(5)} = \frac{20 + (-4.0045134375) - 2 \times 2.0015109375}{4} = 2.998116171875
x_3^{(5)} = \frac{3 - 2 \times (-4.0045134375) + 3 \times 2.998116171875}{10} = 2.0003375390625
\mathbf{x}^{(5)} = (-4.0045134375, 2.998116171875, 2.0003375390625)$，$\|\mathbf{x}^{(5)} - \mathbf{x}^{(4)}\|_{\infty} = 0.0127509375 > 10^{-3}
第 6 次迭代：
x_1^{(6)} = \frac{-12 - 2 \times 2.998116171875 - 2.0003375390625}{5} = -3.9993139765625
x_2^{(6)} = \frac{20 + (-3.9993139765625) - 2 \times 2.0003375390625}{4} = 3.000002736328125
x_3^{(6)} = \frac{3 - 2 \times (-3.9993139765625) + 3 \times 3.000002736328125}{10} = 1.9998636162109375
\mathbf{x}^{(6)} = (-3.9993139765625, 3.000002736328125, 1.9998636162109375)$，$\|\mathbf{x}^{(6)} - \mathbf{x}^{(5)}\|_{\infty} = 0.0051994609375 > 10^{-3}
第 7 次迭代：
x_1^{(7)} = \frac{-12 - 2 \times 3.000002736328125 - 1.9998636162109375}{5} = -3.9999738177734375
x_2^{(7)} = \frac{20 + (-3.9999738177734375) - 2 \times 1.9998636162109375}{4} = 3.000074737451171875
x_3^{(7)} = \frac{3 - 2 \times (-3.9999738177734375) + 3 \times 3.000074737451171875}{10} = 1.9999171847900390625
\mathbf{x}^{(7)} = (-3.9999738177734375, 3.000074737451171875, 1.9999171847900390625)$，$\|\mathbf{x}^{(7)} - \mathbf{x}^{(6)}\|_{\infty} = 0.0006598412109375 < 10^{-3}
迭代在 k=7 时满足终止条件。
故一共迭代了 7 次，方程组的解为：\mathbf{x} \approx (-4.000, 3.000, 2.000)。
直接解方程验证可知原方程的精确解为 \mathbf{x} = (-4, 3, 2)^T，证明迭代无误。