[論文レビュー] Newtonを二つに割ると自由水面が解ける — Casulli–Zanolli (2012)

濡れたセルと乾いたセルが一つの系に共存しても、単調収束が保証される理由。

$V(\eta) + T\eta = \mathbf{b}$ 。たった一行の式です。初めて見たとき、Newton 1回で済むと思いました。ところが井戸の縁——水深が 0 を跨ぐセル——で Jacobian が静かに零になります。同じセルが、ある時刻には水中、次の時刻には干上がる(wetting/drying)。そこですべての Newton 系コードは一度倒れます。Casulli と Zanolli が 2012 年に行ったのは、Newton そのものを二つに割ることでした。外側の反復は下から上へ、内側の反復は上から下へ——両側から正解を挟み込み、単調収束を保証します。

論文情報#

著者: Vincenzo Casulli, Paola Zanolli
雑誌: Journal of Computational and Applied Mathematics, vol. 236 (2012), pp. 3937–3947
DOI: 10.1016/j.cam.2012.04.025
キーワード: mildly nonlinear systems / wetting and drying / free-surface hydrodynamics / Richards equation / confined–unconfined aquifers

Q1. なぜ自由水面で Newton はしばしば死ぬのか#

問題はおとなしい見た目をしています。

\mathbf{V}(\eta) + T\eta = \mathbf{b}

ここで $\eta \in \mathbb{R}^N$ はセル毎の piezometric head、 $T$ は対称半正定 sparse 行列(離散拡散ステンシル)、 $\mathbf{V}$ は対角非線形関数——セル毎体積 $V_i(\eta_i) = \int_{-\infty}^{\eta_i} a_i(z)\,dz$ です。関数 $a_i$ は bounded variation(不連続も可)であればよい。代表的な自由水面の例は三つです。

wetting/drying: $V(\eta) = \max(0, \eta)$ 。乾いたセルでは $V'=0$ 。
confined–unconfined aquifer: $V(\eta) = \max(0, \min(c-(-h), \eta-(-h)))$ 。折れ目が二つ。
Richards 方程式の mixed form: $V(\eta)$ は土壌依存で滑らかだが、ほぼ平坦な区間が広い。

三つに共通するのは、 $V'(\eta)$ が区間ごと消えることです。標準 Newton は $J = V'(\eta^{(k)}) + T$ の逆をとりますが、 $V'=0$ かつ $T$ の対応する行が弱結合(あるいは $T$ の零空間に触れる)なら、 $J$ は特異になります。初期値が真値に十分近くないと発散する——これが 30 年来、自由水面コードを悩ませてきた持病です。

Q2. $\mathbf{V}(\eta)$ を二つに切る — Jordan 分解#

最初の仕掛けは関数の分解です。 $a_i$ が有界変動であるとは、二つの非減少関数の差で書けるということ(Jordan 分解)。

a_i(\eta) = p_i(\eta) - q_i(\eta), \quad p_i, q_i \ge 0,\ \text{非減少}

自動的に $V_i$ も分解されます。

V_i(\eta) = V_{1,i}(\eta) - V_{2,i}(\eta), \quad V_{j,i}(\eta) = \int_{-\infty}^\eta p_i\ \text{or}\ q_i

$V_1, V_2$ はいずれも非減少かつ有界。wetting/drying の例では $p(\eta)=1\,[\eta \ge 0]$ 、 $q\equiv 0$ 。confined–unconfined の例では $p(\eta)=1\,[\eta \ge -h]$ 、 $q(\eta)=1\,[\eta > c]$ 。二つに割れば各片は単調——この一文がすべての出発点です。

Q3. 内と外、二度の線形化#

次の系を解きます。

V_1(\eta) - V_2(\eta) + T\eta = \mathbf{b}

外側(outer)反復は $V_2$ を $\eta^{n-1}$ のまわりで線形化します。

V_1(\eta^n) + (T - Q^{n-1})\eta^n = \mathbf{b} + V_2^{n-1} - Q^{n-1}\eta^{n-1}

ここで $Q^{n-1} = \mathrm{diag}(q_i(\eta_i^{n-1}))$ 。出発は $\eta^0 \le \ell$ (下の柵)。これもまだ非線形なので、内側(inner)反復で $V_1$ を線形化します。

(T + P^{n,m-1} - Q^{n-1})\eta^{n,m} = P^{n,m-1}\eta^{n,m-1} - V_1^{n,m-1} + \mathbf{d}^{n-1}

出発は $\eta^{n,0} \ge u$ (上の柵)。 $\ell$ と $u$ は $T+P-Q$ が Stieltjes 行列(対称 M-行列)になるよう設計します。その結果:

内側は上から単調減少( $\eta^{n,m+1} \le \eta^{n,m}$ )
外側は下から単調増加( $\eta^{n+1} \ge \eta^n$ )
両者が解を挟み込み、一点へ収束

収束は浮動小数の議論ではなく不等式で証明されます。「Newton は真値の近くでしか収束しない」という長年の弱点が消えます。下のシミュレーションを直接いじってみてください。

T = 0.50b = 0.90

step 1/4 · outer n=0 · eta=-0.4000 · |r|=1.10e+0

Pink = outer iterate (rises from below); green = inner iterate (falls from above). They squeeze toward the true root of V(eta)+T*eta=b with V(eta)=max(0,min(1,eta)).

T を 0.05 まで下げると $T+P-Q$ は特異に近づきますが、反復は 5 ステップ以内で終わります。b を 1.5 以上に上げると、解が $V$ の平坦部に落ちて——Newton なら永久振動する場所で——内側は一発で止まります。

Q4. 1D 井戸にポンプを差す — Python 実装#

論文 §6 の confined–unconfined aquifer を 1D に縮めて、ゼロから書きます。トイ問題は回転放物面の底/天井+中央のポンプ。日数が進むと phreatic 域が拡がり、やがて乾燥域が現れます。標準的 CFD のトイ問題(Sod, Karman, blast)には当てはまらないので、論文の設定をそのまま使います。

import numpy as np
 
# セル毎体積: V_i(eta) = a0 * clamp(eta, -h_i, c_i) (-h_i だけ平行移動)
def cell_volume(eta, h, c, a0=0.3):
    return a0 * (np.clip(eta, -h, c) - (-h))
 
# Jordan 分解の P, Q (対角 Jacobian)
def p_diag(eta, h, a0=0.3):
    # V1 = a0 * max(0, eta - (-h)) の微分
    return np.where(eta >= -h, a0, 0.0)
 
def q_diag(eta, c, a0=0.3):
    # V2 = a0 * max(0, eta - c) の微分
    return np.where(eta > c, a0, 0.0)
 
def nested_newton_step(eta_prev, h, c, T, dt, phi, tol=1e-6):
    """1 時間ステップの単調 nested Newton 解法。
 
    V(eta_new) + T eta_new = V(eta_prev) + dt*phi を解く。
    ここで T は -dt*Laplacian から来る SPD 行列。
    """
    N = len(eta_prev)
    V_prev = cell_volume(eta_prev, h, c)
    rhs0 = V_prev + dt * phi
 
    # 外側初期値: 下の柵より十分下 (全セル乾燥)
    eta = np.full(N, -h.max() - 1.0)
    for n_out in range(20):
        Q = np.diag(q_diag(eta, c))
        V2 = np.where(eta > c, 0.3 * (eta - c), 0.0)
        d = rhs0 + V2 - Q @ eta
 
        # 内側初期値: 上の柵より十分上 (全セル加圧)
        eta_in = np.full(N, c.max() + 1.0)
        for m_in in range(30):
            P = np.diag(p_diag(eta_in, h))
            V1 = np.where(eta_in > -h, 0.3 * (eta_in - (-h)), 0.0)
            f = P @ eta_in - V1 + d
            A = T + P - Q                     # Stieltjes (対称 M-行列)
            eta_in = np.linalg.solve(A, f)    # 1D: 密; 2D は PCG
            r_in = V1 + (T - Q) @ eta_in - d
            if np.max(np.abs(r_in)) < tol:
                break
        eta = eta_in
        V_cur = cell_volume(eta, h, c)
        r_out = V_cur + T @ eta - rhs0
        if np.max(np.abs(r_out)) < tol:
            return eta, n_out + 1, m_in + 1
    return eta, 20, m_in + 1
 
# 1D 井戸 (放物面断面)
N = 41
L = 1000.0
x = np.linspace(0, L, N)
h_arr = 10 * np.maximum(0, 1 - ((x - L/2) / (L/2))**2)
c_arr = h_arr.copy()
 
# 拡散行列 T = dt*kappa/dx^2 * Laplacian (zero-flux 境界)
dx = L / (N - 1)
dt = 86400.0     # 1日
kappa = 1.0
alpha = dt * kappa / dx**2
T = -alpha * np.eye(N, k=-1) + 2*alpha*np.eye(N) - alpha*np.eye(N, k=1)
T[0, 0] = alpha; T[0, 1] = -alpha
T[-1, -1] = alpha; T[-1, -2] = -alpha
 
# 初期: 天井に届く加圧状態
eta = c_arr.copy()
phi = np.zeros(N); phi[N//2] = -8e-4   # 中央のポンプ (sink)
 
for day in range(10):
    eta, n_out, m_in = nested_newton_step(eta, h_arr, c_arr, T, dt, phi)
    pre = int(np.sum(eta >= c_arr - 1e-6))
    dry = int(np.sum(eta <= -h_arr + 1e-6))
    print(f"day {day+1}: pre={pre} dry={dry} outer={n_out} inner_last={m_in}")

実行例:

day 1: pre=39 dry=0 outer=5 inner_last=5
day 2: pre=37 dry=0 outer=4 inner_last=4
...
day 9: pre=0 dry=8 outer=1 inner_last=5
day 10: pre=0 dry=12 outer=1 inner_last=5

**外側は平均 3 回、内側は約 5 回で終わります。**論文 Table 1 と定性的に一致します。最も興味深いのは、加圧→phreatic の遷移期に内外どちらも一発で止まる場面がよくあること(論文 Remark 6: $u \le \tilde\eta^{1,1} \le \ell$ なら一発で根に到達)。

Q5. 2D 断面を実際に進化させてみる#

同じアルゴリズムを 1D 井戸断面で日次に動かすデモです。ポンプ強度はスライダー、+1 day で 1 日進めます。

pump flux phi = 0.0008 (m/s per cell)

day = 0 · pressurized cells = 58 · phreatic = 0 · dry = 2 · total water V = 3998.9 m

Cyan = pressurized (confined) region; teal = phreatic; black = dry. Yellow line is piezometric head eta(x). Pump near the center steadily drains the aquifer; the wet/dry boundary advances inward as days pass.

ポンプを 0.0025 以上にして 10 日回すと、中央に黒い乾燥域が現れて広がります。湿セル集合は毎日変わりますが、求解側は知らんふりです。単純 Newton なら、湿/乾ステンシルが変わるたびに Jacobian を作り直す必要があります。nested ループは同じコードでそのまま解き続けます。

このアルゴリズムができないこと — 批判の一段#

万能ではありません。落とし穴をいくつか。

$P-Q = 0$ となるセルの扱い。乾燥セルで $T$ の対応する行も消えると、行列は T2(既約)条件を失います。論文は乾燥セルを系から除外することを薦めますが、有効セル集合のブックキープは実装上は厄介です。
速くはない。1 ステップあたり外側×内側で平均 15〜25 回の線形求解が走ります。CG 付き implicit Euler が収束する場面では明らかに速いことがあります。代償は発散ゼロの保証です。
圧縮性 Euler に直接は使えない。本手法は対角な $V$ で有界変動という縛りに限られます。Riemann ソルバの非線形性はこのようには分解できません。
内側ソルバ自体は別問題。各内側系は SPD ですが、どう前処理するかは別の技術領域です。論文は PCG を薦めますが、大規模格子ではマルチグリッドが上位です。

再現可能性スコア#

数学(§2–§5): 自己完結。鉛筆と紙で不等式を全部追えます。A.
アルゴリズム箱(1, 2): pseudocode が完備。A.
数値例(§6): 格子幅、時間ステップ、全係数が明示。30 行の Python で定性的再現に成功。A.
実務応用: OpenFOAM には first-class な wetting/drying モジュールがなく、meshFlow/overFlowSolver 系に組み込む余地があります。海洋・河口モデル(SCHISM、TELEMAC)では類似アイデアが既に採用されています。

次に読むべき論文#

Brugnano & Casulli (2008/2009) — piecewise linear system の有限ステップ終了定理。本論文の直接の祖先。
Casulli (2009) — 同著者の high-resolution wetting/drying アルゴリズム。
Casulli & Zanolli (2010) — Richards 方程式 mixed form への nested Newton 適用(滑らかな $a(\eta)$ 版)。