跳转至

Chap 6: Direct Methods for Solving Linear Systems⚓︎

3223 个字 87 行代码 预计阅读时间 17 分钟

目标:求解 \(A \bm{x} = \bm{b}\)

Linear Systems for Equations⚓︎

Gaussian Elimination⚓︎

高斯消元法(Gaussian elimination) 的基本思路:

  • 首先将矩阵 \(A\) 归约成一个上三角(upper-triangular) 矩阵
  • 然后通过回代(backward-substitution) 来求解未知量

先来看消元(elimination) 的实现:先令 \(A^{(1)} = A = (a_{ij}^{(1)})_{n \times n}, \bm{b}^{(1)} = \bm{b} = \begin{bmatrix}b_1^{(1)} \\ \vdots \\ b_n^{(1)}\end{bmatrix}\)

  • 1 步:

    • 如果 \(a_{11}^{(1)} \ne 0\),计算 \(m_{i1} = \dfrac{a_{i1}^{(1)}}{a_{11}^{(1)}}, (i = 2, \dots, n)\)
    • 那么增广矩阵 (augmented matrix) 的第 \(i\) \(\text{row}_i\) 为:\(m_{i1} \times \text{row}_1\),得到

      \[ \left[ \begin{array}{cccc|c} a_{11}^{(1)} & a_{12}^{(1)} & \cdots & a_{1n}^{(1)} & b_1^{(1)} \\ 0 & & A^{(2)} & & \bm{b}^{(2)} \\ \end{array} \right] \]

      其中 \(\begin{cases}a_{ij}^{(2)} = a_{ij}^{(1)} - m_{i1} a_{1j}^{1} \\ b_i^{(2)} = b_i^{(1)} - m_{i1} b_1^{(1)}\end{cases}, (i, j = 2, \dots, n)\)

  • k 步:

    • 如果 \(a_{kk}^{(k)} \ne 0\),计算 \(m_{ik} = \dfrac{a_{ik}^{(k)}}{a_{kk}^{(k)}}, (i = k+1, \dots, n)\)
    • \(\begin{cases}a_{ij}^{(k+1)} = a_{ij}^{(k)} - m_{ik} a_{kj}^{k} \\ b_i^{(k+1)} = b_i^{(k)} - m_{ik} b_k^{(k)}\end{cases}, (i, j = k+1, \dots, n)\)
  • n-1 步后:

    \[ \begin{bmatrix} a_{11}^{(1)} & a_{12}^{(1)} & \dots & a_{1n}^{(1)} \\ & a_{22}^{(2)} & \dots & a_{2n}^{(2)} \\ & & \dots & \vdots \\ & & & a_{nn}^{(n)}\end{bmatrix} \begin{bmatrix}x_1 \\ x_2 \\ \vdots \\ x_n\end{bmatrix} = \begin{bmatrix}b_1^{(1)} \\ b_2^{(2)} \\ \vdots \\ b_n^{(n)}\end{bmatrix} \]

接下来看回代

  • \(x_n = \dfrac{b_n^{(n)}}{a_{nn}^{(n)}}\)
  • \(x_i = \dfrac{b_i^{(i)} - \sum\limits_{j=i+1}^n a_{ij}^{(i)} x_j}{a_{ii}^{(i)}}, (i = n - 1, \dots, 1)\)
  • 我们必须找到最小的整数 \(k \ge i\) \(a_{ki}^{(i)} \ne 0\),然后交换第 \(k\) 行和第 \(i\)
代码实现

求解 \(n \times n\) 线性方程组:

\[ \begin{align} E_1: & a_{11} x_1 + a_{12} x_2 + \dots + a_{1n} x_n = a_{1, n+1} \notag \\ E_2: & a_{21} x_1 + a_{22} x_2 + \dots + a_{2n} x_n = a_{2, n+1} \notag \\ \vdots \quad & \quad \vdots \quad \quad \quad\ \vdots \quad \quad \quad \quad \quad \quad \vdots \quad \quad \quad \vdots \notag \\ E_n: & a_{n1} x_1 + a_{n2} x_2 + \dots + a_{nn} x_n = a_{n, n+1} \notag \\ \end{align} \]
  • 输入:未知量和方程的数量 \(n\);增广矩阵 \(A = (a_{ij})\),其中 \(1 \le n \le n, 1 \le j \le n+1\)
  • 输出:解 \(x_1, x_2, \dots, x_n\),或者线性方程组没有唯一解的信息
Step 1  for i =1, ..., n - 1 do Steps 2-4:
    Step 2  Let p be the smallest integer with i <= p <= n and a_pi != 0;
            if no integer p can be found
                then Output('no unique solution exists');
                STOP;
    Step 3  if p != i then perform (E_p) <-> (E_i);
    Step 4  for j = i + 1, ..., n do Step 5 and 6:
        Step 5  Set m_ji = a_ji / a_ii
        Step 6  Perform (E_j - m_ji * E_i) -> (E_i)
    Step 7  if a_nn = 0 then Output('no unique solution exists');
                STOP;
    Step 8  Set x_n = a_{n,n+1} / a_nn    // start backward substitution
    Step 9  for i = n - 1, ..., 1 set x_i = [a_{i, n+1} - sum(j=i+1, n, a_ij * x_j)] / a_ii;
    Step 10 Output(x_1, ..., x_n);
            STOP;    // success

Amount of Computation⚓︎

现在我们来统计一下计算量(仅考虑乘法 / 除法

  • 消元:\(\sum\limits_{k=1}^{n-1} (n - k)(n - k + 2) = \dfrac{n^3}{3} + \dfrac{n^2}{2} - \dfrac{5}{6}n\)
  • 回代:\(1 + \sum\limits_{i=1}^{n-1}(n - i + 1) = \dfrac{n^2}{2} + \dfrac{n}{2}\)

所以对于很大的 \(n\),乘法和除法的总数大约为 \(\textcolor{red}{\dfrac{n^3}{3}}\)。也就是说,高斯消元法的时间复杂度为 \(O(n^3)\)

Pivoting Strategies⚓︎

除非有特殊说明,以下的 \(k\) 指的是第 \(k\) 次高斯消元。

一般高斯消元法的问题

在高斯消元的过程中,如果其中一个主元(pivot) \(a_{kk}^{(k)} = 0\),那么就需要进行行交换 \((E_k) \leftrightarrow (E_p)\),其中 \(p\) 是最小的满足 \(p > k\) \(a_{pk}^{(k)} \ne 0\) 的整数。但为了减小舍入误差,即使主元并不等于 0 的时候也要做行交换。

如果 \(a_{kk}^{(k)}\) 相比 \(a_{jk}^{(k)}\) 较小的话,那么乘数 \(m_{jk} = \dfrac{a_{jk}^{(k)}}{a_{kk}^{(k)}}\) 会大于 1,导致误差的积累。并且在回代的时候,\(x_k\) 的值因为分母 \(a_{kk}^{(k)}\) 的值过小,其计算误差也会被放大。

所以,我们需要选取合适的主元以减小误差。下面给出一些选取主元的策略。

Partial Pivoting⚓︎

部分主元法(partial pivoting)(或称为最大列主元法 (maximal column pivoting):找到最小的 \(p\),使得 \(|a_{pk}^{(k)}| = \max\limits_{k \le i \le n} |a_{ik}^{(k)}|\),然后交换第 \(p\) 行和第 \(k\) 行。

代码实现

求解 \(n \times n\) 线性方程组:

\[ \begin{align} E_1: & a_{11} x_1 + a_{12} x_2 + \dots + a_{1n} x_n = a_{1, n+1} \notag \\ E_2: & a_{21} x_1 + a_{22} x_2 + \dots + a_{2n} x_n = a_{2, n+1} \notag \\ \vdots \quad & \quad \vdots \quad \quad \quad\ \vdots \quad \quad \quad \quad \quad \quad \vdots \quad \quad \quad \vdots \notag \\ E_n: & a_{n1} x_1 + a_{n2} x_2 + \dots + a_{nn} x_n = a_{n, n+1} \notag \\ \end{align} \]
  • 输入:未知量和方程的数量 \(n\);增广矩阵 \(A = (a_{ij})\),其中 \(1 \le n \le n, 1 \le j \le n+1\)
  • 输出:解 \(x_1, x_2, \dots, x_n\),或者线性方程组没有唯一解的信息
Step 1  for i = 1, ..., n set NROW(i) = i;    // initialize row pointer
Step 2  for i = 1, ..., n - 1 do Steps 3-6:    // elimination process
    Step 3  Let p be the smallest integer with i <= p <= n and
            |a(NROW(p), i)| = max_{i <= j <= n}|a(NROW(j), i)|;
    Step 4  if a(NROW(p), i) = 0 then Output('no unique solution exists');
                STOP;
    Step 5  if NROW(i) != NROW(p) then set NCOPY = NROW(i);    // simulated row interchange
                                          NROW(i) = NROW(p);
                                          NROW(p) = NCOPY;
    Step 6  for j = i + 1, ..., n do Steps 7-8:
        Step 7  Set m(NROW(j), i) = a(NROW(j), i) / a(NROW(i), i);
        Step 8  Perform (E_NROW(j) - m(NROW(j), i) * E_NROW(i)) -> (E_NROW(j));
    Step 9  if a(NROW(n), n) = 0 then Output('no unique solution exists');
                STOP;
    // start backward substitution
    Step 10 Set x_n = a(NROW(n), n + 1) / a(NROW(n), n);
    Step 11 for i = n - 1, ..., 1
                set x_i = (a(NROW(i), n + 1) - sum(j=i+1, n, a(NROW(i), j) * x_j)) / a(NROW(i), i);
    Step 12 Output(x_1, ..., x_n);    // procedure completed successfully
            STOP;
例子

求解线性方程组 \(\begin{cases}30.00 x_1 + 594100 x_2 = 591700 \\ 5.291 x_1 - 6.130 x_2 = 46.78\end{cases}\),舍入精度为 4 位。

见教材 \(P_{363}\)

Scaled Partial Pivoting⚓︎

缩放部分主元法(scaled partial pivoting)(或称为缩放列主元法 (scaled-column pivoting):将一行中最大的元素放在主元的位置上。

  • 1 步:定义每一行的缩放因子 \(s_i = \max\limits_{1 \le j \le n} |a_{ij}|\)
  • 2 (对于第 \(k\) 次高斯消元)找到最小的 \(p \ge k\),使得 \(\dfrac{|a_{pk}^{(k)}|}{s_p} = \max\limits_{k \le i \le n} \dfrac{|a_{ik}^{(k)}|}{s_i}\),然后交换第 \(p\) 行和第 \(k\)

缩放因子只计算一次(在高斯消元前,以确保计算效率。

代码实现

和部分主元法的实现相比,区别在于前 3 步,后续步骤都是一样的,因此下面只列出前 3 步:

Step 1  for i = 1, ..., n set s_i = max_{1 <= j <= n}(|a_ij|);
                          if s_i = 0 then Output('no unique solution exists');
                                          STOP;
Step 2  for i = 1, ..., n - 1 do Steps 3-6:    // elimination process
    Step 3  Let p be the smallest integer with i <= p <= n and
            |a(NROW(p), i)| / s(NROW(p)) = max_{i <= j <= n}(|a(NROW(j), i)| / s(NROW(j)));

Complete Pivoting⚓︎

完全主元法(complete pivoting)(或称为最大主元法 (maximal pivoting):搜索所有的元素 \(a_{ij} (i, j = k, \dots, n)\),找出其中数值最大的元素。通过互换(interchange) 行和列,使得该元素来到主元的位置上。

Amount of Computation⚓︎

  • 部分主元法:需要 \(O(n^2)\) 次额外的比较
  • 缩放部分主元法:需要 \(O(n^2)\) 次额外的比较,以及 \(O(n^2)\) 除法
  • 完全主元法:需要 \(O(\dfrac{n^3}{3})\) 次额外的比较

如果新的缩放因子在行交换的时候才被确定,那么缩放部分主元法需要 \(O(\dfrac{n^3}{3})\) 次额外的比较,以及 \(O(n^2)\) 除法

Matrix Factorization⚓︎

高斯消元法是一种简单粗暴的方法,但效率不是很高。因此下面介绍一种基于高斯消元法实现的改进方法——矩阵分解(matrix factorization)。它的计算过程如下:

  • 1 步:

    • \(m_{i1} = \dfrac{a_{i1}}{a_{11}} (a_{11} \ne 0)\)
    • \(L_1 = \begin{bmatrix}1 & & & \\ -m_{21} & 1 & & \\ \vdots & & \ddots \\ -m_{n1} & & &1\end{bmatrix}\)(这就是第一高斯变换矩阵(first Gaussian transformation matrix),那么 \(L_1 [A^{(1)} \quad \bm{b}^{(1)}] = \begin{bmatrix}a_{11}^{(1)} & \dots a_{1n}^{(1)} & b_1^{(1)} \\ 0 & A^{(2)} & \bm{b}^{(2)}\end{bmatrix}\)
  • k 步:

    • k 高斯变换矩阵(kth Gaussian transformation matrix):\(L_k = \begin{bmatrix}1 & & & & \\ & \ddots & & & \\ & & 1 & & \\ & & -m_{k+1, k} & & \\ & & \vdots &\ddots & \\ & & -m_{n, k} & & 1\end{bmatrix}\)(空的地方都是 0
  • n-1 步:

    \[ L_{n-1}L_{n-2} \dots L_1 [A \quad \bm{b}] = \begin{bmatrix}a_{11}^{(1)} & a_{12}^{(1)} & \dots & a_{1n}^{(1)} & b_1^{(1)}\\ & a_{22}^{(2)} & \dots & a_{2n}^{(2)} & b_2^{(2)} \\ & & \dots & \vdots & \vdots\\ & & & a_{nn}^{(n)} & b_n^{(n)}\end{bmatrix} \]

定理

若高斯消元法能够在不使用行互换的基础上求解线性方程组 \(A \bm{x} = \bm{b}\),那么矩阵 \(A\) 可以被因式分解为一个下三角矩阵 \(L\) 和上三角矩阵 \(U\) 的乘积,即 \(A = LU\)

其中:

\[ U = \begin{bmatrix}a_{11}^{(1)} & a_{12}^{(1)} & \dots & a_{1n}^{(1)} \\ 0 & a_{22}^{(2)} & \dots & a_{2n}^{(2)} \\ \vdots & \vdots & \ddots & \vdots \\ 0 & \dots & \dots & a_{nn}^{(n)}\end{bmatrix} \quad L = \begin{bmatrix}1 & 0 & \dots & 0 \\ m_{21} & 1 & \dots & 0 \\ \vdots & \vdots & \ddots & \vdots \\ m_{n1} & \dots & m_{n, n-1} & 0\end{bmatrix} \quad m_{ji} = \dfrac{a_{ji}^{(i)}}{a_{ii}^{(i)}} \]

如果矩阵 \(L\) unitary(也就是说主对角线元素都是 1,那么得到的矩阵分解是唯一的

唯一性证明

如果分解不是唯一的,那么存在 \(L_1, L_2, U_1, U_2\),使得 \(A = L_1 U_1 = L_2 U_2\),因此:

\[ \underbrace{U_1U_2^{-1}}_{\text{Upper-triangular}} = L_1^{-1} L_2 U_2 U_2^{-1} = \underbrace{L_1^{-1} L_2}_{\substack{\text{Lower-triangular} \\ \text{with diagnoal entries } \textcolor{red}{1}}} = I \]

如果 \(U\) 也是 unitary 的,那么这种分解就称为Crout 分解Crout 分解能够通过对 \(A^T\) \(LU\) 分解。也就是说,找到 \(A^T = LU\),那么 \(A = U^T L^T\) 就是 \(A\) Crout 分解。

\(LU\) 分解的代码实现

\(n \times n\) 的矩阵 \(A = (a_{ij})\) 分解为下三角矩阵 \(L = (l_{ij})\) 和上三角矩阵 \(U = (u_{ij})\),也就是说 \(A = LU\),其中 \(L\) \(U\) 的主对角线元素均为 1

  • 输入:维度 \(n\)\(A\) 的元素 \(a_{ij}, 1 \le i, j \le n\)\(L\) 的对角元素 \(l_{11} = \dots = l_{nn} = 1\) \(U\) 的对角元素 \(u_{11} = \dots = u_{nn} = 1\)
  • 输出:\(L\) 的项 \(l_{ij}, 1 \le j \le i, 1 \le i \le n\),以及 \(U\) 的项 \(u_{ij}, i \le j \le n, 1 \le i \le n\)
Step 1  Select l_11 and u_11 satisfying l_11 * u_11 = a_11;
        if l_11 * u_11 = 0 then Output('Factorization impossible');
            Stop;
Step 2  for j = 2, ..., n set u_1j = a_1j / l_11;    // first row of U
                            l_j1 = a_j1 / u_11;    // first column of L
Step 3  for i = 2, ..., n - 1 do Steps 4 and 5:
    Step 4  Select l_ii and u_ii safisfying l_ii * u_ii = a_ii - sum(k=1, i-1, l_ik * u_ki);
            if l_ii * u_ii = 0 then Output('Factorization impossible');
                Stop;
    Step 5  for j = i + 1, ..., n:
                set u_ij = 1 / l_ii * (a_ij - sum(k=1, i-1, l_ik * u_kj));    // ith row of U
                set l_ji = 1 / u_ii * (a_ji - sum(k=1, i-1, l_jk * u_ki));    // ith column of L
Step 6  Select l_nn and u_nn satisfying l_nn * u_nn = a_nn - sum(k=1, n-1, l_nk * u_kn);
// if l_nn * u_nn = 0, then A = LU but A is singular
Step 7  Output(l_ij for j = 1, ..., i and i = 1, ..., n);
        Output(u_ij for j = i, ..., n and i = 1, ..., n);
        STOP;

这样我们只得到了 \(L\) \(U\),还没有解出这个线性方程组(矩阵,所以还要继续往下转换。

因为 \(A \bm{x} = LU \bm{x} = \bm{b}\),所以我们令 \(\bm{y} = U\bm{x}\)

  • 先解 \(L \bm{y} = \bm{b}\)
    • \(y_1 = \dfrac{b_1}{l_{11}}\)
    • \(y_i = \dfrac{1}{l_{ii}} \Big[b_i - \sum\limits_{j=1}^{i-1} u_{ij} y_j\Big]\)
  • 再解 \(U \bm{x} = \bm{y}\)
    • \(x_n = \dfrac{y_n}{u_{nn}}\)
    • \(x_i = \dfrac{1}{u_{ii}} \Big[y_i - \sum\limits_{j=i+1}^n u_{ij} x_j\Big]\)

Special Types of Matrices⚓︎

Strictly Diagonally Dominant Matrix⚓︎

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

\[ |a_{ii}| > \sum\limits_{\substack{j = 1 \\ j \ne i}}^n |a_{ij}| \]

定理

  • 严格对角占优矩阵 \(A\) 非奇异的(nonsigular)(行列式不为 0,且存在逆矩阵)
  • 在这种矩阵上使用高斯消元法无需行或列的互换
  • 并且计算将相对于舍入误差的增长保持稳定
证明
  • \(A\) 是非奇异的——反证法证明
  • 高斯消元法无需行或列的互换——归纳法证明:通过高斯消元法得到的每一个矩阵 \(A^{(2)}, A^{(3)}, \dots, A^{(n)}\) 都是严格对角占优的
  • 略过

Choleski's Method for Positive Definite Matrix⚓︎

定义

对于一个矩阵 \(A\),如果它是对称的,且 \(\forall \bm{x} \ne \bm{0}, \bm{x}^T A \bm{x} > 0\) 成立,那么称该矩阵是正定(positive definite) 矩阵。

定理

正定矩阵 \(A\) 的性质:

a. \(A\) 是非奇异的

b. \(a_{ii} > 0, i = 1, 2, \dots n\)

c. \(\max_{1 \le k, j \le n} |a_{kj}| \le \max_{1 \le i \le n} |a_{ii}|\)

d. \((a_{ij})^2 < a_{ii} a_{jj}, i \ne j\)

PPT 上还有这些性质:

  • \(A^{-1}\) 也是正定的
  • \(A\) 的每个前导主子矩阵(leading principal submatrices) \(A_k\) 的行列式 (determinant) 都是正的
前导主子矩阵

矩阵 \(A\) 的前导主子矩阵为 \(A_k = \begin{bmatrix}a_{11} & a_{12} & \dots & a_{1k} \\ a_{21} & a_{22} & \dots & a_{2k} \\ \vdots & \vdots & \ddots & \vdots \\ a_{k1} & a_{k2} & \dots \\ a_{kk}\end{bmatrix}, 1 \le k \le n\)

证明


我们将 \(A = LU\) 中的 \(U\) 进一步拆分成对角矩阵 \(D\) 和单位上三角矩阵 \(\widetilde{U}\)

可以推导出:\(A\) 是对称矩阵(\(A = A^T \rightarrow LU = LD \widetilde{U} = \widetilde{U^T} DL^T\)\(\Rightarrow L = \widetilde{U}^T \Rightarrow A = LDL^T\)。这样我们得到了另一种矩阵分解—— \(LDL^T\) 分解:

\(LDL^T\) 分解的代码实现

\(n \times n\) 的矩阵 \(A = (a_{ij})\) 分解为 \(LDL^T\) 的形式,其中 \(L\) 是下三角矩阵,对角线元素均为 1\(D\) 为对角矩阵,对角线上的元素均为正数。

  • 输入:维度 \(n\)\(A\) 的元素 \(a_{ij}, 1 \le i, j \le n\)
  • 输出:\(L\) 的项 \(l_{ij}, 1 \le j \le i, 1 \le i \le n\),以及 \(D\) 的项 \(u_{i}, 1 \le i \le n\)
Step 1  for i = 1, ..., n do Steps 2-4:
    Step 2  for j = 1, ..., i - 1, set v_j = l_ij * d_j;
    Step 3  Set d_i = a_ii - sum(j=1, i-1, l_ij * v_j);
    Step 4  for j = i + 1, ..., n, set l_ji = (a_ji - sum(k=1, i-1, l_jk * v_k)) / d_i;

\(D^{\frac{1}{2}} = \begin{bmatrix}\sqrt{u_{11}} & & & \\ & \sqrt{u_{22}} & & \\ & & & & \\ & & & \sqrt{u_{nn}}\end{bmatrix}\)\(\widetilde{L} = LD^{\frac{1}{2}}\) 仍然是一个上三角矩阵,因此 \(A = \widetilde{L} \widetilde{L}^T\)

综上,若 \(A\) 是正定矩阵,那么:

  • \(L\) 是一个对角线元素均为 1 的下三角矩阵,并且 \(D\) 是一个对角项均为正数的对角矩阵时,\(A\) 可被分解为 \(LDL^T\)
  • \(L\) 是一个对角线上均为非零元素的下三角矩阵时,\(A\) 可被分解为 \(LL^T\)
算法:Choleski

目标:将规模为 \(n \times n\) 的对称的正定矩阵 \(A\) 分解为 \(LL^T\),其中 \(L\) 是下三角矩阵。

输入:\(n\) 维矩阵 \(A\),其元素为 \(a_{ij}, 1 \le i, j \le n\)

输出:矩阵 \(L\),其元素为 \(l_{ij}, 1 \le j \le i, 1 \le i \le n\)

Step 1  set l_11 = sqrt(a_11);
Step 2  for j = 2, ..., n, set l_j1 = a_j1 / l_11;
Step 3  for i = 2, ..., n - 1 do steps 4 and 5
    Step 4  set l_ii = sqrt(a_ii - sum(pow(l_ik, 2), 1, i - 1))
    // LDL^T is faster, but must be modified to solve Ax = b

    Step 5  for j = i + 1, ..., n, set l_ji = (a_ji - sum(l_jk * l_ik, 1, i - 1)) / l_ii;
Step 6  set l_nn = sqrt(a_nn - sum(pow(l_nk, 2), 1, n - 1))
Step 7  output (l_ij for j = 1, ..., i and i = 1, ..., n);
Stop.

Crout Reduction for Tridiagonal Linear System⚓︎

带状矩阵

对于一个 \(n \times n\) 的矩阵,如果有整数 \(p, q\),满足 \(1 < p, q < n\),当 \(i + p \le j\) \(j + q \le i\) 时,有 \(a_{ij} = 0\),那么称该矩阵为带状矩阵(band matrix),其带宽(bandwidth) \(w = p + q - 1\)

\(p = q = 2\) 时,\(w = 3\),此时的矩阵称为三对角矩阵(tridiagonal matrix),其形式如下:

\[ A = \begin{bmatrix}a_{11} & a_{12} & 0 & \dots & \dots & 0 \\ a_{21} & a_{22} & a_{23} & \dots & \dots & \vdots \\ 0 & a_{32} & a_{33} & \dots & \dots & \vdots \\ \vdots & \ddots & \ddots & \ddots & \ddots & a_{n-1, n} \\ 0 & \dots & \dots & 0 & a_{n, n-1} & a_{nn}\end{bmatrix} \]

对于上述形式的线性方程组(\(A \bm{x} = \bm{f}\),我们采用一种特殊的 \(LU\) 分解,称为 Crout 分解。下面给出具体求解步骤:

  1. 寻找矩阵 \(A\) Crout 分解,\(L, U\) 分别为:

    \[ L = \begin{bmatrix}l_{11} & 0 & \dots & \dots & 0 \\ l_{21} & l_{22} & \ddots & & \vdots \\ 0 & \ddots & \ddots & \ddots & 0 \\ 0 & \dots & 0 & l_{n, n-1} & l_{nn}\end{bmatrix} \quad U = \begin{bmatrix}1 & u_{12} & 0 & \dots & 0 \\ 0 & 1 & \ddots & \ddots & \vdots \\ \vdots & & \ddots & \ddots & u_{n-1, n} \\ 0 & \dots & \dots & 0 & 1\end{bmatrix} \]
  2. 求解 \(L\bm{y} = \bm{f} \Rightarrow y_1 = \dfrac{f_1}{l_{11}}, y_i = \dfrac{(f_i - l_{i, i-1} y_{i-1})}{l_{ii}} (i = 2, \dots, n)\)

  3. 求解 \(U\bm{x} = \bm{y} \Rightarrow x_n = y_n, x_i = y_i - u_{i,i+1} x_{i+1}\)

定理

如果 \(A\) 是三对角线矩阵,且是对角线占优的,并满足 \(|b_1| > |c_1| > 0, |b_n| > |a_n| > 0, a_i \ne 0, c_i \ne 0\),那么 \(A\) 是非奇异的,对应的线性方程组有解。

  • 如果 \(A\) 是严格对角占优的,那么没有必要让所有的 \(a_i, b_i, c_i\) 都是非零的
  • 该方法是稳定的,因为所有从计算过程中获得的值会被约束在原有元素的范围内
  • 计算量为 \(O(n)\)
代码实现

求解 \(n \times n\) 的线性方程组:

\[ \begin{align} E_1: a_{11} x_1 + a_{12} x_2 \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad & = a_{1, n+1} \notag \\ E_2: a_{21} x_1 + a_{22} x_2 + a_{23} x_3 \quad \quad \quad \quad \quad \quad \quad \quad \quad \ \ & = a_{2, n+1} \notag \\ \vdots \quad \quad \quad \quad \quad \quad \quad \quad \quad \vdots \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \ \ & \ \ \vdots \notag \\ E_{n-1}: \quad a_{n-1, n-2} x_{n-2} + a_{n-1, n-1} x_{n-1} + a_{n-1, n} x_n & = a_{n-1, n+1} \notag \\ E_{n}: \quad \quad \quad \quad \quad \quad \quad \quad \quad \quad \ a_{n, n-1} x_{n-1} + a_{n, n} x_n & = a_{n, n+1} \notag \\ \end{align} \]

假设这个线性方程组有唯一解。

  • 输入:维度 \(n\)\(A\) 的元素
  • 输出:解 \(x_1, \dots, x_n\)
Step 1  Set l_11 = a_11;
            u_12 = a_12 / l_11;
            z_1 = a_{1,n+1} / l_11;
Step 2  for i = 2, ..., n-1 set l_{i, i-1} = a_{i, i-1};    // ith row of L
                                l_ii = a_ii - l_{i, i-1} * u_{i-1, i};
                                u_{i, i+1} = a_{i, i+1} / l_ii;    // (i+1)th column of U
                                z_i = (a_{i, n+1} - l_{i, i-1} * z_{i-1}) / l_ii;
Step 3  Set l_{n, n-1} = a_{n, n-1};    // nth row of L
            l_nn = a_nn - l_{n, n-1} * u_{n-1, n};
            z_n = (a_{n, n+1} - l_{n, n-1} * z_{n-1}) / l_nn;
// Step 4 and 5 solve Ux = z
Step 4  x_n = z_n;
Step 5  for i = n-1, ..., 1 set x_i = z_i - u_{i, i+1} * x_{i+1};
Step 6  Output(x_1, ..., x_n);
        STOP;

对应的作业练习📝

评论区

如果大家有什么问题或想法,欢迎在下方留言~