有限要素法（FEM）の数値解法

有限要素法は、偏微分方程式を弱形式（積分形式）に変換し、要素分割された空間上の基底関数で場を近似して連立方程式として解く枠組みである。複雑形状・材料不均一・異方性・多物理連成を同一の形式で扱える点が核である。

参考ドキュメント

• S. C. Brenner, L. R. Scott, The Mathematical Theory of Finite Element Methods (Springer) https://link.springer.com/book/10.1007/978-1-4757-3658-8
• PETSc Manual: KSP (Linear System Solvers) https://petsc.org/release/manual/ksp/
• 有限要素法入門（講義資料、日本語PDF） https://nkl.cc.u-tokyo.ac.jp/FEM/01-FEMintro.pdf

1. FEMが「数値解法」として強い理由

1. 形状の自由度
   三角形・四面体などの非構造メッシュにより、欠陥・界面・多結晶粒界のような複雑境界を自然に表現できる。

2. 物理法則と整合する離散化
   保存則やエネルギー最小化原理に基づく弱形式は、境界条件（特にフラックス型）を一貫して組み込める。

3. 連成・異方性の取り込みが局所積分に落ちる
   材料テンソル（弾性定数 $C_{ijkl}$、誘電率・透磁率テンソル、拡散テンソルなど）を要素積分に組み込むだけで、複雑な材料応答を同一実装で扱える。

4. 大規模化の焦点が明確
   最終的に「疎行列の線形／非線形方程式を速く・頑健に解く」へ帰着するため、Krylov法・前処理・領域分割などのHPC技術と直結する。

2. 強形式から弱形式へ

代表としてPoisson型を例示する。

$$-\mathrm{\nabla }\cdot (\kappa \mathrm{\nabla }u)=f\text{in }\mathrm{\Omega }$$

境界を $\partial \mathrm{\Omega }={\mathrm{\Gamma }}_D\cup {\mathrm{\Gamma }}_N$ と分け、 Dirichlet境界で $u=u_D$、Neumann境界で $\kappa \mathrm{\nabla }u\cdot n=g$ とする。

試験関数 $v$ を掛けて領域積分し、部分積分（Greenの公式）を行う：

$${\int }_{\mathrm{\Omega }}\kappa \mathrm{\nabla }u\cdot \mathrm{\nabla }vd\mathrm{\Omega }={\int }_{\mathrm{\Omega }}fvd\mathrm{\Omega }+{\int }_{{\mathrm{\Gamma }}_N}gvd\mathrm{\Gamma },(\mathrm{\forall }v\in V_0)$$

ここで $V_0$ は $v{|}_{{\mathrm{\Gamma }}_D}=0$ を満たす関数空間である。

この形を、抽象化して

$$a(u,v)=\ell (v)$$

と書く。ここで $a(\cdot ,\cdot )$ は双一次形式、$\ell (\cdot )$ は線形汎関数である。以後の議論はこの抽象形で統一できる。

3. Galerkin近似

有限次元部分空間 $V_h\subset V$ を選び、基底 $\{{\varphi }_i{\}}_{i=1}^N$ により

$$u_h(x)=\sum _{i=1}^N{U}_i{\varphi }_i(x)$$

と近似する。Galerkin法では試験関数にも同じ空間を用い、$v_h={\varphi }_j$ と置く：

$$a(u_h,{\varphi }_j)=\ell ({\varphi }_j)(j=1,\dots ,N)$$

すると線形系

$$AU=b$$

が得られる。ここで

$$A_{ji}=a({\varphi }_i,{\varphi }_j),b_j=\ell ({\varphi }_j)$$

である。Poisson型では $A$ がいわゆる剛性行列（stiffness matrix）に相当するが、本稿では「離散化行列 $A$」として統一して扱う。

4. メッシュ・要素・形状関数

4.1 要素分割と参照要素

領域 $\mathrm{\Omega }$ を要素 $e$ の集合に分割する。各要素上で参照要素 $\hat{e}$ を導入し、座標写像 $x=x(\xi )$（アイソパラメトリック写像）で物理要素へ移す。

勾配変換はヤコビアン $J=\partial x/\partial \xi$ を用いて

$${\mathrm{\nabla }}_x\varphi =J^{-T}{\mathrm{\nabla }}_{\xi }\hat{\varphi }$$

積分は

$${\int }_e(\cdot )d\mathrm{\Omega }={\int }_{\hat{e}}(\cdot )|detJ|d\xi$$

で評価する。

4.2 代表的な連続要素（Lagrange要素）

次数 $p$ を上げると同じメッシュでも近似能力が上がるが、求積・条件数・ソルバ設計に影響する。

要素	次数	自由度の増え方	特徴
$P1$（一次三角形）	1	少	実装が簡単、メッシュ品質の影響が目立つ
$P2$（二次三角形）	2	中	勾配場が改善、局所特異性にも強くなる
$Q1$（一次四角形）	1	少	構造格子に相性、歪みが大きいと精度低下
$Q2$（二次四角形）	2	中	高精度だが求積とソルバが重要になる

5. 数値積分（ガウス求積）と要素行列

一般に、要素積分は解析的に閉じないため求積点 $\{{\xi }_q,w_q\}$ を用いて

$${\int }_{e}g(x)d\mathrm{\Omega }\approx \sum _q{w}_{q}g(x({\xi }_q))|detJ({\xi }_q)|$$

で評価する。

次数 $p$ の要素では、基底関数の積や勾配の積が現れるため、必要求積次数が上がる。不十分な求積は人工的な剛性低下や不安定化を招き得るため、要素次数と求積をセットで設計する。

6. 局所寄与の散布

要素 $e$ ごとの局所行列 $A^{(e)}$ と右辺 $b^{(e)}$ を作り、大域系へ足し込む。 概念的には

$$A=\sum _e{P}_e^T{A}^{(e)}{P}_e,b=\sum _e{P}_e^T{b}^{(e)}$$

と表せる（$P_e$ は要素自由度から大域自由度への写像行列）。

大規模化では、以下が支配的になる。

• 疎行列の保存形式（CSRなど）
• アセンブリの並列化（領域分割、通信）
• 行列・ベクトル演算（特に疎行列ベクトル積）の効率

行列を明示的に保持しない行列フリー（matrix-free）手法では、$A^{(e)}$ の作用を要素ごとに計算して加算し、メモリ帯域の制約を緩和することがある。

7. 境界条件の組み込み

7.1 Dirichlet境界（本質境界条件）

代表的手法：

• 自由度固定（代入）：$u=u_D$ を満たす節点自由度を固定し、未知数から外す
• ペナルティ法：$\alpha {\int }_{{\mathrm{\Gamma }}_D}(u-u_D)vd\mathrm{\Gamma }$ を追加（$\alpha$過大で条件数悪化）
• Nitsche法：弱く課して整合性と安定性を両立させる

7.2 Neumann境界（自然境界条件）

弱形式の右辺に境界積分として自然に現れる：

$${\int }_{{\mathrm{\Gamma }}_N}gvd\mathrm{\Gamma }$$

熱流束、拡散フラックス、表面荷重などで頻出である。

8. 時間発展問題：質量行列と時間積分

半離散化（空間のみFEM）すると

$$M\dot{U}+AU=f(t)$$

の形が典型である。ここで

$$M_{ij}={\int }_{\mathrm{\Omega }}\rho {\varphi }_i{\varphi }_{j}d\mathrm{\Omega }$$

は質量行列である。

時間積分の代表：

• 後退Euler（一次、強安定）$$(M+\mathrm{\Delta }tA)U^{n+1}=M{U}^n+\mathrm{\Delta }t{f}^{n+1}$$
• Crank–Nicolson（二次、拡散で標準）$$(M+\frac{\mathrm{\Delta }t}{2}A)U^{n+1}=(M-\frac{\mathrm{\Delta }t}{2}A)U^n+\frac{\mathrm{\Delta }t}{2}(f^{n+1}+f^n)$$

拡散支配や硬い系では陰解法が選ばれやすく、その場合「毎ステップ線形系を解く」ことが計算コストを決める。

9. 線形ソルバ：Krylov法と前処理

FEMの大規模計算では、疎線形系 $AU=b$ の反復解法が中心になる。

9.1 Krylov法

• $A$ が対称正定値（SPD）に近い：CG
• 非対称・鞍点・複素：GMRES、BiCGStab など

反復回数は条件数と前処理で大きく変わる。

9.2 前処理

前処理 $P^{-1}$ を導入し、$P^{-1}AU=P^{-1}b$ を解く。代表：

• ILU/IC：中規模で使いやすいが並列性に限界
• 多重格子（MG/AMG）：楕円型に強力でスケールしやすい
• 領域分割（Schwarz、Additive/Multiplicative）：並列化と相性が良い
• ブロック前処理：連成・鞍点に対してSchur補完を使う

多相材料や層状系では係数コントラストが強く条件数が悪化しやすい。前処理は物理構造（スケール、連成、拘束）に沿わせると効くことが多い。

10. 非線形問題：Newton法とNewton–Krylov

非線形残差 $R(U)=0$ を解く場合、Newton法は

$$J(U^k)\delta U=-R(U^k),U^{k+1}=U^k+\delta U$$

である。ここで $J=\partial R/\partial U$ はヤコビアンである。

重要点：

• 各Newton反復で大規模線形系が現れる
• 線形ソルバの許容誤差がNewton収束に影響する（内側反復と外側反復の調整）
• ヤコビアンを組まずに作用だけ使う行列フリーNewton–Krylovも有力である

弾塑性、相場（Allen–Cahn、Cahn–Hilliard）、磁気弾性のような連成非線形で頻出である。

11. 混合要素・鞍点問題と安定性（inf-sup）

速度–圧力、変位–圧力、電磁場の拘束などでは

$$\left(\begin{array}{cc}A& B^T\\ B& -C\end{array}\right)\left(\begin{array}{c}x\\ y\end{array}\right)=\left(\begin{array}{c}f\\ g\end{array}\right)$$

の鞍点構造になりやすい。不適切な離散化はスプリアスモードやロッキングを生む。

代表的対策：

• 安定な要素対を選ぶ（例：Taylor–Hood要素など）
• 安定化項を導入する（PSPG、SUPGなど）
• ブロック前処理でSchur補完を近似して解く

12. 電磁場における要素選択：$H(\mathrm{curl})$ と $H(\mathrm{div})$

Maxwell方程式では回転（curl）や発散（div）が本質である。 スカラーのLagrange要素だけでベクトル場を近似すると、非物理なモードが出ることがある。

代表要素：

• Nédélec要素：$H(\mathrm{curl})$ に整合（渦の表現）
• Raviart–Thomas要素：$H(\mathrm{div})$ に整合（フラックス連続性）

磁束保存 $\mathrm{\nabla }\cdot B=0$ の数値的破れは、物理的にも計算的にも致命的になり得るため、関数空間と要素選択を方程式の構造に合わせるのが要点である。

13. 精度・収束・検証

13.1 誤差の基本像

一般に、要素サイズ $h$ と次数 $p$ により誤差が支配される（$h$-FEM、$p$-FEM、$hp$-FEM）。 Céaの補題により、Galerkin解は最良近似に比例する形で評価される：

$$\parallel u-u_h{\parallel }_V\le \frac{C}{\alpha }\underset{v_h\in V_h}{inf}\parallel u-v_h{\parallel }_V$$

ここで $\alpha$ は強圧性、$C$ は連続性に関わる定数である。

13.2 検証メニュー

• 格子収束：$h$ を半分にして誤差が理論次数で減るか
• 製造解法（MMS）：既知の解析関数 $u$ を代入して $f$ を作り、実装全体を検証
• 保存則／散逸則：質量・エネルギー・磁束などが期待通りか
• 境界条件テスト：Dirichlet/Neumann/周期の実装差が理解通りか

14. 計算量とボトルネック

大規模化で効く観点：

• 自由度数 $N$ に対する計算量の見積もり（反復回数と1反復コスト）
• メモリ帯域の制約（疎行列・ベクトル演算）
• 前処理の並列スケーラビリティ（AMG、領域分割）
• アセンブリと通信の重なり（overlap）
• 行列フリー化と高次要素（演算強度を上げて帯域制約を緩和）

「良い離散化」だけでは速くならず、「良い前処理」だけでも精度が保証されない。FEMでは両者を一体で設計する必要がある。

まとめ

• FEMは、弱形式 $a(u,v)=\ell (v)$ を有限次元空間で近似して疎連立方程式へ落とし込む数値解法である。離散化行列 $A$ は双一次形式の選択と要素設計により決まり、複雑形状・異方性・連成を同一枠組みで扱える。
• 大規模計算では、疎線形系（および非線形問題の反復内に現れる線形系）の解法が支配的である。Krylov法と前処理（多重格子、領域分割、ブロック前処理）を物理構造に沿って設計し、検証（格子収束・MMS・保存則）で数値品質を担保することが要点である。