はじめに
今日は初めて有限要素法について勉強したことを書きます。新しいことを学ぶのは量子計算の話以来です。「勉強しながらまとめていこう」とブログを初めた最初のきっかけを思い出し、とてもワクワクしております。当然のことながらこの話題が今後も続くかは謎ですし、この記事の内容も正しいのかも分かりませんし、そもそもCAE使ったこともないです。
しかし噂によると有限要素法できると1000万円くらいもらえるらしいです。燃えるぜー！皆頑張ってね（不正確：要確認）。
マジかよ！！！今から有限要素法の勉強するわ！！！！（遅い） https://t.co/IoiuzPKUvj
— HELLO CYBERNETICS (@ML_deep) 2019年2月15日
まあ、勉強のきっかけは何でも良いのですが、やはり科学技術、数学として改めて似たような考え方が至るところで出てくるなと感じました。
有限要素法の下地をひとっ走り
まず有限要素法とはなんぞやというところまで超特急で行きましょう。 行けなかったのでガラーキン法までｗ
基本的には微分方程式の数値解法であるのですが、何やら弱形式だとかガラーキン法だとか、「ムムっ？？？」という言葉と共に、「これを解けばええで」的な話で終わってしまうようなケースも見られます。「ああ、解けば良いのね、確かに解けば出るらしい、しかし心が分からん。」というのでは全く面白くありません。
微分方程式を解きたいときの気持ちになって、順々に方法を探り、そして有限要素法と呼ばれる状態まで至りましょう。。行けなかったのでガラーキン法まで。
考える微分方程式
まず2階の微分方程式を考えます。「なんで2階の微分方程式なんだ、3階の微分方程式じゃないのか！？というか一般にN階の微分方程式考えとけや！」と思うかもしれませんが、2階の微分方程式と言えばなぜか物理の問題の至るところで現れます。きっとコイツをカバーしとけば色々な問題をよく近似できるのでしょう（要確認。んー、例えばラプラシアンだとかも二階微分になりますが、なんだろう3階微分以降を特別に固有名詞を付けて扱うようなケースを私は知らない。加速度とは言うけど、それ以降の加加速度とかは普通は言わない。実世界は2階の微分でよく表せられるということだろうか？）
兎にも角にも、位置を表す変数 $x$ と位置 $x$ での変位なり熱なり、電場なり、まあ何かしらの知りたい値 $u(x)$ とします。そして外からの既知の外力なり何なりを $f(x)$ としましょう。
$$ \frac{d^2u(x)}{dx^2} - f(x) = 0 $$
ここで、$f(x)$が既知なのでいい感じの関数であれば簡単に手で解けてしまいます。 が、それは簡単な場合の話です。そういうことが困難であるということにしましょう。
解の候補を絞る
さて微分方程式を解く定石の1つとして代入法があります。線形微分方程式は適当な解の線形結合で一般解を得ることができるので $\exp(\lambda x)$ を解とおいてしまい、適当な解を見つけて後で足し合わせようという方法です。ところでもっと一般的に
$$ u^*(x) = \sum_i a_i \exp(\lambda_i x) $$
と予め置いてしまえば、後は微分方程式にコイツをブチ込んで方程式を満たすような $a_i, \lambda_i$ を求めるという方法も考えらます。この時、本当に $u^*(x)$ という表現方法で解を表しうるのかが問題になりますが、これはいわゆるフーリエ級数展開であり、いくらでも精度よく関数近似ができることが知られていますので、きっとなんとかなるのでしょう…（もとよりフーリエ級数展開も熱伝導方程式を解くために生まれたのでした。似たような方法に回路方程式を解くための演算子法があり、後にラプラス変換として理論的にまとめられるようです。そしてラプラス変換とフーリエ変換はキレイにつながっていたと、いう話は興味深いのですが関係ないので置いておきます）。
我々はここで更に一般的に、級数展開の仕方を雑に下記の用に置くことにします。
$$ u^*(x) = \sum_i a_i \phi_i(x) $$
フーリエ級数展開はそれはもう都合の良い方法ではあるのですが、今回は一先ず色々な基底関数 $\phi_i (x)$ の線形結合で解の候補 $ u^*(x)$ を表現できるのではないかということにします。基底関数がどのようなものであるかは今は定めません。
解の候補が満たすべき素朴な条件：最小二乗法
解の候補 $u^*(x)$が本当に微分方程式の解になっているのであれば、
$$ \begin{align} \frac{d^2u^*(x)}{dx^2} - f(x) &= 0 \\ \Leftrightarrow \ \ \ \frac{d^2(\sum_i a_i \phi_i(x))}{dx^2} - f(x) &= 0 \\ \end{align} $$
を満たします（何を言っている、元々この方程式を満たすような解を探しているのだから当たり前だ）。しかし、実は解の候補 $u^*(x)$ は上記の方程式を完全には満足できる見込みが無いとしましょう（原因はいろいろある。 $a_i$ の選び方が不味いのかもしれないし、仮に $a_i$ を上手に選んだとしても、準備している基底関数 $\phi_i(x)$ の選び方が不味いのであれば、そもそも表現不可能かもしれない）。もしも方程式をしっかり満たしていない、つまり左辺がちゃんと $0$ になっていないのであれば、その値、すなわち $0$ からのズレを
$$ l({\bf a}, \Phi, x) = \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) $$
と表現できます。ここで${\bf a} = \{a_1, a_2, \cdots \}$ で $\Phi = \{\phi_1, \phi_2, \cdots \}$ です（単に係数・基底関数をまとめて表記しただけである）。これを残差と呼びます。残差が $0$ ならば元の方程式そのものです。
ここで $x$ を引数に取っているのは、 $x$ の値に応じて関数のズレ具合が違うであろうと考えられるからです。例えば $x = 1$ では残差が $5$ くらいで他の $x = 2$ のところでは残差が $-2$ くらいかもしれません。これらのズレの総和を評価したい場合に、一番最初の思いつくのは下記の式ではないでしょうか。
$$ \begin{align} L({\bf a}, \Phi) &= \int_x l({\bf a}, \Phi, x)dx \\ &= \int_x \left| \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) \right|^2 dx \end{align} $$
これで残差が $x$ の至るところで $0$ ならば、上記の積分も $0$ になっており、方程式の条件をしっかり満たすということになります。今はそういえば、方程式を完全に満たしている見込みはとりあえず無いという設定なので、残差の二乗を積算した二乗誤差を最小とするような $u^*(x)=\sum_i a_i \phi_i(x)$を求めるという方針が考えられるでしょう。これを最小二乗法と言います。
解の候補が満たすべき面白い条件：重み付き残差法
ところで、例えば $x = 3$ の点でズドーンと大きくズレているのと、 いろんなところで少しずつずれているのと、気持ちの面ではどちらが良いのでしょうか…？仮にそのズレ具合の総和が$100$であること、これらは同じズレ具合と言えるでしょうか。 少なくとも二乗誤差 $L({\bf a}, \Phi)$ で評価するのであれば、同じという扱いになります。
要するに、残差の現れ方にもしっかり着目したほうが良いんじゃないのか？という話です。そこで、残差
$$ l({\bf a}, \Phi, x) = \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) $$
を下記のように書き換えます。
$$ l({\bf a}, \Phi, x, w) = \left[ \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) \right] w(x) $$
これは先程まで考えていた残差に対して、場所 $x$ 毎に異なる値を示す重み $w(x)$ を掛けた物になっています。これで、場所 $x$ 毎に残差の絶対値を調整できるため、 ズレてほしくない場所程、 $w(x)$ の値を大きく取ってあげれば良いことになります。
そして、上記の修正された残差を $x$ 全体で総和を取った式
$$ \begin{align} L({\bf a}, \Phi, w) &= \int_x l({\bf a}, \Phi, x, w)dx \\ &=\int_x \left[ \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) \right] w(x) dx \end{align} $$
を重み付き残差と呼びます（なんとなく、修正された残差を重みつき残差と呼んで、積分されたものを重み付き残差の総和と言いたくなるが、ここらへんは定義の問題であって、心がわかっていれば呼び名はどうでも良い）。
こうして、$w(x)$ で自分たちの重要度を反映した重みを考慮した損失関数らしきものを作ることが出来ました。
欲張りな重み付き残差法
ところで $w(x)$ をどのように決めれば良いのかは自明ではありません。本当ならば本来の微分方程式をキッチリと満たしてあげたいところですが、残念ながら完全に一致するような形で求まらないので、ズレの評価を適切にやってあげたいという話でした。形式的に重み付き残差を定義したものの、さて適切な評価の仕方なんぞやはりわからないのです。
そこで、重み付き残差法では下記のような戦略を取ります。まず重み $w(x)$ に対しても添え字 $j$ を付与します。
$$ L({\bf a}, \Phi, w_j) =\int_x \left[ \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) \right] w_j(x) dx $$
としておいて、あとはいろんな重み $w_j(x)$ を片っ端から準備してみましょう。
$$ \begin{align} L({\bf a}, \Phi, w_1) &=\int_x \left[ \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) \right] w_1(x) dx \\ L({\bf a}, \Phi, w_2) &=\int_x \left[ \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) \right] w_2(x) dx \\ L({\bf a}, \Phi, w_3) &=\int_x \left[ \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) \right] w_3(x) dx \\ & \vdots \\ \end{align} $$
ズレ具合の評価の仕方が $w(x)$ で決まるなら、いろんな $w(x)$ を試してみちゃえという話です。そして色々な $w(x)$ で試して、その結果、どのような重みの付け方に対しても相応に残差が小さくなるのであれば、それは元の微分方程式を上手く近似できていると考えるのです。こうすることで、重み $w(x)$ の付け方に偏りがあったとしても、それを平均的にいい感じに評価できることになります。
重み付き残差の準備の仕方を決める：ガラーキン法
さて、どのような重み $w_j(x)$ を準備してあげましょうか？という話にもちろんなってくるわけですが、ここで $w_j(x) = \phi_j(x)$ としてしまいます。 $\phi_j(x)$ というのは微分方程式を満たす解 $u(x)$ を近似しようとしたときに準備した下記の近似関数
$$ u^*(x) = \sum_i a_i \phi_i(x) $$
の基底関数です。すなわち、ある1つ目の重み付き残差$L({\bf a}, \Phi, x, w_1)$ には $w_1 = \phi_1$ を採用して
$$ L({\bf a}, \Phi, \phi_1) =\int_x \left[ \frac{d^2 (\sum_i a_i \phi_i(x))}{dx^2} - f(x) \right] \phi_1(x) dx $$
とします。2つ目には$w_2 = \phi_2$、3つ目には…4つ目には…として $\Phi$ に入っている基底関数全てを使って1つずつ重み付き残差を準備していきます（なにやら「重みの準備が面倒だから、手元に既に準備した基底関数でいいや」と投げやりになってんじゃないの…？と思うところですが、まあココらへんは、まだ私もよく分からん…ｗ）。
結局、基底関数の個数を $M$ 個と決めておき $\Phi = \{\phi_1, \cdots, \phi_M \}$ としておけば、重み付き残差も $M$ 個考えることができ、
$$ \begin{align} L({\bf a}, \Phi, \phi_1) &=\int_x \left[ \frac{d^2 (\sum_i^M a_i \phi_i(x))}{dx^2} - f(x) \right] \phi_1(x) dx \\ L({\bf a}, \Phi, \phi_2) &=\int_x \left[ \frac{d^2 (\sum_i^M a_i \phi_i(x))}{dx^2} - f(x) \right] \phi_2(x) dx \\ & \vdots \\ L({\bf a}, \Phi, \phi_M) &=\int_x \left[ \frac{d^2 (\sum_i^M a_i \phi_i(x))}{dx^2} - f(x) \right] \phi_M(x) dx \end{align} $$
を同時に評価してやるのがガラーキン法だと言えます。
一先ず今回はココまでで終了！ 後は線形代数のコトバで式を整理してあげたほうがキレイではあるのですが、本題の部分は一先ず良いと思われます。有限要素法までの流れとしては、二階微分が中々に鬱陶しいので、部分積分によって1階の微分までしか出ない形式に整えます。さらに、各々の基底関数を局在した関数にして、基底関数1つ1つがある $x$ の領域だけをそれぞれ担当する形にするようです。そうすることで、上記のように重み付き残差を複数考えた意義も出てきます（各々の重み付き残差が、基底関数が値を持つ $x$ の領域の近似精度を担当することとなる）。
