§1
求解线性方程组 $A\mathbf{x}=\mathbf{b}$(其中 $A$ 是 $n\times n$ 非奇异矩阵,$\mathbf{b}$ 已知)是科学计算中最基本的问题之一。当 $n$ 很大或系数来自微分方程的离散化时,手工消元不可行,必须依靠数值方法。数值解法的两大分支是:
- 直接法:通过有限步代数运算(消元、分解)得到(舍入误差意义下的)精确解,典型如高斯消去法、LU 分解;
- 迭代法:从初猜 $\mathbf{x}^{(0)}$ 出发,构造序列 $\mathbf{x}^{(k)}$ 逐步逼近真解,典型如 Jacobi、Gauss-Seidel 迭代。
本章核心公式即线性系统
$$A\mathbf{x}=\mathbf{b},\qquad A\in\mathbb{R}^{n\times n},\ \mathbf{x},\mathbf{b}\in\mathbb{R}^n.$$
一、高斯消去法
高斯消去法分两步:消元把 $A$ 化为上三角矩阵 $U$,回代解出 $\mathbf{x}$。
设第 $k$ 步已把前 $k-1$ 列消为下三角零,当前主元为 $a_{kk}^{(k)}$。对 $i=k+1,\dots,n$ 计算乘数
$$l_{ik}=\frac{a_{ik}^{(k)}}{a_{kk}^{(k)}},$$
并用第 $k$ 行消去第 $i$ 行第 $k$ 列:
$$a_{ij}^{(k+1)}=a_{ij}^{(k)}-l_{ik}\,a_{kj}^{(k)},\qquad j=k,\dots,n+1,$$
(第 $n+1$ 列即增广的常数项)。消元共 $n-1$ 步,得到上三角系统
$$U\mathbf{x}=\mathbf{c}.$$
回代从最后一行起:
$$x_n=\frac{c_n}{u_{nn}},\qquad x_i=\frac{1}{u_{ii}}\Big(c_i-\sum_{j=i+1}^{n}u_{ij}x_j\Big),\ i=n-1,\dots,1.$$
运算量:消元约需 $\frac{2}{3}n^3$ 次浮点乘除法,回代约 $n^2$ 次,故总复杂度 $O(n^3/3)$。
二、部分主元策略
若主元 $a_{kk}$ 绝对值很小,除数接近于 $0$,会严重放大舍入误差甚至导致结果崩溃。为此采用部分主元(partial pivoting):在第 $k$ 步,于第 $k$ 列、第 $k$ 行及以下寻找绝对值最大的元素,将其所在行与第 $k$ 行交换,再消元。即选
$$r=\arg\max_{i\ge k}|a_{ik}^{(k)}|,\quad \text{交换第 }k\text{ 行与第 }r\text{ 行}.$$
部分主元几乎总是把主元相对大小控制在合理范围,是现代直接法的默认选择。另一种全主元还要在子矩阵中搜索,代价更高,实践中较少使用。
三、LU 分解
高斯消去本质上把 $A$ 分解为
$$A=LU,$$
其中 $L$ 是单位下三角矩阵(对角线为 $1$,下方存乘数 $l_{ik}$),$U$ 是上三角矩阵。一旦得到分解,对任意右端 $\mathbf{b}$ 只需两步三角求解:
1. 前代解 $L\mathbf{y}=\mathbf{b}$;
2. 回代解 $U\mathbf{x}=\mathbf{y}$。
这比每次重新消元高效得多——当需要用同一个 $A$ 求解多个不同的 $\mathbf{b}$(如在不同荷载、不同时间步下)时,一次分解、多次前代/回代,总计算量与一次高斯消去相当(仍约 $O(n^3/3)$ 分解 + $O(n^2)$ 每次求解)。带部分主元的版本写作 $PA=LU$,其中 $P$ 为置换矩阵。
四、病态与条件数
即便算法精确,若系数矩阵本身"接近奇异",解的舍入误差也会被放大。衡量敏感度的量是条件数
$$\kappa(A)=\|A\|\,\|A^{-1}\|$$
(常用 2-范数 $\kappa_2(A)=\sigma_{\max}/\sigma_{\min}$)。条件数越大,问题越病态:输入(或浮点舍入)的微小扰动可使输出解大幅偏离。
经典病态例子是 Hilbert 矩阵 $H_n$,其元素
$$H_{ij}=\frac{1}{i+j-1},\qquad i,j=1,\dots,n.$$
$H_n$ 的条件数随阶数 $n$ 急剧增长:
| $n$ | $\kappa_2(H_n)$(约) |
|---|---|
| 2 | $1.9\times10^{1}$ |
| 4 | $1.5\times10^{4}$ |
| 6 | $1.5\times10^{7}$ |
| 8 | $1.5\times10^{10}$ |
| 10 | $1.6\times10^{13}$ |
到 $n=10$ 时,即使双精度(约 $10^{-16}$ 相对精度)也几乎无法可靠求解。
五、迭代法:Jacobi 与 Gauss-Seidel
把 $A$ 拆成 $A=D+L+U$($D$ 对角、$L$ 严格下三角、$U$ 严格上三角)。Jacobi 迭代用上一步的全部旧值并行更新每个分量:
$$x_i^{(k+1)}=\frac{1}{a_{ii}}\Big(b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}\Big),\qquad i=1,\dots,n.$$
Gauss-Seidel 迭代在更新 $x_i$ 时立即采用已经算出的新值 $x_1^{(k+1)},\dots,x_{i-1}^{(k+1)}$:
$$x_i^{(k+1)}=\frac{1}{a_{ii}}\Big(b_i-\sum_{ji}a_{ij}x_j^{(k)}\Big).$$
收敛性:迭代可统一写成 $\mathbf{x}^{(k+1)}=T\mathbf{x}^{(k)}+\mathbf{c}$。收敛的充要条件是迭代矩阵谱半径
$$\rho(T)<1.$$
对 Jacobi 法 $T=-D^{-1}(L+U)$;Gauss-Seidel 亦有其对应的 $T$。一个重要充分条件:若 $A$ 严格对角占优(即 $|a_{ii}|>\sum_{j\ne i}|a_{ij}|$ 对所有 $i$ 成立),则两种迭代都收敛,且 Gauss-Seidel 通常比 Jacobi 收敛更快。
六、直接法与迭代法的对比
| 方面 | 直接法(高斯 / LU) | 迭代法(Jacobi / GS) |
|---|---|---|
| 适用规模 | 中小规模、稠密矩阵 | 大规模、稀疏矩阵 |
| 是否需存储全矩阵 | 是 | 只需存非零元,内存省 |
| 计算量 | 固定 $O(n^3)$ | 取决于收敛速度,难以预估 |
| 对病态的敏感 | 高(受条件数限制) | 收敛可能极慢 |
| 并行性 | 弱 | Jacobi 强,GS 较弱 |
| 典型用途 | 一般线性系统、需多右端 | 偏微分方程离散后的大型稀疏系统 |
七、例题与解答
例题1(高斯消元). 解方程组
$$\begin{cases}2x+y=3,\\ x+3y=4.\end{cases}$$
其精确解为 $(x,y)=(1,1)$。
解: 增广矩阵
$$\left[\begin{array}{cc|c}2&1&3\\1&3&4\end{array}\right].$$
第 1 步主元 $a_{11}=2$,乘数 $l_{21}=a_{21}/a_{11}=1/2=0.5$。第 2 行减去 $0.5\times$ 第 1 行:
$$\left[\begin{array}{cc|c}2&1&3\\0&2.5&2.5\end{array}\right].$$
回代:$y=2.5/2.5=1$,$x=(3-1\times1)/2=1$。故 $(x,y)=(1,1)$。✓
对应的 LU 分解为 $L=\begin{bmatrix}1&0\\0.5&1\end{bmatrix}$,$U=\begin{bmatrix}2&1\\0&2.5\end{bmatrix}$,验证 $LU=\begin{bmatrix}2&1\\1&3\end{bmatrix}=A$。
例题2(Jacobi 迭代). 用 Jacobi 法迭代上述系统,初值 $(x^{(0)},y^{(0)})=(0,0)$。
迭代格式:
$$x^{(k+1)}=\frac{3-y^{(k)}}{2},\qquad y^{(k+1)}=\frac{4-x^{(k)}}{3}.$$
- 第 1 步:$x^{(1)}=(3-0)/2=1.5$,$y^{(1)}=(4-0)/3=1.333$;
- 第 2 步:$x^{(2)}=(3-1.333)/2=0.833$,$y^{(2)}=(4-1.5)/3=0.833$;
- 第 3 步:$x^{(3)}=(3-0.833)/2=1.083$,$y^{(3)}=(4-0.833)/3=1.056$。
可见序列逐步逼近真解 $(1,1)$。该系统的迭代矩阵谱半径 $\rho(T)\approx0.2<1$,故收敛。✓
练习
1. 用高斯消元解 $\begin{cases}3x+2y=8\\ x+4y=9\end{cases}$,并写出回代过程。
2. 对例题 1 的系统,写出 Gauss-Seidel 迭代的前两步(初值取 $(0,0)$)。
3. 设 $A=\begin{bmatrix}4&1\\1&3\end{bmatrix}$,验证 $A$ 严格对角占优,并据此说明 Jacobi 与 Gauss-Seidel 迭代均收敛。
4. 为什么对 Hilbert 矩阵 $H_{10}$ 直接求解会严重失真?用条件数解释。
参考答案与提示
1. 增广矩阵消元:主元 $3$,乘数 $1/3$,得 $\begin{bmatrix}3&2&8\\0&10/3&19/3\end{bmatrix}$。回代 $y=(19/3)/(10/3)=1.9$,$x=(8-2\times1.9)/3=1.4$。解 $(1.4,1.9)$。
2. Gauss-Seidel 用新值:$x^{(1)}=(3-0)/2=1.5$,$y^{(1)}=(4-1.5)/3=0.833$;$x^{(2)}=(3-0.833)/2=1.083$,$y^{(2)}=(4-1.083)/3=0.972$。
3. 对角元 $|4|>|1|$,$|3|>|1|$,故严格对角占优;由定理两种迭代均收敛。
4. 由第四节表,$\kappa_2(H_{10})\approx1.6\times10^{13}$,条件数极大,舍入误差被放大约 $10^{13}$ 倍,远超双精度有效位数,故解严重失真。
本章小结
- 线性方程组数值解法分直接法(高斯消元、LU)与迭代法(Jacobi、Gauss-Seidel)。
- 高斯消去将 $A$ 化为上三角后回代,运算量约 $O(n^3/3)$。
- 部分主元通过换行避免小主元放大舍入误差。
- LU 分解 $A=LU$ 适合多右端问题,一次分解多次求解。
- 病态程度由条件数 $\kappa(A)$ 度量;Hilbert 矩阵是经典病态例子。
- 迭代收敛充要条件为谱半径 $\rho(T)<1$;严格对角占优保证收敛。
互动演示
拖动下方控件观察动态过程。
Jacobi 迭代解 2x2 线性方程组:迭代点逐步收敛到两直线交点(精确解 (1,1))。