この記事を読んで何ができるようになるのか
この記事は、続編の実装編の理論的背景を理解するものです。有限要素法で問題を解くにあたり、どんな作業が必要なのかを解説します。
実装編では、2次元のポアソン方程式を有限要素法で解くためのシンプルなPythonプログラムを紹介します。
続編
https://qiita.com/fried_nasubi/items/b1da34aa4f4bb5efa7cf
解説の方針
この記事では、二次元ポアソン方程式に対する有限要素法（FEM）をフルスクラッチで実装するための理論の理解を目標とします。
有限要素法には多くの数式や概念が登場しますが、本ノートでは理論の厳密さよりも「実際に何を計算するのか」を重視して説明します。
最終的な目標は、三角形メッシュの頂点座標から大域剛性行列を構築し、連立一次方程式を解いて近似解を得ることです。
本ノートでは次の流れで説明を進めます。
- ポアソン方程式から弱形式を導出する
- 計算領域を三角形へ分割する
- 基底関数を導入し、解を基底関数の線形結合として近似する
- 試験関数を選択し、偏微分方程式を連立一次方程式へ変換する
- 各三角形について局所剛性行列を計算する
- 局所剛性行列を組み合わせて大域剛性行列を構築する
- 境界条件を適用し、解くべき連立一次方程式を得る
- 連立一次方程式を解き、各節点における近似解を求める
特に重要なのは、有限要素法が「偏微分方程式を解く手法」であると同時に、「巨大な連立一次方程式を構築する手法」でもあるという点です。
本ノートでは、この連立一次方程式がどのように作られるのかを追いながら説明していきます。
弱形式の導入
対象とする方程式は、二次元ポアソン方程式
-\Delta u=f
です。
有限要素法では、この方程式をそのまま扱うのではなく、まず弱形式へ変換します。
両辺に任意の関数 $v$ を掛けて領域全体で積分すると
-\int_\Omega dS\,v\Delta u
=
\int_\Omega dS\,fv
となります。この関数$v$は試験関数といいます。
ここで Green の第一恒等式を用いると
\int_\Omega dS\,(\nabla u)\cdot(\nabla v)
=
\int_\Omega dS\,fv
+
\int_\Gamma dl\,\mathbf{n}\cdot(v\nabla u)
を得ます。
これがポアソン方程式の弱形式です。
有限要素法では、この式を出発点として計算を進めます。
以降では、計算領域を三角形へ分割し、弱形式を連立一次方程式へ変換していきます。
計算準備その１ 計算領域の分割
弱形式の左辺第一項
\int_\Omega dS\,(\nabla u)\cdot(\nabla v)
について考えます。
有限要素法では、まず計算領域 $\Omega$ を多数の三角形に分割します。
複雑な形状の領域であっても、十分に細かい三角形の集まりとして近似することができます。そのため、有限要素法では領域を三角形へ分割して扱うことが一般的です。
領域を三角形に分割すると、領域全体での積分は各三角形上での積分の和として表現できます。
\int_\Omega dS\,(\nabla u)\cdot(\nabla v)
=
\sum_i
\int_{\Omega_i}
dS\,(\nabla u)\cdot(\nabla v)
ここで、$\Omega_i$ は第 $i$ 番目の三角形を表します。
有限要素法では、この各三角形上の積分を計算し、それらを組み合わせることで最終的な連立方程式を構成していきます。
計算準備その２ 基底関数の導入
ポイント：基底関数は、微分方程式を離散化するための手段
有限要素法では、各三角形の各頂点に対応し、基底関数と呼ばれる関数を導入します。
第 $i$ 番目の三角形の3つの頂点座標を
(x_{i,1},y_{i,1}),
\quad
(x_{i,2},y_{i,2}),
\quad
(x_{i,3},y_{i,3})
とします。
また、次の量を定義します。
2A_i=
\det
\begin{pmatrix}
x_{i,2}-x_{i,1} & y_{i,2}-y_{i,1}\\
x_{i,3}-x_{i,1} & y_{i,3}-y_{i,1}
\end{pmatrix}
頂点の番号の付け方によっては $A_i$ は負になることがあります。しかし、その絶対値 $|A_i|$ は三角形の面積に等しくなります。
第 $i$ 番目の三角形の第 $j$ 頂点に対応する基底関数を $\phi_{i,j}$ と表現すると、
\phi_{i,1}=
\left\{
\begin{array}{ll}
\displaystyle
\frac{1}{2A_i}
\left\{
(y_{i,2}-y_{i,3})x
+
(x_{i,3}-x_{i,2})y
+
(x_{i,2}y_{i,3}-x_{i,3}y_{i,2})
\right\}
&
\Omega_i \text{内}
\\[6pt]
0
&
\Omega_i \text{外}
\end{array}
\right.
\phi_{i,2}=
\left\{
\begin{array}{ll}
\displaystyle
\frac{1}{2A_i}
\left\{
(y_{i,3}-y_{i,1})x
+
(x_{i,1}-x_{i,3})y
+
(x_{i,3}y_{i,1}-x_{i,1}y_{i,3})
\right\}
&
\Omega_i \text{内}
\\[6pt]
0
&
\Omega_i \text{外}
\end{array}
\right.
\phi_{i,3}=
\left\{
\begin{array}{ll}
\displaystyle
\frac{1}{2A_i}
\left\{
(y_{i,1}-y_{i,2})x
+
(x_{i,2}-x_{i,1})y
+
(x_{i,1}y_{i,2}-x_{i,2}y_{i,1})
\right\}
&
\Omega_i \text{内}
\\[6pt]
0
&
\Omega_i \text{外}
\end{array}
\right.
となります。
これらの関数には次の特徴があります。
- 三角形内部では一次関数である
- 三角形外部では常に0である
- 対応する頂点で1をとる
- 残りの2頂点では0をとる
例えば $\phi_{i,1}$ については、
\phi_{i,1}(x_{i,1},y_{i,1})=1
\phi_{i,1}(x_{i,2},y_{i,2})=0
\phi_{i,1}(x_{i,3},y_{i,3})=0
が成り立ちます。
同様に、$\phi_{i,2}$ は第2頂点で1をとり、$\phi_{i,3}$ は第3頂点で1をとります。
このような性質を持つため、三角形内部の値を頂点上の値から線形補間することができます。
記法の簡略化
今後の計算を簡潔にするため、
b_{i,1}=y_{i,2}-y_{i,3}
c_{i,1}=x_{i,3}-x_{i,2}
b_{i,2}=y_{i,3}-y_{i,1}
c_{i,2}=x_{i,1}-x_{i,3}
b_{i,3}=y_{i,1}-y_{i,2}
c_{i,3}=x_{i,2}-x_{i,1}
と定義します。
すると、例えば $\phi_{i,1}$ は
\phi_{i,1}
=
\frac{1}{2A_i}
\left\{
b_{i,1}x
+
c_{i,1}y
+
(x_{i,2}y_{i,3}-x_{i,3}y_{i,2})
\right\}
と書くことができます。
基底関数による近似
コンピュータで微分方程式を解く際は、何らかの形で関数を離散化する必要があります。
例えば1次元では関数の離散的な近似の1つに「折れ線グラフ化」が考えられます。つまり、本来は滑らかな曲線であるものを短い直線で近似してしまうことです。
さきほど導入した基底関数で二次元の関数近似することは、ちょうどこの概念を2次元に拡張したものです。滑らかな局面を細かい三角形タイルのつなぎ合わせで表現したようなものです。（この三角形タイルを$xy$平面に投影したものが$\Omega_i$です）
関数 $u(x,y)$ は、基底関数$\phi_{ij}$と、各グリッド上に定義した値$u_{i,j}$を用いて以下のように書けます。
\sum_i \sum_{j=1}^{3}
u_{i,j}\phi_{i,j}(x,y) \simeq u(x,y)
分割領域内での計算
弱形式の左辺第一項は
\int_\Omega dS\,(\nabla u)\cdot(\nabla v)
=
\sum_i
\int_{\Omega_i}
dS\,(\nabla u)\cdot(\nabla v)
という形で、個々の三角形において
\int_{\Omega_i}
dS\,(\nabla u)\cdot(\nabla v)
を計算し、足し合わせたものであるといえます。
前節で導入した基底関数により、これを簡単に計算することができるようになりました。
試験関数$v$には任意の関数を代入することができるので、ここで、例えば
v(x,y) = \phi_{k,j}
を代入するとどうなるかを考えてみます。$\phi_{k,j}$は第$k$番目の三角形内部以外では0となるので、
\begin{align}
\int_\Omega dS\,(\nabla u)\cdot(\nabla v)
&=
\sum_{i}\int_{\Omega_i}
dS\,(\nabla u)\cdot(\nabla \phi_{k,j}) \\
&=
\int_{\Omega_k}
dS\,(\nabla u)\cdot(\nabla \phi_{k,j})
\end{align}
となります。
領域$\Omega_k$内部で非ゼロとなる基底関数は${\phi_{k,1}, \phi_{k2}, \phi_{k3} }$しかありません。したがって、先ほどの$u(x,y)$の近似式を代入すると次のようになります。
\int_{\Omega_k}
dS\,(\nabla u)\cdot(\nabla v)
=
\sum_{l=1}^{3}\int_{\Omega_k}
dS\,(\nabla u_{k,l}\phi_{k,l})\cdot(\nabla \phi_{k,j})
$u_{k,l}$は定数のため、積分の外に出すことができます。よって、
\sum_{l}u_{k.l}\int_{\Omega_k}
dS\,(\nabla \phi_{k,l})\cdot(\nabla \phi_{k,j}) =
C_{k,l,1}\,u_{k,1} +
C_{k,l,2}\,u_{k,2} +
C_{k,l,3}\,u_{k,3}
のように${u_{k,1}, u_{k2}, u_{k3} }$の1次関数の形の式を得られます。
ただし
C_{k,l,j} = \int_{\Omega_k}
dS\,(\nabla \phi_{k,l})\cdot(\nabla \phi_{k,j})
です。
ここで弱形式の全体を思い出してみます。
\int_\Omega dS\,(\nabla u)\cdot(\nabla v) = \int_\Omega dS\, fv
\,+\int_\Gamma dl \, \mathbf{n}\cdot(v\nabla u)
左辺はついさっき計算しました。
右辺第一項については$v = \phi_{k,j}$とした影響で、こちらも$\Omega_k$内部の積分だけが残ります。
例えば$f$が定数の場合は簡単に計算できます。$\phi_{k,j}$は三角錐の形をしているので、
\int_{\Omega_k} dS \, \phi_{k,j}
= |A_k|/3 \:\:\:\:\:(|A_k|は \Omega_kの面積 )
となります。
$f$が定数でない場合は何かしら近似をすることになりますが、いずれにしても
\int_\Omega dS\, fv = \int_{\Omega_k} dS\, f \phi_{k,j}
=f_{k,j}
といった「第$k$番目の三角形の第$j$頂点」と結びついた定数になります。
右辺第二項は
\int_\Gamma dl \, \mathbf{n}\cdot(\phi_{k,j}\nabla u)
という式で、計算領域の境界$\Gamma$上での積分です。よって第$k$三角形の第1頂点が$\Gamma$上にない場合は、この項は0になります。
第k三角形の第j頂点がΓ上にある場合
$\Gamma$は境界なので、何かしらの境界条件を与えています。境界条件の与え方によって処理が異なります。
①ノイマン境界条件を課している場合
ノイマン境界条件は、$\Gamma$上での
\mathbf{n}\cdot \nabla u
の値を数値として与えるものなので、この値は既知です。ただし、この項は線積分になっているため、境界条件は節点上だけでなく節点間にまで補間する必要があります。しかしあまり難しい計算ではないです。具体的な計算は実装編でお見せします）
②ディリクレ条件を課している場合
あとで説明しますが、この項の影響はなくなるので未知のままでいいです。
このように、$\phi_{k,j}$の特徴（どんな境界条件が課されているのか、いないのか）によって計算方法が異なりますが、とりあえずこの値を
$
g_{k,j}
$
とでも置いておきましょう。
以上のとおり、弱形式の試験関数として$\phi_{k,j}$を代入すると、
C_{k,j,1}u_{k,1} +
C_{k,j,2}u_{k,2} +
C_{k,j,3}u_{k,3}
= f_{k,j} +g_{k,j}
という式が得られました。$v$として別の基底関数を選べば、異なる$u_{i,j}$を含む式を得られます。このようにして連立方程式を立てることで未知数を解きます。
ここで有限要素法の本質が見えてきます。
もともとの問題は、
-\Delta u=f
という偏微分方程式の解を求める問題でした。
しかし有限要素法では、解そのものを求める代わりに、
「各節点における解の値」
を未知数とする連立一次方程式へ問題を変換します。
つまり有限要素法とは、
偏微分方程式を連立一次方程式へ変換し、その解として節点上の近似解を求める手法
と考えることができます。
剛性行列
前節のとおり、未知数$u_{i,j}$に関する連立方程式を立てる見通しが立ちました。しかし、$u_{i,j}$は各三角形の頂点に対して定義した値なので、本当に求めたい節点上の値としては重複が発生しています。
ここで、計算内容をイメージしやすくするため、計算領域と、その内部の分割のイメージ図を掲載します。
図には、$\Omega$を分割する小三角形とその通し番号、節点とその通し番号を記入しています。さらに各節点を「第$i$番目の三角形の第$j$頂点」としてみたときの番号を「$i$-$j$」のように記入しています。
この図では「第1三角形の第3頂点」、「第2三角形の第3頂点」と「第3三角形の第1頂点」などが重複しています。
三角形の頂点は重複していても、節点そのものは重複していません。未知数の定義位置を「第$i$番目の三角形の第$j$頂点」ではなく「第$k$節点」としてしまえば未知数に重複はなくなります。よって、以降は未知数のインデックスを$u_k$のように表現することにします。
例えば、先の例で出したものは「第4節点」と表現できるため、
u_4 = u_{1,3} = u_{2,3} = u_{3,1}
だといえます。
問題を解くうえでの次の目標は、重複のない未知数${ u_k }$を用いた連立方程式を作ることです。
試験関数として$\phi_{1,1}, \phi_{1,2}, \phi_{1,3}$を選ぶと、順に
\begin{align}
u_{1,1} + 2u_{1,2} + u_{1,3} &= 1 \\
-u_{1,1} + 3u_{1,2} + 2u_{1,3} &= 1 \\
-2u_{1,1} + u_{1,2} + 2u_{1,3} &= 2
\end{align}
試験関数として$\phi_{2,1}, \phi_{2,2}, \phi_{2,3}$を選ぶと、順に
\begin{align}
2u_{2,1} + -u_{2,2} + 3u_{2,3} &= 1 \\
2u_{2,1} + u_{2,2} + u_{2,3} &= 2 \\
u_{2,1} + 3u_{2,2} + 2u_{2,3} &= 1
\end{align}
という式を得たとします。次に「三角形の頂点に結び付けた変数」から「節点に結び付けた変数」への変換を行います。図より、$u_{1,3} = u_4$などと変換することで先の式は
\begin{align}
u_{1} + 2u_{2} + u_{4} &= 1 \\
-u_{1} + 3u_{2} + 2u_{4} &= 1 \\
-2u_{1} + u_{2} + 2u_{4} &= 2
\end{align}
\begin{align}
-u_{1} + 2u_{3} + 3u_{4} &= 1 \\
u_{1} + 2u_{3} + u_{4} &= 2 \\
3u_{1} + u_{3} + 2u_{4} &= 1
\end{align}
のようになります。
有限要素法の計算過程では
「ある三角形の頂点に結び付いた基底関数$\phi_1,\phi_2,\phi_3$を用いて導いた3つの連立方程式」を、「すべて足し合わせる」という操作をします。
ここでは、第１三角形と、第2三角形の連立方程式を組み合わせる過程を紹介します。
まず、必要に応じ、連立方程式の順番を入れ替えます。
一組の連立方程式を解くのに式の順番は関係ありませんが、有限要素法では、個別の三角形ごとに用意した複数の連立方程式を組み合わせて解くため、あとの処理の都合上、このような処理をしておきます。
さて、第１三角形に関する連立方程式は、順に第１、第２、第４節点に関する基底関数を試験関数として使ったものです。これは（$1\rightarrow 2\rightarrow 4$）と順になっているため並べ替えは不要です。
一方、第２三角形に関する連立方程式は、順に第３、第１、第４節点に関する式です。こちらは式が節点の番号順になるように並び替えます。
\begin{align}
u_{1} + 2u_{3} + u_{4} &= 2 \\
-u_{1} + 2u_{3} + 3u_{4} &= 1 \\
3u_{1} + u_{3} + 2u_{4} &= 1
\end{align}
次に、それぞれを変数$u_1, u_2, u_3, u_4$を変数として持つ方程式に拡張します。
\begin{pmatrix}
1 &2 &0 &1 \\
-1 &3 &0 &2 \\
0 &0 &0 &0 \\
-2 &1 &0 &2
\end{pmatrix}
\begin{pmatrix}
u_1 \\
u_2 \\
u_3 \\
u_4
\end{pmatrix}
=
\begin{pmatrix}
1 \\
1 \\
0 \\
2
\end{pmatrix}
\begin{pmatrix}
1 &0 &2 &1 \\
0 &0 &0 &0 \\
-1 &0 &2 &3 \\
3 &0 &1 &2
\end{pmatrix}
\begin{pmatrix}
u_1 \\
u_2 \\
u_3 \\
u_4
\end{pmatrix}
=
\begin{pmatrix}
2 \\
0 \\
1 \\
1
\end{pmatrix}
この式の両辺を足し合わせると、
\begin{pmatrix}
2 &2 &2 &2 \\
-1 &3 &0 &2 \\
-1 &0 &2 &3 \\
1 &1 &1 &4
\end{pmatrix}
\begin{pmatrix}
u_1 \\
u_2 \\
u_3 \\
u_4
\end{pmatrix}
=
\begin{pmatrix}
3 \\
1 \\
1 \\
3
\end{pmatrix}
を得ます。この「両辺同士を足し合わせること」が、弱形式の
\int_\Omega dS\,(\nabla u)\cdot(\nabla v)
=
\sum_{i}\int_{\Omega_i}
dS\,(\nabla u)\cdot(\nabla \phi_{k,1}) \\
の$\sum_i$に対応します。
今回は例として二つの三角形に関する「足し合わせ」を示しましたが、実際の計算ではすべての三角形に関して足し合わせを行います。
このように
- 個別の三角形において、3変数の連立方程式を立てる
- 関係ない変数も含む式に拡張する
- それらを足し合わせて、全変数を含む連立方程式を立てる
という手順を取ります。
1つの三角形に対して計算した、3つの変数のみを含む行列形式
\begin{pmatrix}
1 &2 &1 \\
-1 &3 &2 \\
-2 &1 &2
\end{pmatrix}
\begin{pmatrix}
u_1 \\
u_2 \\
u_4
\end{pmatrix}
=
\begin{pmatrix}
1 \\
1 \\
2
\end{pmatrix}
の左辺に出てくる行列を局所剛性行列と呼びます。
そして、すべての三角形について足し合わせた式、
K\mathbf{u} = \mathbf{f}, \:\:\:\:\:\:\:\:
\mathbf{u} =
\begin{pmatrix}
u_1 \\
u_2 \\
\vdots \\
u_n
\end{pmatrix}
に出てくる行列$K$を大域剛性行列といいます。
注：試験関数について
この記事では、基底関数として三角形の頂点に対応した$\phi_{i,j}$を導入し、一貫してこの記法を用います。これは実装との対応を分かりやすくするためです。一方、文献によっては、節点ごとに重複のない基底関数$\phi_i$を定義して説明する流儀もあります。
ディリクレ条件の適用
節点上の$u$の値はすべてが未知ではなく、ディリクレ条件を与えた境界点では既知になっています。したがって、
\mathbf{u} =
\begin{pmatrix}
u_1 \\
c \\
u_3 \\
\vdots \\
u_n
\end{pmatrix}
のように既知節点が入っていて、行列形式はこのような形になります。
\begin{pmatrix}
K_{1,1} & K_{1,2} & \cdots &K_{1,n}\\
K_{2,1} & K_{2,2} & \cdots & K_{2,n}\\
\vdots & \vdots & \ddots & \vdots\\
K_{n,1} & K_{n,2} & \cdots & K_{n,n}
\end{pmatrix}
\begin{pmatrix}
u_1 \\
c \\
u_3\\
\vdots \\
u_n
\end{pmatrix}
=
\begin{pmatrix}
f_1 \\
f_2 \\
f_3 \\
\vdots \\
f_n
\end{pmatrix}
これは連立方程式の問題なので、既知変数があればそれだけ式の数（行列の行）を減らせます。
この例の場合、第2列に既知数が入っているので、両辺から第2行を消し、左辺の行列の第2列に$c$を掛けて右辺に移項します。
（このような処理でいい理由はあとで簡単な例で確認しましょう）
よって、
\begin{pmatrix}
K_{1,1} & K_{1,2} & \cdots &K_{1,n}\\
K_{3,1} & K_{3,2} & \cdots & K_{3,n}\\
K_{4,1} & K_{4,2} & \cdots & K_{4,n}\\
\vdots & \vdots & \ddots & \vdots\\
K_{n,1} & K_{n,2} & \cdots & K_{n,n}
\end{pmatrix}
\begin{pmatrix}
u_1 \\
u_3 \\
u_4\\
\vdots \\
u_n
\end{pmatrix}
=
\begin{pmatrix}
f_1 \\
f_3 \\
f_4 \\
\vdots \\
f_n
\end{pmatrix}
-c
\begin{pmatrix}
K_{2,1} \\
K_{2,3} \\
K_{2,4} \\
\vdots \\
K_{2,n}
\end{pmatrix}
のようになります。
前後の式を見比べてみると$f_2$という値が消えています。
ところで、この消去した2行目の式が何だったか考えてみます。
で説明したように、大域剛性行列の第２行には、第２節点に関する基底関数を試験関数として作った式が足しあわされています。
図から第２節点を含む三角形を確認すると、$f_2$には
\int_\Gamma dl \, \mathbf{n}\cdot(\phi_{1,2}\nabla u)
\int_\Gamma dl \, \mathbf{n}\cdot(\phi_{3,2}\nabla u)
\int_\Gamma dl \, \mathbf{n}\cdot(\phi_{4,1}\nabla u)
が含まれていることがわかります。これは「分割領域内での計算」の節の「②ディリクレ条件を課している場合」で未知のまま放置していた項です。式の並べ替えをしておいたことで、これらの未知項を「あとで削除する行」にすべて押し込めたわけです。
この作業を$\mathbf{u}$からすべての既知数が取り除かれるまで繰り返します。
既知変数削除の考え方
式変形を検算できるように$4\times4$の行列に関する式で試してみましょう。
\begin{pmatrix}
1 & 2 &3 &4\\
3 & 5 & 1 &7 \\
2 & 7 &3 & 3\\
1 & 2 &4 & 2
\end{pmatrix}
\begin{pmatrix}
u_0 \\
c \\
u_2\\
u_3
\end{pmatrix}
=
\begin{pmatrix}
1 \\
2 \\
3 \\
4
\end{pmatrix}
この式から既知変数$c$を除くには以下のように式変形します。
既知数が入っているのは第「２」行のため、
両辺から第2行をを消し、左辺の第「２」列に$c$を掛け右辺から引きます。
\begin{pmatrix}
1 & 3 &4\\
2 & 3 & 3\\
1 & 4 & 2
\end{pmatrix}
\begin{pmatrix}
u_0 \\
u_2\\
u_3
\end{pmatrix}
=
\begin{pmatrix}
1 \\
3 \\
4
\end{pmatrix}
-c
\begin{pmatrix}
2 \\
7 \\
2
\end{pmatrix}
解くべき式が完成した
なんだかんだでいろいろな作業をしてきましたが、その結果
-\Delta u = f
というポアソン方程式を、
K\mathbf{u} = \mathbf{f}, \:\:\:
\mathbf{u}=
\begin{pmatrix}
u_1 \\
u_3 \\
u_4\\
\vdots \\
u_n
\end{pmatrix}
という形に変形することができました。
この形式さえ作れてしまえば、本当に連立方程式を解くだけなので、例えばPythonなら
u = np.linalg.solve(K, F)
としてしまえば節点上の$u$が求まります。
以上で、有限要素法によってポアソン方程式を連立一次方程式へ変換する流れを確認しました。
重要なのは、偏微分方程式を直接解くのではなく、節点上での解の値を未知数とする行列方程式へ帰着できることです。また、その行列は各三角形ごとの局所的な計算を組み合わせることで構築できることも分かりました。
次は実装編として、ここで導いた式を実際のプログラムへ落とし込み、局所剛性行列や大域剛性行列をどのように計算するのかを見ていきます。
次回記事
https://qiita.com/fried_nasubi/items/b1da34aa4f4bb5efa7cf