MMSソース項の自動導出（シンボリック計算）

MMSとソース項の理論基礎

コード検証における位置づけ

支配方程式を \( L(u) = 0 \) と書くと、任意に選んだ製造解 \( u_m \) は一般にこれを満たしません。そこで残差をソース項として定義し、方程式に加えます。

$$ s(\mathbf{x}, t) := L(u_m), \qquad L(u) = s $$

改造後の方程式に対して \( u_m \) は構成上の厳密解です。ソルバーにソース項 \( s \) と（製造解から導いた）境界条件・初期条件を与えて解き、数値解と \( u_m \) の誤差がメッシュ細分化で理論通りの次数で減ることを確認します。

観測次数による合否判定

細分比2のメッシュ系列で誤差ノルム \( E_h \) を測ると、観測収束次数は

$$ p_{obs} = \frac{\ln(E_{2h}/E_h)}{\ln 2} $$

で得られます。2次精度の離散化なら \( p_{obs} \to 2 \) に漸近するはずで、届かなければ離散化・境界実装・ソース項処理のどこかに欠陥があるという明確な合否基準になります。これがMMSの強みで、「答えがそれらしい」という主観的確認を「次数が出る/出ない」の客観的判定に置き換えられます。

製造解の選定基準

製造解は物理的に意味がある必要はまったくありませんが、検証能力を最大化する選定基準があります。

• 十分に滑らか — 形式次数の収束を観測するため \( C^\infty \) 級（三角関数・指数関数の合成が定番）
• 方程式の全項を働かせる — どの偏微分項もゼロにならない解を選ぶ。定数や1次式では拡散項のバグを検出できない
• 各方向に非自明な変化 — 多次元では方向ごとに異なる波数を与え、方向の取り違えバグを検出
• 物理制約の考慮 — 乱流量や密度など正値であるべき変数は正の製造解にする（ソルバーの正値性制限が働かないように）

シンボリック計算によるソース項導出

なぜ自動導出が必須か

1次元熱伝導程度なら手計算でもソース項を導けます。しかし圧縮性Navier-Stokes方程式に3次元の製造解を代入すると、ソース項は数百項規模の式になり、手作業の微分は誤りなしには不可能です。ソース項自体が間違っていれば「コードのバグ」と「ソース項のバグ」が区別できず検証が成立しないため、数式処理系（CAS）による自動導出はMMSの実質的な前提条件です。SymPy・Mathematica・Mapleのいずれでも手順は同じです。

1. 製造解 \( u_m(\mathbf{x},t) \) をシンボリック式として定義
2. 支配方程式の微分演算子をそのまま適用（\( \partial_t \)、\( \nabla\cdot \)、\( \nabla^2 \) …）
3. 得られた式を整理（simplify）し、コード生成機能でソルバー言語へ出力

導出の具体例——非定常熱伝導

熱伝導方程式 \( \partial T/\partial t - \alpha \nabla^2 T = s \) に製造解 \( T_m = \sin(\pi x)\cos(\pi y)\,e^{-t} \) を選ぶと、ソース項は代入するだけで

$$ s = \frac{\partial T_m}{\partial t} - \alpha \nabla^2 T_m = \left(2\pi^2 \alpha - 1\right)\sin(\pi x)\cos(\pi y)\,e^{-t} $$

と閉形式で得られます。境界条件も同じ製造解から導きます：Dirichlet境界には \( T_m \) の境界値を、Neumann境界には \( -k\,\partial T_m/\partial n \) を与えます。境界条件を製造解と整合させ忘れるのがMMS実装の最頻出ミスです。

コード生成の実務——桁落ちと共通部分式

シンボリック式をそのまま出力すると、実行時間と数値精度の両方で問題が出ることがあります。実務上の定石は3つです。①共通部分式除去（CSE）で三角関数の重複評価をまとめる、②大きな項同士の引き算になっている形は simplify で整理してから出力する（桁落ち回避）、③生成コードを2つ以上の独立な方法（別のCAS、または有限差分による数値微分）と突き合わせてソース項自体を検証する。特に③は「検証ツールの検証」として省略しがちですが、乱数点数点での数値微分との一致確認だけでも導出ミスの大半を捕捉できます。

実務での適用手順

SymPyによる導出とコード出力

無償で完結する標準構成がSymPyです。導出からC/Fortranコード出力までの骨格を示します。

import sympy as sp

x, y, t, alpha = sp.symbols("x y t alpha")
T_m = sp.sin(sp.pi * x) * sp.cos(sp.pi * y) * sp.exp(-t)   # 製造解

L = sp.diff(T_m, t) - alpha * (sp.diff(T_m, x, 2) + sp.diff(T_m, y, 2))
s = sp.simplify(L)                                          # ソース項

print(sp.ccode(s))          # C言語式として出力 (UDF/user subroutine 用)
print(sp.fcode(s))          # Fortran式として出力
s_num = sp.lambdify((x, y, t, alpha), s)                    # Python検証用関数

Navier-Stokes級の長い式では sp.cse(s) で共通部分式に分解してから出力すると、生成コードが桁違いに短く速くなります。

検証ランの標準構成

1. メッシュ系列 — 最低3水準、可能なら4水準。一様細分（細分比2）が次数評価を単純にする
2. 誤差ノルム — \( L^2 \) ノルムを主、\( L^\infty \) ノルムを併記（局所的な次数低下は \( L^\infty \) に先に現れる）
3. 時間と空間の分離 — 空間次数を見るときは時間刻みを十分小さく固定（または定常MMS）、時間次数を見るときはその逆
4. 収束判定の厳格化 — 反復残差は離散化誤差より2桁以上小さく落とす。さもないと測っているのが反復誤差になる
5. 次数プロット — 横軸 \( \log h \)・縦軸 \( \log E \) の傾きを理論次数の参照線と重ねて報告

ソース項の数値積分に注意

有限要素法では、ソース項は要素ごとに数値積分されて荷重ベクトルに入ります。製造解に高い波数を選ぶとソース項が急峻になり、既定の積分次数では積分誤差が離散化誤差を上回って観測次数を汚すことがあります。対策は、積分次数を1〜2段上げて結果が変わらないことを確認するか、製造解の波数を下げることです。有限体積法でもセル中心値×体積の近似がセル内変化の大きいソース項で同じ問題を起こします。

ツール対応とソース項の注入方法

導出側のツール

ソルバー側への注入手段

先端研究の動向

複雑物理への拡張

MMSの研究前線は「どこまで複雑な物理で次数検証を成立させるか」にあります。乱流モデル（Spalart-Allmaras・k-ω等）の輸送方程式に対する製造解は正値性と生成・散逸項のバランス設計が難しく、専用の製造解セットが論文・MASAで整備されてきました。多相流の界面追跡、化学反応系の硬い（stiff）ソース項、移動境界・ALE定式化など、非滑らかさや強い非線形を含む系へのMMS適用が現在も活発な研究対象です。

非滑らかな問題と次数低下の理論

衝撃波を含む圧縮性流れでは、解の不連続により高次スキームでも大域次数が1次に落ちることが理論的に知られています。この場合のMMSは「滑らかな領域で次数が出るか」「不連続の捕捉幅が理論通りか」を分けて検証する形に進化しており、不連続を跨ぐ誤差ノルムの扱い（重み付きノルム、領域分割評価）が研究されています。次数が出ないことが仕様通りのケースを知らないと、正常なコードを不合格にしてしまいます。

CI（継続的インテグレーション）への組込み

近年の数値計算コード開発では、MMS次数テストを自動回帰テストとしてCIに組み込むのが標準的な実践になっています。コミットごとに粗いメッシュ2水準で観測次数を計算し、許容帯（例：理論次数±0.2）を外れたらビルドを失敗させる仕組みです。シンボリック導出からコード生成までをスクリプト化してあれば、方程式や離散化の変更に追従するコストも小さく、「いつの間にか次数が落ちていた」型の退行を機械的に防げます。

トラブル対応

観測次数が出ないときの診断表

次数低下の犯人はたいてい境界

MMSは「コードが正しい」ことは証明しませんが、「この次数で収束する」という反証可能な性質を系統的に検査できる、コード検証の最有力手法です。導出の自動化まで整えれば維持コストは小さく、回帰テストとして資産になります。関連記事：MMSの概要、MMS収束次数のトラブルシュート、メッシュ収束性検証。