有限要素法(FEM)の概略をメモ

重み付き残差法の大まかな流れをメモする。
ここでは
$$ \frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f=0 $$を有限要素法で解くことを考える
これは、任意の重み関数$ w $をかけて領域積分した値が0に等しいという以下の式に等価である。
$$ \int_{\Omega}w\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f \right)dS=0 $$
このまま($w$が任意)では難しいので$w=w_i$　$i=1,2\cdots ,n $について満たす解を求める
領域を少領域に分け、それぞれの領域で（例えば）1次式で解を
$$ u\approx \alpha_1+\alpha_2x+\alpha_3y $$と近似する。ただし、領域ごとに係数は異なる。このままでは扱いにくいので領域に接する節点における物理量で補間して解を表すことを考える。
小領域が三角形とすると、小領域$ e $、接点$j=(x_j,y_j)$における物理量は
$$ u^e_j=\alpha_1+\alpha_2x^e_j+\alpha_3y^e_j $$
なので、
$$
\begin{pmatrix}
u^e_1 \\
u^e_2 \\
u^e_3 \\
\end{pmatrix}
=
\begin{pmatrix}
1&x^e_1&y^e_1 \\
1&x^e_2&y^e_2 \\
1&x^e_3&y^e_3 \\
\end{pmatrix}
\begin{pmatrix}
\alpha_1\\
\alpha_2\\
\alpha_3\\
\end{pmatrix}
$$
となる。
これを解くと、$\alpha_i$が求めることができ、解（の近似）が求められる。つまり、$u^e_j$が求まれば、解が求まることがわかる。
扱いやすいように形を整えるために、
$$N^e_j(x,y)=a^e_j+b^e_jx+c^e_jy $$を用いて以下のように表すことが出来る。ただし、$a^e_j,b^e_j,c^e_j$は$x^e_i,y^e_i$を用いて表され、既知である。
$$ u\approx=\alpha_1+\alpha_2x+\alpha_3y=\sum_{j=1}^3N^e_j(x,y)u^e_j$$
$N^e_j(x,y)$は、内挿関数や形状関数と呼ばれ、節点の物理量で補間するための重み（割合）であると考えることが出来る。
接点における重みを考えると、
$$ N^e_j(x_i,y_i)=\delta_{ij} $$
であり、任意の点$(x,y) $における解において各接点の割合の和は１になることから、
$$\sum_{j=1}^3N^e_j(x,y)=1$$
となる。
ここで一度原点に戻ると、少領域$e$で
$$ \int_{\Omega_e}w_i\left(\frac{\partial^2 u}{\partial x^2}+\frac{\partial^2 u}{\partial y^2}+f \right)dS=0 $$
となる$u$を求めたかったわけだが、$u=\sum_{j=1}^3N^e_j(x,y)u^e_j$を代入すると
$$ \int_{\Omega_e}w_i\left(\frac{\partial^2 \sum_{j=1}^3N^e_j(x,y)u^e_j}{\partial x^2}+\frac{\partial^2 \sum_{j=1}^3N^e_j(x,y)u^e_j}{\partial y^2}+f \right)dS=0 $$
これは以下のように変形できる。
$$
\sum_{j=1}^3 \int_{\Omega_e}w_i\left(\frac{\partial^2 N^e_j}{\partial x^2}+\frac{\partial^2 N^e_j}{\partial y^2} \right)dSu^e_j+\int_{\Omega_e}w_ifdS=0
$$

ここからガラーキン法(Galerkin法)を使う。
$$ w_i=-N^e_i(x,y) $$とする。ただし、小領域外では0であり、マイナスは式変形しやすくするために便宜上つけたものである。
すると、
$$
\sum_{j=1}^3 \int_{\Omega_e}\left(-N^e_i(x,y)\frac{\partial^2 N^e_j}{\partial x^2}+-N^e_i(x,y)\frac{\partial^2 N^e_j}{\partial y^2} \right)dSu^e_j=\int_{\Omega_e}N^e_ifdS
$$
これを式1とする。

ここで、次に行う式変形のために2次元のガウスの発散定理を思い出す。
$$ \int \mbox{div}{\bf{F}} dS=\int {\bf{F}}\cdot {\bf{n}}dC $$
$ {\bf{F}}=\left(N_i\frac{\partial N_j}{\partial x},N_i\frac{\partial N_j}{\partial y}\right)$を代入すると
$$\int\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+N_i\frac{\partial^2 N_j}{\partial x^2}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}+N_i\frac{\partial^2 N_j}{\partial y^2}dS=\int\left(N_i\frac{\partial N_j}{\partial x}n_x+N_i\frac{\partial N_j}{\partial y}n_y\right)dC$$
これを用いると式1は（面倒なので$e$省略）
$$
\sum_{j=1}^3 \int_{\Omega_e}\left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}\right)dSu^e_j=\int_{\Omega_e}N_ifdS+\int_{\partial\Omega_e}N_i\frac{\partial u}{\partial {\bf{n}}}dC
$$
ただし最後の項は法線方向微分である。
$$ k^e_{ij}=\int_{\Omega_e}\left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+\frac{\partial N_i}{\partial y}\frac{\partial N_j}{\partial y}\right)dS$$とすると、$N^e_j(x,y)=a^e_j+b^e_jx+c^e_jy $より、
$$ k^e_{ij}=\int_{\Omega_e}\left(b_ib_j+c_ic_j\right)dS$$なので計算できる。
$$ \frac{\partial u}{\partial {\bf{n}}} $$は境界条件を設定するのに使える。
$K_e=(k^e_{ij})$とすると、$K_e$は対称行列であり要素剛性行列と呼ばれる。
各要素の$K_e$を上手く合わせれば1次連立方程式に帰着され、$ u^e_j $を求めることが出来る。