1次元Poisson方程式に対する有限要素法

Last updated at 2024-12-02Posted at 2024-08-28

初めに

何度か，有限差分法（FDM: Finite Difference Method）による数値計算を行ってきた（e.g., shear-driven cavity flow, von Kármán vortex street）．同手法は非常に強力であるものの，形状表現に限界を感じる（乱暴な言い方をすれば，有限差分法には，もうそろそろ，飽きてきた）．

一方，過去に，地盤沈下（1次元のPoisson方程式）に対して最小二乗法やGalerkin法を考えたことがある．Galerkin法の考えを引き継いだ有限要素法（FEM: Finite Element Method）には，その幾何学的柔軟性に魅力を感じる．先ずは，1次元のPoisson方程式に対する有限要素法を学ぶ．勿論，1次元では幾何学的柔軟性を味わえないので，後々，2次元・3次元に拡張していきたいと思っているが，いきなり難しいことをするのは難しい．

以下が勉強になる:

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$は有界領域，$\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 
    &:= \left\{ 
        u \in H^1 \left( \Omega \right) 
        \mid 
        \left. u \right|_{\Gamma_D} = g_D
    \right\} \\
    V 
    &:= \left\{ 
        v \in H^1 \left( \Omega \right)
        \mid
        \left. v \right|_{\Gamma_D} = 0 
    \right\}
\end{align}

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

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

Greenの第一恒等式（Gauss-Greenの式，付録参照）から，

\begin{align}
    - \int_{\Omega} \nabla^2 u v dx
    &= \int_{\Omega} \nabla u \cdot \nabla v dx - \int_{\Gamma} \nabla u \cdot \boldsymbol{n} v ds \\
    &= \int_{\Omega} \nabla u \cdot \nabla v dx - \int_{\Gamma_N} g_N v ds
\end{align}

より，

\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 f, g \rangle 
    &= \int_{\Omega} \nabla f \cdot \nabla g dx \\
    \left(  f, g \right) 
    &= \int_{\Omega} f g dx \\
    \left[  f, g \right] 
    &= \int_{\Gamma_N} f g ds \\
\end{align}

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

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

弱形式に対する解は弱解（weak solution）と呼ばれていて，この弱解を元の問題の解として採用する，という方針である．ところで，解に対する微分の要求が緩められている（微分が1階分だけ，試験関数に押し付けられている）．

(2024/09追記)
最近になって，やっと，試行関数を取り出す空間を

\begin{align}
    V_{g_D} 
    &:= \left\{ 
        w \in H^1 \left( \Omega \right) 
        \mid 
        \left. w \right|_{\Gamma_D} = g_D
    \right\} \\
\end{align}

と書いておいた方が，ずっとすっきりすることが分かった．こうしておくと，試験関数を取り出す空間が

\begin{align}
    V_{0} 
    &:= \left\{ 
        w \in H^1 \left( \Omega \right) 
        \mid 
        \left. w \right|_{\Gamma_D} = 0
    \right\} \\
\end{align}

と書けて，非常にすっきりする（文字が増え辛いし，共通点と相違点が明確になる）．そして，多くの教科書はこのような書き方をしている（$V(g_D)$と$V(0)$など）．勉強を始めた当時の私は，試行関数，試験関数を$u, \ v$と書くから，それぞれが属する空間は$U, \ V$と書くのが明解だろう，などと考えていたのだが，多くの教科書には，当然，意図があるよね，，

しかし，全部書き直す勇気は無いので，以降の議論は相変わらず，$u \in U, \ v \in V$と書く．

Galerkin法

弱解の有限次元近似を作り，これを近似解として採用する．

先ず，関数空間$U, V$の有限次元近似$\hat{U}, \hat{V}$を用意する．これには，

\begin{alignat}{2}
    \hat{g}_D &\simeq g_D \quad &&\text{on} \quad \Gamma_D \\
    \psi_i &= 0 \quad &&\text{on} \quad \Gamma_D
\end{alignat}

なる$\hat{g}_D, \psi_i$を用いて（ただし，$\psi_i \in V \ (i = 1, 2, \ldots, N)$は1次独立），

\begin{align}
    \hat{U} 
    &= \left\{
        \hat{g}_D + \sum_{i=1}^{N} \beta_i \psi_i 
        \mid \boldsymbol{\beta} \in \mathbb{R}^{N}
    \right\} \\
    \hat{V} 
    &= \left\{ 
        \sum_{i=1}^{N} \alpha_i \psi_i 
        \mid \boldsymbol{\alpha} \in \mathbb{R}^{N}
    \right\}
\end{align}

を定める¹．ここから取り出す試行関数$\hat{u} \in \hat{U}$と試験関数$\hat{v} \in \hat{V}$を用いて弱形式を近似する．

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

因みに，

$\hat{U}, \hat{V}$を同じ基底$\lbrace \psi_i \rbrace$で張る方法はGalerkin法
$\hat{U}, \hat{V}$を異なる基底$\lbrace \psi_i \rbrace, \ \lbrace \phi_i \rbrace$で張る方法はPetrov-Galerkin法

と呼ばれる．特に，区別を強調したいとき，前者をBubnov-Galerkin法と呼ぶ場合がある．

さて，弱形式は，任意の$\hat{v} \in \hat{V}$に対して成り立つから，特に$\psi_i$に対しても成り立つ．したがって，

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

を得る．ここで，

\begin{align}
    \hat{u} = \hat{g}_D + \sum_{i=1}^{N} \beta_i \psi_i
\end{align}

を代入すると，

\begin{align}
    \langle \hat{g}_D, \psi_i \rangle + \sum_{j=1}^{N} \beta_j \langle \psi_j, \psi_i \rangle 
    &= \left( f, \psi_i \right) + \left[ g_N, \psi_i \right] \\
    \sum_{j=1}^{N} \beta_j \langle \psi_j, \psi_i \rangle 
    &= \left( f, \psi_i \right) + \left[ g_N, \psi_i \right] - \langle \hat{g}_D, \psi_i \rangle
\end{align}

これが任意の$i = 1, 2, \ldots, N$に対して成り立つから，これは$N$元の連立方程式である．即ち，

\begin{align}
    \boldsymbol{A} \boldsymbol{\beta} &= \boldsymbol{b}
\end{align}

と書ける．ここで，

\begin{align}
    A_{ij} 
    &= \langle \psi_j, \psi_i \rangle \\
    b_i 
    &= \left( f, \psi_i \right) 
        + \left[ g_N, \psi_i \right] 
        - \langle \hat{g}_D, \psi_i \rangle \\
\end{align}

である．

== 雑談ここから ==
以下，菊池1999より．

古典的なRitz-Galerkin法の枠組みでは，基底関数に，固有関数，直交多項式，三角関数などを用いることが多い．こうした関数は本質的に1次元的であり，長方形，直方体には拡張し易いが，境界の形状が複雑になると，Dirichlet条件の近似が難しくなり，近似関数の構成が困難になる．一方，Neumann条件は弱形式に取り込めるので，取り扱いが簡単である．

例えば，菊池1999の例3.1では，1次元のPoisson方程式に対して$\psi_i = \sin (i \pi x) \ (1 \le i \le N)$を用意している．例3.2では，2次元の正方形領域上でのPoisson方程式に対して，$\psi_{ij} = \sin (i \pi x) \sin (j \pi y) \ (1 \le i, j \le N)$を用意している（スペクトル法っぽい？理解が間違っているかもしれないが）．

Neumann条件の扱いの簡便さはあるが，古典的なRitz-Galerkin法の枠組みでは，その幾何学的柔軟性は，凡そ有限差分法と同じ程度なのだろうか？
== 雑談ここまで ==

有限要素法

以降，簡単のため，1次元の問題

\begin{alignat}{2}
    - \frac{d^2}{dx^2} u 
    &= f 
    \quad &&\text{in} \quad x \in (0, 1) \\
    u 
    &= g_D 
    \quad &&\text{on} \quad x = 0 \\
    \frac{d}{dx} u 
    &= g_N 
    \quad &&\text{on} \quad x = 1
\end{alignat}

を考える．関数空間$U, \ V$を

\begin{align}
    U
    &= \left\{ w \in H^1 (\Omega) \mid w(0) = g_D \right\} \\
    V
    &= \left\{ w \in H^1 (\Omega) \mid w(0) = 0 \right\}
\end{align}

と定め，

\begin{align}
    \langle u, v \rangle 
    &= \int_{\Omega} \nabla u \cdot \nabla v dx 
     = \int_{0}^{1} u^{\prime} v^{\prime} dx \\
    \left( f, v \right) 
    &= \int_{\Omega} f v dx 
     = \int_{0}^{1} f v dx \\
\end{align}

とおくと，弱形式

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

を満たす弱解$u \in U$を探すことが，我々のやることである．

有限要素法は，試行関数$\hat{u}$と試験関数$\hat{v}$に区分多項式を用いるGalarkin法である．$U, V$の有限次元近似$\hat{U}, \hat{V}$を定めれば，一つのGalarkin解が得られる．特に，区分的多項式を用いて$\hat{U}, \hat{V}$を定めることで有限要素解を定義する．

領域を有限個の要素に分割し（$\bar{\Omega} \simeq \hat{\Omega} := \bigcup_{k=1}^{K} e_k$²），各要素上でGalerkin法を適用する．1次元の問題では，$\bar{\Omega} = [0, 1]$を$K$個の要素（finite element，或いは単にelement）に分割する．各要素$e_k$は，

\begin{align}
    e_k 
    = \left[ x_{k-1}, x_{k} \right] \quad (k = 1, 2, \ldots, K)
\end{align}

とする．即ち，$K$個の要素（element）と$K+1$個の節点（node）がある．

各要素$e_k$上で1次多項式を用意し，それらを連続に接続する．そのような関数を区分1次多項式と呼び，その空間を$\hat{W}$と書く．試行関数$\hat{u}$と試験関数$\hat{v}$として，区分1次多項式を用いる．即ち，$\hat{U} \subset \hat{W}$，$\hat{V} \subset \hat{W}$を満たすよう定める．

$\hat{W}$の基底関数$\lbrace \psi_i \rbrace_{i=0}^{K}$を，次のように定める: $\psi_i$は区分1次多項式で，節点$x_i$で$\psi_i = 1$，節点$x_j \ (j \neq i)$で$\psi_i = 0$という値を取る．つまり，$\psi_i (x_j) = \delta_{ij} \ (j = 0, 1, \ldots, K)$．

具体的な形を書けば，

\begin{align}
    \psi_i (x)
    =
    \begin{cases}
        \frac{x - x_{i-1}}{x_{i} - x_{i-1}} & (x \in [x_{i-1}, x_i]) \\
        \frac{x_{i+1} - x}{x_{i+1} - x_{i}} & (x \in [x_i, x_{i+1}]) \\
        0 & (\mathrm{otherwise})
    \end{cases}
\end{align}

である（が，覚えておく必要は無く，ハット型の関数の集まりというイメージが大切らしい）．

さて，試行関数の空間$\hat{U}$と試験関数の空間$\hat{V}$を以下のように定める．

\begin{align}
    \hat{U} 
    &= \left\{ \hat{w} \in \hat{W} \mid \hat{w} (0) = g_D \right\} \\
    \hat{V} 
    &= \left\{ \hat{w} \in \hat{W} \mid \hat{w} (0) = 0 \right\}
\end{align}

基底関数$\psi_i$を用いて書くと，

\begin{align}
    \hat{U} &= \left\{ 
        \hat{g}_D \psi_0 + \sum_{i=1}^{K} \beta_i \psi_i 
        \mid \boldsymbol{\beta} \in \mathbb{R}^{K}
    \right\} \\
    \hat{V} &= \left\{ \sum_{j=1}^{K} \alpha_j \psi_j \mid \boldsymbol{\alpha} \in \mathbb{R}^{K} \right\}
\end{align}

以上で定まるGalerkin解$\hat{u} \in \hat{U}$は，

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

を満たす．この$\hat{u}$を，区分的1次要素（P1要素³）を用いた有限要素解と呼ぶ．

ここからの方針は，各要素$e_k$毎に要素係数行列$\boldsymbol{A}_k$と要素右辺ベクトル$\boldsymbol{b}_k$を計算し，それらを組み合わせて全体系での連立一次方程式を作り，解く．

長さ座標

要素$e_k = [ x_{k-1}, x_k ]$について，以下の1次関数$L_0, L_1$を用意する．

\begin{align}
    L_0 (x) &= \frac{x_k - x}{h_k} \\
    L_1 (x) &= \frac{x - x_{k-1}}{h_k}
\end{align}

ここで，$h_k = x_k - x_{k-1}$は要素$e_k$の長さである．$L_0, L_1$を要素$e_k$の長さ座標，或いは形状関数と呼ぶ．これは，

\begin{align}
    L_0 (x_{k-1}) &= 1, \quad L_0 (x_k) = 0 \\
    L_1 (x_{k-1}) &= 0, \quad L_1 (x_k) = 1 \\
    L_0 (x) + L_1 (x) &= 1 \quad (x \in e_k)
\end{align}

を満たす．少し前に登場した$\hat{W}$の基底関数と，いま導入した長さ座標（形状関数）を図示しておく．

基底関数（P1要素のもの），$\psi$	長さ座標（形状関数），$L_0, L_1$

$\hat{w} \in \hat{W}$を$L_0, L_1$を用いて次のように展開する．

\begin{align}
    \hat{w} (x)
    &= w (x_{k-1}) L_0 (x) + w (x_k) L_1 (x) \\
    &= w_{k-1} L_0 (x) + w_k L_1 (x) \\
    &= \sum_{j=0}^{1} w_{k+j-1} L_j (x) \quad (x \in e_k)
\end{align}

ここで，$w_k = w (x_k)$とした．

弱形式に登場する

\begin{align}
    \langle u, v \rangle 
    &= \int_{\Omega} \nabla u \cdot \nabla v dx 
     = \int_{0}^{1} u^{\prime} v^{\prime} dx \\
    \left( f, v \right) 
    &= \int_{\Omega} f v dx 
     = \int_{0}^{1} f v dx \\
\end{align}

に対して，各要素$e_k$で

\begin{align}
    \langle u, v \rangle_{e_k} 
    &= \int_{e_k} \nabla u \cdot \nabla v dx 
     = \int_{x_{k-1}}^{x_{k}} u^{\prime} v^{\prime} dx \\
    \left( f, v \right)_{e_k} 
    &= \int_{e_k} f v dx 
     = \int_{x_{k-1}}^{x_{k}} f v dx \\
\end{align}

という記号を導入すると，先に記した弱形式

\begin{align}
    \langle \hat{u}, \hat{v} \rangle 
    &= \left( f, \hat{v} \right) 
        + g_N \hat{v} (1) \quad (\hat{v} \in \hat{V})
\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} 
        + g_N \hat{v} (1) \quad (\hat{v} \in \hat{V})
\end{align}

と書ける（$\int_{\Omega} = \sum_{k=1}^{K} \int_{e_k}$）．

さて，

\begin{align}
    \hat{W} \ni \hat{w}
    = \sum_{j=0}^{1} w_{k+j-1} L_j
\end{align}

を思い出しながら$\langle \hat{u}, \hat{v} \rangle_{e_k}$を具体的に書くと，

\begin{align}
    \langle \hat{u}, \hat{v} \rangle_{e_k}
    &= \langle \sum_{j=0}^{1} u_{k+j-1} L_j, \sum_{i=0}^{1} v_{k+i-1} L_i \rangle_{e_k} \\
    &= \sum_{j=0}^{1} \sum_{i=0}^{1} u_{k+j-1} v_{k+i-1} \langle L_j, L_i \rangle_{e_k} \\
    &= \sum_{i=0}^{1} \sum_{j=0}^{1} v_{k+i-1} A_{ij}^{(k)} u_{k+j-1}
\end{align}

ただし，$A_{ij}^{(k)} := \langle L_j, L_i \rangle_{e_k}$とした．

一方，$\left( f, \hat{v} \right)_{e_k}$は，

\begin{align}
    \left( f, \hat{v} \right)_{e_k}
    &= \left( f, \sum_{j=0}^{1} v_{k+j-1} L_j \right)_{e_k} \\
    &= \sum_{j=0}^{1} v_{k+j-1} \left( f, L_j \right)_{e_k} \\
    &= \sum_{j=0}^{1} v_{k+j-1} f_{j}^{(k)}
\end{align}

ただし，$f_{j}^{(k)} := \left( f, L_j \right)_{e_k}$とした．

また，$g_N \hat{v} (1)$については，$\hat{v} (1) = v_K$より，

\begin{align}
    g_N \hat{v} (1)
    &= g_N v_K
\end{align}

を得る．

要素係数行列，要素自由項ベクトル

要素$e_k = [ x_{k-1}, x_k ]$において，

\begin{align}
    \boldsymbol{u}_{k} &= (u_{k-1}, u_k)^{\top} \\
    \boldsymbol{v}_{k} &= (v_{k-1}, v_k)^{\top} \\
    \boldsymbol{A}_{k} &=
    \begin{pmatrix}
        A_{00}^{(k)} & A_{01}^{(k)} \\
        A_{10}^{(k)} & A_{11}^{(k)}
    \end{pmatrix} \\
    \boldsymbol{f}_k &= (f_{0}^{(k)}, f_1^{(k)})^{\top} \\
    \boldsymbol{g}_K &= \left( 0, g_N \right)^{\top}
\end{align}

とすると，$k = 1, 2, ..., 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 \\
    g_N \hat{v} (1)
    &= \boldsymbol{v}_{K}^{\top} \boldsymbol{g}_{K}
\end{align}

を得る．

$\boldsymbol{u}_k, \ \boldsymbol{v}_k$を要素節点パラメータ・ベクトル
$\boldsymbol{A}_{k}$を要素係数行列
$\boldsymbol{f}_{k}$を要素自由項ベクトル

と呼ぶ．

いま，要素係数行列$\boldsymbol{A}_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_{x_{k-1}}^{x_k} L^{\prime}_j \cdot L^{\prime}_i dx \\
    &= \int_{x_{k-1}}^{x_k} \frac{\varepsilon}{h_k^2} dx \\
    &= \frac{\varepsilon}{h_k^2} \left. x \right|_{x_{k-1}}^{x_k} \\
    &= \frac{\varepsilon}{h_k} \quad (\varepsilon = 1 \text{ if } i = j, \ \varepsilon = -1 \text{ if } i \neq j)
\end{align}

即ち，

\begin{align}
    \boldsymbol{A}_{k} 
    &= \frac{1}{h_k}
    \begin{pmatrix}
         1 & -1 \\
        -1 &  1
    \end{pmatrix}
     = \frac{1}{x_{k} - x_{k-1}}
    \begin{pmatrix}
         1 & -1 \\
        -1 &  1
    \end{pmatrix} \\
\end{align}

$\boldsymbol{f}_k$も計算してみる．

\begin{align}
    f_{j}^{(k)}
    &= \left( f, L_j \right)_{e_k} \\
    &= \int_{x_{k-1}}^{x_k} f (x) L_j (x) dx \quad (j = 0, 1)
\end{align}

こちらは，何らかの求積法を用いて計算しておく．例えば，

\begin{align}
    f (x) \simeq f (x_{k-1}) L_0 (x) + f (x_{k}) L_1 (x) \quad (x \in e_k)
\end{align}

という線形補間を利用し，

\begin{align}
    f_{j}^{(k)}
    &= \left( f, L_j \right)_{e_k} \\
    &\simeq f (x_{k-1}) \left( L_0, L_j \right)_{e_k} 
            + f (x_{k}) \left( L_1, L_j \right)_{e_k}
\end{align}

として，$x = x_{k-1} + t (x_{k} - x_{k-1}) \ (0 \le t \le 1)$と変数変換すれば，

\begin{align}
    \left( L_0, L_0 \right)_{e_k}
     = \left( L_1, L_1 \right)_{e_k}
    &= h_k \int_{0}^{1} t^2 dt
     = \frac{h_k}{3} \\
    \left( L_0, L_1 \right)_{e_k}
     = \left( L_1, L_0 \right)_{e_k}
    &= h_k \int_{0}^{1} t (1 - t) dt
     = \frac{h_k}{6}
\end{align}

故に，

\begin{align}
    \boldsymbol{f}_k 
    &\simeq \frac{h_k}{6} (2 f (x_{k-1}) + f (x_k), f (x_{k-1}) + 2 f (x_k))^{\top} \\
\end{align}

以上より，弱形式$\langle \hat{u}, \hat{v} \rangle = \left( f, \hat{v} \right) + \left[ g_N, \hat{v} \right]$は，以下のように改められる．

\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}

或いは，

\begin{align}
    \sum_{k=1}^{K} \boldsymbol{v}_{k}^{\top} \boldsymbol{A}_k \boldsymbol{u}_{k} &= \sum_{k=1}^{K} \boldsymbol{v}_{k}^{\top} \boldsymbol{f}_k + \boldsymbol{v}_{K}^{\top} \boldsymbol{g}_{K}
\end{align}

である．

直接剛性法

ところで，$\boldsymbol{u}_k, \boldsymbol{v}_k, \boldsymbol{A}_k, \boldsymbol{f}_k, \boldsymbol{g}_K$は要素$e_k$での話であった．

いま，1次元の問題を考えているので，この要素は2つの節点$x_{k-1}, x_k$を持ち，これら（$\boldsymbol{u}_k, \boldsymbol{v}_k, \boldsymbol{A}_k, \boldsymbol{f}_k, \boldsymbol{g}_K$）は2次元のベクトルや2x2の行列である．

全節点での自由度に拡張する．まず，$\boldsymbol{u}, \boldsymbol{v}$を

\begin{align}
    \boldsymbol{u} &= (u_0, u_1, \ldots, u_K)^{\top} \\
    \boldsymbol{v} &= (v_0, v_1, \ldots, v_K)^{\top}
\end{align}

と定める．次に，

\begin{align}
    \boldsymbol{A}_k^{*}
    &= \begin{pmatrix}
        0 & \cdots & 0 & 0 & \cdots & 0 \\
        \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
        0 & \cdots & A_{00}^{(k)} & A_{01}^{(k)} & \cdots & 0 \\
        0 & \cdots & A_{10}^{(k)} & A_{11}^{(k)} & \cdots & 0 \\
        \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\
        0 & \cdots & 0 & 0 & \cdots & 0
    \end{pmatrix}, \quad
    \boldsymbol{f}_k^{*}
    = \begin{pmatrix}
        0 \\
        \vdots \\
        f_{0}^{(k)} \\
        f_{1}^{(k)} \\
        \vdots \\
        0
    \end{pmatrix}, \quad
    \boldsymbol{g}_K^{*} 
    = \begin{pmatrix}
        0 \\
        \vdots \\
        0 \\
        0 \\
        \vdots \\
        g_N
    \end{pmatrix}
\end{align}

とすると，

\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} \\
    \left( f, \hat{v} \right)_{e_k}
    &= \boldsymbol{v}_{k}^{\top} \boldsymbol{f}_k
    = \boldsymbol{v}^{\top} \boldsymbol{f}_k^{*} \\
    \left[ g_N, \hat{v} \right]
    &= \boldsymbol{v}_{K}^{\top} \boldsymbol{g}_{K}
    = \boldsymbol{v}^{\top} \boldsymbol{g}_K^{*}
\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^{*} \quad (\boldsymbol{v} \in \check{V})
\end{align}

と書ける．ここで，

\begin{align}
    \check{V} := \left\{ \boldsymbol{v} \in \mathbb{R}^{K+1} \mid v_0 = 0 \right\}
\end{align}

である．この式は，任意の$\boldsymbol{v} \in \check{V}$に対して

\begin{align}
    \boldsymbol{v}^{\top} \left( \sum_{k=1}^{K} \boldsymbol{A}_k^{*} \right) \boldsymbol{u} 
    &= \boldsymbol{v}^{\top} \left( \sum_{k=1}^{K} \boldsymbol{f}_k^{*} 
        + \boldsymbol{g}_K^{*} \right) \quad (\boldsymbol{v} \in \check{V}) \\
\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 \quad (\boldsymbol{v} \in \check{V})
\end{align}

$\boldsymbol{v} \in \check{V}, \ \check{V} = \lbrace \boldsymbol{v} \in \mathbb{R}^{K+1} \mid v_0 = 0 \rbrace$より，

\begin{align}
    \left(
        \boldsymbol{A}^{*} \boldsymbol{u} - \boldsymbol{b}^{*}\text{の最初の成分以外}
    \right)
    &= 0
\end{align}

が要求される．即ち，

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

ここで，

\begin{align}
    \boldsymbol{A}^{**}
    &:= \boldsymbol{A}^{*} \text{の第0行を除いた} K \times (K+1) \text{行列}
    \\
    \boldsymbol{b}^{**}
    &:= \boldsymbol{b}^{*} \text{の第0成分を除いた} K \text{次元縦ベクトル}
\end{align}

となっている．$\boldsymbol{u} = (u_0, u_1, \ldots, u_K)^{\top}$であり，$u_0$には$u_0 = g_D$なるDirichlet境界条件が課されているので，ここを取り除けば$K \times K$の正方行列になる．

\begin{align}
    \boldsymbol{u}^{*}
    &:= \boldsymbol{u} \text{の第0成分を除いた} K \text{次元縦ベクトル} \\
    \boldsymbol{A}
    &:= \boldsymbol{A}^{**} \text{の第0列を除いた} K \text{次正方行列} \\
\end{align}

を定義する．即ち，

\begin{align}
    \boldsymbol{u} 
    &= (u_0, \boldsymbol{u}^{*})^{\top} 
     = (g_D, \boldsymbol{u}^{*})^{\top} \\
    \boldsymbol{A}^{**}
    &= \left(
        \begin{array}{c|ccc} 
            A_{10}^{*} & A_{11}^{*} & \cdots & A_{1K}^{*} \\
            \vdots     & \vdots     & \ddots & \vdots     \\
            A_{K0}^{*} & A_{K1}^{*} & \cdots & A_{1K}^{*} \\
        \end{array} 
    \right)
     = \left(
        \begin{array}{c|ccc} 
            A_{10}^{*} &  &   & \\
            \vdots     &  & A & \\
            A_{K0}^{*} &  &   & \\
        \end{array} 
    \right)
\end{align}

上式を改めて見ると，

\begin{align}
    \boldsymbol{A}^{**} \boldsymbol{u} 
    &= \left(
        \begin{array}{c|ccc} 
            A_{10}^{*} &  &   & \\
            \vdots     &  & A & \\
            A_{K0}^{*} &  &   & \\
        \end{array} 
    \right)
    \left(
        \begin{array}{c} 
            g_D \\
            \hline
            \boldsymbol{u}^{*} \\
        \end{array} 
    \right) \\
    &= g_D \left(
        \begin{array}{c} 
            A_{10}^{*} \\
            \vdots     \\
            A_{K0}^{*} \\
        \end{array} 
    \right)
    + A \boldsymbol{u}^{*}
\end{align}

より，

\begin{align}
    \boldsymbol{A}^{**} \boldsymbol{u} 
    &= \boldsymbol{b}^{**} \\
    A \boldsymbol{u}^{*}
    &= \boldsymbol{b} 
    := \boldsymbol{b}^{**} 
    - g_D \left(
        \begin{array}{c} 
            A_{10}^{*} \\
            \vdots     \\
            A_{K0}^{*} \\
        \end{array} 
    \right)
\end{align}

を解けば良い．

以下が重要である:

係数行列はDirichlet条件に対応する節点番号の行と列を抜いたもの
Dirichlet条件の情報は右辺のベクトルに組み込む

Results

簡単な問題

左端にDirichlet条件が，右端にNeumann条件が課された1次元のPoisson方程式を考える．

\begin{alignat}{2}
    - \frac{d^2}{dx^2} u 
    &= f 
    \quad &&\text{in} \quad x \in (0, 1) \\
    u 
    &= g_D 
    \quad &&\text{on} \quad x = 0 \\
    \frac{d}{dx} u 
    &= g_N 
    \quad &&\text{on} \quad x = 1
\end{alignat}

$f$が定数関数（例えば，$f = 1$）ならば，

\begin{align}
    \frac{d^2}{dx^2} u 
    &= - 1 \\
    \frac{d}{dx} u 
    &= - x + c_1 \\
    u 
    &= - \frac{1}{2} x^2 + c_1 x + c_2 \\
\end{align}

境界条件から，

\begin{align}
    c_1
    &= g_N + 1 \\
    c_2
    &= g_D \\
\end{align}

よって，

\begin{align}
    u 
    &= - \frac{1}{2} x^2 + (g_N + 1) x + g_D \\
\end{align}

$(g_D, g_N) = (0, 0)$の場合を，要素数$K$を変化させながら解く．ソルバーには，昔自作したCG法を使った（numpy.linalg.solveでも同じ結果であった）．

$K=2$	$K=4$	$K=8$

$K = 4$に固定し，$(g_D, g_N)$を変化させながら解く．

$(g_D, g_N) = (0, -1)$	$(g_D, g_N) = (0, -0.4)$	$(g_D, g_N) = (1, 1)$

正しそうだ．

Galerkin法を学んだときに見たように，Neumann条件は厳密に満たされていない．が，"柔らかい"意味で満たされている．また，Neumann条件を厳密に満たしていなくても，近似解は良質と言えるだろう．

（比較的？）簡単でない問題

同じ問題を考える．

\begin{alignat}{2}
    - \frac{d^2}{dx^2} u &= f 
    \quad &&\text{in} \quad x \in (0, 1) \\
    u &= g_D 
    \quad &&\text{on} \quad x = 0 \\
    \frac{d}{dx} u &= g_N 
    \quad &&\text{on} \quad x = 1
\end{alignat}

創成解の方法（MMS: Method of Manufactured Solutions, see, e.g., Roache2002）を用いて，解$u$を用意する．以下のように定めておく．

\begin{equation}
    u = \pi + x \sin(2 \pi x)
\end{equation}

よって，

\begin{align}
    f &= 4 \pi \left( \pi x \sin(2 \pi x) - \cos(2 \pi x) \right) \\
    g_D &= \pi \\
    g_N &= 2 \pi \\
\end{align}

を得る．$u$と$f$の概形を描画すると，以下のようになる（描画すること自体に意味を感じないが，取り敢えずイメージを掴んでおく）．

$u$	$f$

先程と同様に，要素数$K$を変化させながら解く．

$K=8$	$K=16$	$K=32$	Error

先程と異なり，節点で必ず厳密解と合致する，とはいかないが，要素数の増加に伴い，合致するようになってくる．若干不安になりかけたが，"続・有限要素法による流れのシミュレーション"第2章にも，"解が3次までの多項式の場合は，1次要素では節点上で真なる解に一致する"という旨の記述があり，つまり3次までの多項式で表せない場合，節点上で解に一致しなくても正しいのだろう．また，解が滑らかなとき，1次要素の誤差は要素の最大長さ$h$の2乗に比例する，という記述があるが，手元の計算結果と合致する．

因みに，同書の問題（の一つ）は，

\begin{align}
    u 
    &= -1 + \frac{5}{3} x + \sin \left( \frac{5}{3} \pi x \right) \\
    f 
    &= \frac{25}{9} \pi^2 \sin \left( \frac{5}{3} \pi x \right) \\
    g_D 
    &= -1 \\
    g_N 
    &= \frac{5}{3} + \frac{5}{6} \pi \\
\end{align}

であり，手元の計算で同書と殆ど同じ計算結果を得る．

$K=4$	$K=8$	$K=16$	Error

独学で作ったコードが教科書と同じ挙動をするのは喜ばしい．

ところで，要素長を小さくすることを，$h$-refinementと呼ぶ．基底の次数を大きくすることを，$p$-refinementと呼ぶ．$h$-refinementと$p$-refinementを同時に使う場合，$hp$-refinementという呼び方をする（例えば，松尾+2017）．

$h$- / $p$-refinementに関する説明は，こちらが分かり易い．この仲間に，自由度は変えないまま，節点を移動させることで，メッシュの濃度を動的に変化させる$r$-refinementというものもあるらしい（AMR（adaptive mesh refinement）は自由度が保存されているか？それとも変化しているか？良く知らない，，）．

終わりに

昔から，何度も勉強しては理解できなくてを繰り返していたことだが，やっと，少しだけ，分かった気になった．たった1次元で，こんなに難しく感じるのは，気のせいだろうか．

コードを自作したが，今時はFreeFEMやFEniCSなどを使うのが普通なんだろう（過去，有限要素法の講義を受けたことがあるが，そこでもFreeFEMを使った）．

Appendix

Greenの恒等式（Green's identities）

Greenの第一恒等式

弱形式を得る途中に用いたGreenの第一恒等式を記しておく．

まず，Gaussの発散の定理を記す．
$$
\begin{equation}
\int_{\Omega} \nabla \cdot \boldsymbol{A} dx
= \int_{\Gamma} \boldsymbol{A} \cdot \boldsymbol{n} ds
\end{equation}
$$
いま，$\boldsymbol{A}$が，スカラー関数$u, v$を用いて$\boldsymbol{A} = v \nabla u$と表されるとする．

このとき，

\begin{align}
    \nabla \cdot \boldsymbol{A} 
    &= \nabla \cdot (v \nabla u) \\
    &= \nabla v \cdot \nabla u + v \nabla^2 u
\end{align}

および

\begin{align}
    \boldsymbol{A} \cdot \boldsymbol{n}
    &= (v \nabla u) \cdot \boldsymbol{n} \\
    &= v \nabla u \cdot \boldsymbol{n}
\end{align}

より，元のGaussの発散の定理から，

\begin{align}
    \int_{\Omega} \nabla v \cdot \nabla u dx + \int_{\Omega} v \nabla^2 u dx
    &= \int_{\Gamma} v \nabla u \cdot \boldsymbol{n} ds
\end{align}

を得る．

この式は，教科書や文脈によって，Gauss-Greenの式などと呼ばれる場合もある（個人的には，わざわざ新しい名前を与えなくても，Gaussの発散の定理と呼んで良いのではないか，とも感じるが，まあしかし，GaussもGreenも偉大だ）．

Greenの第二恒等式

本稿には登場しないが，勉強の記録として，Greenの第二定理を記しておく．

Gaussの発散の定理において，ベクトル場$\boldsymbol{A}$を$\boldsymbol{A} = u \nabla v - v \nabla u$で与える．

このとき，

\begin{align}
    \nabla \cdot \boldsymbol{A} 
    &= \nabla \cdot (u \nabla v - v \nabla u) \\
    &= \nabla u \cdot \nabla v + u \nabla^2 v - \nabla v \cdot \nabla u - v \nabla^2 u \\
    &= u \nabla^2 v - v \nabla^2 u
\end{align}

および

\begin{align}
    \boldsymbol{A} \cdot \boldsymbol{n}
    &= (u \nabla v - v \nabla u) \cdot \boldsymbol{n} \\
\end{align}

より，元のGaussの発散の定理から，

\begin{align}
    \int_{\Omega}
        \left( u \nabla^2 v - v \nabla^2 u \right)
    dx
    &= \int_{\Gamma} 
        \left( u \nabla v - v \nabla u \right) \cdot \boldsymbol{n} 
    ds \\
    &= \int_{\Gamma} 
        \left(
            u \frac{\partial v}{\partial \boldsymbol{n}} 
            - v \frac{\partial u}{\partial \boldsymbol{n}}
        \right)
    ds \\
\end{align}

を得る．

Greenの第三恒等式

こちらも，本稿には登場しない．記そうか迷ったが，勉強している内に大変面白そうで，かつ大きい内容でありそうだと思ったので，議論の発散を避けて別に纏めることにする．こちらとこちらが詳しい．

最良近似 / 誤差最小の原理

最良近似 / 誤差最小の原理．どちらが正しい名前か分からないが．

Poisson方程式の弱解$u \in U$は，

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

を満たす．一方，Ritz-Galerkin解$\hat{u} \in \hat{U}$は，

\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}

を満たす．仮定より，$\hat{U} \subset U, \ \hat{V} \subset V$である．$\hat{V} \subset V$より，

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

が成り立つ．これより，上記2式の差分を取り，

\begin{align}
    \langle \hat{u} - u, \hat{v} \rangle 
    &= 0 \quad (\hat{v} \in \hat{V})
\end{align}

を得る．ここで，適当な$\hat{w} \in \hat{U}$を取れば，$\hat{u} - \hat{w} \in \hat{V}$なので，

\begin{align}
    \langle \hat{u} - u, \hat{u} - \hat{w} \rangle 
    &= 0 \\
    \langle \hat{u} - u, \hat{u} - u + u - \hat{w} \rangle 
    &= 0 \\
    \therefore
    \langle \hat{u} - u, \hat{u} - u \rangle 
    &= \langle \hat{u} - u, \hat{w} - u \rangle \\
\end{align}

以下のノルムを定義し，

\begin{align}
    || u ||
    &= \left( u, u \right)^{1/2} 
     = \left( \int_{\Omega} u^2 dx \right)^{1/2} \\
    ||| u |||
    &= \langle u, u \rangle^{1/2} 
     = \left( \int_{\Omega} 
         \left( \frac{\partial u}{\partial x_1} \right)^2
         + \left( \frac{\partial u}{\partial x_2} \right)^2
         + \left( \frac{\partial u}{\partial x_3} \right)^2
     dx \right)^{1/2} \\
\end{align}

右辺をSchwarzの不等式（$|(u, v)| \le ||u|| \ ||v||$）で押さえつける．

\begin{align}
    ||| \hat{u} - u |||^2
    &\le  ||| \hat{u} - u ||| \ ||| \hat{w} - u ||| \\
    \therefore
    ||| \hat{u} - u |||
    &\le ||| \hat{w} - u ||| \\
\end{align}

左辺はRitz-Galerkin解と弱解との誤差であり，右辺は$\hat{U}$に属する任意の関数と弱解との誤差である．即ち，Ritz-Galerkin解$\hat{u}$は，$\hat{U}$から選ぶ如何なる近似解$\hat{w}$より，弱解$u$に近い．

\begin{align}
    ||| \hat{u} - u |||
    &= \min_{\hat{w} \in \hat{U}} ||| \hat{w} - u ||| \\
\end{align}

とも書ける．

関数空間

Jiao+2024がすっきり纏まっている．

数学が絶望的に苦手なので自信が無いが，自分でも文字に起こして理解の深化を試みる（そういう訳で，間違いを含む可能性が大いにある）．

(2024/12追記)
Evans2010を読め．

Lebesgue空間

Lebesgue空間（$L^p (\Omega)$空間）を次のように定める．

\begin{align}
    L^p (\Omega)
    := \lbrace
        u: \Omega \rightarrow \mathbb{R}
        \mid || u ||_{p} < \infty
    \rbrace
\end{align}

但し，

\begin{align}
    || u ||_{p}
    := \left(
        \int_{\Omega}
            | u |^{p}
        dx
    \right)^{1/p} \quad (1 \le p < \infty)
\end{align}

を$L^p$ノルム，或いは$p$-ノルムなどと呼ぶ（横着して$p = \infty$の場合を省いている）．$L^p$ノルムが有限の値で押さえられるとき，これを（$\Omega$上で）$p$乗可積分な関数と呼ぶ．$L^p (\Omega)$空間は，$p$乗可積分な関数を全て集めてきた集合のことを指す．

多分，厳密には，測度（measure）とか可測（measurable）とかの概念を理解していないといけない．昔，大学院で受講した有限要素法の講義（数学科の講義を受講した）では，

\begin{align}
    L^p (\Omega)
    := \lbrace
        u: \Omega \rightarrow \mathbb{R}
        \mid \text{measurable, } || u ||_{p} < \infty
    \rbrace
\end{align}

但し，

\begin{align}
    || u ||_{p}
    := \left\{
        \begin{array}{ll}
            \left( \int_{\Omega} | u |^{p} dx \right)^{1/p} & (1 \le p < \infty) \\
            \text{ess.sup} \lbrace | u | \mid x \in \Omega \rbrace & (p = \infty) \\
        \end{array}
    \right.
\end{align}

とされていた（$p = \infty$の場合も含めたノルムの定義）．

Sobolev空間

Sobolev空間を次のように定める．

\begin{align}
    W^{k, p} (\Omega)
    := \lbrace
        u: \Omega \rightarrow \mathbb{R}
        \mid D^{\alpha} u \in L^p (\Omega), \ |\alpha| \le k
    \rbrace
\end{align}

但し，$k \in \mathbb{N}_0 := \mathbb{N} \cup \lbrace 0 \rbrace, 1 \le p < \infty$（横着して$p = \infty$の場合を省いている）．特に，$p = 2$のとき，$W^{k, p} (\Omega)$を$H^{k} (\Omega)$と書く．

\begin{align}
    H^{k} (\Omega)
    := W^{k, 2} (\Omega)
    := \lbrace
        u: \Omega \rightarrow \mathbb{R}
        \mid D^{\alpha} u \in L^2 (\Omega), \ |\alpha| \le k
    \rbrace
\end{align}

$W^{k, p} (\Omega)$のことをSobolev空間と呼んだり，$H^{k} (\Omega)$のことをSobolev空間と呼んだりする．或いは，$H^{k} (\Omega)$のことを$k$階のSobolev空間と呼び，区別するようだ．

なお，$\alpha = (\alpha_1, ..., \alpha_d) \in \mathbb{N}_0^d$は多重指数（multi-index）であり，$D^{\alpha}$は以下で定める偏微分を表す．

\begin{align}
    D^{\alpha}
    &:= \frac{\partial^{|\alpha|}}
            {\partial x_1^{\alpha_1} \cdots \partial x_d^{\alpha_d}} \\
\end{align}

ここで，$|\alpha| := \sum_{j=1}^{d} \alpha_j$を$\alpha$の長さと呼ぶ．こちらが分かり易い．

例えば，$|\alpha| = 1$のとき，これを満たすのは$\alpha = (1, 0, 0), (0, 1, 0), (0, 0, 1)$であるから，

\begin{align}
    D^{\alpha}
    &= \frac{\partial}{\partial x_1}, 
        \frac{\partial}{\partial x_2}, 
        \frac{\partial}{\partial x_3} \\
\end{align}

$|\alpha| = 2$のとき，

\begin{align}
    D^{\alpha}
    &= \frac{\partial^2}{\partial x_1^2}, 
        \frac{\partial^2}{\partial x_2^2}, 
        \frac{\partial^2}{\partial x_3^2},
        \frac{\partial^2}{\partial x_1 \partial x_2}, 
        \frac{\partial^2}{\partial x_2 \partial x_3}, 
        \frac{\partial^2}{\partial x_3 \partial x_1}
\end{align}

有限要素法の誤差解析では，$H^1$ノルムを見る場合がある．その定義を記しておきたい．準備として，$L^p$ノルム（特に，$p \in [ 0, \infty )$）を，

\begin{align}
    || u ||_{L^{p}(\Omega)}
    := \left(
        \int_{\Omega}
            | u |^{p}
        dx
    \right)^{1/p}
\end{align}

と書くことにする（記号を改めただけ）．このとき，$H^1$ノルムは，

\begin{align}
    || u ||_{H^{1}(\Omega)}
    &:= \left(
        \sum_{|\alpha| \le 1}^{}
            || D^{\alpha} u ||_{L^{2}(\Omega)}^2
    \right)^{1/2} \\
    &= \left(
        || D^{(0, 0, 0)} u ||_{L^{2}(\Omega)}^2
        + || D^{(1, 0, 0)} u ||_{L^{2}(\Omega)}^2
        + || D^{(0, 1, 0)} u ||_{L^{2}(\Omega)}^2
        + || D^{(0, 0, 1)} u ||_{L^{2}(\Omega)}^2
    \right)^{1/2} \\
    &= \left(
        || u ||_{L^{2}(\Omega)}^2
        + \sum_{i=1}^{d}
            \left\| \frac{\partial}{\partial x_i} u \right\|_{L^{2}(\Omega)}^2
    \right)^{1/2} \\
    &= \left(
        \int_{\Omega} \left| u \right|^{2} dx
        + \sum_{i=1}^{d}
            \int_{\Omega} \left| \frac{\partial}{\partial x_i} u \right|^{2} dx
    \right)^{1/2} \\
    &= \left(
        \int_{\Omega} \left| u \right|^{2} dx
        + \int_{\Omega} \sum_{i=1}^{d} \left| \frac{\partial}{\partial x_i} u \right|^{2} dx
    \right)^{1/2} \\
\end{align}

と書ける．冗長だろうが，人生で一度は具体的に書き出してみたかった．一般的には，

\begin{align}
    || u ||_{H^{1}(\Omega)}
    &:= \left(
        || u ||_{L^{2}(\Omega)}^2
        + \sum_{i=1}^{d} 
            \left\| 
                \frac{\partial u}{\partial x_i} 
            \right\|_{L^{2}(\Omega)}^2
    \right)^{1/2} \\
\end{align}

や

\begin{align}
    || u ||_{H^{1}(\Omega)}
    &:= \left(
        || u ||_{L^{2}(\Omega)}^2
        + || \nabla u ||_{L^{2}(\Omega)}^2
    \right)^{1/2} \\
    &= \left(
        \int_{\Omega} \left| u \right|^{2} dx
        + \int_{\Omega} \left| \nabla u \right|^{2} dx
    \right)^{1/2} \\
\end{align}

などと書かれる印象がある．ここで，

\begin{align}
    \left| \nabla u \right|^{2}
    = \sum_{i=1}^{d} \left| \frac{\partial u}{\partial x_i} \right|^2
    = \left| \frac{\partial u}{\partial x_1} \right|^2
        + \left| \frac{\partial u}{\partial x_2} \right|^2
        + \left| \frac{\partial u}{\partial x_3} \right|^2
\end{align}

は，勾配場のノルムの2乗である．

Karniadakis&Sherwin2005（のch2.2.1.3やch2.2.3）では，試験関数を$\hat{u} = \hat{g}_D + \hat{u}_H$と定めている．ここで，$\hat{g}_D$はDirichlet条件を満たす関数（Karniadakis&Sherwin2005は，これを"lifting function"とか"extending function"とか呼んでいる，厳密にこの言葉は登場しないが），$\hat{u}_H$はDirichlet境界で斉次な関数である．結局，本文と同じ意味だが，境界条件に影響を与えない関数$\hat{u}_H$を用意する，という意味ではこの記述の方が明解だと思う． ↩
$\Omega$が曲線や曲面で閉じられている場合，これを多角形や多面体で表現することになる．したがって，形状を全く厳密に表現することは難しく，区分的な直線や平面で近似している，という意味で近似の記号を用いた． ↩
P1という名前は，Polynomial of degree 1を指すらしい．その他，P2要素というものもあり，これはPolynomial of degree 2: 区分的に2次の多項式を指す．これらは，三角形（2D），四面体要素（3D）の話である．一方，Q1，Q2要素もある．これはQuadrilateral（四角形，六面体）の要素内で，1次や2次の多項式で以て試行関数を構成するものである．他に，P1b要素（Polynomial of degree 1 with bubble function）もある．こうした記法の起こりは，Ciarlet+1978らしい．詳しくは，例えば，中山司: 流れ解析のための有限要素法入門や日本計算工学会流れの有限要素法研究委員会: 続・有限要素法による流れのシミュレーションを参照． ↩

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up