はじめに
前回・前々回と有限差分法(FDM)の「基礎編」「実践編」と称して書いていった。
次は有限要素法(Finite Element Method, FEM)について執筆していく。
今回も基礎編・実践編の２つの記事に分けて説明していく。
(雑になりすぎたので書き直す可能性大)
復習
簡単に復習しておこう。
電場を解析するためにはポアソン方程式:
$$
\nabla \cdot \nabla u = f
$$
を解けばいいのだった。
さらに、これを一意の解とするには境界条件が必要である。
境界条件にはディリクレ境界条件とノイマン境界条件の２つが有名である。
ディリクレ境界条件は境界上の$u$の値を、ノイマン境界条件は境界上の$u$の垂直方向の微分を指定するものだ。
このポアソン方程式は一部の特殊な時を除き、殆どの場合数学的に解くことが出来ない。
そのため、近似解を求める手法が重要になる。その中の一つが前回・前々回と説明した有限差分法(FDM)であった。
今回はより一般に利用されている有限要素法について説明する。
有限要素法(FEM)
FDMと比べて、FEMのアルゴリズムは結構複雑だ。
あまり一般化しすぎた話をするとイメージが沸かないと思うので、今回の説明は(二次元)ポアソン方程式に主眼をおいて話していこう。
とりあえずイメージを掴む、というのは何に置いても大事なことだ。
と、その前にFEMの利点について話しておこう。
FEMはFDMと比べて複雑な形状にも容易に対応できるのである。
FDMでは前回のような台形ですら結構境界の取り扱いが面倒だった(実は記事にしてない工夫も結構あったりする……)。これは、$x$軸, $y$軸と平行に格子点(ノード)を置く必要があるからだ(平行じゃないノードだと、前回のノイマン境界条件で説明したような処理が全領域に渡って必要になる)。
一方で、FEMはそれぞれを有限要素に分割するのだが、この分割方法にFDMのような面倒なことを考えず、適当に格子点をおいても問題ない。しかも、境界付近の取り扱いも容易である。
格子分割というものをプログラムで実装するのは思った数十倍も面倒である。というのも、アルゴリズム自体の面倒さだけじゃなく、丸め誤差などがかなり大きい影響を与えてしまう。
そこがある程度自由にできるなら、それに越したことはないだろう。
前回使ったような規則通りの格子を構造格子と呼び、一方で不規則な格子を非構造格子と呼ぶということもここに記しておこう。
格子分割
上記で説明したとおり、FEMは有限要素に分割する。
この有限要素とは、二次元なら微小面積、三次元なら微小体積のことを言う。下図に簡単な例を作ってみた。
図1格子分割イメージ
一般に、有限要素法では上図のように三角形に分割することが多い。
一応四角形などの要素に分割することも可能で、そのほうが精度(の次数)が高くなるが最初は三角形の分割について解説していこう。理解すれば、四角形などで精度が高くなるのも当然と思えるようになると思う。
なお、この記事ではこの要素をセルと呼ぶことにする。
この記事だけでそう呼ぶっていうだけで、一般的な言葉というわけではないので(普通はElement?)、他のところで注釈なしで使うのはやめておこう。
導出
ポテンシャルエネルギー最小の原理からの導出
殆どのFEMの解説書では、重み付き残差法を基礎において解説している。
ただ、今回は『数値電界計算の基礎と応用』に倣って、ポテンシャルエネルギー最小の原理から説明してから、最後に重み付き残差法からの解釈を書くことにする。
さらに、最初はポアソン方程式$\nabla \cdot \nabla u = f$ではなく、右辺を0にしたラプラス方程式$\nabla \cdot \nabla u = 0$で説明しよう。
ポテンシャルエネルギー最小の原理は変分法というものから導出される。
ちょっと変分法について書いていくと長くなってしまってダレてしまいかねないので、本当はとても大事な理論なのだが、今回は省く。
かなり簡単で曖昧な説明をすると、微分が変数の変化量を求めるものなのに対して、変分は「関数の変化量」を求めるものである。
とりあえず、次のロジックが成り立つということだけ説明なしに納得していってほしい。
- ポテンシャルエネルギー最小の原理
\begin{align}
\frac{\partial}{\partial x}\frac{\partial h}{\partial \frac{\partial u}{\partial x}} +
\frac{\partial}{\partial y}\frac{\partial h}{\partial \frac{\partial u}{\partial y}} -
\frac{\partial h}{\partial u} = 0 \tag{1}
\end{align}
をときたいときは、
$$
X(u) = \iint h(x,y,\frac{\partial u}{\partial x},\frac{\partial u}{\partial y})dxdy \tag{2}
$$
が最小値となる$u$を求めれば良い。
ここで、
$$
h =
\frac{1}{2}\left(\frac{\partial u}{\partial x}\right)^2 +
\frac{1}{2}\left( \frac{\partial u}{\partial y} \right)^2
$$
とすると、式(1)は、次のようにラプラス方程式になってくれる。
$$
\frac{\partial^2 u}{\partial x^2} +
\frac{\partial^2 u}{\partial y^2} = 0
$$
式(1)が偏微分の偏微分になっているのでちょっと複雑にみえるが、普通に$h' = 1/2 x^2+1/2y^2$のような式を$x,y,u$で偏微分してるのと同じだ。
同様に式(2)に代入してみると、次のようになる。
$$
X(u) = \iint \frac{1}{2}\left(\frac{\partial u}{\partial x}\right)^2 +
\frac{1}{2}\left( \frac{\partial u}{\partial y} \right)^2 dxdy
$$
以上の議論をまとめよう。
ラプラス方程式
$$
\frac{\partial^2 u}{\partial x^2} +
\frac{\partial^2 u}{\partial y^2} = 0
$$
を解くには、
$$
X(u) = \iint \frac{1}{2}\left(\frac{\partial u}{\partial x}\right)^2 +
\frac{1}{2}\left( \frac{\partial u}{\partial y} \right)^2 dxdy
$$
が最小値となる$u$を求めればいい。
……え？ 逆に複雑になってる気がするって？ 筆者もそう思うが、ここは目を瞑って進めていこう。
形状関数
ちょっと話が切り替わる。
FDMが各点の値を近似するのに対して、FEMは式を近似するという違いがある。
といっても、式のすべてを近似するなんて難しそうなことは出来ないので、式の係数を変えて近似していくイメージを持ってくれればいい。
最も簡単なのは一次式だ。有限要素法でもほとんどこれが使われる。
まずは、一次元で例示してみよう。
一次元の一次式は次のように表せることは高校で習っただろう。
$$
u = \alpha_1 x + \alpha_2
$$
ここで、一次式を通る２つの点がわかれば$\alpha_1,\alpha_2$が一意に決まる、ということは当然わかると思う。
この考えを二次元に拡張しよう。
二次元での一次式は一次元と同様に次のように表せる。
$$
u = \alpha_1 x + \alpha_2 y + \alpha_3 \tag{3}
$$
当然これは、３つの点がわかれば$\alpha_1 ,\alpha_2 , \alpha_3$が一意に決まるということになる。
つまり、三角形の各頂点の値がわかれば良いのだ。
イメージ図も下に用意した。
図2：１次関数
最初の格子分割に絡めて話そう。
$u$は未知数だが、$x$, $y$については格子分割は自分で定義しているわけだし既知の値のはずだ。
三角形の頂点の座標を$(x_1, y_1), (x_2, y_2), (x_3, y_3)$とすると、次のようになる。
\begin{align}
u_1 = x_1 \alpha_1 + y_1 \alpha_2 + \alpha_3 \\
u_2 = x_2 \alpha_1 + y_2 \alpha_2 + \alpha_3 \\
u_3 = x_3 \alpha_1 + y_3 \alpha_2 + \alpha_3
\end{align}
行列に書き直すと、
\begin{align}
\left[
\begin{matrix}
u_1 \\
u_2 \\
u_3
\end{matrix}
\right] =
\left[
\begin{matrix}
x_1 & y_1 & 1 \\
x_2 & y_2 & 1 \\
x_3 & y_3 & 1
\end{matrix}
\right]
\left[
\begin{matrix}
\alpha_1 \\
\alpha_2 \\
\alpha_3
\end{matrix}
\right]
\end{align}
これによって、$\alpha_1, \alpha_2, \alpha_3$が求まるはずだ。
逆行列をとったものも書き下しておく。
\begin{align}
\left[
\begin{matrix}
\alpha_1 \\
\alpha_2 \\
\alpha_3
\end{matrix}
\right]
=
\left[
\begin{matrix}
x_1 & y_1 & 1 \\
x_2 & y_2 & 1 \\
x_3 & y_3 & 1
\end{matrix}
\right]^{-1}
\left[
\begin{matrix}
u_1 \\
u_2 \\
u_3
\end{matrix}
\right]
\end{align}
上の式を解いて、もう一度最初の式(3)に代入すると、何らかの係数$N_1, N_2, N_3$を用いて次のような一次式(関数)を作ることができる。
$$
u = N_1 u_1 + N_2 u_2 + N_3 u_3 \tag{4}
$$
この$u$を形状関数という。
$N_1, N_2, N_3$は地道に計算すればできるが、ちょっと面倒なので、次のようになるということだけ書き留めておく。
\begin{align}
N_1 &= \frac{1}{2S}(a_1x + b_1y + c_1) \\
a_1 &= y_2 - y_3 \\
b_1 &= x_3 - x_2 \\
c_1 &= x_2 y_3 - x_3 y_2 \\
\end{align}
ここで、$S$は注目している三角形の面積である。
他の$N_2, N_3$も同様な結果になる。また、これらの$N$はすべて$x, y$の関数である、ということを覚えておこう。
有限要素法はこのような形状関数を各セル内で求める手法なのである。
なんとなくイメージが湧いてきただろうか？
定式化
話を戻そう。
$$
X(u) = \iint \frac{1}{2}\left(\frac{\partial u}{\partial x}\right)^2 +
\frac{1}{2}\left( \frac{\partial u}{\partial y} \right)^2 dxdy
$$
が最小となる$u$を求めるのであった。
これはつまり、$X(u)$の$u$に関する微分が$0$となる$u$を求めればいいということだ。
……といっても、このまま$u$で微分しろと言われても困る。
とりあえず、上の形状関数を代入してみた。
$$
X(u) = \iint \frac{1}{2}\left(\frac{\partial }{\partial x}(N_1 u_1 + N_2 u_2 + N_3 u_3)\right)^2 +
\frac{1}{2}\left( \frac{\partial }{\partial y}(N_1 u_1 + N_2 u_2 + N_3 u_3) \right)^2 dxdy
$$
$N$が$x,y$の関数であることに注目すれば、
$$
X = \iint \frac{1}{8S^2}\left[\left(a_1u_1 + a_2u_2 + a_3u_3\right)^2 +
\left( b_1u_1 + b_2u_2 + b_3u_3\right)^2 \right] dxdy \tag{5}
$$
とすることができる。式がどんどん具体的になってきた。
積分領域を、注目しているセル上に渡る積分としよう。
上式は既に$x,y$によらない。よって、積分値は積分領域の面積となる。つまり、
$$
X = \frac{1}{8S}\left[ \left(a_1u_1 + a_2u_2 + a_3u_3\right)^2 +
\left( b_1u_1 + b_2u_2 + b_3u_3\right)^2 \right]
$$
となる。
$u$で微分するのは無理だが、$u_1$で微分することは可能だ。
$$
\frac{\partial X}{\partial u_1} = \frac{1}{4S}\left[ a_1\left(a_1u_1 + a_2u_2 + a_3u_3\right) +
b_1\left( b_1u_1 + b_2u_2 + b_3u_3\right) \right]
$$
すべての$u_i$で上式がゼロになってくれば、少なくとも注目するセル上の$\frac{\partial X}{\partial u} = 0$が(ある程度)成り立っていると言えそうだ。
この考えを一つのセルだけでなく、計算領域全体に広げて考えてみると、今までの議論で$X$は一つのセル上の関数であり、セルごとに違う関数のはずだ。そのため、セル番号$j$に対応させて$X_j$と名付けなおすことにする。
すると、全体のポテンシャルエネルギー$\mathcal{X}$は以下のように表せる。
$$
\mathcal{X} = \sum_j X_j \text{ (}X_j\text{はセル$j$外では0)}
$$
ここで、最小値を求める問題というのは、すべての$u_i$上で
$$
\frac{\partial \mathcal{X}}{\partial u_i} = 0 \tag{6}
$$
が成り立てば良いと言い換えることができる。
なんとなく差分法と同じように、連立方程式を建てることができそうだな～っていうのが見えてきた。
さて、上記の$u_i$は下図のように複数のセルでも使われているはずだ。
図3：一つのノード周りのセル
上図の赤点を$u_1$とすると、色のついたセルのすべてが$u_1$の点を頂点としているため、形状関数内に$u_1$が存在することになる。
逆に言うと、これ以外のセル上の形状関数には$u_1$は存在していないといえるため、$u_1$での微分値は0になる。
つまり、上図の例だと、以下の式のように$X_1, X_2, X_3, X_4, X_5$以外の項はすべて消えてしまうことになる。
$$
\text{式(6)} \iff \frac{\partial X_1}{\partial u_1} + \frac{\partial X_2}{\partial u_1} + \frac{\partial X_3}{\partial u_1} + \frac{\partial X_4}{\partial u_1} + \frac{\partial X_5}{\partial u_1} = 0
$$
これを、すべての$u_i$で立ててしまえば連立方程式が完成だ。
連立方程式の未知数は$u_i$, 式の数はノードの数となる。
境界条件
FEMの境界条件の定式化はFDMと比べてとても簡単だ。
ディリクレ境界条件
これは、FDMと同様に、$u = g$と与えてしまえばいい。
ノイマン境界条件
ノイマン境界条件もFDMと比べてとても簡単で、
$$
\frac{\partial u}{\partial n} = g
$$
を境界に接するセルの形状関数を式(4)のように、微分して導出すればいいだけだ。
このノイマン境界条件の簡単さも、FEMがよく使われている理由の一つだろう。
Galerkin法による導出
以上で、ラプラス方程式の有限要素法のプログラムを作るだけの理論は説明できた。
ただ、上記の説明は変分法から導出されるものだ。ポテンシャルエネルギー最小の原理という謎の理論が突然出てきて納得行かない人もいるだろう。
そこで、重み付き残差法から説明してみよう。普通の数値計算の解説書としてはこちらから導出するほうが一般的だ。
というのも、ポテンシャルエネルギー最小化からの手法が物理的な背景から導出されているが、Galerkin法のほうがより一般的な微分方程式に対応できるからだ。
といっても、導出法が違うだけで、結局は同じ式になる。
今まではラプラス方程式を扱っていたが、より一般的にポアソン方程式$\nabla \cdot \nabla u = f$を解いてみよう。
重み付き残差法
まず、Galerkin法の基礎概念、重み付き残差法について説明しよう。
一般的な微分方程式は以下のように表せる。
$$
Lu = f
$$
$L$は偏微分の演算子である。例えば、ポアソン方程式なら、
$$
L = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}
$$
である。ここで、移項すれば以下のようになるはずだ。
$$
Lu^* - f = 0
$$
$u$は厳密解であるということを明記するために、$u^* $と置き直した。
さて、ここで$u^* $の近似解を$u$とすると、次のようになる。
$$
Lu - f = \epsilon \tag{7}
$$
$\epsilon$は誤差である。この$\epsilon$、定数に見えるかもしれないが、関数であることに注意しよう。
この$u$を有限個の何らかの関数(厳密には線形独立な関数群)$u_i$の足し算(線形結合)で表す。
式で表したほうがわかりやすいかもしれない。
$$
u = \sum_{i=1}^N C_i u_i
$$
この$C_i$は係数である。
当然だが、$u$は近似解なので、式(7)の$\epsilon$は0にならない。
ここで式(7)に$u_i$を代入して$\epsilon$が(定数関数)0になるべく近くなるような$C_i$を選べばいいのである。
ただ、上記で説明したように、$\epsilon$は関数であるため、この係数$C_i$を選ぶというのはなんとなく難しそうな感じがする。
そこで、何らかの重み関数$\psi_j$を用意して、
$$\iint (Lu - f)\psi_j \,dxdy = 0 \tag{8}$$
とする。これで左辺は『関数の積分値』となり、関数ではなく何らかの定数になる。これで$C_i$を選択しやすくなった。
ここでいろいろな重み関数で成り立つような$C_i$を見つけられれば、$Lu - f \approx 0$になりそうだ。
この重み関数$\psi_j$の選び方で、選点法、最小二乗法、モーメント法などに分岐する。
今回説明するのはGalerkin法である。
Galerkin法
Galerkin法では、この関数$\psi_j$を上記の$u_i$それ自身とする手法だ。
その前に、イメージしやすいように$u_i$をどう置くか？と言う話を先にしておこう。
これも色々あるのだが、注目するノード上ので隆起するような一次関数としておくのが一番一般的だ。一次元の例を下に図に表してみた。なお、この$u_i$を基底関数という。
図4：一次元基底関数
二次元の場合は、下図のように$u_i$はそれぞれのノードの周りのセルを線形にせり上がったような関数となる。なんとなく形状関数に似てると感じられると思う。
図5：二次元基底関数
どうせ係数$C_i$があるのだから、頂点の値は何でもいいが、まあ1とするのが普通だろう。
上記の二次元の一次関数は、ノード周りの形状関数を参考にして、
$$
u_i = \sum_j N_{i,j} \tag{9}
$$
のように表すことができる。
形状関数の式のせり上げた$u_j$を1その他の２つの項を0と置いたものだ。
混乱する書き方になってしまったので注記しておくが、形状関数の項で説明した$u_i$は定数であり、式(9)の$u_i$は関数である。変数名を変えればよかったかもしれない……。
もう一つ、弱形式について説明する。
まず、これまでの議論から式(8)は次のように変形することができる。
$$
\iint \left[\left( \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2} \right)\sum_{j=1}^N C_{i,j}u_j - f \right]u_i \,dxdy = 0
$$
ここで、以下の公式を利用する。
$$
\iint (\nabla \cdot \nabla u )vdxdy = -\iint \nabla u \cdot \nabla vdxdy + \int v\nabla u \cdot d\boldsymbol{S}
$$
第二項は、領域の外側に向かう積分である。これは、部分積分を利用した、
$$
\int \frac{\partial^2 u}{\partial x^2}v dx = - \left[ \frac{\partial u}{\partial x} v \right] + \int \frac{\partial u}{\partial x}\frac{\partial v}{\partial x}dx
$$
からわかる。これはグリーンの恒等式と名付けられている。
積分領域を計算領域全体として、
$$
\sum_{i=1}^N \left(\iint_S
\left[C_{i,j}
\left(
\frac{\partial u_i }{\partial x}
\frac{\partial u_j }{\partial x}+
\frac{\partial u_i }{\partial y}
\frac{\partial u_j}{\partial y}
\right)
+ f u_i
\right] dxdy
- \int_\Gamma u_i\nabla u_j \cdot d\boldsymbol{S} \right)= 0$$
という式が導出できる。
基底関数を見ればわかるが、ほとんどが0であり、当然その微分値も0になる。
そのため、左辺の$u_i\nabla u_j$は境界上のノード以外は0になる。
また、$u_i$と$u_j$が近い位置にない場合(より正確には一致・もしくは隣にある場合は)0になる。
これで定式化すると、上のポテンシャルエネルギー最小の原理からの導出と全く同じ式になる。
詳細な証明はしないが、式の形が似ていることと、複数の線形関数に近似している、ということからなんとなくで納得してほしい。納得できない方は自分で調べてみてね。
まとめ
今回は有限要素法の基礎理論について説明した。
定式化まとめ
定式化の手順をまとめよう。
変分法による導出
ポテンシャルエネルギー最小の原理から、ラプラス方程式
$$\nabla \cdot \nabla u = 0$$
は、次の式が最小値となる$u$を求める、という問題と同値である。
$$X(u)=\frac{1}{2}\iint \left( \frac{\partial u}{\partial x} \right)^2 + \left( \frac{\partial u}{\partial y} \right)^2 dxdy$$
上式の最小値を求めるには、$u$で微分して値が0となればいいのだが、関数$u$で微分するのは難しいので、計算領域内の離散的な値$u_i$の微分値を求め、全てが0になるようにする。
それぞれの$u_i$は三角形セル$j$内の一次関数(形状関数)で表せば、上式に代入して微分することができる。これによって、以下のような式にできる。
\begin{align}
\sum_j \frac{\partial X_j}{\partial u_i} &= 0 \\
\frac{\partial X_j}{\partial u_i} &= \frac{1}{4S}\left[ a_1\left(a_1u_1 + a_2u_2 + a_3u_3\right) +
b_1\left( b_1u_1 + b_2u_2 + b_3u_3\right) \right]
\end{align}
すべての$u_i$でこの式を立てて、連立方程式を解けば、解が求められる。
重み付き残差法による導出
一般の微分方程式
$$
Lu^* = f
$$
に対して、厳密解$u^*$を近似解$u$に変えると、
$$
Lu - f = \epsilon
$$
となる。
この$u$は線形独立な関数群(一般には一次関数)の線形結合によって以下のように表される。
$$
u = \sum_i^N C_iu_i
$$
$\epsilon$を0に近づけるために、重み関数$\psi$を利用して、
$$
\iint \psi_j \epsilon dxdy= 0
$$
となる$C_i$を求める。
ここで、$\psi_j$を$u_j$と置くのがガラーキン法である。
あとは、代入すれば、連立方程式を建てることができる。
あとがき
証明していない公式なども使ったりして、あまりスッキリしない説明だったかもしれないが……、それはいつの日か証明する記事を書くかもしれない。
また、なるべく変数の文字数を少なくしようと頑張った結果、逆にわかりにくくなってるかもしれない。
わかりにくい表現などがあったらコメントを下さい。
次回は有限要素法の実装をしていこう。