同じ PDE、3つの道#

PDE を計算機に載せるとき、先に決まるのは未知数をどこに置くかです。式が同じでも、点値か、基底の係数か、セル平均かで手法は別物になります。

1D の移流拡散を一本だけ置きます。

$$\frac{\partial u}{\partial t} + a\frac{\partial u}{\partial x} = \nu\frac{\partial^2 u}{\partial x^2}$$

ドロップダウンで FDM、FEM、FVM を切り替えると、同じ格子の上で自由度が節点、基底係数、セル中心のあいだを移ります。

FDM は点が並び、FVM はセルの真ん中に値を一つ置きます。FEM は節点の係数と、要素をまたぐ基底です。

---

1. 有限差分法（Finite Difference Method, FDM）#

コア・アイディア

“

微分を 差分（difference） で直接近似する。

最も直感的なアプローチです。格子点（node）における関数値をテイラー展開（Taylor expansion）し、導関数を近似します。

数学的な出発点：テイラー展開

点 $x_i$ におけるテイラー展開は以下の通りです。

$$u(x_i + \Delta x) = u(x_i) + \Delta x \frac{\partial u}{\partial x}\bigg|_i + \frac{\Delta x^2}{2}\frac{\partial^2 u}{\partial x^2}\bigg|_i + \cdots$$

ここから差分近似を導き出します。

前進差分（Forward Difference）：

$$\frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1} - u_i}{\Delta x} + O(\Delta x)$$

中心差分（Central Difference）：

$$\frac{\partial u}{\partial x}\bigg|_i \approx \frac{u_{i+1} - u_{i-1}}{2\Delta x} + O(\Delta x^2)$$

二階導関数：

$$\frac{\partial^2 u}{\partial x^2}\bigg|_i \approx \frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2} + O(\Delta x^2)$$

移流拡散方程式への適用

中心差分を使用すると、以下のようになります。

$$\frac{du_i}{dt} + a\frac{u_{i+1} - u_{i-1}}{2\Delta x} = \nu\frac{u_{i+1} - 2u_i + u_{i-1}}{\Delta x^2}$$

これは格子点の未知数 $u_i$ に関する連立常微分方程式（ODE）システムとなります。

メリットとデメリット

メリット：

• 概念が単純で実装が容易
• 構造格子（structured grid）において非常に効率的
• 高次精度の達成が容易（コンパクトスキーム、スペクトルライクスキームなど）
• 直交格子において行列構造が綺麗（帯行列）

デメリット：

• 非構造格子（unstructured grid）への適用が困難 — これが致命的な欠点です
• 複雑な形状に対する格子生成が難しい
• 保存則を自動的には満足しない

代表的な適用分野

• DNS/LES（直交格子ベースの乱流シミュレーション）
• 気象・海洋モデル（構造格子）
• 地震波伝播シミュレーション
• コンパクトスキームベースの高精度計算

---

2. 有限要素法（Finite Element Method, FEM）#

コア・アイディア

“

解を 基底関数（basis function）の線形結合 で近似し、重み付き残差（weighted residual） を最小化する。

FEM は微分方程式を直接解くのではありません。代わりに 弱形式（weak formulation） と呼ばれる積分形式に変換します。

数学的な出発点：弱形式（Weak Form）#

元の PDE（強形式、strong form）にテスト関数 $w$ を掛けて積分します。

$$\int_\Omega w \left(\frac{\partial u}{\partial t} + a\frac{\partial u}{\partial x} - \nu\frac{\partial^2 u}{\partial x^2}\right) dx = 0$$

拡散項に部分積分（integration by parts）を適用すると、以下のようになります。

$$\int_\Omega w \frac{\partial u}{\partial t}\,dx + \int_\Omega w\, a\frac{\partial u}{\partial x}\,dx + \int_\Omega \nu\frac{\partial w}{\partial x}\frac{\partial u}{\partial x}\,dx = \left[\nu w \frac{\partial u}{\partial x}\right]_\Gamma$$

これが 弱形式 です。以下の重要な変化に注目してください。

• 元々は二階導関数が必要でしたが、部分積分後は 一階導関数のみ で済みます。
• 解の連続性に対する要求が緩和されます（$C^1 \to C^0$）。
• 境界条件が自然に含まれます（ノイマン境界条件 = 右辺）。

ガラーキン近似（Galerkin Approximation）#

解を基底関数 $\phi_j$ の線形結合で表現します。

$$u^h(x, t) = \sum_{j=1}^{N} U_j(t)\,\phi_j(x)$$

ガラーキン法では、テスト関数と基底関数を同じ空間から選びます（$w = \phi_i$）。

$$\sum_j \left(\int_\Omega \phi_i \phi_j\,dx\right) \frac{dU_j}{dt} + \sum_j \left(\int_\Omega \phi_i\, a\, \phi_j'\,dx + \int_\Omega \nu\, \phi_i'\, \phi_j'\,dx\right) U_j = 0$$

行列形式にすると以下の通りです。

$$\mathbf{M}\frac{d\mathbf{U}}{dt} + \mathbf{K}\mathbf{U} = \mathbf{f}$$

• $\mathbf{M}$：質量行列（mass matrix）
• $\mathbf{K}$：剛性行列（stiffness matrix）
• $\mathbf{f}$：境界・ソースベクトル

基底関数の選択

最も多く使われるのは、ラグランジュ多項式（Lagrange polynomial） ベースの要素です。

メリットとデメリット

メリット：

• 非構造格子と相性が良い — 三角形や四面体メッシュで複雑な形状を処理できる
• 数学的に厳密な誤差評価（a priori / a posteriori 誤差評価）が可能
• 適応格子（adaptive mesh refinement）との相性が非常に良い
• p-refinement（多項式の次数を上げる）と h-refinement（格子を細分化する）を柔軟に選択できる

デメリット：

• 保存則を直接的には満足しない（標準的なガラーキン法）
• 移流支配問題で不安定（振動）が生じる — SUPG や GLS などの安定化が必要
• 質量行列の逆行列が必要（または質点化、mass lumping）
• FDM/FVM に比べて実装の複雑度が高い

代表的な適用分野

• 固体力学・構造解析（FEM の元々の本拠地）
• 電磁界解析
• 熱伝導
• バイオメカニクス
• 地盤・地質工学

---

3. 有限体積法（Finite Volume Method, FVM）#

コア・アイディア

“

保存則を コントロール・ボリューム（control volume） に対して積分し、ボリューム境界を通過する フラックス（flux） のバランスを合わせる。

FVM は物理的な保存則を最も直接的に反映しています。

数学的な出発点：積分型保存則

移流拡散方程式を積分型で書くと以下のようになります。

$$\frac{d}{dt}\int_{V_i} u\,dx + \oint_{\partial V_i} \left(au - \nu\frac{\partial u}{\partial x}\right) \cdot n\, dS = 0$$

コントロール・ボリューム $V_i = [x_{i-1/2},\, x_{i+1/2}]$ に対して、以下の通りとなります。

$$\frac{d\bar{u}_i}{dt} = -\frac{1}{\Delta x}\left(\hat{F}_{i+1/2} - \hat{F}_{i-1/2}\right)$$

ここで、$\bar{u}_i$ はセル平均値、$\hat{F}_{i+1/2}$ はセル境界での 数値フラックス です。

数値フラックスの決定

FVM の核心は、セル境界のフラックスをいかに計算するかという点にあります。

移流項 — リーマン・ソルバー（Riemann solver）やアップウィンド（upwind）ベース：

$$\hat{F}^{adv}_{i+1/2} = \text{RiemannSolver}(\bar{u}_i,\, \bar{u}_{i+1})$$

拡散項 — 中心差分：

$$\hat{F}^{diff}_{i+1/2} = -\nu\frac{\bar{u}_{i+1} - \bar{u}_i}{\Delta x}$$

再構築（Reconstruction）#

セル平均値だけでは一次精度に留まります。高次精度を達成するためには、セル内部の分布を 再構築（reconstruction） する必要があります。

MUSCL（2次精度）：

$$u_{i+1/2}^L = \bar{u}_i + \frac{1}{2}\psi(r_i)\,(\bar{u}_i - \bar{u}_{i-1})$$

ここで、$\psi$ は勾配制限関数（slope limiter）、$r_i$ は連続する勾配の比です。

WENO（5次精度）：

3つの2次多項式の候補を非線形に加重平均することで、5次精度を達成します。

$$u_{i+1/2}^L = \sum_{k=0}^{2} \omega_k\, q_k(x_{i+1/2})$$

滑らかさの指標（smoothness indicator） $\beta_k$ に応じて重み $\omega_k$ が決定され、不連続付近では自動的に滑らかなステンシルに重みを置きます。

面の拡散は、直交間隔で割った $-\nu(\bar{u}_{i+1}-\bar{u}_i)/\Delta x$ です。法線が格子線と揃うときだけ、その差が法線勾配になります。面が傾くと、隣の平均の差だけでは向きが合いません。その穴を埋めるのが 非直交拡散の補正 です。補正を外すと、保存形でも拡散の向きがずれます。

前に書いた前進差分は $O(\Delta x)$、中心差分と二階導関数ステンシルは $O(\Delta x^2)$ でした。前進はテイラーの $\Delta x^2/2$ を残します。中心は $u_{i+1}$ と $u_{i-1}$ を引いてその項を消します。間隔を半分にすると、前進の誤差はおおよそ半分、$O(\Delta x^2)$ 側はおおよそ 1/4 です。セル平均だけでは再構築なしに一次のままです。高次は面へ持ち上げるときに付きます。

MUSCL の $1/2$ は、平均から面まで線形に上げる係数です。$\psi(r_i)=1$ なら傾きを全部使い、$0$ なら面の値は $\bar{u}_i$ に戻って再び一次です。WENO が足す候補は $k=0,1,2$ の三つだけです。各 $q_k$ は二次多項式で、$\omega_k$ は $\beta_k$ が大きいステンシルの重みを削ります。滑らかな区間で三つが五次にまとまる、という文は加重和の次数です。$\psi$ が 0 に落ちると、五次の式はあってもその面は一次に下がります。三角形の節点は P1 が 3、P2 が 6、P3 が 10 です。次数を 1 から 3 に上げると、要素の未知数は 3 を超えて 10 になります。弱形式の $\mathbf{M}$ は $\int \phi_i \phi_j\,dx$、$\mathbf{K}$ は $\phi_i \phi_j'$ と $\nu\,\phi_i'\phi_j'$ を束ねます。$\nu$ が正なら拡散ブロックはほぼ対称です。移流の $a$ は非対称を残します。

メリットとデメリット

メリット：

• 保存則が離散レベルで厳密に満足される — CFD において決定的なメリット
• 衝撃波や不連続面の捉え方に優れる
• 非構造格子への適用が可能
• 物理的な直感と直結するフラックスベースの考え方

デメリット：

• 高次精度の達成が FEM よりも難しい（再構築が複雑になる）
• 拡散方程式においては FEM よりも非効率な場合がある
• 非構造格子における高次再構築の実装が非常に煩雑

代表的な適用分野

• CFD 全般（OpenFOAM, ANSYS Fluent, SU2 など）
• 圧縮性・非圧縮性流れ
• 多相流、燃焼、反応流
• 気象・気候モデルの力学コアの一部

---

主な比較

数学的な出発点

「何を」離散化するか

これが最も本質的な違いです。

• FDM：微分演算子（$\partial/\partial x$）を離散化
• FEM：解空間（solution space）を離散化
• FVM：積分型方程式のフラックスを離散化

未知数の位置

• FDM：格子点（node）における点値（point value）
• FEM：節点（node）における点値（基底関数の係数）
• FVM：セル（cell）の平均値（cell average）

総合比較

---

境界が曖昧になる現代の手法

近年では、3つの手法のメリットを組み合わせたハイブリッドな手法が活発に研究されています。

不連続ガラーキン法（Discontinuous Galerkin, DG）#

FEM の基底関数 ＋ FVM のフラックスの概念を組み合わせた手法です。各セル内部では多項式で解を近似（FEM）し、セル境界ではリーマン・ソルバーでフラックスを計算（FVM）します。

$$\int_{V_i} w \frac{\partial u^h}{\partial t}\,dx - \int_{V_i} \frac{\partial w}{\partial x} F(u^h)\,dx + \hat{F}_{i+1/2} w_{i+1/2}^- - \hat{F}_{i-1/2} w_{i-1/2}^+ = 0$$

保存性 ＋ 高次精度 ＋ 非構造格子 のすべてを満足します。ただし、計算コストが高く、衝撃波付近ではリミッター（limiter）が必要です。

スペクトル差分法（Spectral Difference） / フラックス再構築法（Flux Reconstruction）#

DG と似ていますが、積分を行わずに微分形式で解くことで効率を高めた手法です。GPU 加速との相性が非常に良く、次世代の CFD ソルバーとして注目されています。

粒子法（Meshless Methods, SPH など）#

格子そのものをなくすアプローチ。粒子（particle）ベースで近似します。自由表面流や大変形問題に強みがありますが、精度や一貫性の問題があります。

---

保存が要るなら差分ステンシルは畳む

直方体の DNS で格子も直交なら、中心差分の $O(\Delta x^2)$ で足ります。三点が二階導関数を作り、行列は帯のままです。形状が折れ、非構造になると、そのステンシルは隣を定義できません。そこで弱形式が幾何を飲み込みます。ただし標準ガラーキンの $\mathbf{M}$ と $\mathbf{K}$ は、セルごとのフラックス相殺を約束しません。

衝撃波が通り、積分値が数字のまま残らなければならないなら FVM です。面フラックス $\hat{F}_{i+1/2}$ は左セルの損失であり右セルの利得なので、内部面は総和で消えます。高次と保存を一つのスキームに置くなら、不連続ガラーキン のようにセル内は多項式、セル面はリーマンフラックスです。リミッターを外すと、その混合も衝撃の前で振動します。