您现在的位置：程式碼之家 >> 編程語言 > >> C語言 >> C >> 關於C

手把手實現高斯消去法：從數學原理到 C 語言代碼

編輯：關於C

問題背景

求解線性方程組 Ax = b，高斯消去法通過初等行變換將增廣矩陣化為上三角矩陣，再回代求解。本文記錄從零編寫 4×4 方程組求解器的完整思路與代碼。

增廣矩陣示例

方程組：
x₁ + x₂     - 4x₄ =  1
-x₁ + x₂ + x₃ + 3x₄ = -2
x₁ + 3x₂ + 5x₃ - 4x₄ = -4
     x₂ + 2x₃ - x₄ = -2

增廣矩陣 [A|b] (4×5)：
[ 1   1   0  -4 |  1 ]
[-1   1   1   3 | -2 ]
[ 1   3   5  -4 | -4 ]
[ 0   1   2  -1 | -2 ]

算法兩階段

前向消元：對第 i 列（i = 0..n-2），將第 i 行作為主元行，消去下方各行第 i 列元素
回代求解：從最後一行向上，利用已知變量逐步求解

完整 C 代碼

#include <stdio.h>

#define ROW 4
#define COL 5  // 增廣矩陣列數 = 變量數 + 1

float a[ROW][COL] = {
    { 1,  1,  0, -4,  1},
    {-1,  1,  1,  3, -2},
    { 1,  3,  5, -4, -4},
    { 0,  1,  2, -1, -2}
};

void print_matrix() {
    for (int i = 0; i < ROW; i++) {
        for (int j = 0; j < COL; j++)
            printf("%-8.4f", a[i][j]);
        printf("\n");
    }
    printf("-----------------------\n");
}

int main() {
    print_matrix();  // 初始矩陣

    /* 前向消元 */
    for (int i = 0; i < ROW - 1; i++) {          // 主元行
        for (int j = i + 1; j < ROW; j++) {      // 被消去行
            float factor = a[j][i] / a[i][i];    // 消元因子
            for (int k = 0; k < COL; k++)        // 整行運算
                a[j][k] -= factor * a[i][k];
        }
    }
    print_matrix();  // 上三角矩陣

    /* 回代求解 */
    float x[ROW] = {0};
    for (int w = ROW - 1; w >= 0; w--) {
        float sum = a[w][COL - 1];               // 右端常數項
        for (int c = 0; c < COL - 1; c++)
            sum -= a[w][c] * x[c];               // 減去已知項貢獻
        x[w] = sum / a[w][w];                    // 求解當前變量
        printf("x%d = %-8.4f\n", w + 1, x[w]);
    }
    return 0;
}

運行輸出

初始矩陣：
1.0000  1.0000  0.0000 -4.0000  1.0000
-1.0000 1.0000  1.0000 3.0000  -2.0000
1.0000  3.0000  5.0000 -4.0000 -4.0000
0.0000  1.0000  2.0000 -1.0000 -2.0000
-----------------------
上三角矩陣：
1.0000  1.0000  0.0000 -4.0000  1.0000
0.0000  2.0000  1.0000 -1.0000 -1.0000
0.0000  0.0000  4.0000 -1.0000 -3.0000
0.0000  0.0000  0.0000 -0.2500 -1.2500
-----------------------
x4 = 5.0000
x3 = 0.5000
x2 = -1.7500
x1 = 3.2500

關鍵技術點

1. 硬編碼先行，再泛化

作者先寫「硬編碼」版回代：

x[w] = (a[w][4] - (a[w][0]*x[0] + a[w][1]*x[1] + a[w][2]*x[2] + a[w][3]*x[3])) / a[w][w];

驗證邏輯正確後，再改為通用循環：

float sum = a[w][COL-1];
for (int c = 0; c < COL-1; c++) sum -= a[w][c] * x[c];
x[w] = sum / a[w][w];

這是「遞增式開發」的典型實踐：先跑通具體案例，再抽象泛化。

2. 列主元選擇

當前代碼假設主元 a[i][i] != 0。工程級實現需加入部分列主元：在消第 i 列前，尋找 |a[k][i]| 最大的行 k（k >= i），與第 i 行交換，提高數值穩定性。

3. 複雜度分析

時間：O(n³)（三重循環）
空間：O(n²) 存儲矩陣
對於 n > 1000，建議使用 LU 分解、迭代法（共軛梯度）或專業庫（Eigen、MKL、OpenBLAS）

驗證解的正確性

/* 代入原方程驗證 */
printf("\n驗證：\n");
printf("Eq1: %.4f + %.4f - 4*%.4f = %.4f (期望 1)\n", x[0], x[1], x[3], x[0]+x[1]-4*x[3]);
printf("Eq2: -%.4f + %.4f + %.4f + 3*%.4f = %.4f (期望 -2)\n", x[0], x[1], x[2], x[3], -x[0]+x[1]+x[2]+3*x[3]);
printf("Eq3: %.4f + 3*%.4f + 5*%.4f - 4*%.4f = %.4f (期望 -4)\n", x[0], x[1], x[2], x[3], x[0]+3*x[1]+5*x[2]-4*x[3]);
printf("Eq4: %.4f + 2*%.4f - %.4f = %.4f (期望 -2)\n", x[1], x[2], x[3], x[1]+2*x[2]-x[3]);