您现在的位置：程式碼之家 >> 編程語言 > >> 更多程式語言 >> 程式設計綜合問答

C/C++ 大型線性方程組求解：高斯消去法與 LU 分解實作指南 (矩陣維度 > 10x10)

編輯：程式設計綜合問答

問題定義

求解任意大小的線性方程組 Ax = b，其中 A 為 n × n 方陣，n > 10。需在 C/C++ 中實作數值穩定、效能良好的演算法。

推薦演算法：部分樞軸高斯消去法

針對中小型稠密矩陣 (n < 500)，部分樞軸高斯消去法 是最標準、易實作且數值穩定的選擇。核心步驟：

前向消去：將增廣矩陣 [A|b] 轉為上三角形式，每列選擇當前欄絕對值最大的行作為樞軸列交換。
回代：從最後一行往回求解未知數。

C++ 完整實作範例 (含記憶體管理與錯誤處理)

#include 
#include 
#include 
#include 

using Matrix = std::vector>;
using Vector = std::vector;

/**
 * 部分樞軸高斯消去法求解 Ax = b
 * @param A 係數矩陣 (n x n)，函數內部會被修改
 * @param b 常數向量 (n)，函數內部會被修改
 * @return 解向量 x
 * @throws std::runtime_error 矩陣奇異或接近奇異
 */
Vector gaussianElimination(Matrix& A, Vector& b) {
    const size_t n = A.size();
    if (n == 0 || A[0].size() != n || b.size() != n)
        throw std::invalid_argument("矩陣維度不匹配");

    const double EPS = 1e-12;

    // 前向消去
    for (size_t col = 0; col < n; ++col) {
        // 1. 部分樞軸：尋找當前欄絕對值最大的行
        size_t pivotRow = col;
        double maxVal = std::abs(A[col][col]);
        for (size_t row = col + 1; row < n; ++row) {
            if (std::abs(A[row][col]) > maxVal) {
                maxVal = std::abs(A[row][col]);
                pivotRow = row;
            }
        }

        // 2. 檢查奇異性
        if (maxVal < EPS)
            throw std::runtime_error("矩陣奇異或接近奇異，無法求解");

        // 3. 交換行 (係數矩陣與常數向量同步交換)
        if (pivotRow != col) {
            std::swap(A[pivotRow], A[col]);
            std::swap(b[pivotRow], b[col]);
        }

        // 4. 消去下方行
        for (size_t row = col + 1; row < n; ++row) {
            double factor = A[row][col] / A[col][col];
            for (size_t c = col; c < n; ++c)
                A[row][c] -= factor * A[col][c];
            b[row] -= factor * b[col];
        }
    }

    // 回代
    Vector x(n);
    for (int i = n - 1; i >= 0; --i) {
        double sum = b[i];
        for (size_t j = i + 1; j < n; ++j)
            sum -= A[i][j] * x[j];
        x[i] = sum / A[i][i];
    }
    return x;
}

// 測試主程式
int main() {
    // 範例：3x3 系統，實際使用時替換為 n>10 的資料
    Matrix A = {{2, 1, -1}, {-3, -1, 2}, {-2, 1, 2}};
    Vector b = {8, -11, -3};

    try {
        Vector x = gaussianElimination(A, b);
        std::cout << "解向量 x:\n";
        for (size_t i = 0; i < x.size(); ++i)
            std::cout << "x[" << i << "] = " << x[i] << '\n';
    } catch (const std::exception& e) {
        std::cerr << "錯誤: " << e.what() << '\n';
    }
    return 0;
}