初めに

1次元の続き．

有限要素法の幾何学的柔軟性を味わいたい．

Methods

Poisson方程式

以下のPoisson方程式を考える．

\begin{alignat}{2}
    - \nabla^2 u 
    &= f \quad &&\text{in} \quad \Omega \\
    u 
    &= g_D \quad &&\text{on} \quad \Gamma_D \\
    \boldsymbol{n} \cdot \nabla u 
    &= g_N \quad &&\text{on} \quad \Gamma_N
\end{alignat}

ここで，$\Omega \subset \mathbb{R}^d$は有界領域（$d = 2$），$\Gamma = \partial \Omega$は$\Omega$の境界，$\Gamma_D$はDirichlet境界，$\Gamma_N$はNeumann境界である．また，$\Gamma = \bar{\Gamma}_D \cup \bar{\Gamma}_N$，$\Gamma_D \cap \Gamma_N = \emptyset$，$\Gamma_D \neq \emptyset$，$\boldsymbol{n}$は$\Gamma$上の外向き単位法線ベクトルである．なお，$f: \Omega \to \mathbb{R}$，$g_D: \Gamma_D \to \mathbb{R}$，$g_N: \Gamma_N \to \mathbb{R}$は既知の関数とする．

弱形式化

関数空間

\begin{align}
    U 
    &= \lbrace 
        u \in H^1 \left( \Omega \right) 
        \mid \left. u \right|_{\Gamma_D} = g_D
    \rbrace \\
    V 
    &= \lbrace 
        v \in H^1 \left( \Omega \right) 
        \mid \left. v \right|_{\Gamma_D} = 0 
    \rbrace
\end{align}

を定める．任意の$v \in V$をPoission方程式の両辺に乗じて，$\Omega$上で積分する．

\begin{align}
    - \int_{\Omega} \nabla^2 u v \; dx
    &= \int_{\Omega} f v \; dx \\
\end{align}

Greenの第一恒等式から，

\begin{align}
    \int_{\Omega} \nabla u \cdot \nabla v \; dx
    &= \int_{\Omega} f v \; dx + \int_{\Gamma_N} g_N v \; ds
\end{align}

を得る．これを，元のPoission方程式に対する弱形式（weak form）と呼ぶ．なお，

\begin{align}
    \langle u, v \rangle 
    &= \int_{\Omega} \nabla u \cdot \nabla v \; dx \\
    \left( u, v \right) 
    &= \int_{\Omega} u v \; dx \\
    \left[ u, v \right] 
    &= \int_{\Gamma_N} u v \; ds \\
\end{align}

という記号を導入しておくと，弱形式は次のように書ける．

\begin{align}
    \langle u, v \rangle 
    &= \left( f, v \right) 
        + \left[ g_N, v \right] 
\end{align}

弱形式に対する解は弱解（weak solution）と呼ばれていて，この弱解を元の問題の解として採用する，という方針である．

有限要素法

領域を有限個の要素に分割し（$\bar{\Omega} \simeq \hat{\Omega} := \bigcup_{k=1}^{K} e_k$），各要素上でGalerkin近似を考える．2次元問題で最も広く利用される，三角形1次要素（P1要素）を念頭に考える．要素分割の注意として，

要素間の隙間
要素同士の重なり
ある頂点が他の要素の辺上に存在

を認めないこととする．節点には，$\lbrace P_i \rbrace_{i=1}^{M}$という番号を与える（2次元なので，既に$K + 1 = M$などという単純な関係は無い）．

基底関数には，1次元のときの自然な拡張として，$\hat{\Omega}$上で連続で，$e_k$上で$x$と$y$の1次多項式を区分的1次多項式と呼び，その全体を$\hat{W}$で表す．2変数の1次関数$f = c_0 + c_1 x + c_2 y$は平面であるから，$\hat{W}$は三角形を連続に繋ぐ折れ面である．

試行関数の空間$\hat{U}$，試験関数の空間$\hat{V}$は，次のように選ぶ（$\hat{U} \subset \hat{W}, \ \hat{V} \subset \hat{W}$，両者を同じ基底で張る (Bubnov–) Galerkin型の方法）．

\begin{align}
    \hat{U} 
    &= \lbrace 
        \hat{u} \in \hat{W}
        \mid \left. \hat{u} \right|_{\Gamma_D} = \hat{g}_D
    \rbrace \\
    \hat{V} 
    &= \lbrace 
        \hat{v} \in \hat{W}
        \mid \left. \hat{v} \right|_{\Gamma_D} = 0
    \rbrace
\end{align}

弱形式

\begin{align}
    \langle u, v \rangle 
    &= \left( f, v \right) 
        + \left[ g_N, v \right] \quad (v \in V)
\end{align}

を有限次元近似（区分的1次近似）して，

\begin{align}
    \langle \hat{u}, \hat{v} \rangle 
    &= \left( f, \hat{v} \right) 
        + \left[ g_N, \hat{v} \right] \quad (\hat{v} \in \hat{V})
\end{align}

を考える．

面積

要素$e_k$に属する節点を，半時計回りに$N^k_i = (x^k_i, y^k_i) \ (i = 0, 1, 2)$と名付ける．この要素の面積$|e_k|$は，上付き文字を省略して，

\begin{align}
    |e_k|
    &= \frac{1}{2} 
    \det \left(
        \begin{array}{cc}
            x_1 - x_0 & x_2 - x_0 \\
            y_1 - y_0 & y_2 - y_0 \\
        \end{array}
    \right) \\
    &= \frac{1}{2} \left(
        (x_1 - x_0) (y_2 - y_0) - (x_2 - x_0) (y_1 - y_0)
    \right)
\end{align}

相対位置ベクトルが為す平行四辺形の面積を利用している．

面積座標（重心座標）

要素$e_k$に属する節点を，半時計回りに$N^k_i = (x^k_i, y^k_i) \ (i = 0, 1, 2)$と呼んでいる（$k$番目の要素についての話であることに気を付けながら，見易さのため，肩の数字を省略する）．

任意の1次関数

\begin{align}
    \hat{u} (x, y)
    &= \alpha_0 + \alpha_1 x + \alpha_2 y \quad (x, y \in e_k)
\end{align}

は，3節点$N_i$における値$u^i = \hat{u} (N_i) \ (i = 0, 1, 2)$が定まれば，定まる．特に，節点$N_i$で1, その他の節点で0となる1次関数$L_i \ (i = 0, 1, 2)$を導入する．つまり，

\begin{align}
    L_i (N_j)
    &= L_i (x_j, y_j)
     = \delta_{ij} \quad (i, j = 0, 1, 2)
\end{align}

$(L_0, L_1, L_2)$を用いると，任意の1次関数$\hat{u}$は，

\begin{align}
    \hat{u} (x, y)
    &= \sum_{i=0}^{2} u^i L_i (x, y) \quad (x, y \in e_k)
\end{align}

と書ける．即ち，$(L_0, L_1, L_2)$は1次関数の基底になる．

任意の点$P (x, y) \in e_k$に対して，3つの実数の組$(L_0(P), L_1(P), L_2(P))$を，点$P$の面積座標（area coordinate），或いは重心座標（barycentric coordinate）と呼ぶ．多分，数学ではこういう呼び方をしていて，工学では形状関数（shape function）と呼ぶ方が馴染みがある（と思う）．

任意の点$P (x, y) \in e_k$に対して，

\begin{align}
    L_0 (P) + L_1 (P) + L_2 (P) 
    &= 1
\end{align}

が成り立つ．

面積座標における積の積分には，$i, j, k \in \mathbb{Z}_{\ge 0}$に対して

\begin{align}
    \int_{e_l}
        L_0^i L_1^j L_2^k
    \; dx
    &= \frac{2 |e_l| i! j! k!}{(2 + i + j + k)!}
\end{align}

が成り立つことが知られている（全然良く理解していないが，取り敢えず認めよう）．

面積座標の係数

$L_i \ (i = 0, 1, 2)$は，節点$N_i$で1, その他の節点で0となる．その係数を

\begin{align}
    L_i
    &= a_i + b_i x + c_i y
\end{align}

とすると，

\begin{align}
    \delta_{ij}
    &= L_i (N_j)
     = L_i (x_j, y_j) \\
    &= a_i + b_i x_j + c_i y_j \\
    &= \left( 
        \begin{array}{ccc}
            1 & x_j & y_j
        \end{array}
    \right) 
    \left( 
        \begin{array}{c}
            a_i \\
            b_i \\
            c_i
        \end{array}
    \right) \\
\end{align}

より（Kroneckerのdeltaを回すとidentityになるから），

\begin{align}
    I
    &= \left( 
        \begin{array}{ccc}
            1 & x_0 & y_0 \\
            1 & x_1 & y_1 \\
            1 & x_2 & y_2 \\
        \end{array}
    \right) 
    \left( 
        \begin{array}{ccc}
            a_0 & a_1 & a_2 \\
            b_0 & b_1 & b_2 \\
            c_0 & c_1 & c_2 \\
        \end{array}
    \right) \\
\end{align}

であるから，

\begin{align}
    B
    &:= \left( 
        \begin{array}{ccc}
            1 & x_0 & y_0 \\
            1 & x_1 & y_1 \\
            1 & x_2 & y_2 \\
        \end{array}
    \right) 
\end{align}

とすると，

\begin{align}
    \left( 
        \begin{array}{ccc}
            a_0 & a_1 & a_2 \\
            b_0 & b_1 & b_2 \\
            c_0 & c_1 & c_2 \\
        \end{array}
    \right) 
    &= B^{-1}
\end{align}

である．$B^{-1}$を具体的に計算しておく．

\begin{align}
    B^{-1}
    &= \left( 
        \begin{array}{ccc}
            1 & x_0 & y_0 \\
            1 & x_1 & y_1 \\
            1 & x_2 & y_2 \\
        \end{array}
    \right)^{-1}
     = \frac{1}{\det B} \tilde{B} \\
\end{align}

基本変形と行列式の関係から，

\begin{align}
    \det B
    &= \det
    \left( 
        \begin{array}{ccc}
            1 & x_0 & y_0 \\
            1 & x_1 & y_1 \\
            1 & x_2 & y_2 \\
        \end{array}
    \right) \\
    &= \det
    \left( 
        \begin{array}{ccc}
            1 & x_0       & y_0       \\
            0 & x_1 - x_0 & y_1 - y_0 \\
            0 & x_2 - x_0 & y_2 - y_0 \\
        \end{array}
    \right) \\
    &= \det
    \left( 
        \begin{array}{cc}
            x_1 - x_0 & y_1 - y_0 \\
            x_2 - x_0 & y_2 - y_0 \\
        \end{array}
    \right) \\
    &= 2 |e_k| \\
    &= (x_1 - x_0) (y_2 - y_0) - (x_2 - x_0) (y_1 - y_0)
\end{align}

$B$の余因子行列$\tilde{B}$は，

\begin{align}
    \tilde{B}
    &= \left( 
        \begin{array}{ccc}
            (-1)^{1+1} |\tilde{B}_{11}| & (-1)^{1+2} |\tilde{B}_{12}| & (-1)^{1+3} |\tilde{B}_{13}| \\
            (-1)^{2+1} |\tilde{B}_{21}| & (-1)^{2+2} |\tilde{B}_{22}| & (-1)^{2+3} |\tilde{B}_{23}| \\
            (-1)^{3+1} |\tilde{B}_{31}| & (-1)^{3+2} |\tilde{B}_{32}| & (-1)^{3+3} |\tilde{B}_{33}| \\
        \end{array}
    \right)^{\top} \\
    &= \left( 
        \begin{array}{ccc}
            x_1 y_2 - x_2 y_1 & y_1 - y_2 & x_2 - x_1 \\
            x_2 y_0 - x_0 y_2 & y_2 - y_0 & x_0 - x_2 \\
            x_0 y_1 - x_1 y_0 & y_0 - y_1 & x_1 - x_0 \\
        \end{array}
    \right)^{\top} \\
    &= \left( 
        \begin{array}{ccc}
            x_1 y_2 - x_2 y_1 & x_2 y_0 - x_0 y_2 & x_0 y_1 - x_1 y_0 \\
            y_1 - y_2         & y_2 - y_0         & y_0 - y_1 \\
            x_2 - x_1         & x_0 - x_2         & x_1 - x_0 \\
        \end{array}
    \right) \\
\end{align}

即ち，

\begin{align}
    B^{-1}
    &= \left( 
        \begin{array}{ccc}
            a_0 & a_1 & a_2 \\
            b_0 & b_1 & b_2 \\
            c_0 & c_1 & c_2 \\
        \end{array}
    \right) \\
    &= \frac{1}{2 |e_k|} \left( 
        \begin{array}{ccc}
            x_1 y_2 - x_2 y_1 & x_2 y_0 - x_0 y_2 & x_0 y_1 - x_1 y_0 \\
            y_1 - y_2         & y_2 - y_0         & y_0 - y_1 \\
            x_2 - x_1         & x_0 - x_2         & x_1 - x_0 \\
        \end{array}
    \right) \\
\end{align}

要素係数行列

弱形式（の有限次元近似）を再掲する．

\begin{align}
    \langle \hat{u}, \hat{v} \rangle 
    &= \left( f, \hat{v} \right) 
        + \left[ g_N, \hat{v} \right] \quad (\hat{v} \in \hat{V})
\end{align}

記号

\begin{align}
    \langle \hat{u}, \hat{v} \rangle_{e_k}
    &= \int_{e_k} \nabla \hat{u} \cdot \nabla \hat{v} \; dx \\
    \left( f, \hat{v} \right)_{e_k} 
    &= \int_{e_k} f \hat{v} \; dx \\
\end{align}

を導入し，積分は，要素毎に計算しておき，後から和を取る方針（$\int_{\Omega} = \sum_{k} \int_{e_k}$）を取る．即ち，

\begin{align}
    \sum_{k=1}^{K} \langle \hat{u}, \hat{v} \rangle_{e_k} 
    &= \sum_{k=1}^{K} \left( f, \hat{v} \right)_{e_k} 
        + \left[ g_N, \hat{v} \right] \quad (\hat{v} \in \hat{V})
\end{align}

を考える．ここで，$u^j := \hat{u} (N_j), \ v^j := \hat{v} (N_j) \ (j = 0, 1, 2)$を用いて，

\begin{align}
    \hat{u}
    = \sum_{j=0}^{2} u^j L_j , \quad
    \hat{v}
    = \sum_{i=0}^{2} v^i L_i \quad (\text{on } e_k) 
\end{align}

という補間を与えると，

\begin{align}
    \langle \hat{u}, \hat{v} \rangle_{e_k}
    &= \sum_{j=0}^{2} \sum_{i=0}^{2} v^i A_{ij}^{(k)} u^j \\
    \left( f, \hat{v} \right)_{e_k}
    &= \sum_{i=0}^{2} v^i f_i^{(k)} \\
\end{align}

但し，

\begin{align}
    A_{ij}^{(k)}
    &:= \langle L_j, L_i \rangle_{e_k} \\
    f_i^{(k)}
    &:= \left( f, L_i \right)_{e_k}
\end{align}

いま，

\begin{align}
    \boldsymbol{u}_k
    &:= \left( u^0, u^1, u^2 \right)^{\top} \\
    \boldsymbol{v}_k
    &:= \left( v^0, v^1, v^2 \right)^{\top} \\
    \boldsymbol{f}_k
    &:= \left( f_0^{(k)}, f_1^{(k)}, f_2^{(k)} \right)^{\top} \\
    \boldsymbol{A}_k
    &:= \left(
        \begin{array}{ccc}
            A_{00}^{(k)} & A_{01}^{(k)} & A_{02}^{(k)} \\
            A_{10}^{(k)} & A_{11}^{(k)} & A_{12}^{(k)} \\
            A_{20}^{(k)} & A_{21}^{(k)} & A_{22}^{(k)} \\
        \end{array}
    \right)
\end{align}

とおくと，$\boldsymbol{A}_k$は対称で，

\begin{align}
    \langle \hat{u}, \hat{v} \rangle_{e_k}
    &= \boldsymbol{v}_k^{\top} \boldsymbol{A}_k \boldsymbol{u}_k \\
    \left( f, \hat{v} \right)_{e_k}
    &= \boldsymbol{v}_k^{\top} \boldsymbol{f}_k \\
\end{align}

と書ける．ここで，$A_{ij}^{(k)} = \langle L_j, L_i \rangle_{e_k}$は，

\begin{align}
    A_{ij}^{(k)}
    &= \langle L_j, L_i \rangle_{e_k} \\
    &= \int_{e_k} \nabla L_j \cdot \nabla L_i dx \\
    &= \int_{e_k} \nabla (a_j + b_j x + c_j y) \cdot \nabla (a_i + b_i x + c_i y) dx \\
    &= \int_{e_k} 
        \left( 
            \begin{array}{c}
                 b_j \\ c_j 
            \end{array}
        \right)
        \cdot
        \left( 
            \begin{array}{c}
                 b_i \\ c_i
            \end{array}
        \right)
    dx \\
    &= \left( b_j b_i + c_j c_i \right) |e_k|
\end{align}

一方，$f_i^{(k)} = \left( f, L_i \right)_{e_k}$については，$f^j := f (N_j) \ (j = 0, 1, 2)$を用いて，

\begin{align}
    f
    = \sum_{j=0}^{2} f^j L_j
\end{align}

という補間を与え，

\begin{align}
    f_i^{(k)}
    &= \left( f, L_i \right)_{e_k} \\
    &= \left( \sum_{j=0}^{2} f^j L_j, L_i \right)_{e_k} \\
    &= \sum_{j=0}^{2} f^j \left( L_j, L_i \right)_{e_k} \\
    &= \sum_{j=0}^{2} f^j \int_{e_k} L_j L_i dx \\
\end{align}

面積座標における便利な公式

\begin{align}
    \int_{e_k}
        L_0^p L_1^q L_2^r
    dx
    &= \frac{2 \; |e_k| \; p! \; q! \; r!}{(2 + p + q + r)!}
\end{align}

から，

\begin{align}
    \int_{e_k} L_0^2 \; dx
    &= \frac{2 \; |e_k| \; 2! \; 0! \; 0!}{(2 + 2 + 0 + 0)!}
     = \frac{|e_k|}{6} \\
    \int_{e_k} L_0 L_1 \; dx
    &= \frac{2 \; |e_k| \; 1! 1\; ! \; 0!}{(2 + 1 + 1 + 0)!} 
     = \frac{|e_k|}{12} \\
    \int_{e_k} L_1 L_2 \; dx
    &= \frac{2 \; |e_k| \; 0! \; 1! \; 1!}{(2 + 0 + 1 + 1)!} 
     = \frac{|e_k|}{12} \\
    \int_{e_k} L_0 L_2 \; dx
    &= \frac{2 \; |e_k| \; 1! \; 0! \; 1!}{(2 + 1 + 0 + 1)!} 
     = \frac{|e_k|}{12} \\
\end{align}

より，

\begin{align}
    f_i^{(k)}
    &= \sum_{j=0}^{2} f^j \int_{e_k} L_j L_i dx \\
    f_0^{(k)}
    &= f^0 \int_{e_k} L_0 L_0 dx 
        + f^1 \int_{e_k} L_1 L_0 dx
        + f^2 \int_{e_k} L_2 L_0 dx \\
    &= \frac{|e_k|}{12} \left( 2 f^0 + f^1 + f^2 \right) \\
    f_1^{(k)}
    &= \frac{|e_k|}{12} \left( f^0 + 2 f^1 + f^2 \right) \\
    f_2^{(k)}
    &= \frac{|e_k|}{12} \left( f^0 + f^1 + 2 f^2 \right) \\
\end{align}

を得る．

直接剛性法

1次元のときと同様に，系全体に拡大する．全節点には，$\lbrace P_i \rbrace_{i=1}^{M}$という番号を与えていることを思い出す．

$u_i := \hat{u} (P_i), v_i := \hat{v} (P_i) \ (i = 1, ..., M)$として，

\begin{align}
    \boldsymbol{u}
    &:= \left( u_1, u_2, ..., u_M \right)^{\top} \\
    \boldsymbol{v}
    &:= \left( v_1, v_2, ..., v_M \right)^{\top} \\
\end{align}

とする．

要素$e_k \ (k = 1, ..., K)$に属する節点$N^k_0, N^k_1, N^k_2$に対し，

\begin{align}
    N^k_0 = P_{i^k_0}, \
    N^k_1 = P_{i^k_1}, \
    N^k_2 = P_{i^k_2}
\end{align}

となる整数$i^k_0, i^k_1, i^k_2$を取る．これらを全体節点番号と呼ぶ．

続いて，$\boldsymbol{f}_{k}^{*} \in \mathbb{R}^M$を，第$i^k_0, i^k_1, i^k_2$成分が$f^{(k)}_0, f^{(k)}_1, f^{(k)}_2$であり，それ以外の成分がゼロのベクトルとすると，

\begin{align}
    \left( f, \hat{v} \right)_{e_k}
    = \boldsymbol{v}_k^{\top} \boldsymbol{f}_k
    = \boldsymbol{v}^{\top} \boldsymbol{f}_k^{*} \quad (k = 1, ..., K)
\end{align}

同様に，$\boldsymbol{g}_{k}^{*} \in \mathbb{R}^M$をNeumann条件が課された節点でのみ値を持ち，それ以外の成分がゼロのベクトルとすると，

\begin{align}
    \left[ g_N, \hat{v} \right] 
    = \boldsymbol{v}^{\top} \boldsymbol{g}_k^{*} \quad (k = 1, ..., K)
\end{align}

さらに，$\boldsymbol{A}_{k}^{*} \in \mathbb{R}^{M \times M}$をゼロで埋められた行列とする．

但し，第$(i^k_m, i^k_n)$成分には$A_{mn}^{(k)}$を代入する．これを用いると，

\begin{align}
    \langle \hat{u}, \hat{v} \rangle_{e_k}
    = \boldsymbol{v}_k^{\top} \boldsymbol{A}_k \boldsymbol{u}_k
    = \boldsymbol{v}^{\top} \boldsymbol{A}_k^{*} \boldsymbol{u} \quad (k = 1, ..., K)
\end{align}

以上より，弱形式

\begin{align}
    \sum_{k=1}^{K} \langle \hat{u}, \hat{v} \rangle_{e_k} 
    &= \sum_{k=1}^{K} \left( f, \hat{v} \right)_{e_k} 
        + \left[ g_N, \hat{v} \right] 
\end{align}

は，

\begin{align}
    \sum_{k=1}^{K} \boldsymbol{v}^{\top} \boldsymbol{A}_k^{*} \boldsymbol{u}
    &= \sum_{k=1}^{K} \boldsymbol{v}^{\top} \boldsymbol{f}_k^{*}
        + \boldsymbol{v}^{\top} \boldsymbol{g}_k^{*} \\
    \boldsymbol{v}^{\top} 
    \left(
        \sum_{k=1}^{K} \boldsymbol{A}_k^{*} \boldsymbol{u}
        - \sum_{k=1}^{K} \boldsymbol{f}_k^{*}
        - \boldsymbol{g}_k^{*}
    \right)
    &= 0 \\
\end{align}

に同値となる．

\begin{align}
    \boldsymbol{A}^{*}
    &= \sum_{k=1}^{K} \boldsymbol{A}_k^{*} \\
    \boldsymbol{b}^{*}
    &= \sum_{k=1}^{K} \boldsymbol{f}_k^{*} + \boldsymbol{g}_k^{*} \\
\end{align}

とおくと，

\begin{align}
    \boldsymbol{v}^{\top} 
    \left(
        \boldsymbol{A}^{*} \boldsymbol{u}
        - \boldsymbol{b}^{*}
    \right)
    &= 0 \\
\end{align}

である．$\boldsymbol{v} \in \check{V} = \lbrace \boldsymbol{v} \in \mathbb{R}^{M} \mid v_i = 0 \text{ for } i \text{ s.t. } P_i \in \Gamma_D \rbrace$の任意の元であるから，

\begin{align}
    \boldsymbol{A}^{*} \boldsymbol{u}
    - \boldsymbol{b}^{*}
    \in \breve{V}
    &= \lbrace
        \boldsymbol{w} \in \mathbb{R}^{M} 
        \mid w_i = 0 \text{ for } i \text{ s.t. } P_i \notin \Gamma_D
    \rbrace
\end{align}

即ち，

\begin{align}
    \boldsymbol{A}^{**} \boldsymbol{u}
    &= \boldsymbol{b}^{**}
\end{align}

ここで，

\begin{align}
    \boldsymbol{A}^{**}
    &:= \boldsymbol{A}^{*} \text{の第$i$行（$P_i \in \Gamma_D$なる$i$）を除いた行列} \\
    \boldsymbol{b}^{**}
    &:= \boldsymbol{b}^{*} \text{の第$i$成分（$P_i \in \Gamma_D$なる$i$）を除いた縦ベクトル} \\
\end{align}

である．いま，$\boldsymbol{u} \in \mathbb{R}^M$に対して，Dirichlet条件$u_i = g_D (P_i) \ (P_i \in \Gamma_D)$を課す．

\begin{align}
    \boldsymbol{A}^{**}
    &= \left(
        \begin{array}{cccc}
            \boldsymbol{a}_1 & \boldsymbol{a}_2 & ... & \boldsymbol{a}_M \\
        \end{array}
    \right)
\end{align}

と書くと，

\begin{align}
    \sum_{i=1}^{M} \boldsymbol{a}_i u_i
    &= \boldsymbol{b}^{**} \\
    \sum_{P_i \in \Gamma_D} \boldsymbol{a}_i u_i
    + \sum_{P_i \notin \Gamma_D} \boldsymbol{a}_i u_i
    &= \boldsymbol{b}^{**} \\
    \sum_{P_i \notin \Gamma_D} \boldsymbol{a}_i u_i
    &= \boldsymbol{b}^{**} 
    - \sum_{P_i \in \Gamma_D} \boldsymbol{a}_i u_i \\
    \boldsymbol{A} \boldsymbol{u}^{*}
    &= \boldsymbol{b}
\end{align}

を解けば良い．但し，

\begin{align}
    \boldsymbol{A}
    &:= \boldsymbol{A}^{**} \text{の第$i$列（$P_i \in \Gamma_D$なる$i$）を除いた正方行列} \\
    \boldsymbol{u}^{*}
    &:= \boldsymbol{u} \text{の第$i$成分（$P_i \in \Gamma_D$なる$i$）を除いた縦ベクトル} \\
    \boldsymbol{b}
    &:= \boldsymbol{b}^{**} 
    - \sum_{P_i \in \Gamma_D} \boldsymbol{a}_i u_i \\
\end{align}

である．

Results

Square

幾何学的柔軟性は見え辛いが，先ずは最も素朴に．

斉次Dirichlet条件

Dirichlet問題

\begin{alignat}{2}
    - \nabla^2 u &= f \quad &&\text{in} \quad \Omega \\
    u &= g_D \quad &&\text{on} \quad \Gamma_D \\
\end{alignat}

を考え，$\Omega = (-1, 1)^2, \ \Gamma_D = \partial \Omega$とし，$f = 1, \ g_D = 0$とする（斉次Dirichlet境界値問題）．ネット上の拾い物だが，こちらに変数分離による解が載っている．

\begin{align}
    u
    &= \frac{1 - x^2}{2}
        - \frac{16}{\pi^3} \sum_{k = 1, \ k \in \mathrm{odd}}^{\infty}
            \frac{\sin(k \pi (1 + x) / 2)}{k^3 \sinh(k \pi)}
            \left(
                \sinh(k \pi (1 + y) / 2) + \sinh(k \pi (1 - y) / 2)
            \right) \\
    &\simeq \frac{1 - x^2}{2}
        - \frac{16}{\pi^3} \sum_{k = 1, \ k \in \mathrm{odd}}^{201}
            \frac{\sin(k \pi (1 + x) / 2)}{k^3 \sinh(k \pi)}
            \left(
                \sinh(k \pi (1 + y) / 2) + \sinh(k \pi (1 - y) / 2)
            \right) \\
\end{align}

これを有限和で近似しておき，解像度を徐々に増加させながら解くことで，数値解の収束を確認する．scipy.spatial.Delaunayを利用してメッシュを作成し，ソルバーには昔作ったCG法を使った（numpy.allcloseを使って確認したところ，numpy.linalg.solveでも同じ結果であった）．CSCとかCSRとかの工夫はしていない．後で触れるため，有限差分法による解を併記する（今は触れない）．

$M: \text{# node}$	$M=25$	$M=81$	$M=289$	$u (0, 0)$	$\mathrm{Error}$
$u_{\text{FDM}}$
$u_{\text{FEM}}$

1次元のときと同じく，誤差の最大値が2次の速度で収束する．また，全体的な誤差は1次の速度で収束する（1次元のときは見えなかったが，これが良く聞く$\mathcal{O}(h^p)$の収束のことだろうか？と思ったが，Karniadakis&Sherwin2005を見ると，一般には，誤差の$H^1$ノルムが$\sim \mathcal{O}(h^p)$で減少するらしい）．ところで，後で登場するGmshと比べると，隅付近のメッシュの品質に疑問があるが，まあいいか．

斉次Dirichlet条件 + 斉次Neumann条件

混合境界値問題

\begin{alignat}{2}
    - \nabla^2 u &= f \quad &&\text{in} \quad \Omega \\
    u &= g_D \quad &&\text{on} \quad \Gamma_D \\
    \boldsymbol{n} \cdot \nabla u &= g_N \quad &&\text{on} \quad \Gamma_N
\end{alignat}

を考え，

\begin{align}
    \Omega
    &= (0, 1)^2 \\
    \Gamma_D
    &= \lbrace (x, 0) \cup (0, y) \mid x \in (0, 1), y \in (0, 1) \rbrace \\
    \Gamma_N
    &= \lbrace (x, 1) \cup (1, y) \mid x \in (0, 1), y \in (0, 1) \rbrace \\
\end{align}

とし，$f = 1, \ g_D = g_N = 0$とする．要素数$K$を徐々に増加させながら解く．再度，scipy.spatial.Delaunayを利用してメッシュを作成し，ソルバーには昔作ったCG法を使った．CSCとかCSRとかの工夫はしていない．

	$K=8$	$K=32$	$K=128$	$K=512$
$2D$
$3D$

こちらを見ると，恐らく正しそうだ．

また，線形ソルバー（反復法）の比較もしてみたが，昔の実験と同じような結果を得た．

	Residual	Solution
$K=512$
$K=2,048$
$K=8,192$

Residualのカラムで，$M_{(\cdot)}$は反復行列，$\varrho (\cdot)$はスペクトル半径である（タイトルの$M$は節点数，$K$は要素数，$\kappa_{(\cdot)}$は条件数）．SOR法の緩和パラメータは，えいやっと$\omega = 1.8$を選んだ¹．最適 / 準最適ではないだろうが，少なくとも悪い数字ではないはずだ（経験的に，大きい問題ほど，$\omega \to 2$を攻めると効果がある²）．CG法については，戸川1977の

共役勾配法によれば有限回（$n$元の方程式の場合，高々$n$回）で厳密解を得ることができるが，実際には，$n$よりももっとはるかに少ない回数で十分な精度の解が得られる場合が少なくない．

という記述の確からしさが確認できる．自由度が小さいときは何れの方法も収束するが，自由度が大きくなるとGauss-Seidel法などは（現実的な時間内に）収束し切れず，劣悪な解しか得られない（が，スムージングされていることは強調すべきであろう, cf. マルチグリッド法との関係）．

斉次Dirichlet条件 + 非斉次Neumann条件 + non-constant source

昔，大学院で有限要素法の講義を受講したとき，以下のような問題をFreeFEMで解いた事がある（cf. github）．

\begin{alignat}{2}
    - \nabla^2 u 
    &= - \sin \left( 2 \pi \left( x + y \right) \right) \quad &&\text{in} \quad \Omega = (0, 1)^2 \\
    u 
    &= 0      \quad &&\text{on} \quad y=0 \\
    \nabla u  \cdot \boldsymbol{n}
    &= 0.05 \quad &&\text{on} \quad y=1,\ x=0,\ x=1
\end{alignat}

以下はFreeFEMによる有限要素法の解と自作コードの有限差分法の解．同じことを，自作の有限要素法でできるか？当然，できなければならない（勿論，離散化スキームやメッシャーの違いなどに起因する数値解の若干の揺らぎはあるだろうが，この程度の簡単な問題で，各手法における数値解の間に大きな乖離が存在するとは考え辛い）．

先ずは，non-constantなsourceの導入（$h$は解像度を意味する離散化パラメータ）．

$h$	$\sim 0.12$	$\sim 0.06$	$\sim 0.03$
$u_{\text{FDM}}$
$u_{\text{FEM}}$

次いで，斉次Neumann条件への書き換え．

$h$	$\sim 0.12$	$\sim 0.06$	$\sim 0.03$
$u_{\text{FDM}}$
$u_{\text{FEM}}$

最後に，非斉次Neumann条件の導入．

$h$	$\sim 0.12$	$\sim 0.06$	$\sim 0.03$
$u_{\text{FDM}}$
$u_{\text{FEM}}$

FDM，FEM両者は同じ自由度（節点数）である．さらに，全ての節点は全く同じ座標に置いている．但し，FEMにはP1要素を使っているので，要素数としてはFEMが多い．また，FDMにおけるLaplace作用素は2次精度の近似，Neumann境界の取り扱いは1次精度の近似である（外側にghost cellを置いて2次精度にするなどの工夫はしていない）．加えて，線形ソルバーが違っていて，FDMはpoint Jacobi，FEMはConjugate Gradientとしている．全くフェアとは言えないが，収束判定の閾値は同じ値にしている（$\le 10^{-6}$）．また，可視化も同じ方法（matplotlib.pyplot.tricontourf）．

上記のような違いはあるが，FreeFEM，自作FDM，自作FEMの何れの方法においても，殆ど同等の数値解を得たと言って良いだろう．本稿のスコープとしては，自作FEMがFreeFEMと同じような解を得たのは喜ばしい．ところで，FreeFEMの解と自作FDM，自作FEMの$h \sim 0.03$は殆ど同じくらいの節点数を持つ（1軸を32分割，要素の品質は違うから，要素数も多分少し違う）．

なお，解像度の意味では，FDMよりFEMの方が収束速度が速いだろうか？誤差解析とかをきちんと勉強すれば，きちんと比較できる気がする．菊池先生の本（菊池1994，菊池1999）のどちらかに載っていたように記憶している（多分，前者だろうね）．FDMの誤差解析は存在するのだろうか？余り聞いたことが無い（今の時代，FDMをやっている人も少ないし）．

興味深いこととして，最も粗い解像度（$h \sim 0.12$）のときのFDM solとFEM solを比較すると，FEMの方がNeumann境界付近での解が良好であるように見える（特に右上のピーク）．FDMではNeumann条件を1次精度で陽に課した一方，有限要素法ではこれを弱形式に取り込み陽に課した訳ではない．にも関わらず，FEMの方がパワフルに見える（見えるだけかもしれない）．これも，きちんと誤差解析を勉強すれば議論できそうな気がするが，未だ十分な知識を持ち合わせていない．

非斉次Dirichlet条件

非斉次なDirichlet条件を課す問題を探していたところ，Laplace方程式だが，以下を見つけた（Millan+2015）．

$\Omega = (0, 1)^2$において，

\begin{alignat}{2}
    \nabla^2 u &= 0 \quad &&\text{in} \quad \Omega \\
    u &= 0      \quad &&\text{on} \quad x=0,\ x=1,\ y=0 \\
    u &= x(1-x) \quad &&\text{on} \quad y=1
\end{alignat}

に対して，

\begin{align}
    u
    &= \frac{8}{\pi^3}
        \sum_{n=1, \ 3, \ 5, \ ...}^{\infty}
        \frac{1}{n^3 \sinh(n \pi)}
        \sin  (n \pi x)
        \sinh (n \pi y)
\end{align}

が解．$N = \lbrace 1, 21, 51, 101, 201, 401, 801 \rbrace$までの有限級数を作ってみたところ，$N = \lbrace 401, 801 \rbrace$では解が吹き飛んだ．3乗とかの速度で数値が大きくなることによる数値誤差の影響だろう．話題がすっかり変わるが，尤度関数の代わりに対数尤度関数の計算を使うのも，こういう理由で都合が良いのだろう．

$N = 21$	$N = 201$	$N = 401$	$u_{\text{max}}$	$\| u_{\text{max}} - 0.25 \|$

安全のため，$N = 101$までの有限級数で近似し，これを厳密解として取り扱う．

上記の問題から引き続き，FDMとFEMで解く．今回は厳密解が分かっているので，FDM，FEM両者で誤差の収束が議論できるのが嬉しい．

$h$	$\sim 0.12$	$\sim 0.06$	$\sim 0.03$	$\text{Error}$
$u_{\text{FDM}}$
$u_{\text{FEM}}$

相変わらず，FDMとFEM両者は同じだけの節点を持ち，節点の座標はFDM，FEMで共通している．FDMは正方形のセル，FEMはそれを二分するP1要素で解いている．FDMの係数行列が対称でないから，両者にscipy.sparse.linalg.spsolve（直接法）を使った．非対称な系のための反復法（e.g., bicgstab, gmres）も幾つか試したのだけれど，途中から爆発していき，うまく収束の傾向を見られなかった（V字のような収束の傾向）．scipyを使っただけなので理由が良く分からなかったのが残念だが，取り敢えず，ソルバーに起因する問題と解像度に起因する誤差の収束は切り分けるべきだろうと思ったので，一貫して安定して解けた直接法を使った．スコープをFDMとFEMの比較に絞れば，かなりフェアにできただろうかと思っている（或いはバイアスされているかもね）．可視化には，FDMでは四角形セルを用いたことを強調するためにplot_surfaceを，FEMでは三角形要素を用いたことを強調するためにplot_trisurfを使った（可視化の都合だが，いまは分かり易さが大切だと思った）．

FDMは，$L^2$，$L^\infty$の両者が1次程度の速度で収束した（因みに，最初に取り上げた素朴な問題（$- \nabla^2 u = 1 \text{ in } \Omega, u = 0 \text{ on } \Gamma$）に対しても同様の調査を行ったが，やはりFDMでは，$L^2$，$L^\infty$共に1次収束であった）．FEMは，他の問題でも見た通り，$L^2$が1次程度，$L^\infty$が2次程度の速度で収束した．係数行列は殆ど同じなのだけれど，こうして少しでも違いがあるとは面白い（同じになると予想していた）．が，$L^2$の意味では同じくらいだし，大差無いと言ってしまっても良いだろうか？ただ，個人的には，有限差分法の収束を見ることができて嬉しい．だって，局所的な離散化は2次精度（Laplacianは2次精度中心差分近似した）なのに，大域的には1次精度でしか収束しないだなんて，どっひゃー面白いじゃないか（どうでも良いことだが，陽的Euler法でも同様の性質を見たことを思い出す）．

以上の観測結果から大雑把なことを言ってしまうと，FDMとFEMの違いは$L^\infty$の収束速度であるように感じる（勿論，このような簡単な問題に限った話）．混合境界値問題のときに見た，FEMのピークを捉える上手さは，ここから来ているのだろうか．観測ベースなので数学的でないのが残念だが，現在の私の直感的な理解はこんなところだ．やはり，きちんと数学を学ばなければそろそろdead endだ．

Pentagon

アメリカ国防総省の本庁舎．頂点が1つ増えただけだが，有限差分法では既に取り扱いが難しい（有限差分法でも工夫すれば不可能ではない？cf. こちら）．外側の節点の座標はこちらで計算し，メッシュの作成にはGmsh（Geuzaine&Remacle2009）を用いた．$f = 1$とする．

斉次Dirichlet条件

Gmshでは，目安となるメッシュサイズ$h$を指定する（target mesh sizeと呼ばれる）．必ず達成できる訳ではない．図中には，$h$に加えて節点数$M$，要素数$K$を併記した．

	$h=1.6$	$h=0.8$	$h=0.4$	$h=0.2$	$h=0.1$
$2D$
$3D$

2Dの図における背景のグレーの線は，メッシュを意味する．特に，低解像度の解から，要素内の補間がpolynomial of 1st order（折れ面）で形成されているのが確認できる（あれ，可視化上の都合かな，，まあ都合の良い解釈で良いか）．

ところで，scipy.sparse.linalg.spsolveを利用し，且つ係数行列をCSR形式で渡しておくと，問題が大きいときにかなり有利であることが分かった（単にscipy.sparse.linalg.spsolveを使うだけだとpay-offは少ない³）．なお，CSRとCSCは余り差が無かった，対称な所為だろうか，，（きちんと読んでいないが，類似の議論があったのでメモしておく: こちら）

pseudocode

import numpy as np
import scipy as sp
A, b = assemble(A_k, b_k, ...)   # A_k, b_k: element-wise matrix, vector
A, b = set_boundary_condition(A, b, ...)
u_np = np.linalg.solve(A, b)
u_sp = sp.linalg.solve(A, b, assume_a="sym")
u_ss = sp.sparse.linalg.spsolve(A, b)
u_sc = sp.sparse.linalg.spsolve(sp.sparse.csr_matrix(A), b)

target mesh size $h$ を変えながら64回解き，要した時間の平均値と（Bessel補正された）標本標準偏差を記録した．マーカー付き実線は平均値，shadeは1倍の標本標準偏差．

横軸はメッシュサイズなので，右側ほど問題規模が小さく，左側ほど問題規模が大きい（後から気づいたが，sp.sparse.linalg.spsolveは直接法らしい，おっとっと，，）．

plt.spyとかplt.imshowとかを使えば，Aがスパースであることが確認できる．問題規模が大きいほど，sparsityがある．

$h=1.6$	$h=0.4$	$h=0.1$

ここでは，plt.imshowを使った．青い箇所が$<0$，赤い箇所が$>0$，白い箇所が$=0$（背景の白色と同化しているが）．

斉次Dirichlet条件 + 斉次Neumann条件 (1)

左側の境界に斉次Dirichlet条件を課す．その他は斉次Neumann条件．

$h=1.6$	$h=0.8$	$h=0.4$	$h=0.2$	$h=0.1$

斉次Dirichlet条件 + 斉次Neumann条件 (2)

$y \le 0$なる境界ノードに斉次Dirichlet条件を課す．その他は斉次Neumann条件．

$h=1.6$	$h=0.8$	$h=0.4$	$h=0.2$	$h=0.1$

左端に$y = 0$なる境界ノードが有ったり無かったりで様子が違うが，まあご愛嬌かな．

Icosagon (20-gon)

正20角形．本当は円がやりたかったが，Gmshでの円の定義，というか，適切なメッシュ切りが良く理解できなかったので，誤魔化す．勉強し直して戻って来たい．

斉次Dirichlet条件

$h=0.4$	$h=0.2$	$h=0.1$

斉次Dirichlet条件 + 斉次Neumann条件

$h=0.4$	$h=0.2$	$h=0.1$

L-shape

L字型．$f = 1$とする．

斉次Dirichlet条件

	$h=0.8$	$h=0.4$	$h=0.2$	$h=0.1$
$2D$
$3D$

$h=0.8$は劣悪だが，やはりP1要素の補間を感じることができる．高解像度化するにつれ，数値解が収束してゆき，こちらと合致してくるので，正しくコードが書けていると思う．

ところで，Gmshでは，特定の節点周辺のみ高解像度化することができる．Gmshに特別な話ではなく，一般的な有限要素法でも行われる．非一様なメッシュと呼ばれる．

	$h_f=h_c/2$	$h_f=h_c/4$	$h_f=h_c/8$	$h_f=h_c/16$
$h_c=0.8$
$h_c=0.4$
$h_c=0.2$
$h_c=0.1$

$h_f$は凹な箇所付近での細かい（fine）メッシュのサイズ，$h_c$はそれ以外の箇所での粗い（coarse）メッシュのサイズである（タイトルの$h$は$h_c$に対応している）．凹な箇所付近だけ高解像度化できている．自由度の増加を適切に抑えつつ，興味がある領域のみ高解像度化できる（局所的な高解像度化）．

斉次Dirichlet条件 + 斉次Neumann条件

	$h=0.8$	$h=0.4$	$h=0.2$	$h=0.1$
$2D$
$3D$

参照できる解が無いが，等値線がperpendicularに突き刺さり，正しそうだ．

終わりに

ムズいなぁ，，

実は，手法の議論は，殆ど桂田先生の資料を焼き増した形になってしまった，，余り良くない（桂田先生の資料は非常に明快なのだが，我が物顔で色々言うのは違う，という意味）．数値実験を色々試し，その中で自分が気になったことに触れる，という形で差分を作った（つもりだ）．

しかし，コードを自作して，当初の予定通り有限要素法の幾何学的柔軟性を体験できたり，予定外だったけどP1の補間の様子を体験できたり，局所的高解像度化の様子を体験できたりしたのは，喜ばしいことだよね．

それと，Neumann条件を厳密に満たすように近似解を構成している訳ではないのに，全体として近似の精度が良い関数を構成できる，というのは面白い．これはどちらかと言うと，弱形式の枠組み自体の特徴だろうか（cf. 最小二乗法 vs. Galerkin法）．Karniadakis&Sherwin2005のch2.2.1.1には，"in the Galerkin formulation Dirichlet boundary conditions have to be specified explicitly whereas Neumann conditions are dealt with implicitly as part of the formulation."という記述があり，この文は個人的に分かり易いと感じた．

最適な緩和パラメータは，Jacobi法の反復行列のスペクトル半径を用いて計算できるそうだが，わざわざJacobi法の反復行列を作り，固有値計算までやらなければならないのは骨が折れる（大きい問題では特に）． ↩
同様のことが，こちらでも言及されている． ↩
実は，scipy.sparse.linalg.spsolveには，"The square matrix A will be converted into CSC or CSR form"と書いてあるから，特別何もしなくても速くなることを期待したのだけれど（利用しているscipy==1.13.1のドキュメンテーションを見ている）．CSC/CSRに変更しないままspsolveに渡すと，SparseEfficiencyWarning: spsolve requires A be CSC or CSR matrix formatと警告される． ↩