高斯列主元消去法

1. 高斯消去法

高斯消去法是一種求解線性方程組 $A\mathbf{x}=\mathbf{b}$ 的方法，通過逐步化簡增廣矩陣，將其變為上三角矩陣，從而方便求解未知數。

線性方程組的一般形式為：
$$?\begin{array}{l}{a}_{11}{x}_1+a_{12}{x}_2+?+a_{1n}{x}_n=b_1\\ a_{21}{x}_1+a_{22}{x}_2+?+a_{2n}{x}_n=b_2\\ ?\\ a_{n1}{x}_1+a_{n2}{x}_2+?+a_{nn}{x}_n=b_n\end{array}$$
其增廣矩陣形式為：
$$\left[\begin{array}{ccccc}{a}_{11}& a_{12}& ?& a_{1n}& b_1\\ a_{21}& a_{22}& ?& a_{2n}& b_2\\ ?& ?& ?& ?& ?\\ a_{n1}& a_{n2}& ?& 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}?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)

• 主程序從文件讀取增廣矩陣數據。
• 調用 column_elimination 函數求解線性方程組。
• 輸出解向量。