陰解法とトーマス法 — 1次元拡散を無条件安定で解くひと筋のトリック

陽解法の時間刻みは格子幅の2乗にロックされる。一次元でその錠を外す代金は、トーマス法の掃引が1回である。

夜11時に回した計算が、朝にはNaNでした。格子を半分に切った直後で、時間刻みはそのままです。CFL(Courant–Friedrichs–Lewy)の1行を抜いた請求が、クラスタ1時間分でした。

格子幅が半分なら、陽解拡散の時間刻みは4分の1です。1024点でも、爆発しないためだけに数十万ステップを踏みます。抜け道は短いです。右辺を未来の時刻に移し、その三重対角をトーマス法で解きます。

rのスライダーを0.05から5まで動かすと、同じガウシアンパルスが陽解と陰解でどこから割れるか、1枚の画面に重なります。

CFL number r = D·dt / dx²: 0.45Explicit: stable (r ≤ 0.5)Implicit: unconditionally stable

0.5より下では2本の曲線がほぼ重なります。0.5を越えた瞬間にオレンジ(陽解)がのこぎり状に震え、1を越えると画面の外へ飛び、シアン(陰解)だけが線として残ります。

CFLが細格子を縛る#

拡散方程式は次の通りです。

\frac{\partial q}{\partial t} = D \frac{\partial^2 q}{\partial x^2}

ここで $q$ は濃度(または温度)、 $D$ は拡散係数です。最も素朴な陽解法(forward Euler + 中心差分)は

q_i^{n+1} = q_i^{n} + r\,(q_{i+1}^{n} - 2 q_i^{n} + q_{i-1}^{n}),\quad r = \frac{D\,\Delta t}{\Delta x^2}

安定であるためには $r \le 1/2$ が必要です。von Neumann解析を一行でまとめると、最短波長(波長 $= 2\Delta x$ )の増幅率が $|1 - 4r|$ で、これを $\le 1$ にするには $r \le 1/2$ が強制されます。

$1024 \to 2048$ で $\Delta t$ は $1/4$ です。もう一段細かくして $4096$ にすると、1024点基準の $\Delta t$ は $1/16$ 、ステップ数は16倍です。点数を8倍にすると、同じ時刻まで64倍のステップになります。請求書の本体はセル数ではなく、この2乗です。

$r = 0.6$ なら最短波長の増幅率は $|1 - 4r| = 1.4$ です。60ステップは、その率が60乗で効く場所です。出力の $1.33\times 10^{12}$ には、初期パルスがそのモードに載せていた分まで入っています。 $r = 0.4$ では率が $0.6$ なので、陽解のピークは $4.21\times 10^{-1}$ に留まります。

未来値を右辺へ — 陰解法の単純な仕掛け

直す行は1本です。右辺の $q^n$ を $q^{n+1}$ に書き換えます。

q_i^{n+1} = q_i^{n} + r\,(q_{i+1}^{n+1} - 2 q_i^{n+1} + q_{i-1}^{n+1})

未来値が右辺にあるので、セルを順に解くことはできません。すべてのセルが同時に結び付き、行列方程式になります。

A\,\mathbf{q}^{n+1} = \mathbf{q}^{n}

同じ解析で増幅率は $1/(1 + 4r)$ です。正の $r$ では常に1より小さいです。 $\Delta t$ を大きくしても発散しません。無条件安定(unconditionally stable)です。

代償は精度です。極端に大きい $\Delta t$ は安定でも正確ではありません。それでも発散しないという保証は、計算を一晩置く側には小さくありません。

同じ $r = 0.4$ 、60ステップのピークは陽解 $4.21\times 10^{-1}$ 、陰解 $0.4456$ です。安定な範囲なのに、小数第3位で既に割れています。時刻は同じで、差は打ち切り誤差です。

三重対角行列とトーマス法のひと撫で

行列 $A$ を書くと、1次元ではセル $i$ が $i-1$ と $i+1$ にしか結合しません。

A = \begin{pmatrix} 1+r & -r & & & \\ -r & 1+2r & -r & & \\ & -r & 1+2r & -r & \\ & & \ddots & \ddots & \ddots \\ & & & -r & 1+r \end{pmatrix}

これは 三重対角(tridiagonal) 行列です。3本の対角しかゼロではありません。一般のガウス消去法は $\mathcal{O}(N^3)$ ですが、三重対角なら $\mathcal{O}(N)$ で終わります。このアルゴリズムをトーマス法と呼びます。1949年にLlewellyn ThomasがIBMの社内メモにまとめたものです。

トーマス法は2回掃きます。

前進消去 — 副対角を順に消しながら、新しい対角 $b'$ と右辺 $d'$ を計算します。
後退代入 — 最終行から逆に値を代入します。

各段階はセルあたり4〜5回の演算で済み、メモリも3つの配列で十分です。行列を丸ごと持つ必要はありません。

境界の対角だけが $1+r$ で、内部は $1+2r$ です。ゼロフラックスではゴーストセルを境界値と等しく置くので、外側の隣が消えます。対角の $2r$ から $r$ が1つ落ち、その隣へ向いていた $-r$ も一緒に消えます。だからコードは b[0] = 1 + r です。内部の式を境界に写すと、総量が漏れます。

NumPyで陽解法と陰解法を比較#

import numpy as np
 
def explicit_diffusion_step(q, r):
    """前進オイラー + 中心差分。"""
    next_q = q.copy()
    next_q[1:-1] = q[1:-1] + r * (q[2:] - 2 * q[1:-1] + q[:-2])
    return next_q
 
def thomas_solver(a, b, c, d):
    """三重対角系 A x = d を O(N) で解く。
    a: 副対角(a[0]は未使用)
    b: 主対角
    c: 上対角(c[-1]は未使用)
    d: 右辺
    """
    n = len(d)
    cp = np.zeros(n)
    dp = np.zeros(n)
    cp[0] = c[0] / b[0]
    dp[0] = d[0] / b[0]
    for i in range(1, n):
        m = b[i] - a[i] * cp[i - 1]
        cp[i] = c[i] / m
        dp[i] = (d[i] - a[i] * dp[i - 1]) / m
    x = np.zeros(n)
    x[-1] = dp[-1]
    for i in range(n - 2, -1, -1):
        x[i] = dp[i] - cp[i] * x[i + 1]
    return x
 
def implicit_diffusion_step(q, r):
    """後退オイラー — トーマス法で (I - r·L) q^{n+1} = q^n を解く。"""
    n = len(q)
    a = np.full(n, -r)
    b = np.full(n, 1 + 2 * r)
    c = np.full(n, -r)
    b[0] = 1 + r        # ゼロフラックス境界
    b[-1] = 1 + r
    return thomas_solver(a, b, c, q.copy())
 
# ガウシアンパルスで比較
N = 80
x = np.linspace(0, 1, N)
q0 = np.exp(-((x - 0.5) ** 2) / 0.005)
D, dx = 1.0, 1.0 / N
r_values = [0.4, 0.6, 2.0]
for r in r_values:
    dt = r * dx**2 / D
    qe, qi = q0.copy(), q0.copy()
    for _ in range(60):
        qe = explicit_diffusion_step(qe, r)
        qi = implicit_diffusion_step(qi, r)
    print(f"r={r:.2f} | explicit max={np.max(np.abs(qe)):.2e} | implicit max={np.max(np.abs(qi)):.4f}")
# r=0.40 | explicit max=4.21e-01 | implicit max=0.4456
# r=0.60 | explicit max=1.33e+12 | implicit max=0.3873
# r=2.00 | explicit max=2.78e+38 | implicit max=0.2517

$r=0.6$ で陽解法は60ステップで $10^{12}$ まで爆発しています。同じ時間刻みで陰解法は滑らかに広がり、ピークは $0.3873$ です。 $r=2$ では陽解が $2.78\times 10^{38}$ 、陰解のピークは $0.2517$ まで下がります。

出力の3行は同じ時刻ではない#

$\Delta t = r\,\Delta x^2 / D$ で $N = 80$ 、 $D = 1$ なら、60ステップの物理時間は $r$ に比例します。 $r = 2$ の行は $r = 0.4$ の行より5倍先の時刻です。ピーク $0.4456$ と $0.2517$ を並べて、陰解のほうが死んでいると読むのは誤りです。時刻を揃えるには $r$ とステップ数の積を同じにします。 $r = 2$ の12ステップが $r = 0.4$ の60ステップと同じ時刻ですが、出力はその揃えをしていません。

陰解法が常に正解ではない理由

無条件安定だからといって陰解法が万能ではありません。コストは3方向に付きます。

1. 精度 vs 安定性 — 安定性は振幅が爆発しない保証であり、答えが正しい保証ではありません。後退オイラーは時間に関して $\mathcal{O}(\Delta t)$ です。大きい $\Delta t$ で精度も要るなら Crank–Nicolson( $\mathcal{O}(\Delta t^2)$ 、これも無条件安定)に切り替えるのが妥当です。

Crank–Nicolsonの最短波長の増幅率は $(1-2r)/(1+2r)$ です。 $r = 2$ では $-0.6$ です。大きさは1より小さく爆発はしませんが、符号が負なので位相が毎ステップ反転します。 $r$ をさらに上げると $-1$ に張り付きます。

2. 非線形性 — 拡散係数が解に依存する場合( $D = D(q)$ )、行列が毎ステップ変わり、PicardやNewton反復がトーマスの上に必要です。

3. 多次元 — 2D/3Dではもう三重対角ではありません。5点・7点ステンシルは対角から離れた位置にバンドを作ります。ADI(Alternating Direction Implicit)で1D三重対角を2回・3回に分けるか、JacobiからBiCGSTABまで続く反復法へ移ります。

3点目は、圧縮性から非圧縮性へ移るときにそのまま再登場します。圧力ポアソン方程式 $\nabla^2 P = -\rho \nabla\!\cdot(\mathbf{u}\!\cdot\!\nabla\mathbf{u})$ が同じ形の大規模スパース行列だからです。陰解拡散を解いたことがある人は、圧力も同じ道具で解きます。

大きいrを正解と書くと失敗する#

$r = 2$ の陰解ピーク $0.2517$ は、正解に近づいた記録ではありません。物理時間は5倍で、時間誤差は $\Delta t$ に比例します。曲線を残す前に、そのステップの物理時間を一緒に書かないと、比較は嘘になります。

格子が2次元になると、この三重対角コードはそのまま使えません。行列を組まず掃引だけするほうが安く見えるとき、陽解が壊れたあと行列なしで陰的に進むLU-SGSはその安定の保証をどこまで持っていけるでしょうか。掃引の向きを誤ると、その保証も一緒に崩れるのか。そこはまだ開いています。