高斯列主元消去法——python实现

高斯列主元消去法

1. 高斯消去法

高斯消去法是一种求解线性方程组 $A\mathbf{x}=\mathbf{b}$ 的方法，通过逐步化简增广矩阵，将其变为上三角矩阵，从而方便求解未知数。

线性方程组的一般形式为：
$$\{\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+\cdots +a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+\cdots +a_{2n}{x}_n=b_2\\ ⋮\\ a_{n1}{x}_1+a_{n2}{x}_2+\cdots +a_{nn}{x}_n=b_n\end{array}$$
其增广矩阵形式为：
$$\left[\begin{array}{ccccc}{a}_{11}& a_{12}& \cdots & a_{1n}& b_1\\ a_{21}& a_{22}& \cdots & a_{2n}& b_2\\ ⋮& ⋮& \ddots & ⋮& ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}& b_n\end{array}\right]$$

2. 列主元消去法

在高斯消去法中，可能会遇到主元（当前行的对角线元素）过小，导致计算误差放大的问题。为了避免这种情况，引入列主元消去法。
列主元消去法的核心思想是：在每一步消元之前，选择当前列中绝对值最大的元素作为主元，并将含有该主元的行与当前行交换。通过行交换将其移到主对角线位置。这样做有两个主要目的：

• 避免小主元造成的数值不稳定
• 减少舍入误差的积累

3. 算法步骤

高斯列主元消去法的算法分为两个主要阶段：前向消元和回代求解。

(1) 前向消元

• 步骤 1：从左到右逐列进行消元。

  • 在当前列中找到绝对值最大的元素（列主元）。
  • 将包含列主元的行与当前行交换。
  • 如果当前列主元为零，说明矩阵奇异（无唯一解），算法终止。
  • 用当前行的主元将当前列下方的所有元素变为零，即将矩阵变为上三角矩阵。

• 数学表示：

  • 对于第 $j$ 列，找到主元 $a_{k,j}$，满足
    $$k=\arg \underset{i\ge j}{\max }|a_{i,j}|$$
  • 交换第 $j$ 行与第 $k$ 行。
  • 消元操作：
    计算乘数 $m_{ij}=a_{ij}/a_{jj}$
    对第 $k$ 行执行消元操作：$a_{i,j}=a_{i,j}-m_{ij}\cdot a_{j,j}$

(2) 回代求解

• 步骤 2：一旦矩阵变为上三角矩阵，从最后一行开始逐行求解未知数。
  • 对于第 $k$行（从下往上），计算
    $$x_k=\frac{b_k-{\sum }_{j=k+1}^n{a}_{k,j}{x}_j}{a_{k,k}}$$

4. 数值稳定性分析

• 主元选择：选择列中绝对值最大的元素作为主元，避免除以接近零的小数
• 误差控制：减少舍入误差的传播，提高计算精度
• 适用性：对于绝大多数非奇异矩阵都能稳定求解

5. 时间复杂度

• 前向消元：$O(n^3)$
• 回代过程：$O(n^2)$
• 总体复杂度：$O(n^3)$

二、代码实现与讲解

import numpy as np

def column_elimination(A):
    """
    使用高斯列主元消去法求解线性方程组
    :param A: 增广矩阵（numpy数组），最后一列为常数向量
    :return: 解向量（numpy数组），或None（如果奇异）
    """
    n = len(A)  # 获取矩阵行数（即方程个数）

    # 前向消元过程
    for j in range(n):
        # 步骤1: 寻找列主元（当前列中绝对值最大的元素）
        max_row = j  # 初始化主元行为当前行
        for k in range(j + 1, n):
            # 比较当前列中各元素的绝对值
            if abs(A[k, j]) > abs(A[max_row, j]):
                max_row = k  # 更新主元行
        
        # 步骤2: 交换当前行和主元行
        A[[j, max_row]] = A[[max_row, j]]
        
        # 步骤3: 检查主元是否为零（奇异矩阵）
        if abs(A[j, j]) < 1e-15:  # 设置一个小的阈值避免浮点误差
            return None  # 矩阵奇异，无唯一解
        
        # 步骤4: 消元操作
        for i in range(j + 1, n):
            # 计算乘数因子
            mij= A[i, j] / A[j, j]
            # 对第i行进行消元：A[i, j:] 表示从第j列开始到最后一列
            # 这里使用了向量化操作，提高计算效率
            A[i, j:] -= mij* A[j, j:]

    # 回代过程
    x = np.zeros(n)  # 初始化解向量
    # 从最后一行开始向上回代
    for k in range(n - 1, -1, -1):
        # 计算已知解的部分和：A[i, i+1:n] 与 x[i+1:n] 的点积
        # 然后从常数项中减去这个和
        # 最后除以主对角线元素
        x[k] = (A[k, -1] - np.dot(A[k, k + 1:n], x[k + 1:n])) / A[k, k]
    
    return x

# 从文件读取矩阵数据
def read_file(filename):
    """
    从文本文件读取矩阵数据
    :param filename: 文件名
    :return: NumPy数组表示的矩阵
    """
    data = []  # 初始化数据列表
    with open(filename, 'r') as file:
        for line in file:
            # 处理行中的多个空格分隔
            # 分割每行数据并转换为浮点数
            row = [float(num) for num in line.split()]
            data.append(row)
    # 将列表转换为NumPy数组
    return np.array(data, dtype=float)

# 主程序
if __name__ == "__main__":
    # 从文件读取数据
    filename = 'equation2.txt'  # 数据文件名
    data = read_file(filename)  # 读取数据
    
    print("增广矩阵：")
    print(data)
    
    # 执行高斯列主元消去法
    # 使用data.copy()避免修改原始数据
    solution = column_elimination(data.copy())
    
    if solution is None:
        print("\n矩阵奇异，无唯一解")
    else:
        print("\n方程组的解：")
        print(solution)

1. 函数定义：column_elimination(A)

def column_elimination(A):
    """
    使用高斯列主元消去法求解线性方程组
    :param A: 增广矩阵（numpy数组）
    :return: 解向量（numpy数组），或None（如果奇异）
    """

• 函数接受一个增广矩阵 ( A )，并通过高斯列主元消去法求解线性方程组。
• 如果矩阵奇异（无唯一解），函数返回 None。

2. 前向消元过程

# len(A)————矩阵的行数
n = len(A)

# 前向消元过程
for j in range(n):
    # 寻找列主元（当前列中绝对值最大的元素）
    max_row = j
    for k in range(j + 1, n):
        # A[i, j] 表示 A 的第 i 行第 j 列的元素
        if abs(A[k, j]) > abs(A[max_row, j]):
            max_row = k

• 在每一步消元中，先找到当前列 $j$ 中绝对值最大的元素作为主元，并记录其所在行 max_row。

    # 交换当前行和主元行
    A[[j, max_row]] = A[[max_row, j]]

• 将包含主元的行与当前行进行交换，确保当前列的主元在对角线位置。

    # 检查主元是否为零（奇异矩阵）
    if abs(A[j, j]) < 1e-15:
        return None

• 如果主元的绝对值过小（小于阈值 $10^{-15}$），认为矩阵奇异（可能无解或有无穷多解），直接返回 None。

    # 消元操作（对当前行以下的所有行进行消元）
    # 通过消元将矩阵变为上三角矩阵
    for i in range(j + 1, n):
        mij = A[i, j] / A[j, j]
        A[i, j:] -= mij * A[j, j:]

• 对当前列 $j$ 下方的所有行进行消元操作。
• 计算乘子 $m_{ij}=\frac{a_{i,j}}{a_{j,j}}$，并将当前行的 $j$ 列及其右侧的元素减去主元行的对应元素乘以乘子$m_{ij}$，使当前列的下方元素变为零。

3. 回代求解过程

# 回代过程
x = np.zeros(n)   #x——[0. 0. 0.]

• 初始化解向量 ( x ) 为零向量，长度为矩阵的行数 ( n )。

# 从最后一行向上逐行求解，计算每个未知数的值
for k in range(n - 1, -1, -1):
    x[k] = (A[k, -1] - np.dot(A[k, k + 1:n], x[k + 1:n])) / A[k, k]

• 从最后一行开始，逐行计算未知数 $x_k$。
• 根据公式
  $$x_k=\frac{b_k-{\sum }_{j=k+1}^n{a}_{k,j}{x}_j}{a_{k,k}}$$
  计算每个未知数的值。
• 使用 np.dot 计算当前行右侧未知数的线性组合。

return x

• 返回解向量 ( x )。

4. 从文件读取矩阵数据

def read_file(filename):
    data = []
    with open(filename, 'r') as file:
        for line in file:
            # 处理行中的多个空格分隔
            row = [float(num) for num in line.split()]
            data.append(row)
    return np.array(data, dtype=float)

• 从文件中逐行读取数据，并将其转换为浮点数存储为二维列表。
• 最后，将数据转换为 numpy 数组以方便后续计算。

5. 主程序

if __name__ == "__main__":
    # 从文件读取数据
    filename = 'equation2.txt'
    data = read_file(filename)
    print(data)

    # 执行高斯列主元消去法
    # 深拷贝data.copy()————传入增广矩阵的副本，以避免修改原矩阵
    solution = column_elimination(data.copy())

    print("\n方程组的解：")
    print(solution)

• 主程序从文件读取增广矩阵数据。
• 调用 column_elimination 函数求解线性方程组。
• 输出解向量。