📜  数学 |线性方程组(1)

📅  最后修改于: 2023-12-03 15:39:57.684000             🧑  作者: Mango

数学 | 线性方程组

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

简介

线性方程组是一组形如 $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实现。在实际应用中,需要根据具体问题选择合适的求解方法,以达到最优的求解效果。