› Richardson補外と観測収束次数シミュレーター 戻る
数値解析

Richardson補外と観測収束次数シミュレーター

3つの刻み幅で得た数値解から、観測収束次数 p と Richardson補外による真値推定 u_ext をリアルタイムで計算します。CFD・FEMのメッシュ収束（V&V）、ODE解法の精度評価、中点則・台形則・RK法の次数チェックに使える数値解析ツールです。

パラメータ設定
粗ステップ h₃
最も粗い刻み幅（h₃ = 最大）
中ステップ h₂
中間の刻み幅（通常 h₃/2）
細ステップ h₁
最も細い刻み幅（h₁ = 最小）
粗解 u₃
h₃ で得られた数値解
中解 u₂
h₂ で得られた数値解
細解 u₁
h₁ で得られた数値解
計算結果
—
ステップ比 r = h₂/h₁
—
観測収束次数 p
—
Richardson補外値 u_ext
—
細解の打切り誤差
—
漸近誤差比 e₁/e₂
—
次数の判定
—
log-log プロット — 誤差 vs ステップ h（傾き = p）

青点：3刻み h₁,h₂,h₃ の誤差。黄線：3点を通す回帰直線（傾き = 観測収束次数 p）。緑印：Richardson補外値（h→0）。

解 vs ステップ h（補外値を含む）
誤差 vs ステップ h（両対数軸）
理論・主要公式

$$u_{exact} \approx u_1 + \frac{u_1-u_2}{r^p-1},\qquad p=\frac{\ln\bigl(|u_2-u_3|/|u_1-u_2|\bigr)}{\ln r}$$

r は刻み幅の比 (h₂/h₁)、p は観測収束次数。u₁,u₂,u₃ は細・中・粗の数値解。3解が漸近領域にあるとき、u_ext は元の解より1〜2桁小さい誤差で真値を近似する。

$$e_h \approx C\,h^{p},\qquad \frac{e_1}{e_2}=\frac{1}{r^{p}}$$

離散誤差 e_h は刻み h の p乗で減少する。漸近誤差比 e₁/e₂ が 1/r^p に近ければ漸近領域に入っており、観測次数は信頼できる。

$$\mathrm{GCI}=\frac{F_s\,|u_1-u_2|}{(r^p-1)\,|u_1|}$$

Grid Convergence Index（Roache 1998）。安全係数 F_s=1.25（3メッシュ）で、Richardson補外誤差の信頼バンドを与える。CFD の V&V 報告で標準。

Richardson補外と収束次数

🙋
「Richardson補外」って初めて聞きました。なんで3つも刻み幅で計算する必要があるんですか？1つで十分じゃないんですか？
🎓
いい疑問だね。1つの刻み h で計算した数値解 u_h には、未知の誤差 e_h が含まれている。これを「どれくらい外れているか」「どの方向に外れているか」を知るには、最低でも2つの刻みでの解の差を見ないといけない。そして「収束次数 p を未知数として推定する」となると、3点必要になるんだ。例えば h=0.1 で u=0.8125 と出ても、「真値は 0.81 か 0.82 か 0.85 か」は1点だけでは判断できない。3つあれば「ああ、これは 0.811 に向かって p=2 で収束しているな」と逆算できる。
🙋
なるほど。じゃあ「観測収束次数」と「理論収束次数」って違うんですか？教科書には「中点則は2次精度」って書いてあるけど、それと別物？
🎓
そこが大事なところ。「理論次数」は紙の上で誘導した値（中点則=2、RK4=4 など）で、「観測次数」は実際にコードを走らせて出てきた値だ。両者が一致するのは、(1)プログラムにバグがない、(2)刻みが十分細かい（漸近領域）、(3)解が滑らか、の3条件が揃った時だけ。例えば「2次精度のはず」のコードで観測次数が 0.8 しか出なかったら、ほぼ確実に境界条件が1次に落ちている。Richardson補外は「自分のコードが本当に主張通りの精度を持っているか」を検証する診断ツールなんだ。
🙋
デフォルト値だと p=1.983 で「2次」って出ますね。これは中点則とかRK2とかですか？
🎓
そう、まさにそれ。この例は ∫₀¹ 1/(1+x²) dx = π/4 ≈ 0.7854 を中点則で計算した値だよ。中点則は理論的に2次精度（誤差 ∝ h²）で、観測値も 1.983 とほぼ2に乗った。Richardson補外値 u_ext=0.8111 も h=0.1 の解 0.8125 より真値に近い…と言いたいところだけど、この問題は実は対称性のおかげで中点則が「超収束」して、真値より少し離れた極限に向かっている。これは「観測次数だけ見て安心せず、収束先 u_ext が物理的に正しいか別途確かめる」べき教訓でもあるんだ。
🙋
CFDで「メッシュ収束をやれ」って言われるけど、実務だとどう使うんですか？
🎓
手順は3ステップだ。(1)同じ物理条件で、節点数が約2倍ずつ違う3メッシュ（粗・中・細）を作る、(2)代表量（抗力係数、圧損、最大温度など）を u₃,u₂,u₁ として記録、(3)このツールに入力して p と GCI を出す。学会・査読論文だと「観測次数 p が理論次数の±0.5以内、GCI が許容誤差以下」を達成して初めて「メッシュ非依存」と認められる。逆に言うと、シングルメッシュで「いい絵が出ました」だけのCFD報告は、V&V 観点では信頼度ゼロ。AIAA G-077 や ASME V&V20 はこの手順を標準化した文書だよ。
🙋
漸近誤差比 e₁/e₂ が「1/r^p に近いか」で漸近領域を判定するって書いてありますね。これがズレてたらどう解釈すれば？
🎓
e₁/e₂ が期待値（r=2,p=2 なら 0.25）から大きく外れていたら、3点はまだ漸近領域に入っていない。つまり「次数 p の推定そのものが信頼できない」ということだ。対処法は2つあって、(a) もっと細いメッシュで4点目を取って 2:1 比のペアを追加、(b) 解に不連続や強い特異点がないか確認。例えば衝撃波がある流れや、応力集中のある構造解析では、解の正則性で次数が頭打ちになり、p が 1〜1.5 で止まることもよくある。その場合は無理に2次に押し込めず、「観測次数1.2で運用する」と素直に報告するのが正しい姿勢だね。

よくある質問

3つの刻み幅 h₁<h₂<h₃（ステップ比 r=h₂/h₁=h₃/h₂）で得た数値解 u₁,u₂,u₃ から、p = ln(|u₂−u₃|/|u₁−u₂|) / ln(r) で求まります。理論次数 p_theory（例：オイラー法=1、台形則やRK2=2、RK4=4）と一致すれば「コードと格子は漸近領域にある」と判定できます。p が理論値から大きくずれている場合、刻みがまだ粗い・コードにバグがある・境界条件で次数が落ちている、のいずれかを疑います。本ツールは r が等比でないときも同じ式で見積もります。
Richardson補外値 u_ext = u₁ + (u₁−u₂)/(r^p − 1) は、漸近誤差展開の主項を消去した推定値です。3つの解が「漸近領域」にあるとき、u_ext の誤差は元の解の誤差より 1〜2桁小さくなります。判定の目安は漸近誤差比 e₁/e₂ が 1/r^p に近いこと。例えば r=2, p=2 なら e₁/e₂ ≈ 0.25 が期待値で、これから大きく外れる場合は次数を信用してはいけません。Roache (1998) の GCI (Grid Convergence Index) は安全係数 Fs=1.25 を掛けた誤差バンドを使います。
主な原因は4つです。(1) 刻みがまだ粗く漸近領域に入っていない — さらに細かい h で再計算。(2) 境界条件・初期条件が低次精度で実装されている（例：1次の階段近似） — これは離散誤差全体を支配し、観測次数が0.5〜1に落ちます。(3) 解に特異点や不連続がある（衝撃波、応力集中、リエントラント角） — 解の正則性で次数が決まり、滑らかな解の理論次数は得られません。(4) 反復ソルバーの収束が不十分で離散誤差より反復誤差が大きい — 残差を1桁以上小さく取り直します。
はい、まさにこの目的のためのツールです。CFDの V&V では「3メッシュで同じ解析を行い、観測収束次数と Richardson補外値、GCIを報告すること」が AIAA G-077 や ASME V&V20 で要求されます。手順は、(1)粗・中・細の3メッシュで同条件・同モノマー計算、(2)代表量（抗力係数、ピーク応力、Nu数 など）を u₃,u₂,u₁ として入力、(3) p と u_ext と GCI を確認、です。観測次数 p が理論次数（FVM/FEMで通常2）と±0.5以内に収まれば信頼できる結果と判定します。

実世界での応用

CFD の V&V（検証・妥当性確認）：航空機の抗力計算、ターボ機械の効率、自動車外装の空力など、商用CFD（ANSYS Fluent、OpenFOAM、STAR-CCM+）の報告には観測収束次数と GCI の記載がほぼ必須です。AIAA G-077（航空宇宙）、ASME V&V20（流体・熱）、NAFEMS のガイドラインはいずれも「3メッシュ最小」「観測次数の確認」「GCI 報告」を要求しています。シングルメッシュの「いい絵」報告は学会では却下対象です。

有限要素法の応力解析：応力集中部のピーク応力やJ積分は、メッシュ細分化で値が変動しやすい代表量です。粗・中・細の3メッシュで再計算し、Richardson補外で「真の応力」を推定。リエントラント角などの特異点があると観測次数は1以下に落ち込み、補外しても収束しないため「報告値は最も細いメッシュの値とその不確かさ」とすべきです。NAFEMS のベンチマーク問題の多くがこの種類の収束評価を含みます。

ODE 解法の精度評価：新しい時間積分スキーム（Runge-Kutta、Adams-Bashforth、IMEX など）を実装したとき、ステップ Δt を半分にしたときに誤差がきちんと 2^p で減るかチェックします。バグがあると観測次数が「中途半端な値」（例：3.2 のはずが 1.8）になり、デバッグの強力な手がかりになります。教科書の演習問題や論文の数値実験で頻出する標準的な検証手法です。

気象・海洋モデル：大気大循環モデル（GCM）や海洋モデル（NEMO、MITgcm）でも空間・時間解像度の収束評価が行われます。ただしこれらは強い非線形・乱流のため、観測次数が理論次数より低く出るのが普通で、漸近領域に入るのに膨大な計算資源が必要です。実務では「収束しないが、関心の物理量が定常的に変化しなくなる」点で打ち切ることも多く、Richardson補外の限界を学ぶ良い例でもあります。

よくある誤解と注意点

最も多い誤りが、「3点の差が小さければ収束している」と思い込むことです。例えば u₃=1.000, u₂=1.001, u₁=1.0011 と並んでいても、u₃→u₂ で 0.001 動き、u₂→u₁ で 0.0001 しか動かないと、観測次数は ln(10)/ln(2) ≈ 3.3 と計算されます。これを見て「3次精度だ！」と喜ぶのは早計で、もしかすると単に「3点全部が同じ間違った値（バイアス）に張り付いているだけ」かもしれません。Richardson補外は「収束方向と次数」を当てるツールであり、「真値を保証するツール」ではない。物理現象として妥当な値か、別の手法（解析解・実験・別ソルバー）と必ず突き合わせてください。

次に、「等比でないステップで p の式を使う」誤り。本ツールの p の式は、ステップ比 r=h₂/h₁=h₃/h₂ が等しいことを暗黙に仮定しています。CFD で「適応メッシュ」「ハイブリッドリファイン」を使うと、粗・中・細の比率が 2 / 1.7 などにバラけ、観測次数の推定値が大きく狂います。実務では Roache (1998) の不等比メッシュ用の式（反復解法で p を求める）を使うか、可能な限り厳密な r=2 のメッシュペアを作ってください。本ツールでは r ≒ const のときだけ意味のある値を返します。

最後に、「観測次数が高ければ細かいメッシュは不要」という誤解。p=4（RK4 や高次有限要素）が出ても、「現在のメッシュで誤差が許容値以下か」は別問題です。誤差は e ≈ C·h^p と書け、C（誤差定数）が大きいと p が高くてもメッシュは粗いほど大きな誤差になります。GCI の式に E_h = (u₁−u₂)/(r^p−1) として現れる「現在の誤差そのもの」をチェックすべきで、p だけで安心してはいけません。「次数が高い手法は同じ精度を粗いメッシュで達成できる」のは「漸近領域に入った後」の話で、入り口の粗いメッシュでは1次手法より誤差が大きいことすらあります。

使い方ガイド

  1. 粗・中・細の3つの刻み幅h₁、h₂、h₃を入力（例：CFDメッシュの要素数から逆算、刻み幅比r=2推奨）
  2. 各刻み幅で得た数値解u₁、u₂、u₃を入力（例：梁の最大たわみ、圧力抵抗、熱流束）
  3. シミュレーターが観測収束次数p、Richardson補外値u_ext、打切り誤差εを自動計算し、数値手法の収束性を判定

具体的な計算例

4点支持梁のFEM解析で、メッシュ分割数8→16→32要素（刻み幅比r=2）のとき、最大たわみu₁=10.24mm、u₂=10.12mm、u₃=10.08mmが得られた場合：p=（ln|u₁-u₂|/|u₂-u₃|）/ln(r)≈2.0となり、2次精度を確認。Richardson補外値u_ext≈10.07mmが真値推定値。RK2法やCrank-Nicolson法で同様に次数検証可能。

実務での注意点