高斯列主元消去法—

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

高斯列主元消去法

1. 高斯消去法

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

线性方程组的一般形式为：
$\begin{cases} 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 \\ \vdots \\ a_{n1}x_1 + a_{n2}x_2 + \cdots + a_{nn}x_n = b_n \\ \end{cases}$
其增广矩阵形式为：
$\begin{bmatrix} a_{11} & a_{12} & \cdots & a_{1n} & b_1 \\ a_{21} & a_{22} & \cdots & a_{2n} & b_2 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} & b_n \\ \end{bmatrix}$

2. 列主元消去法

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

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

在这里插入图片描述

3. 算法步骤

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

(1) 前向消元

步骤 1：从左到右逐列进行消元。
- 在当前列中找到绝对值最大的元素（列主元）。
- 将包含列主元的行与当前行交换。
- 如果当前列主元为零，说明矩阵奇异（无唯一解），算法终止。
- 用当前行的主元将当前列下方的所有元素变为零，即将矩阵变为上三角矩阵。
数学表示：
- 对于第 $j$ 列，找到主元 $a_{k,j}$ ，满足
  $\arg\max_{i \geq j} |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 npdef 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 = jfor 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)