9.1 數值模擬中的迭代求解器#

概念說明#

工程模擬（結構應力、流場、電磁場、熱傳）的標準流程：偏微分方程 → 網格離散化 → 大型稀疏線性（或非線性）系統 → 求解。百萬網格的系統若用直接法（LU），記憶體以 \(O(n^2)\) 、時間以 \(O(n^3)\) 增長——千萬維問題直接出局。迭代求解器因此是模擬軟體的心臟，工業級組合如下：

Krylov 子空間法：4.3 節 CG 的家族。對稱正定用 CG；非對稱用 GMRES（最小殘差法，每步擴充 Krylov 基、用最小化 \(\|A\mathbf{x} - \mathbf{b}\|\) 選解；需重啟以控制記憶體）或 BiCGSTAB（短遞推、更省記憶體但收斂較不平滑）。
預處理（preconditioning）是成敗關鍵：把 \(\kappa\) 壓小。常用階梯——Jacobi/對角縮放（便宜、弱）、不完全 LU（ILU，中庸）、代數多重網格（AMG，最強，與問題尺寸無關的收斂輪數）。
多重網格（multigrid）：細網格上誤差分頻率——迭代法（如 Gauss-Seidel，又稱 smoother）快壓高頻、慢壓低頻；多重網格把低頻誤差投影到粗網格求解再插值回來，各頻段各擊破。收斂輪數與網格大小無關（ \(O(1)\) 輪），理論上是偏微分方程求解的最優方法。
非線性問題：牛頓-克雷洛夫（Newton-Krylov）：外層牛頓迭代（平方收斂、少輪數），內層每個牛頓步用 Krylov 法解（不必組裝 Jacobian 矩陣——只需 Jacobian-向量積，可用有限差分或自動微分），外平方內迭代，是現代非線性模擬的標準架構。
瞬態問題與熱啟動：時間推進的每一步用上一步的解當初值，迭代次數驟降——「好的初值是最強的預處理」。

共通工程準則：殘差導向停機（ \(\|r\| \le \varepsilon \|b\|\) ，而不是只看步長）；監控收斂歷史的形狀（平坦→預處理不良；鋸齒→非對稱問題用了 CG；暴增→發散，需檢查時間步長或鬆弛因子）。

生活類比#

多重網格像公司治理：基層（細網格）解決日常瑣事（高頻誤差），解不開的大方針（低頻誤差）上呈主管（粗網格），決定後層層下達（插值回細網格）。各層只做自己擅長的尺度，整體效率遠高於基層硬扛一切。

演算法步驟#

GMRES 重啟版（概念）：

給初值 \(\mathbf{x}_0\) 、殘差 \(\mathbf{r}_0 = \mathbf{b} - A\mathbf{x}_0\) 、重啟長度 \(m\) 。
建 Krylov 基 \(\{\mathbf{r}_0, A\mathbf{r}_0, \dots, A^{m-1}\mathbf{r}_0\}\) （Arnoldi 正交化）。
在此子空間解最小化 \(\|A(\mathbf{x}_0 + \mathbf{v}) - \mathbf{b}\|\) 的小問題（Hessenberg 最小平方）。
更新 \(\mathbf{x}\) ；殘差夠小則停，否則以新殘差重啟回到 2。

Python 範例#

用 SciPy 對 2D 拉普拉斯系統比較求解器：

import numpy as np
from scipy.sparse import diags
from scipy.sparse.linalg import spsolve, cg, gmres, spilu, LinearOperator

def laplace_system(n):
    """2D 五點差分，n x n 內點 -> (n^2) x (n^2) SPD 稀疏矩陣"""
    N = n * n
    main = -4.0 * np.ones(N)
    off = np.ones(N - 1); off[n-1::n] = 0      # 斷開跨行連接
    A = diags([main, off, off, np.ones(N - n), np.ones(N - n)],
              [0, 1, -1, n, -n], format="csr")
    b = np.ones(N)                              # 常數源項
    return A.tocsr(), b

n = 80
A, b = laplace_system(n)
print(f"未知數 {A.shape[0]}，非零元素 {A.nnz}")

import time
t = time.perf_counter(); x1 = spsolve(A, b); t1 = time.perf_counter() - t
print(f"直接法  : {t1:.3f}s, 殘差 {np.linalg.norm(A @ x1 - b):.1e}")

t = time.perf_counter()
x2, info = cg(A, b, rtol=1e-10)
t2 = time.perf_counter() - t
print(f"CG      : {t2:.3f}s, 殘差 {np.linalg.norm(A @ x2 - b):.1e}, info={info}")

# ILU 預處理的 CG
iluw = spilu(A.tocsc(), drop_tol=1e-4)
M = LinearOperator(A.shape, matvec=iluw.solve)
t = time.perf_counter()
x3, info = cg(A, b, rtol=1e-10, M=M)
t3 = time.perf_counter() - t
print(f"CG+ILU  : {t3:.3f}s, 殘差 {np.linalg.norm(A @ x3 - b):.1e}, info={info}")

典型結果：直接法最慢且吃記憶體；裸 CG 收斂慢；ILU 預處理後迭代次數大減、總時間最短。

複雜度與收斂性分析#

Krylov 法每輪 \(O(\text{nnz})\) ，收斂輪數由 \(\sqrt{\kappa}\) （CG）或譜分布（GMRES）控制；多重網格 \(O(N)\) 總成本（每輪掃一遍未知數、輪數 \(O(1)\) ），理論最優。記憶體：直接法 \(O(N^{1.5}) \sim O(N^2)\) ，Krylov \(O(N)\) ，GMRESTruncation 後 \(O(mN)\) 。實務調參順序：先選對 Krylov 法（對稱？）、再上預處理（ILU/AMG）、最後調停機準則。非線性模擬的收斂失敗大多不是線性求解器的錯，而是牛頓步太大或時間步長太粗。

練習#

對 \(n = 40, 80, 160\) 的拉普拉斯系統，記錄 CG 與直接法的時間，驗證 CG 總成本近線性增長。
改變 ILU 的 drop_tol（1e-2, 1e-4, 1e-6），畫出「預處理建立時間＋CG 迭代時間」對 drop_tol 的曲線，找出甜蜜點。
實作 9.1 之 V-cycle 多重網格的簡化版（一維問題），與 200 輪 Gauss-Seidel 比較達殘差 \(10^{-8}\) 的時間。