数値解析の基礎 — FEM・FVM・FDMの理論と精度

PDEと数値解法の全体像

CAEのすべての解析は、結局は偏微分方程式（PDE）を数値的に解く作業です。PDEをコンピュータで扱うには、連続体を有限個の点・セル・要素に離散化する必要があります。その手法として FDM・FEM・FVM の3つが主流です。

1.1 PDEの分類

2階線形PDEは係数の性質によって3種類に分類されます：

1.2 各手法の哲学的な違い

有限差分法（FDM）

2.1 テイラー展開と差分近似

テイラー展開：

これから1階微分の近似を導出：

2階微分の2次精度中心差分：

2.2 陽解法の安定性条件 — CFL条件とフォン・ノイマン安定性解析

1次元拡散方程式 $\partial u/\partial t = \alpha \partial^2 u/\partial x^2$ をFTCS（前進時間・中心空間）スキームで離散化：

フォン・ノイマン安定性解析：$u_i^n = G^n e^{ikx_i}$ を代入すると増幅因子 $G$：

安定条件 $|G| \leq 1$ より：

移流方程式 $\partial u/\partial t + c\partial u/\partial x = 0$ のCFL（Courant-Friedrichs-Lewy）条件：

物理的意味：1タイムステップで移流される距離が格子幅以下でなければならない。

2.3 高次差分スキームの一覧

有限要素法（FEM）の数学的基礎

3.1 重み付き残差法とガラーキン法

楕円型PDE の例：$-\nabla \cdot (k\nabla u) = f$ in $\Omega$、境界条件: $u = g$ on $\Gamma_D$。

強形式: $r(u) = -\nabla \cdot (k\nabla u) - f = 0$ がすべての点で成立。

重み付き残差法： 重み関数 $v$ に対して：

弱形式（グリーンの定理で部分積分）：

ガラーキン法： 近似解と重み関数を同じ空間から選ぶ — $u \approx \sum_j u_j N_j(\mathbf{x})$, $v = N_i(\mathbf{x})$：

これが行列形式で $\mathbf{K}\mathbf{u} = \mathbf{f}$ となります。

3.2 形状関数の種類

代表的な形状関数：

• ラグランジュ補間基底： $N_i$ は節点 $j \neq i$ でゼロ、節点 $i$ で1。$C^0$ 連続性。構造・熱解析で標準的。
• ハーミート（Hermite）基底： 変位と勾配を節点値に持つ。$C^1$ 連続性。梁・板解析で使用。
• セレンディピティ要素： 辺中点節点を持つ高次ラグランジュ要素（8節点四辺形など）。

1次元ラグランジュ形状関数（2節点線形要素、$\xi \in [-1,1]$）：

3節点2次要素（$\xi \in [-1,1]$）：

3.3 アイソパラメトリック要素と数値積分

物理座標 $\mathbf{x}$ と参照座標 $\boldsymbol{\xi}$ の変換（アイソパラメトリック変換）：

積分の変換：

ガウス・ルジャンドル求積（$n_g$ 点）：

$2n_g - 1$ 次多項式まで厳密に積分できます（$n_g$ 点ガウス積分の精度）。

有限体積法（FVM）

4.1 保存則の積分形

一般的な保存方程式（移流-拡散）の積分形（各コントロールボリューム CV に適用）：

$\phi$: 一般変数（密度・運動量・エネルギーなど）、$\Gamma$: 拡散係数、$S_\phi$: ソース項。

物理的意味：「CVに蓄積される量 = 流入する移流フラックス - 流出する移流フラックス + 拡散フラックス + 体積ソース」。

4.2 セルセンタ vs セルバーテックス

4.3 界面フラックスの補間スキーム

界面 $f$ でのスカラー値 $\phi_f$ を計算するスキーム比較：

$F_f = \mathbf{v}_f \cdot \mathbf{S}_f$ は界面フラックス。

高解像度スキーム（TVD: Total Variation Diminishing）では、フラックスリミッター $\psi(r)$ を用いて精度と有界性を両立します：

線形連立方程式ソルバー

FEM/FVMで最終的に解くべき大規模疎行列方程式 $\mathbf{A}\mathbf{x} = \mathbf{b}$ のソルバーは、問題のサイズと性質によって選択します。

5.1 直接法

LU分解： $\mathbf{A} = \mathbf{L}\mathbf{U}$（$\mathbf{L}$: 下三角、$\mathbf{U}$: 上三角）。前進/後退代入で解を求める。

コレスキー分解（対称正定値行列用）： $\mathbf{A} = \mathbf{L}\mathbf{L}^T$。LU分解の約半分の計算量。構造解析の静解析でよく使用。

直接法のスケーリング（バンド行列、帯幅 $b$）：

• メモリ: $O(n \cdot b)$
• 計算量: $O(n \cdot b^2)$
• 3次元問題では帯幅 $b \sim \sqrt{n}$ なので計算量 $O(n^2)$ — 大規模問題には反復法が必要

5.2 反復法 — 共役勾配法（CG法）

対称正定値行列 $\mathbf{A}$ に対するCG法のアルゴリズム：

初期値: x₀（通常 x₀ = 0）
残差: r₀ = b - Ax₀
探索方向: p₀ = r₀

for k = 0, 1, 2, ...
  α_k = (r_k^T r_k) / (p_k^T A p_k)   // ステップ幅
  x_{k+1} = x_k + α_k p_k              // 解の更新
  r_{k+1} = r_k - α_k A p_k            // 残差の更新
  β_k = (r_{k+1}^T r_{k+1}) / (r_k^T r_k)
  p_{k+1} = r_{k+1} + β_k p_k          // 探索方向の更新

収束条件：$\|r_k\| / \|b\| < \varepsilon$（通常 $\varepsilon = 10^{-6}$ 〜 $10^{-8}$）

5.3 前処理（プリコンディショニング）

条件数 $\kappa(\mathbf{A})$ が大きい（悪条件の）行列は収束が遅い。前処理行列 $\mathbf{M} \approx \mathbf{A}$ を用いて、等価な系 $\mathbf{M}^{-1}\mathbf{A}\mathbf{x} = \mathbf{M}^{-1}\mathbf{b}$ を解きます。

5.4 疎行列の格納形式

大規模疎行列を効率的に格納するデータ構造：

• CSR（Compressed Sparse Row）： 行ごとに非ゼロ要素の列インデックスと値を格納。行演算（行列-ベクトル積）が高速。FEM・FVMで標準的。
• CSC（Compressed Sparse Column）： 列ごとに格納。LU分解に向く。
• COO（Coordinate List）： (行, 列, 値) のリスト。アセンブリが容易だが演算は遅い。

$n$ 節点・$nnz$ 非ゼロ要素のCSR行列のメモリ：$3 \times nnz + (n+1)$ 整数 + $nnz$ 浮動小数点数。

固有値ソルバー

6.1 べき乗法と逆反復法

べき乗法（Power method）： 最大固有値を求める最もシンプルな手法：

収束率：$|\lambda_2/\lambda_1|$（2番目に大きい固有値と最大固有値の比）。固有値が密集していると収束が非常に遅い。

逆反復法（Inverse iteration）： シフト $\mu$ 付き逆反復で指定固有値 $\lambda \approx \mu$ に最も近いモードを求める：

6.2 ランチョス法

大規模対称固有値問題で最も広く使われる手法。クリロフ部分空間 $\mathcal{K}_m(\mathbf{A}, \mathbf{v}_1) = \text{span}\{\mathbf{v}_1, \mathbf{A}\mathbf{v}_1, \ldots, \mathbf{A}^{m-1}\mathbf{v}_1\}$ に対して三重対角化を行います：

$\mathbf{T}_m$ は $m \times m$ 三重対角行列で、その固有値（ライツ値）が $\mathbf{A}$ の固有値の良い近似になります。$n = 10^6$ 自由度の系でも $m = 100 \sim 300$ 程度のランチョスベクトルで十分な精度が得られます。

6.3 サブスペース反復法とARPACK

実用的な大規模固有値ソルバーの多くはARPACK（Arnoldi PACKage）を基盤とし、非対称行列にも対応した Implicitly Restarted Arnoldi Method (IRAM) を実装しています。

ARPACK を利用した固有値ソルバーの例：

• SLEPc（OpenFOAM, FEniCS など）
• LAPACK/ScaLAPACK（小〜中規模）
• FEAST固有値ソルバー（特定帯域の固有値を高並列で計算）

数値精度・誤差管理

7.1 誤差の種類

7.2 リチャードソン外挿とGCI

メッシュ収束スタディで3つのメッシュ解を使って真解を推定するリチャードソン外挿（Richardson Extrapolation）：

2つのメッシュ（$h_1 < h_2$, 細化比 $r = h_2/h_1$）から：

グリッド収束指数（GCI: Grid Convergence Index）： Roache（1994）が提案した標準的な誤差評価指標：

$e_{12} = (f_2 - f_1)/f_1$（相対誤差）、$F_s = 1.25$（安全係数、Roache推奨）。

GCI < 5% なら「十分な収束が達成されている」と判断するのが一般的な基準です。ASME V&V 20-2009 規格でも採用されています。

7.3 数値的な落とし穴と対策

最適解法の選択ガイド

8.1 問題タイプ別ソルバー選択チャート

8.2 並列計算スケーラビリティ

アムダールの法則： 並列化できる割合 $p$ のとき、$N$ プロセッサでの理論最大速度向上比：

$p = 0.95$（95%並列化可能）でも $N \to \infty$ の理論限界は $1/(1-0.95) = 20$ 倍。

大規模FEM/FVM並列計算での注意点：

• 領域分割（Domain Decomposition）： MeTIS, SCOTCH などでグラフ分割。ハロー交換通信コストを最小化する。
• AMGの並列化： 粗化段階でプロセス間通信が必要。スケーラビリティに限界がある。
• 実効並列効率： 実務では50〜70%程度。通信が支配するとそれ以下になる。

関連インタラクティブツール