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

发表于 分类于

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

直接法

高斯消去法(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} \] 1. 将当前列主元 \(a\_{kk}^{(k)}\) 变为1 2. 将 \(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} \]

直接法求解的特点

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

判断矩阵是否为病态矩阵

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

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

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

迭代法

将待求解方程组改写成其等价形式\(\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迭代收敛。

迭代法的特点

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

迭代法的优点

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

解线性方程组的误差分析

范数

向量范数

注：将向量看作矩阵时，应视为\(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\)