格子を2倍細かくしたら誤差が3倍になった — Jacobi反復の増分収束判定が隠すもの

細かい格子のほうが間違った答えを出した

2003年のある学部CFD課題レポートに、こんな一文があります。「格子間隔が狭いぶん計算量は増えたが、それに見合う精度は得られなかった」。2次元定常熱伝導をJacobi反復(各点の値を上下左右4点の平均で置き換える反復法)で解いた結果です。格子を細かくすれば精度が上がるという常識とは逆です。この記事では、原因が格子ではなく止めるルールにあることを示します。読み終えると、「収束した」という判定を格子と一緒にどう変えるべきかがわかります。

2003年の課題レポートの一行#

問題は単純です。長方形の板の下壁が100 °C、残り三つの壁が0 °Cです。内部の温度はラプラス方程式 $\nabla^2 T = 0$ に従います。5点差分(中心点と上下左右4点で作る2次精度の近似)で離散化すると連立方程式になります。レポートはこれをJacobi、点Gauss–Seidel、線Gauss–Seidel、PSOR、LSORで解きました。

止めるルールはどの方法でも同じでした。1回の反復で変わった量が許容値 $e_s$ を下回ったら止めます。レポートは $e_s = 10^{-4}$ を使いました。そして二つのことを観察しています。格子を細かくすると誤差がかえって大きくなった。SORの $\omega$ は0.8、1.0、1.2のうち1.2が最も速かった。この二つは同じ根から出ています。

増分は誤差ではない

反復 $k$ での解を $T^k$、離散方程式の真の解を $T^*$ とします。欲しいのは誤差 $e_k = T^k - T^*$ です。しかし $T^*$ はわからないので、代わりに増分 $d_k = T^{k+1} - T^k$ を測ります。両者の関係は反復行列のスペクトル半径(最大固有値の絶対値)$\rho$ で決まります。

ここで $\rho$ は最も遅く減る誤差モードの減衰率です。数回反復するとそのモードだけが残るので、この近似はよく当たります。つまり $|d_k| < e_s$ で止めると、残る誤差はおよそ $e_s/(1-\rho)$ です。$\rho$ が1に近いほど、この倍率は大きくなります。

下のシミュレーションで実際に操作してみましょう。最も遅いモードを一つだけ取り出した図です。

$n$ を10から160に上げてみてください。青い増分の曲線はほぼ底に張り付きますが、オレンジの誤差の山はほとんど縮みません。右のグラフで停止した瞬間の黄色い棒がその差、つまり $1/(1-\rho)$ 倍です。

ρは格子が決める

Jacobiの $\rho$ は計算できます。大きさ $L_x \times L_y$ の板に間隔 $h$ の一様格子なら、次のとおりです。

$h$ は格子間隔、$L_x, L_y$ は板の横と縦の長さです。要点は $1-\rho_J \propto h^2$ であることです。$h$ を半分にすると $1-\rho_J$ は1/4になります。同じ $e_s$ で止めれば、残る反復誤差は4倍になります。

一方、離散化誤差(5点差分そのものの誤差)は $h^2$ で減ります。だから全誤差は二つの項の和になります。

最後の式は $L_x=1$、$L_y=2$ の場合です。前の項は減り、後ろの項は増えます。$h$ がある値より小さくなると、格子を細かくするほど答えが悪くなります。2003年のレポートが見たのはこれです。

コード:同じ板、三つの格子、三つの反復法

レポートの21×41格子を正方形セルとみなし、板を1×2としました。格子は21×41、41×81、81×161と一様に細かくします。止めるルールは $\max|d_k| < 10^{-3}$ です。レポートの $e_s$ は相対誤差のようですが、ここでは解析を簡単にするため絶対値を使いました。

参照解を二つ用意します。一つは5点離散式の厳密解 $T_{\text{discrete}}$ で、反復誤差はこれと比べます。もう一つは連続問題のフーリエ級数解で、全誤差はこれと比べます。下壁の両端の角では温度が不連続で誤差が大きいため、全誤差は $y \ge L_y/8$ の領域だけで測ります。SORは赤黒順序(チェス盤のように2色に分けて交互に更新する方法)でベクトル化しました。

import numpy as np
 
LX, LY, TB = 1.0, 2.0, 100.0   # 1 x 2 の板、下壁 100 °C、他の壁 0
 
def exact_plate(nx, ny, terms=399):
    # 連続ラプラス問題の解析解(フーリエ級数)
    x = np.linspace(0, LX, nx); y = np.linspace(0, LY, ny)
    X, Y = np.meshgrid(x, y, indexing="ij")
    T = np.zeros_like(X)
    for n in range(1, terms + 1, 2):
        k = n * np.pi / LX
        decay = np.exp(-k * Y) * (1 - np.exp(-2 * k * (LY - Y))) / (1 - np.exp(-2 * k * LY))
        T += 4 * TB / (n * np.pi) * np.sin(k * X) * decay
    return T
 
def discrete_plate(nx, ny):
    # 5点離散式の厳密解(x はサイン変換、y は漸化式)
    T = np.zeros((nx, ny)); i = np.arange(1, nx - 1); J = ny - 1; j = np.arange(ny)
    for p in range(1, nx - 1):
        s = np.sin(p * np.pi * i / (nx - 1))
        c = 2 / (nx - 1) * np.sum(TB * s)
        lam = 2 - 2 * np.cos(p * np.pi / (nx - 1))
        r = (2 + lam + np.sqrt((2 + lam) ** 2 - 4)) / 2
        T[1:-1, :] += np.outer(s, c * (r ** (j - J) - r ** (J - j)) / (r ** -J - r ** J))
    T[1:-1, 0] = TB
    return T
 
def jacobi_sweep(T, w, masks):
    Tn = T.copy()
    Tn[1:-1, 1:-1] = 0.25 * (T[2:, 1:-1] + T[:-2, 1:-1] + T[1:-1, 2:] + T[1:-1, :-2])
    return Tn
 
def rb_sor_sweep(T, w, masks):
    T = T.copy()
    for m in masks:                       # 赤セル、次に黒セル
        gs = np.zeros_like(T)
        gs[1:-1, 1:-1] = 0.25 * (T[2:, 1:-1] + T[:-2, 1:-1] + T[1:-1, 2:] + T[1:-1, :-2])
        T[m] = (1 - w) * T[m] + w * gs[m]
    return T
 
def omega_opt(nx, ny):
    rho_j = 0.5 * (np.cos(np.pi / (nx - 1)) + np.cos(np.pi / (ny - 1)))
    return rho_j, 2 / (1 + np.sqrt(1 - rho_j ** 2))
 
def solve_until(nx, ny, sweep, w, es):
    T = np.zeros((nx, ny)); T[1:-1, 0] = TB
    I, Jg = np.meshgrid(np.arange(nx), np.arange(ny), indexing="ij")
    inner = (I > 0) & (I < nx - 1) & (Jg > 0) & (Jg < ny - 1)
    masks = [inner & ((I + Jg) % 2 == c) for c in (0, 1)]
    for k in range(1, 100000):
        Tn = sweep(T, w, masks)
        d = np.abs(Tn - T).max(); T = Tn
        if d < es:
            return T, k
 
es = 1e-3                                 # max |T^(k+1) - T^k| < es で停止
print(f"stop rule: max|dT| < {es}   error measured on y >= LY/8 (away from hot corners)")
print(f"{'grid':>8} {'method':>10} {'omega':>6} {'iters':>6} {'es/(1-rho)':>11} {'iter_err':>9} {'disc_err':>9} {'total':>8}")
for nx, ny in [(21, 41), (41, 81), (81, 161)]:
    Td, Ta = discrete_plate(nx, ny), exact_plate(nx, ny)
    band = slice((ny - 1) // 8, None)
    disc = np.abs(Td - Ta)[:, band].max()
    rho_j, w_opt = omega_opt(nx, ny)
    for name, sweep, w, rho in [("Jacobi", jacobi_sweep, 1.0, rho_j),
                                ("SOR", rb_sor_sweep, 1.2, None),
                                ("SOR-opt", rb_sor_sweep, w_opt, w_opt - 1)]:
        T, k = solve_until(nx, ny, sweep, w, es)
        bound = f"{es / (1 - rho):11.3f}" if rho is not None else f"{'-':>11}"
        print(f"{nx:>3}x{ny:<4} {name:>10} {w:6.3f} {k:6d} {bound} "
              f"{np.abs(T - Td).max():9.4f} {disc:9.4f} {np.abs(T - Ta)[:, band].max():8.4f}")

stop rule: max|dT| < 0.001   error measured on y >= LY/8 (away from hot corners)
    grid     method  omega  iters  es/(1-rho)  iter_err  disc_err    total
 21x41       Jacobi  1.000    629       0.130    0.1264    0.1691   0.1525
 21x41          SOR  1.200    256           -    0.0422    0.1691   0.1639
 21x41      SOR-opt  1.780     47       0.005    0.0021    0.1691   0.1689
 41x81       Jacobi  1.000   1831       0.519    0.4922    0.0439   0.4850
 41x81          SOR  1.200    791           -    0.1683    0.0439   0.1616
 41x81      SOR-opt  1.883     90       0.009    0.0033    0.0439   0.0435
 81x161      Jacobi  1.000   4815       2.075    1.8097    0.0112   1.8078
 81x161         SOR  1.200   2260           -    0.6486    0.0112   0.6468
 81x161     SOR-opt  1.940    174       0.017    0.0048    0.0112   0.0107

上のコードをブラウザで実際に動かしてみましょう。左がJacobi、右がSORで、同じ板の上を同時に進みます。

格子ボタンを21×41から81×161に切り替えると、Jacobiのtotal errが0.15から1.81まで増えます。$\omega$ を1.2に下げると、SORも同じ落とし穴にはまります。show error field を押すと、残った誤差が板の中央に大きな山として残っているのが見えます。

1.81という数字が語ること#

出力表のJacobiの行だけを見てください。離散化誤差は0.169 → 0.044 → 0.011と、きっちり1/4ずつ減っています。格子は役目を果たしました。ところが反復誤差は0.126 → 0.492 → 1.810と4倍ずつ増えています。予測値 $e_s/(1-\rho)$ の0.130、0.519、2.075とほぼ一致します。

全誤差はこの二つを合わせた結果です。21×41では離散化誤差が支配的で0.15です。41×81からは反復誤差が支配的になり、0.49、1.81と増えます。最も細かい格子の答えが最も粗い格子より12倍悪いのです。反復回数は629 → 1831 → 4815と増えているにもかかわらず、です。

21×41で全誤差0.1525が離散化誤差0.1691より小さいのは偶然です。反復誤差の符号が離散化誤差を一部打ち消しました。この打ち消しを精度と取り違えてはいけません。

ω=1.2で止めるべきではなかった#

レポートはSORの $\omega$ を1.2までしか試していません。しかしこの問題の最適値はずっと大きいのです。

$\rho_{\text{SOR}}$ は最適 $\omega$ でのSORのスペクトル半径です。$1-\rho$ が $h^2$ ではなく $h$ に比例します。出力では最適 $\omega$ が1.780、1.883、1.940です。反復回数は47、90、174で、Jacobiより13〜28倍少なくなります。

もっと重要なのは、止めたときに残る誤差です。$e_s/(1-\rho)$ は0.005、0.009、0.017で、$1/h$ でしか増えません。そのため全誤差は0.169 → 0.044 → 0.011と離散化誤差にそのまま従います。格子を細かくしたぶん答えが良くなります。$\omega=1.2$ は $\rho$ を少ししか下げないので、反復誤差は0.042 → 0.168 → 0.649と依然4倍ずつ増えます。

最適 $\omega$ の公式は、この記事のような長方形の板と一様格子でのみ厳密です。実際の問題ではレポートの言うとおり試行錯誤が必要です。それでも方向ははっきりしています。1より少し大きい値ではなく、2に近い値から試すべきです。情報の伝わる順序が収束を変える別の例は、LU-SGSのセル番号並べ替えの記事で扱いました。

止めるルールを格子と一緒に変える方法

対策は三つあります。一つ目は、増分ではなく誤差を推定することです。連続する二つの増分の比 $\hat\rho = |d_k|/|d_{k-1}|$ で $\rho$ を推定し、$|d_k|\,\hat\rho/(1-\hat\rho) < e_s$ で判定します。コード1行で済みます。

二つ目は、許容値を格子に合わせることです。反復誤差が離散化誤差より十分小さければよいのです。Jacobiなら $e_s \propto h^4$、最適SORなら $e_s \propto h^3$ で小さくする必要があります。格子を変えるたびに $e_s$ をそのままにすることが、今回の落とし穴の直接の原因です。

三つ目は、判定量を残差に変えることです。残差 $r = b - AT$ は誤差に $A$ を掛けたものです。これも格子によって倍率が変わるので、初期残差に対する相対値を使います。SIMPLE系の外部反復でも同じ問題が起きます。アンダーリラクセーションが収束解を変える記事と合わせて読むと、停止判定と緩和係数が結果に残す痕跡を区別できます。

格子収束の研究では、この点が特に重要です。反復誤差が混ざると、観測される次数が2次から外れます。検証ケースが実際に測っているもので見た次数の汚染と同じ種類です。

ソルバーの「converged」を信じる前に#

ログにconvergedと出たら、もう二つ数字を出力してみてください。最後の二つの増分の比 $\hat\rho$ と $1/(1-\hat\rho)$ です。この倍率が数百を超えるなら、許容値は実際の誤差よりそれだけ楽観的です。格子を一段細かくしてこの倍率が4倍に跳ねたなら、Jacobi系の $h^2$ の落とし穴です。2003年のレポートの観察は間違っていませんでした。ただ、犯人は格子ではなく止めるルールだったのです。