数值计算方法 Chapter5. 解线性方程组的直接法本文已参与「新人创作礼」活动，一起开启掘金创作之路。数值计算方法

本文已参与「新人创作礼」活动，一起开启掘金创作之路。

数值计算方法 Chapter5. 解线性方程组的直接法

0. 问题描述

这一章节考察的就是如何求解线性方程组：

\left\{ \begin{aligned} a_{11}x_1 + a_{12}x_2 + ... + a_{1n}x_n &= b_1 \\ a_{21}x_1 + a_{22}x_2 + ... + a_{2n}x_n &= b_2 \\ ... \\ a_{n1}x_1 + a_{n2}x_2 + ... + a_{nn}x_n &= b_n \end{aligned} \right.

或者可以用矩阵来表达：

\begin{pmatrix} a_{11} & a_{12} & ... & a_{1n} \\ a_{21} & a_{22} & ... & a_{2n} \\ ... \\ a_{n1} & a_{n2} & ... & a_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ ... \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ ... \\ b_n \end{pmatrix}

1. 消元法

1. 三角方程组

首先，我们来考察一些特殊形式的方程：

1. 对角方程组

对角方程组的函数形式如下：

\begin{pmatrix} a_{11} & & & \\ & a_{22} & & \\ & & ... & \\ & & & a_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ ... \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ ... \\ b_n \end{pmatrix}

显然，可以直接写出 $x_i = b_i / a_{ii}$ ，其中， $i \in [1, n]$ 。

2. 下三角方程组

下面，我们考察一下一个稍微复杂一点的情况，即下三角矩阵的情况：

\begin{pmatrix} a_{11} & & & \\ a_{21} & a_{22} & & \\ & & ... & \\ a_{n1} & a_{n2} & ... & a_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ ... \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ ... \\ b_n \end{pmatrix}

此时，问题的难度会多少复杂一点，但是也同样可以直接给出解答如下：

x_i = \frac{b_i - \sum_{j < i}a_{ij}x_j}{a_ii}

我们可以给出python伪代码如下：

def down_triangle(A, b):
    n = len(A)
    x = [0 for _ in range(n)]
    for i in range(n):
        y = b[i]
        for j in range(i):
            y -= A[i][j] * x[j]
        x[i] = y / A[i][i]
    return x

3. 上三角方程组

同样的，对于上三角函数的情况，我们同样有：

\begin{pmatrix} a_{11} & a_{12} & ... & a_{1n} \\ & a_{22} & ... & a_{2n} \\ & & ... & \\ & & & a_{nn} \end{pmatrix} \begin{pmatrix} x_1 \\ x_2 \\ ... \\ x_n \end{pmatrix} = \begin{pmatrix} b_1 \\ b_2 \\ ... \\ b_n \end{pmatrix}

可以解得：

x_i = \frac{b_i - \sum_{j > i}a_{ij}x_j}{a_ii}

同样的，我们可以给出python伪代码如下：

def up_triangle(A, b):
    n = len(A)
    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        y = b[i]
        for j in range(i+1, n):
            y -= A[i][j] * x[j]
        x[i] = y / A[i][i]
    return x

2. Gauss消元法

现在，我们来考察一下一般形式的多元线性方程的解法。

其核心的思路其实还是将其转换成三角矩阵然后进行求解。

我们同样给出python伪代码如下：

def gauss_elimination(A, b):
    n = len(A)
    for i in range(n-1):
        for k in range(i+1, n):
            c = A[k][i] / A[i][i]
            for j in range(i, n):
                A[k][j] -= c * A[i][j]
            b[k] -= c * b[i]
    return up_triangle(A, b)

3. Gauss-Jordan消元法

Gauss-Jordan消元法和上述Gauss消元法本质上是一样的，不过Gauss消元法是将一般矩阵转换成三角矩阵，而Gauss-Jordan消元法是将一般矩阵转换成对角矩阵。

同样的，我们给出python伪代码如下：

def gauss_jordan_elimination(A, b):
    n = len(A)
    for i in range(n-1):
        for k in range(i+1, n):
            c = A[k][i] / A[i][i]
            for j in range(i, n):
                A[k][j] -= c * A[i][j]
            b[k] -= c * b[i]
    
    for j in range(n-1, 0, -1):
        for i in range(j):
            c = A[i][j] / A[j][j]
            A[i][j] = 0
            b[i] -= c * b[j]

    res = [b[i] / A[i][i] for i in range(n)]
    return res

2. 直接分解法

直接分解法和上述消元法其实并没有本质上的不同，不过区别在于，直接分解法的核心思路在于基于三角阵的特异性从而不断地尝试将一般矩阵转换为三角阵的形式然后进行求解。

1. Dolittle分解

Dolittle分解的思路是说将一般 $n$ 阶矩阵 $A$ 转换成 $A = LU$ ，其中， $L$ 是一个对角元素为1的下三角矩阵，而 $U$ 是一个上三角矩阵。

即：

\begin{pmatrix} a_{11} & a_{12} & ... & a_{1n} \\ a_{21} & a_{22} & ... & a_{2n} \\ & & ... & \\ a_{n1} & a_{n2} & ... & a_{nn} \end{pmatrix} = \begin{pmatrix} 1 & & & \\ l_{21} & 1 & & \\ & ... & & \\ l_{n1} & l_{n2} & ... & 1 \end{pmatrix} \begin{pmatrix} u_{11} & u_{12} & ... & u_{1n} \\ & u_{22} & ... & u_{2n} \\ & & ... & \\ & & & u_{nn} \end{pmatrix}

而两个三角矩阵 $L$ 和 $U$ 都是相对比较容易求解的。

我们易有：

\left\{ \begin{aligned} u_{ij} &= a_{ij} - \sum_{k=1}^{i-1}l_{ik}u_{kj} \\ l_{ij} &= (a_{ij} - \sum_{k=1}^{j-1}l_{ik}u_{kj}) / u_{jj} \end{aligned} \right.

此时，对于待求解方程 $Ax = b$ 即可变为：

$LUx = b$

更进一步的，我们可以将其拆分为两个方程：

\left\{ \begin{aligned} Ly &= b \\ Ux &= y \end{aligned} \right.

分别解上述两个方程，即可得到最终的解 $x$ ：

\left\{ \begin{aligned} y_i &= b_i - \sum_{j=1}^{i-1} l_{ij}y_{j} \\ x_i &= (y_{i} - \sum_{j=i+1}^{n} u_{ij}x_j) / u_{ii} \end{aligned} \right.

同样的，我们可以给出python伪代码如下：

def dolittle_solve(A, b):
    n = len(A)
    L = [[0 for _ in range(n)] for _ in range(n)]
    U = [[0 for _ in range(n)] for _ in range(n)]

    for k in range(n):
        L[k][k] = 1
        for j in range(k, n):
            U[k][j] = A[k][j]
            for r in range(k):
                U[k][j] -= L[k][r] * U[r][j]
        for i in range(k+1, n):
            L[i][k] = A[i][k]
            for r in range(k):
                L[i][k] -= L[i][r] * U[r][k]
            L[i][k] /= U[k][k]
    
    y = [0 for _ in range(n)]
    for i in range(n):
        y[i] = b[i]
        for j in range(i):
            y[i] -= L[i][j] * y[j]

    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = y[i]
        for j in range(i+1, n):
            x[i] -= U[i][j] * x[j]
        x[i] /= U[i][i]

    return x

2. Courant分解

Courant分解和Dolittle分解本质上并无区别，依然是将一般的参数矩阵 $A$ 拆分成两个三个矩阵，不过对于Courant分解而言， $L$ 是一个一般的下三角矩阵，而 $U$ 是对角元素全为1的上三角矩阵。

即：

\begin{pmatrix} a_{11} & a_{12} & ... & a_{1n} \\ a_{21} & a_{22} & ... & a_{2n} \\ & & ... & \\ a_{n1} & a_{n2} & ... & a_{nn} \end{pmatrix} = \begin{pmatrix} l_{11} & & & \\ l_{21} & l_{22} & & \\ & ... & & \\ l_{n1} & l_{n2} & ... & l_{nn} \end{pmatrix} \begin{pmatrix} 1 & u_{12} & ... & u_{1n} \\ & 1 & ... & u_{2n} \\ & & ... & \\ & & & 1 \end{pmatrix}

此时，仿上我们很快就能得到 $L$ 和 $U$ 的表达式如下：

\left\{ \begin{aligned} l_{ij} &= a_{ik} - \sum_{r=1}^{k-1} l_{ir}u_{rk} \\ u_{ij} &= (a_{kj} - \sum_{r=1}^{k-1}l_{kr}u_{rj}) / l_{kk} \end{aligned} \right.

同样的，我们可以通过下述函数对最终的解 $x$ 进行求解：

\left\{ \begin{aligned} Ly &= b \\ Ux &= y \end{aligned} \right.

最终结果表达式如下：

\left\{ \begin{aligned} y_i &= (b_i - \sum_{j=1}^{i-1} l_{ij}y_{j}) / l_{ii} \\ x_i &= y_{i} - \sum_{j=i+1}^{n} u_{ij}x_j \end{aligned} \right.

同样给出python伪代码实现如下：

def courant_solve(A, b):
    n = len(A)
    L = [[0 for _ in range(n)] for _ in range(n)]
    U = [[0 for _ in range(n)] for _ in range(n)]

    for k in range(n):
        for i in range(k, n):
            L[i][k] = A[i][k]
            for r in range(k):
                L[i][k] -= L[i][r] * U[r][k]
        U[k][k] = 1
        for j in range(k+1, n):
            U[k][j] = A[k][j]
            for r in range(k):
                U[k][j] -= L[k][r] * U[r][j]
            U[k][j] /= L[k][k]
        
    
    y = [0 for _ in range(n)]
    for i in range(n):
        y[i] = b[i]
        for j in range(i):
            y[i] -= L[i][j] * y[j]
        y[i] /= L[i][i]

    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = y[i]
        for j in range(i+1, n):
            x[i] -= U[i][j] * x[j]

    return x

3. 追赶法

追赶法本质上只是Courant分解的一种特殊情况而已。

它的解法和Courant分解完全一致，唯一不同的是，这里的 $A$ 矩阵是比较特殊的三元矩阵，即：

A = \begin{pmatrix} a_1 & b_1 \\ c_2 & a_2 & b_2 \\ & & ... & \\ & & c_{n-1} & a_{n-1} & b_{n-1} \\ & & & c_{n} & a_{n} \end{pmatrix}

由此，展开矩阵 $L$ 和 $U$ 也可以同步修改为：

L = \begin{pmatrix} u_1 \\ w_2 & u_2 \\ & ... & ... \\ & & w_n & u_n \end{pmatrix} , U = \begin{pmatrix} 1 & v_1\\ & 1 & v_2\\ & & ... \\ & & & 1 \end{pmatrix}

此时，由于矩阵 $L$ 和 $U$ 的特殊性，我们可以直接写出最终的解如下：

\left\{ \begin{aligned} u_i &= a_i - c_iv_{i-1} \\ v_i &= b_i / u_i \\ y_i &= (f_i - c_i y_{i-1}) / u_i \\ x_i &= (y_i - v_i x_{i+1}) \end{aligned} \right.

4. 对称正定矩阵的 $LDL^T$ 分解

同样的， $LDL^T$ 分解和上述Dolittle分解和Courant分解其实本质上也差不多，但是不同的是，这里针对的问题是对称正定矩阵。

因此，对于对称正定矩阵，我们总可以有：

\begin{aligned} A &= LDL^T \\ &= \begin{pmatrix} 1 \\ l_{21} & 1 \\ & ... \\ l_{n1} & l_{n2} & ... & 1 \end{pmatrix} \begin{pmatrix} d_1 \\ & d_2 \\ & & ... \\ & & & d_n \end{pmatrix} \begin{pmatrix} 1 & l_{12} & ... & l_{1n}\\ & 1 & ... & l_{2n}\\ & & ... \\ & & & 1 \end{pmatrix} \end{aligned}

此时，我们可以解得：

\left\{ \begin{aligned} d_k &= a_{kk} - \sum_{r=1}^{k-1}d_rl_{kr}^2\\ l_{ik} &= (a_{ik} - \sum_{r=1}^{k-1}d_rl_{ir}l_{kr}) / d_k \end{aligned} \right.

进而，我们可以将方程 $Ax = LDL^Tx = b$ 拆成三步进行完成：

\left\{ \begin{aligned} Lz &= b \\ Dy &= z \\ L^T x &= y \end{aligned} \right.

解得：

\left\{ \begin{aligned} z_i &= b_i - \sum_{j=1}^{i-1}l_{ij}z_{j} \\ y_i &= z_i / d_i \\ x_i &= y_i - \sum_{j=i+1}^{n}l_{ji}x_j \end{aligned} \right.

同样，我们给出python伪代码如下：

def ldl_solve(A, b):
    
    n = len(A)
    l = [[0 for _ in range(n)] for _ in range(n)]
    d = [0 for _ in range(n)]
    for k in range(n):
        d[k] = A[k][k]
        for r in range(k):
            d[k] -= d[r] * l[k][r] * l[k][r]
        
        l[k][k] = 1
        for i in range(k+1, n):
            l[i][k] = A[i][k]
            for r in range(k):
                l[i][k] -= d[r] * l[i][r] * l[k][r]
            l[i][k] /= d[k]

    z = [0 for _ in range(n)]
    for i in range(n):
        z[i] = b[i]
        for j in range(i):
            z[i] -= l[i][j] * z[j]
    y = [z[i]/d[i] for i in range(n)]
    x = [0 for _ in range(n)]
    for i in range(n-1, -1, -1):
        x[i] = y[i]
        for j in range(i+1, n):
            x[i] -= l[j][i] * x[j]
    return x

数值计算方法 Chapter5. 解线性方程组的直接法