Gauss列主元消去法、平方根法與追趕法：C++ 實現線性方程組求解器

專案背景

本文實現三類線性方程組的數值求解器：一般線性方程組（Gauss 列主元消去法）、對稱正定方程組（平方根法/Cholesky 分解）、三對角線性方程組（追趕法/Thomas 算法）。所有代碼基於 C++，並對比驗證數值結果。

核心數值方法

1. Gauss 列主元消去法

針對一般 n×n 線性系統 A x = b。消元階段每列選擇絕對值最大元素作為主元，避免「大數吃小數」導致的數值不穩定。流程分兩階段：

• 前向消元：將增廣矩陣化為上三角形式
• 回代求解：從最後一行向上逐步求出未知數

偽代碼關鍵步驟：
for k=1..n-1: 選主元 → 行交換 → 消去第k列下方元素

2. 平方根法

專門求解對稱正定矩陣 A。利用 Cholesky 分解 A = L L^T，其中 L 為下三角矩陣。分解公式：
L_kk = √(a_kk - Σ_{j Lkj²)
Lik = (aik - Σj LijLkj) / Lkk (i > k)}

求解分兩步：L y = b（前向代入）→ L^T x = y（回代）。

3. 追趕法

針對三對角矩陣（主對角線及其上下各一條對角線非零）。將係數矩陣分解為 A = L U，其中 L 為單位下雙對角、U 為上雙對角。遞推公式：
α₁ = a₁₁, β₁ = b₁ / α₁
α_i = a_ii - a_i,i-1 a_i-1,i / α_i-1
β_i = (b_i - a_i,i-1 β_i-1) / α_i-1

最後回代：x_n = β_n, x_i = β_i - a_i,i+1 x_i+1 / α_i。

C++ 實現要點

// Gauss 列主元核心片段
for (int k = 1; k <= n; ++k) {
    // 選主元
    int maxRow = k;
    for (int i = k+1; i <= n; ++i)
        if (fabs(a[i][k]) > fabs(a[maxRow][k])) maxRow = i;
    if (maxRow != k) swapRows(k, maxRow);
    // 消元
    for (int i = k+1; i <= n; ++i) {
        double factor = a[i][k] / a[k][k];
        for (int j = k; j <= n+1; ++j) a[i][j] -= factor * a[k][j];
    }
}
// 回代
for (int i = n; i >= 1; --i) {
    double sum = 0;
    for (int j = i+1; j <= n; ++j) sum += a[i][j] * x[j];
    x[i] = (a[i][n+1] - sum) / a[i][i];
}

數值結果與誤差分析

三種方法均能正確求解測試方程組。但浮點運算引入的舍入誤差不可避免：消元後下三角區殘留 ~1e-15 量級非零值。去除「將極小值截為 0」的處理後，解向量出現 1.23e-14 這類本應為 0 的值。這正是數值分析研究誤差控制的核心動機。

總結

• 一般矩陣首選 Gauss 列主元，穩健通用
• 對稱正定矩陣用 平方根法，儲存量減半、速度更快
• 三對角系統專用 追趕法，O(n) 複雜度極高效
• 工程實踐中務必關注條件數與誤差傳播