3.2 高斯消去法：手算與 np.linalg.solve#

本節目標#

• 理解高斯消去法的直覺：把方程組「削成階梯」，再倒著解回去
• 完整手算一個 \(3\times 3\) 例子，並與 np.linalg.solve / scipy.linalg.solve 對照
• 認識增廣矩陣與列運算三律
• 理解部分樞軸（partial pivoting）：為什麼「小主元」危險
• 親手寫一個教學版的高斯消去函數

一、直覺：把方程組削成階梯#

解方程組最原始的策略是 消元：用一條方程去消掉另一條裡的某個未知數。對 \(3\times 3\) 的目標是把形狀

這個右上三角形的形狀叫 列階梯形。為什麼要削成這樣？因為最後一列只剩一個未知數，可以直接解出來；倒數第二列只剩兩個……一路 由下往上倒帶，每列只解一個未知數。

關鍵在於：消元的過程只是「重新包裝」方程組，解集合完全不變。

二、列運算三律#

以下三種操作都不改變解集合，是消去法唯一允許的動作：

1. 交換兩列（方程式的順序本來就無關緊要）
2. 用一個 非零 常數 乘 某一列（兩邊同乘，等式仍成立）
3. 把某一列的 倍數加到 另一列（兩式相減消元）

三、增廣矩陣#

把係數與右邊的 \(\mathbf{b}\) 併在一起寫，就是增廣矩陣 \([A\,|\,\mathbf{b}]\) 。所有列運算都作用在這整個矩陣上 —— 右邊那行是「陪著走的常數」，忘了同步更新是最常見的手算錯誤。

四、完整手算一個 \(3\times 3\)#

例：

增廣矩陣：

第一步：用第 1 列消去下面的 \(x\) 。 \(R_2 \leftarrow R_2 - 2R_1\) ， \(R_3 \leftarrow R_3 - R_1\) ：

第二步：用第 2 列消去第 3 列的 \(y\) 。 \(R_3 \leftarrow R_3 - \left(-\tfrac{1}{3}\right)R_2 = R_3 + \tfrac{1}{3}R_2\) ：

回代：由第 3 列 \(-\tfrac{7}{3}z=-7\) 得 \(z=3\) ；代入第 2 列 \(-3y-3=-9\) 得 \(y=2\) ；代入第 1 列 \(x+2+3=6\) 得 \(x=1\) 。解為 \((1,2,3)\) 。

五、用 NumPy 驗證#

import numpy as np
from scipy.linalg import solve   # 與 np.linalg.solve 幾乎同用法

A = np.array([[1.0,  1.0, 1.0],
              [2.0, -1.0, 1.0],
              [1.0,  2.0, -1.0]])
b = np.array([6.0, 3.0, 2.0])

x_np = np.linalg.solve(A, b)
x_sp = solve(A, b)
print(x_np)   # [1. 2. 3.]
print(x_sp)   # [1. 2. 3.]
print(np.allclose(A @ x_np, b))   # True：代回檢驗

np.linalg.solve 內部做的就是（帶部分樞軸的）高斯消去 —— LU 分解的近親，第 6 章會正式揭曉。手算是為了 看懂，實務計算交給它。

六、小主元的危險與部分樞軸#

消元時用來「砍」下面各列的那個數叫 主元（pivot）。我們的教學流程是「選左上角那格當主元」，但電腦不能這麼天真：

• 若主元是 0，根本沒辦法當除數 → 必須先 換列（列運算三律第 1 條登場）。
• 若主元 很小但不為 0（例如 \(10^{-16}\) ），除法會把誤差 放大到爆炸。

看個例子（在紙上是好例子，在電腦上會翻車）：

不換列直接消去：第一列乘 \(10^{16}\) 加到第二列，第二列變成天文數字，最後算出的答案被浮點誤差嚴重污染。

部分樞軸（partial pivoting） 的策略：每一欄消元前，在「該欄以下」挑 絕對值最大 的那一列當主元列，先交換再消去。這是 np.linalg.solve 預設行為，也是第 9 章討論數值穩定性的第一課。

eps = 1e-16
A = np.array([[eps, 1.0],
              [1.0, 1.0]])
b = np.array([1.0, 2.0])
print(np.linalg.solve(A, b))   # [1. 1.]：樞軸策略給出正確答案

七、手寫教學版高斯消去#

def gauss_solve(A, b, pivoting=True):
    A = A.astype(float).copy()
    b = b.astype(float).copy()
    n = len(b)

    for col in range(n):                       # 前向消去
        if pivoting:
            pr = col + np.argmax(np.abs(A[col:, col]))
            if np.isclose(A[pr, col], 0):
                raise ValueError("矩陣奇異：找不到非零主元")
            if pr != col:
                A[[col, pr]] = A[[pr, col]]
                b[[col, pr]] = b[[pr, col]]
        elif np.isclose(A[col, col], 0):
            raise ValueError("主元為 0，需要換列")

        for row in range(col + 1, n):
            f = A[row, col] / A[col, col]
            A[row, col:] -= f * A[col, col:]
            b[row] -= f * b[col]

    x = np.zeros(n)                            # 回代
    for row in range(n - 1, -1, -1):
        x[row] = (b[row] - A[row, row+1:] @ x[row+1:]) / A[row, row]
    return x


A = np.array([[1.0, 1.0, 1.0], [2.0, -1.0, 1.0], [1.0, 2.0, -1.0]])
b = np.array([6.0, 3.0, 2.0])
print(gauss_solve(A, b))                    # [1. 2. 3.]
print(np.allclose(gauss_solve(A, b), np.linalg.solve(A, b)))  # True

消去法的計算量大約是 \(\tfrac{2}{3}n^3\) 次浮點運算（回代再加 \(n^2\) ）： \(n=1000\) 時約 \(7\times 10^8\) 次 —— 電腦瞬間完成，但這也預告了第 9 章：矩陣越大，演算法的效率與穩定性越重要。

小結與練習#

• 高斯消去 = 前向消去（削成上三角）+ 回代（由下往上解），全程只動用列運算三律。
• 增廣矩陣 \([A\,|\,\mathbf{b}]\) 讓右邊的常數陪著列運算一起走。
• 部分樞軸 = 每欄挑絕對值最大的主元列；小主元會放大浮點誤差，這是 np.linalg.solve 預設做法。
• 實務永遠用 solve；手寫版本是為了理解機器在做什麼。

練習：

1. 用增廣矩陣手算 \(\begin{cases}2x+y-z=1\\ x-y+2z=4\\ 3x+2y+z=7\end{cases}\) ，再用 np.linalg.solve 驗證。
2. 把 gauss_solve 的 pivoting=False 拿去解第六節的 \(10^{-16}\) 例子，與 pivoting=True、np.linalg.solve 三者比較誤差 np.abs(A@x-b).max()。
3. 修改 gauss_solve，讓它印出每一步的增廣矩陣，重播第四節的手算過程。
4. 用 timeit 比較你手寫的 gauss_solve 與 np.linalg.solve 在 \(n=300\) 隨機矩陣上的速度，體會「向量化」的差距（第 9.4 節的主題）。