CFD線形ソルバーの選び方：JacobiからBiCGSTABまでの収束速度比較

1952年、コーネリアス・ランチョス、マグナス・ヘステネス、エドゥアルト・スティフェルの3人は、ある学会で偶然にも独立して同じ結果を発表しました。対称正定値（SPD）行列をわずか $N$ 回の反復で解くアルゴリズムです。これが現在、私たちが**共役勾配法（Conjugate Gradient, CG）**と呼んでいる手法です。興味深いのは、70年経った今日も、CFDソルバーの核心にはこのアイデアがほぼそのまま生き続けている点です。それは、Krylov（クリロフ）部分空間で解を探すというアイデアです。このポストでは、CFDで最も頻繁に使用される4つのソルバー（Jacobi、Gauss-Seidel、CG、BiCGSTAB）の収束速度を、Pythonとインタラクティブなシミュレーションで比較します。最後まで読めば、実務において「この問題にはどのソルバーが適しているか」を即座に判断できるようになるでしょう。

なぜCFDでは線形システムを反復法で解くのか#

ナビエ・ストークス方程式を離散化すると、各ステップで $A\mathbf{x}=\mathbf{b}$ という形式の巨大な線形システムが現れます。ここで、$A$ は係数行列、$\mathbf{x}$ はセル中心の変数ベクトル、$\mathbf{b}$ は右辺のソース項です。3D格子の100万セルの場合、$A$ は $10^6 \times 10^6$ サイズの疎行列になります。LU分解（行列を下三角行列と上三角行列の積に分解する直接法）をこの規模で実行すると、メモリとCPUを膨大に消費します。OpenFOAM、SU2、GerrisなどのオープンソースCFDコードが例外なく**反復法（iterative method）**を採用しているのはそのためです。

反復法の基本アイデアは単純です。初期推定値 $\mathbf{x}_0$ から出発し、残差 $\mathbf{r} = \mathbf{b} - A\mathbf{x}$ が許容値以下になるまで $\mathbf{x}$ を少しずつ修正していきます。問題は、その「少しずつ」をどう定義するかです。この定義こそがソルバーの名前を分ける境界線となります。

Krylov部分空間というアイデア#

初期残差 $\mathbf{r}_0$ に行列 $A$ を繰り返し掛けたベクトルが生成する空間があります。

$\mathcal{K}_m$ が Krylov部分空間、$m$ は次元、$\text{span}$ は線形結合で張られる空間です。Krylov系のソルバーは、解を $\mathbf{x}_m \in \mathbf{x}_0 + \mathcal{K}_m$ の中で探しつつ、残りの残差 $\mathbf{r}_m$ がある最適性条件を満たすように決定します。

解を部分空間に限定することには2つの利点があります。第一に、行列とベクトルの積以外の演算が必要ないため、疎行列構造をそのまま活用できます。第二に、理論的にその部分空間内で最適な解を抽出できます。CGはAノルム誤差を、GMRESは2ノルム残差を、Krylov空間内で正確に最小化します。

4つのソルバーの特性比較#

圧力ポアソン方程式 $-\nabla^2 p = f$ はSPDであるため、CGが定石です。一方、対流と拡散が混ざった運動量方程式は非対称になるため、BiCGSTABまたはGMRESが選択されます。OpenFOAMの PCG/PBiCGStab や、SU2の BCGSTAB は、まさにこの役割分担に従っています。

Pythonで見る収束速度#

1Dポアソン問題を例に挙げましょう。 $-u'' = 1$、境界条件 $u(0)=u(1)=0$ を有限差分で離散化すると、三対角行列が得られます。Jacobi法とCG法を直接比較してみましょう。

import numpy as np
 
def build_poisson_1d(n):
    A = np.zeros((n, n))
    for i in range(n):
        A[i, i] = 2.0
        if i > 0:     A[i, i - 1] = -1.0
        if i < n - 1: A[i, i + 1] = -1.0
    return A, np.ones(n)
 
def jacobi_sweep(A, b, tol=1e-10, max_iter=2000):
    D  = np.diag(A)
    LU = A - np.diag(D)
    x  = np.zeros_like(b)
    b0 = np.linalg.norm(b)
    history = []
    for _ in range(max_iter):
        x = (b - LU @ x) / D
        res = np.linalg.norm(b - A @ x) / b0
        history.append(res)
        if res < tol:
            break
    return x, history
 
def cg_sweep(A, b, tol=1e-10, max_iter=2000):
    x  = np.zeros_like(b)
    r  = b - A @ x
    p  = r.copy()
    rr = r @ r
    b0 = np.linalg.norm(b)
    history = [np.sqrt(rr) / b0]
    for _ in range(max_iter):
        Ap    = A @ p
        alpha = rr / (p @ Ap)
        x    += alpha * p
        r    -= alpha * Ap
        rr_new = r @ r
        history.append(np.sqrt(rr_new) / b0)
        if history[-1] < tol:
            break
        p  = r + (rr_new / rr) * p
        rr = rr_new
    return x, history
 
A, b = build_poisson_1d(n=60)
_, hist_j  = jacobi_sweep(A, b)
_, hist_cg = cg_sweep(A, b)
print(f"Jacobi: {len(hist_j):4d}回, 最終残差 {hist_j[-1]:.2e}")
print(f"CG:     {len(hist_cg):4d}回, 最終残差 {hist_cg[-1]:.2e}")

結果は一目瞭然です。

Jacobi: 2000回, 最終残差 2.11e-02
CG:       60回, 最終残差 1.87e-13

Jacobi法が2000回反復しても $10^{-2}$ に届かない一方で、CG法は $N=60$ 回でマシン精度（machine epsilon）に到達します。理論上、CG法は最大 $N$ 回で厳密解に到達します（浮動小数点誤差を無視した場合）。SPD行列に対して、CG法は魔法のような力を発揮します。

収束の競争を実際に見てみよう#

以下のシミュレーションで、直接操作してみましょう。スライダーで格子サイズ $N$ と最大反復回数を変更すると、4つのソルバーの残差曲線が対数スケールで再描画されます。

$N$ を80まで大きくすると、Jacobi法とGauss-Seidel法の曲線はほぼ水平になります。一方、CG法は階段状に下がり、$N$ 付近で急激に収束します。BiCGSTABは変動しながら下降しますが、この「ジグザグ」はbreak-down（分子または分母が0に爆発し、数値的に破綻すること）のリスクを示す視覚的な信号です。実戦でBiCGSTABが突然発散した場合は、前処理（preconditioner）を強化するか、GMRESに切り替えるのが定石です。

実務ソルバー選択ガイド#

圧力ポアソン方程式、拡散のみの熱伝導: CG + ILU(0)前処理。対称性が保証されているなら、他の選択肢は計算リソースの無駄です。

非対称な運動量方程式、対流-拡散混合: BiCGSTABが基本です。収束が不安定だったり停滞したりする場合は、GMRES($m$)に切り替えます。 $m$ はリスタート周期であり、メモリと速度のトレードオフを決定します。業界の慣行では $m=30$ 前後です。

並列CFD: Gauss-Seidel法は逐次的な性質のため、GPUでは力を発揮できません。ここ10年の標準は、Jacobiベースの平滑化（smoother）と代数マルチグリッド（AMG）の組み合わせです。

SU2のLU-SGS: SU2には、Krylov系の他にLU-SGS（Lower-Upper Symmetric Gauss-Seidel）という定常反復法も搭載されています。圧縮性流体において暗解法（implicit）の時間前進と組み合わせる際、メモリ使用量を大幅に抑えられます。ただし、Krylov系ほど頑健ではありません。

Tip: 実務では「ソルバーの選択」よりも「前処理の選択」の方が、収束速度をはるかに大きく左右します。素のCGを使うよりも、PCG（Preconditioned CG） + AMGの組み合わせの方が100倍速いことも珍しくありません。

3行まとめ#

• Jacobi法とGauss-Seidel法は単純ですが、実際のCFD問題での収束率は $1-O(1/N^2)$ と低調です。問題が大きくなるほど反復回数が2乗で増加します。
• Krylov系（CG、BiCGSTAB、GMRES）は行列ベクトル積のみで、 $N$ 回以内にマシン精度に到達します。SPDにはCG、非対称にはBiCGSTABまたはGMRESが実務上のデフォルトです。
• ソルバーよりも前処理が重要です。ILU(0)、AMG、Jacobiスムーザー、これら3つの組み合わせが10倍から100倍の速度差を生みます。

下のκスライダで各ソルバの限界がひと目で見える。