数学 | 线性方程组

线性方程组是一类重要的数学问题，也是很多领域中常见的问题之一，包括计算机科学。在程序设计中，解决线性方程组是很常见的需求，例如计算机图形学中的矩阵变换、图像处理中的滤波效果、工程计算中的结构分析等等。

简介

线性方程组是一组形如 $Ax=b$ 的方程集合，其中 $A$ 是一个 $n\times n$ 的矩阵，$x$ 是一个 $n\times 1$ 的向量，$b$ 是一个 $n\times 1$ 的向量。在数学上，解线性方程组即为找到一个 $n\times 1$ 的向量 $x$，使得 $Ax=b$ 成立。

在计算机科学中，解线性方程组有多种方法，例如高斯消元法、LU分解法、Jacobi迭代法、Gauss-Seidel迭代法等等。不同的方法各有特点，可以根据实际需求选择使用哪种方法。

高斯消元法

高斯消元法是求解线性方程组的最基本算法之一，其基本思想是通过初等变换将矩阵 $A$ 转化为简化阶梯形矩阵，然后求解得到 $x$ 的值。高斯消元法具体步骤如下：

1. 构造增广矩阵 $[A,b]$，即将矩阵 $A$ 和向量 $b$ 拼接成一个新的矩阵；
2. 从第一行开始逐行处理，对于每一行 $i$，找到该行及以下所有行中第 $i$ 个元素（即主元）不为0的行 $k$，使得 $|a_{ik}|$ 最大，交换第 $i$ 行和第 $k$ 行，此时第 $i$ 行及以下所有行的主元都不为0；
3. 将第 $i$ 行的主元除以 $a_{ii}$，然后用第 $i$ 行逐一减去下方的各行，使得下方各行的第 $i$ 列元素变为0，即做行变换，使得矩阵 $A$ 转化为阶梯形矩阵；
4. 重复步骤2和步骤3，直到所有行都处理完毕，此时矩阵 $A$ 转化为简化阶梯形矩阵；
5. 从最后一行开始，逆向求解出 $x$ 的值，即直接代入式子 $Ax=b$ 中，依次求解得到 $x_n,...,x_1$。

高斯消元法的python实现如下：

import numpy as np

def gauss_elimination(A, b):
    n = len(A)
    # 构造增广矩阵
    augment_A = np.hstack([A, b.reshape(n, 1)])
    # 高斯消元
    for i in range(n):
        # 找到主元不为0的行，并将其移动到上方
        max_row = i
        for j in range(i+1, n):
            if abs(augment_A[j, i]) > abs(augment_A[max_row, i]):
                max_row = j
        augment_A[[i, max_row], :] = augment_A[[max_row, i], :]
        # 消元
        for j in range(i+1, n):
            ratio = augment_A[j, i] / augment_A[i, i]
            augment_A[j, :] -= ratio * augment_A[i, :]
    # 回代求解
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (augment_A[i, -1] - augment_A[i, i+1:-1] @ x[i+1:]) / augment_A[i, i]
    return x

LU分解法

LU分解法是另一种常见的求解线性方程组的方法，其基本思想是将矩阵 $A$ 分解为一个下三角矩阵 $L$ 和一个上三角矩阵 $U$ 的乘积，即 $A=LU$，然后用前向代换和后向代换求解 $L$ 和 $U$ 中的未知数，最终得到 $x$ 的值。LU分解法的python实现如下：

import numpy as np

def lu_decomposition(A, b):
    n = len(A)
    # 构造增广矩阵
    augment_A = np.hstack([A, b.reshape(n, 1)])
    # LU分解
    L = np.eye(n)
    U = np.zeros((n, n))
    for i in range(n):
        # 处理第i列
        for j in range(i, n):
            U[i, j] = augment_A[i, j] - L[i, :i] @ U[:i, j]
        # 处理第i行
        for j in range(i+1, n):
            L[j, i] = (augment_A[j, i] - L[j, :i] @ U[:i, i]) / U[i, i]
    # 前向代换
    y = np.zeros(n)
    for i in range(n):
        y[i] = augment_A[i, -1] - L[i, :i] @ y[:i]
    # 后向代换
    x = np.zeros(n)
    for i in range(n-1, -1, -1):
        x[i] = (y[i] - U[i, i+1:] @ x[i+1:]) / U[i, i]
    return x

Jacobi迭代法

Jacobi迭代法是一种迭代方法，其基本思想是从一个初始值开始，通过迭代的方式逼近真实的解。具体来说，Jacobi迭代法的思想是每次只更新 $x_i$ 的值，保持其余 $x_j$ 的解不变。其迭代公式如下：

$$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum\limits_{j\neq i}a_{ij}x_j^{(k)})$$

其中，$k$ 表示迭代次数，$x_i^{(k)}$ 表示在第 $k$ 次迭代后 $x_i$ 的解。Jacobi迭代法的python实现如下：

import numpy as np

def jacobi_iteration(A, b, max_iter=100, tol=1e-6):
    n = len(A)
    # 初始化
    x = np.zeros(n)
    for i in range(max_iter):
        # 迭代更新
        x_new = np.zeros(n)
        for j in range(n):
            x_new[j] = (b[j] - A[j, :] @ x + A[j, j] * x[j]) / A[j, j]
        # 判断收敛条件
        if np.max(np.abs(x_new - x)) < tol:
            break
        x = x_new
    return x

Gauss-Seidel迭代法

Gauss-Seidel迭代法是Jacobi迭代法的改进版，其基本思想是每更新一个 $x_i$ 的值，就直接将其用于后续的更新计算，而不是等待下一次迭代。其迭代公式如下：

$$x_i^{(k+1)}=\frac{1}{a_{ii}}(b_i-\sum\limits_{j=1}^{i-1}a_{ij}x_j^{(k+1)}-\sum\limits_{j=i+1}^na_{ij}x_j^{(k)})$$

其中，$k$ 表示迭代次数，$x_i^{(k)}$ 表示在第 $k$ 次迭代后 $x_i$ 的解。Gauss-Seidel迭代法的python实现如下：

import numpy as np

def gauss_seidel_iteration(A, b, max_iter=100, tol=1e-6):
    n = len(A)
    # 初始化
    x = np.zeros(n)
    for i in range(max_iter):
        # 迭代更新
        for j in range(n):
            x[j] = (b[j] - A[j, :j] @ x[:j] - A[j, j+1:] @ x[j+1:]) / A[j, j]
        # 判断收敛条件
        if np.max(np.abs(A @ x - b)) < tol:
            break
    return x

总结

线性方程组求解是一个重要的数学问题，也是计算机程序设计中常见的问题之一。本文介绍了高斯消元法、LU分解法、Jacobi迭代法、Gauss-Seidel迭代法等多种线性方程组求解方法，并提供了对应的Python实现。在实际应用中，需要根据具体问题选择合适的求解方法，以达到最优的求解效果。