【数值分析】线性方程组求解(一)

本篇作为数值分析的入门教程，介绍一下线性方程组求解的常用算法。

概述

线性方程组求解算法是解决多个线性方程同时成立的数学方法。

常用的算法包括直接求解的高斯消元法、LU分解法等，适合求解稠密矩阵；而迭代法适用于大规模稀疏矩阵，通过不断逼近逐步得到解。

这些算法在计算机科学、工程计算等领域应用广泛，是解决实际问题的重要工具。

本文将对常用的直接求解的算法进行介绍；迭代求解算法将在下一篇博客更新。

高斯消元

高斯消元（Gaussian Elimination）是一种解决线性方程组的方法，它通过行操作将增广矩阵转换为行简化阶梯形矩阵，进而求解线性方程组。

给定一个线性方程组,将其转为增广矩阵的形式： $\left[ \begin{array}{cccc|c} a_{11} & a_{12} & \cdots & a_{1n} & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{m1} & a_{m2} & \cdots & a_{mn} & b_m \\ \end{array} \right]$

高斯消元的步骤

行交换：如果需要，交换两行以确保主对角线上的元素不为零。
行缩放：将某一行乘以一个非零常数，使得主对角线上的元素为1。
行加法：将一个行的倍数加到另一行上，以消除主对角线以下的元素。处理完后增广矩阵成为上三角矩阵。
回代求解：从最后一行开始向上回代，可以逐步求解出每个未知数 $ x_j $ 的值。

高斯消元的条件

显然，不是所有的方程组都可以高斯消元的。

考虑特例，一个方程的任意倍数组成的方程组就无法高斯消元求解。

归纳法易证，可以高斯消元充要条件为原方阵顺序主子式非0。

列主元高斯消元

当主元素$a_{kk}$相对其他元素很小时，做除法会产生很大的舍入误差。

因此，我们可以交换行，来避免绝对值小的数做除数。

即在消元前，按列选最大的主元，交换到当前行。

这些交换行的操作可以记作初等矩阵P

$PAx = b$

三角分解

高斯消元的每一步记作$L_i$, 显然每步都是下三角矩阵。最终$A$变为上三角矩阵，记作$U$, 则有

\[A = L_{n-1}...L_1A=U\] \[A = LU\] \[L = L_{1} ^{-1}L_{2} ^{-1}...L_{n-1} ^{-1} = \left[ \begin{array}{cccc} 1 & 0 & \cdots & 0 \\ a_{21} & 1 & \cdots \\ \vdots & \vdots & \ddots & \vdots \\ a_{m1} & a_{m2} & \cdots & 1 \\ \end{array} \right]\]

于是问题变成求解

\[Ly = b, Ux = y\]

直接三角分解

直接三角分解需要保证$ A $ 是一个非奇异矩阵，且所有的主对角线元素都不为零。

对于 $ A $ 中的每个元素 $ a_{ij} $：如果 $ i \leq j $（即在或低于主对角线上），计算 $ U $ 的元素：

\[u_{ij} = a_{ij} - \sum_{k=1}^{i-1} l_{ik} u_{kj}\]

如果 $ i > j $（即高于主对角线），计算 $ L $ 的元素：

\[l_{ij} = \frac{1}{u_{jj}} (a_{ij} - \sum_{k=1}^{j-1} l_{ik} u_{kj})\]

对应的 $Ly = b, Ux = y$

对于 $ i = 1, 2, \ldots, n $： $y_i = \frac{1}{l_{ii}} \left( b_i - \sum_{j=1}^{i-1} l_{ij} y_j \right)$

对于 $ i = n, n-1, \ldots, 1 $： $x_i = \frac{1}{u_{ii}} \left( y_i - \sum_{j=i+1}^{n} u_{ij} x_j \right)$

选主元三角分解

和高斯消元一样，可以使用列主元三角分解，$PA = LU$ ，这里不再赘述。

import numpy as np

def lu_decomposition(A):
    n = A.shape[0]
    L = np.zeros((n, n))
    U = np.zeros((n, n))
    
    for i in range(n):
        L[i, i] = 1
        
        for j in range(i, n):
            sum = 0
            for k in range(i):
                sum += L[i, k] * U[k, j]
            U[i, j] = A[i, j] - sum
            
        for j in range(i+1, n):
            sum = 0
            for k in range(i):
                sum += L[j, k] * U[k, i]
            L[j, i] = (A[j, i] - sum) / U[i, i]
    
    return L, U

def solve(L, U, b):
    n = len(b)
    y = np.zeros(n)
    x = np.zeros(n)
    for i in range(n):
        s = sum(L[i][j] * y[j] for j in range(i))
        y[i] = b[i] - s

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

    return x


# 示例
A = np.array([[4, 3], [6, 3]])
L, U = lu_decomposition(A)
print(L)
print(U)

b = [1,1]
x = solve(L, U, b)
print(x)

Cholesky 分解

Cholesky 分解是一种特殊的三角分解，它适用于对称正定矩阵。

对于一个对称正定矩阵 $ A $，Cholesky 分解将其分解为一个下三角矩阵 $ L $ 和其转置 $ L^T $ 的乘积，即：

\[A = LU = LDU_0\]

这里$U_0$对角线全为1，D为对角阵。

$A = A^T = U_0^TDL^T$ $L = U_0^T$ $A = LDL^T$

由于矩阵正定，可以将D开方分到两个L中去，记作 $A = LL^T$

其中 $ L $ 是一个下三角矩阵，其对角线上的元素都是正数。

分解过程

对于矩阵 $ A $ 中的每个元素 $ a_{ij} $：

对角线元素 $ l_{ii} $ 的计算公式为： $l_{ii} = \sqrt{a_{ii} - \sum_{k=1}^{i-1} l_{ik}^2}$

非对角线元素 $ l_{ij} $（其中 $ i > j $）的计算公式为： $l_{ij} = \frac{1}{l_{jj}} \left( a_{ij} - \sum_{k=1}^{j-1} l_{ik} l_{jk} \right)$

给定 $ A = LL^T $，求解 $ Ly = b $ 和 $ L^T x = y $：

对于 $ Ly = b $（前向代入）：对于 $ i = 1, 2, \ldots, n $： $y_i = \frac{1}{l_{ii}} \left( b_i - \sum_{j=1}^{i-1} l_{ij} y_j \right)$

对于 $ L^T x = y $（后向代入）：对于 $ i = n, n-1, \ldots, 1 $： $x_i = \frac{1}{l_{ii}} \left( y_i - \sum_{j=i+1}^{n} l_{ji} x_j \right)$

显然的，由于只需计算一个矩阵L,因此计算量大约是LU分解的一半。

由于分解过程$l_{jk}$数量级不增长，且对角元素恒正，因此计算结果的数值稳定性较好。

import numpy as np

def cholesky_decomposition(A):
    n = len(A)
    L = [[0 for _ in range(n)] for _ in range(n)]
    
    for i in range(n):
        for j in range(i, n):
            sum = 0
            for k in range(j):
                sum += L[i][k] * L[j][k]
            
            if i == j:
                L[i][j] = (A[i][j] - sum) ** 0.5
            else:
                L[j][i] = (A[j][i] - sum) / L[i][i]
                L[i][j] = L[j][i]
    
    return L

# 将下三角矩阵转换为NumPy数组
def lower_triangular_to_numpy(L):
    return np.array([[L[i][j] if j <= i else 0 for j in range(len(L))] for i in range(len(L))])

# 示例
A = np.array([[4, 12], [12, 37]])
L = cholesky_decomposition(A)
L_np = lower_triangular_to_numpy(L)

print(L_np)

追赶法

追赶法，也称为三对角矩阵算法（Tridiagonal Matrix Algorithm, TDMA），是一种用于解决三对角线性方程组的高效算法。

三对角矩阵是指矩阵中只有主对角线、主对角线上方和下方的元素不为零。

追赶法利用了三对角矩阵的特殊结构，通过消元过程求解方程组，其时间复杂度为 $O(n)$，远优于一般的高斯消元法。

追赶法的基本步骤

给定一个三对角线性方程组：

\[\begin{cases} b_1 x_1 + c_1 x_2 = d_1 \\ a_i x_{i-1} + b_i x_i + c_i x_{i+1} = d_i, & i = 2, 3, \ldots, n-1 \\ a_n x_{n-1} + b_n x_n = d_n \end{cases}\]

其中，$a_i, b_i, c_i, d_i$ 是已知系数，$x_i$ 是未知数。

将A进行LU分解，使上三角矩阵主对角线元素为1，上对角线元素为$\beta_i$, 计算$ \beta_i $： $\beta_1 = \frac{c_1}{b_1}$ $\beta_i = \frac{c_i}{b_i - \beta_{i-1} a_i}, \quad i = 2, 3, \ldots, n-1$

解向量的求解和LU分解一致，解 $Ly=f, Ux = y $，只是相邻元素可以直接求解。

前向带入求解： $y_1 = \frac{d_1}{b_1}$ $y_i = \frac{d_i - a_i y_{i-1}}{b_i - \beta_{i-1} a_i}$

后向带入求解： $x_n = y_n$ $x_i = y_i - \beta_i x_{i+1}$

如果三对角矩阵严格对角占优，能有更好的数值稳定性。

def tdma(a, b, c, d):
    n = len(b)
    beta = [0] * n
    y = [0] * n
    x = [0] * n
    
    # 计算beta_i
    beta[0] = c[0] / b[0]
    for i in range(1, n-1):
        beta[i] = c[i] / (b[i] - beta[i-1] * a[i-1])
    
    # 前向代入求解y
    y[0] = d[0] / b[0]
    for i in range(1, n):
        y[i] = (d[i] - a[i-1] * y[i-1]) / (b[i] - beta[i-1] * a[i-1])
    
    # 后向代入求解x
    x[n-1] = y[n-1]
    for i in range(n-2, -1, -1):
        x[i] = y[i] - beta[i] * x[i+1]
    
    return x

# 示例
a = [1,1]  # 下对角线元素
b = [1,2,3]  # 主对角线元素
c = [1,1]  # 上对角线元素
d = [1, 2, 3]  # 常数项

# 解方程组
x = tdma(a, b, c, d)
print(x)