计算方法-[3]线性方程组的数值解法

发表于 2019-05-08 分类于计算方法

计算方法第三章线性方程的数值解法课程内容整理，尚未包含代码部分。

直接法

高斯消去法(Gaussian Elimination)

$$
\mathbf{Ax} = \mathbf{b} \Rightarrow
\begin{bmatrix}
a_{11} & \cdots & a_{1n} \\
& \ddots & \vdots \\
& & a_{nn}^{(n-1)}
\end{bmatrix}\mathbf{x} =
\begin{bmatrix}
b_1 \\
\vdots \\
b_n^{(n-1)}
\end{bmatrix}
$$

消去：

$$
l_{ik} = a_{ik}^{(k)}/a_{kk}^{(k)} \
$$
$$
\begin{cases}
a_{ij}^{(k+1)} = a_{ij}^{(k)} - l_{ik} a_{kj}^{(k)} \\
b_{i}^{(k+1)} = b_{i}^{(k)} - l_{ik} b_{k}^{(k)}
\end{cases}
(k = 1,2，\cdots, n - 1 , \quad i，j = k + 1, k + 2, \cdots, n)
$$

回代：

$$
\begin{cases}
x_n = b_n^{(n)} / a^{(n)}_{nn}\\
x_i = (b_i^{(i)} - \sum\limits^n_{j = i+1})/a_{ii}^{(i)}
\end{cases}
\quad i = n - 1, n - 2, \cdots, 1
$$

复杂度： $O(n^3)$

列主元消去法

主元素：$a_{kk}^{(k)}$ (作为除数的元素)

使用高斯消去法时，若除数很小，除法计算时将产生较大的舍入误差，影响后续消元。

列主元消去法即改进的高斯消去法，在选取主元素时按列选取待消去的矩阵元素中 绝对值 最大的元素与当前行交换，成为新的主元素。即：
$$
if \quad \max_{k \le i\le n}|a_{i_k}| \ne 0 \quad and \quad i_k \ne k \\
then \quad i_k \leftrightarrow k
$$

相比于高斯消去法，列主元消去法更稳定，但行列交换需花费较多时间，计算量高于高斯消去法 $O(n^2/3)$。

高斯-若当消去法（Gauss-Jordan Method)

$$
\mathbf{Ax} = \mathbf{b} \Rightarrow
\begin{bmatrix}
a_{11} & & \\
& \ddots &\\
& & a_{nn}^{(n-1)}
\end{bmatrix}\mathbf{x} =
\begin{bmatrix}
b_1 \\
\vdots \\
b_n^{(n-1)}
\end{bmatrix}
$$

将当前列主元 $a_{kk}^{(k)}$ 变为1
将 $a_{kk}^{(k)}$ 所在列的上下元素变为0

计算量比高斯消去法多$(n^3-3n^2 + 2n)/6$ 次乘除。常用于求逆矩阵。

高斯若当法求逆矩阵

$$
[A|E] \Rightarrow[E| A^{-1}]
$$

LU 分解法

$n$阶矩阵$A$,若其所有顺序主子式均不为0，则可将其分解为一个单位下三角矩阵$L$ 和一个上三角矩阵$U$的乘积，且这种分解唯一。
$$
\mathbf{A} =
\begin{bmatrix}
1 & & & \\
l_{21} & 1 & & \\
\vdots & \vdots & \ddots & \\
l_{n1} & l_{n2} & \cdots & 1
\end{bmatrix}
\begin{bmatrix}
u_{11} & u_{12} & \cdots & u_{1n} \\
& u_{22} & \cdots & u_{2n} \\
& & \ddots & \vdots \\
& & & u_{nn}
\end{bmatrix}
$$

求解过程：
$$
\mathbf{A}^{(k)} = \mathbf{L}_{k-1} \mathbf{A}^{(k-1)}\\
\mathbf{U} = \mathbf{A}^{(n)} = \mathbf{L}_{n-1}\mathbf{A}^{(k-1)}\\
\mathbf{L}_k =
\begin{bmatrix}
1 & & & \\
& 1 & & \\
& -l_{k+1, k} & & \\
& \vdots & \ddots\\
& -l_{n, k} & 0 & 1
\end{bmatrix},
\quad l_{ik} = \frac{a_{ik}^{(k)}}{a_{kk}^{(k)}}
\
$$
回代：
$$
\begin{cases}
y_i = b_i - \sum \limits_{k = 1}^{i-1} k_{ik}y_k, \quad i = 1, 2, \cdots, n \\
x_i = y_i - \sum \limits^{n}_{k = i + 1}u_{ik}x_{k}/u_{ii}, \quad i = n, n-1, \cdots, 1
\end{cases}
$$

直接法求解的特点

经过可预先确定的有限次算术运算 求出精确解；
由于有舍入误差，只能得到近似解；
需要对解进行误差分析；
不适合用于求解病态方程组；
一般适合于解矩阵 $A$ 为低阶稠密矩阵的方程组。

判断矩阵是否为病态矩阵

判断矩阵是否病态，常由经验得出：

行列式很大或很小（如某些行、列近似相关）；
元素间相差大数量级，且无规则；
主元消去过程中出现小主元；
特征值相差大数量级。

使用直接法求解病态矩阵的方程组可能使解严重偏离真实值。

迭代法

将待求解方程组改写成其等价形式$\mathbf{x = Bx+g}$ ，以迭代序列的形式逼近方程的根。

$\mathbf{L}$：下三角矩阵；

$\mathbf{U}$：上三角矩阵；

$\mathbf{D}$：对角矩阵。

雅可比迭代法(Jacobi Iterative Methods)

$$
\mathbf{x}^{(k+1)} = \mathbf{B}_J \mathbf{x}^{(k)} + \mathbf{g}_J \\
\begin{cases}
\mathbf{B}_J = -\mathbf{D}^{-1}(\mathbf{L+U}) \\
\mathbf{g}_J = \mathbf{D}^{-1} \mathbf{b}
\end{cases}
$$

高斯-塞尔德(Gauss-Seidel Iterative Methods)

$$
\mathbf{x}^{(k+1)} = \mathbf{B}_S \mathbf{x}^{(k)} + \mathbf{g}_S\\
\begin{cases}
\mathbf{B}_S = -\mathbf{(D+L)}^{-1}\mathbf{U} \\
\mathbf{g}_S = \mathbf{(D+L)}^{-1} \mathbf{b}
\end{cases}
$$

迭代法的收敛性

迭代法收敛的充要条件是谱半径小于1。
$$
\rho(\mathbf{B}) < 1
$$
谱半径(spectral radius)：$\rho(\mathbf{B}) = \max\limits_{1\le i\le n}|\lambda_i|$，$\lambda\ _i$ 为$\mathbf{B}$的特征值。

特殊的迭代矩阵

若$\mathbf{A}$ 为严格对对角占优矩阵，则方程组的Jacobi和 Gauss-Seidel 迭代均收敛。

严格对角占优矩阵 (strictly diagonally dominant matrix)

$\mathbf{A}$的每一行对角线元素绝对值都严格大于同行其他元素绝对值之和：
$$
|a_{ii}|>\sum\limits^n_{j = 1\\j\ne i}|a_{ij}| \quad (i = 1, 2, \cdots, n)
$$

若$\mathbf{A}$ 不可约，即不能通过出初等行变换变为上三角矩阵的形式，且按行弱对角占优，则方程组的 Jacobi 和 Gauss-Seidel迭代均收敛。

弱对角占优；
$$
|a_{ii}|\ge\sum\limits^n_{j = 1\\j\ne i}|a_{ij}| \quad (i = 1, 2, \cdots , n)
$$
且至少有一个$i$使得大于号成立，则$\mathbf{A}$为弱对角占优矩阵。

若$\mathbf{A}$ 为对称正定矩阵，则 Gauss-Seidel迭代收敛。

迭代法的特点

构造出合适的迭代格式可加快收敛速度；
不容易构造能收敛到解的迭代函数；

迭代法的优点

舍入误差影响小；
对于高阶方程组计算量小于直接法；
适合于求解系数矩阵问题

解线性方程组的误差分析

范数

向量范数

注：将向量看作矩阵时，应视为$n\times1$的列向量。

$\mathbf{R}^n$ 空间的向量范数$||\cdot||$ 对任意$\mathbf{x, y} \in \mathbf{R}^n$ 满足：

正定性(positive definite)： $||\mathbf{x}||\ge 0; \quad ||\mathbf{x}|| = 0 \Leftrightarrow \mathbf{x} = 0$
齐次性(homogeneous)：$||\alpha \mathbf{x}|| = |\alpha| \cdot ||\mathbf{x}||$
三角不等式(triangle inequality)：$||\mathbf{x +y}|| \le||\mathbf{x}|| + ||\mathbf{y}||$

常用向量范数：

$||x||_1 = \sum \limits_{i = 1}^n |x_i|$
$||x||_2 = \sqrt{\sum \limits_{i = 1}^n | x_i| ^2}$
$||x||_p = (\sum\limits_{i = 1}^n |x_i|^p)^{1/p}$
$||x||_{\infty} = \max\limits_{1\le i \le n}|x_i|$

矩阵范数

$\mathbf{R}^{m\times n}$ 空间的矩阵范数$||\cdot||$ 对任意 $\mathbf{A,B} \in \mathbf{R}^{m\times n}$ 满足：

正定性：$||\mathbf{A}|| \ge 0; \quad ||\mathbf{A}|| = 0 \Leftrightarrow \mathbf{A} = 0$
齐次性：$||\alpha \mathbf{A}|| = |\alpha| \cdot ||\mathbf{A}||$
三角不等式：$||\mathbf{A+B} || \le || \mathbf{A}|| + || \mathbf{B}||$
相容(consistent)（当$m = n$ 时）：$||\mathbf{AB}||\le ||\mathbf{A}|| \cdot || \mathbf{B}||$

常用矩阵范数：

行和范数： $||\mathbf{A}||_{\infty} = \max \limits_{1\le i \le n} \sum \limits_{j=1}^n |a_{ij}$
列和范数：$||\mathbf{A}||_1 = \max \limits_{1\le j\le n} \sum \limits_{i = 1}^n |a_{ij}|$
谱范数(spectral norm)：$||\mathbf{A}||_2 = \sqrt{\rho(\mathbf{A^TA)}}$

条件数

矩阵$\mathbf{A}$的条件数记为$cond(\mathbf{A})$，越大则矩阵越病态，越难得到准确解。
$$
cond(\mathbf{A}) = ||\mathbf{A}|| \cdot ||A^{-1}||
$$
性质：

若$\mathbf{A}$ 对称，则$cond(\mathbf{A})_2 = \frac{max|\lambda|}{min|\lambda|}$
若$\mathbf{A}$ 可逆，则$cond(\mathbf{A})_p \ge 1$
若$\mathbf{A}$ 可逆， $\alpha \in R$，则$cond(\alpha \mathbf{A}) = cond(\mathbf{A})$
若$\mathbf{A}$ 正交，则$cond(\mathbf{A})_2 = 1$
若$\mathbf{A}$ 可逆， $\mathbf{R}$ 正交，则 $cond(\mathbf{RA})_2 = cond(\mathbf{AR})_2 = cond(\mathbf{A})_2$