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

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

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

以下が勉強になる:

* [菊池1994: 有限要素法の数理 数学的基礎と誤差解析](http://www.baifukan.co.jp/cgi-bin/db/fromjapanbook.pl?ISBN=4-563-03388-X)
* [菊池1999: 有限要素法概説 [新訂版] 理工学における基礎と応用](https://www.saiensu.co.jp/search/?isbn=978-4-7819-0911-0&y=1999)（個人的には，こちらの方が分かり易かった）
* [桂田祐史: 有限要素法への入門](http://nalab.mind.meiji.ac.jp/~mk/lecture/applied-complex-function-2019/members/fem.pdf)
* [Pythonで1次元有限要素法（Poisson方程式）](https://qiita.com/Sunset_Yuhi/items/4c4ddc25609a7619cce0)
* [日本数値流体力学会有限要素法研究委員会: 有限要素法による流れのシミュレーション](https://www.kinokuniya.co.jp/f/dsg-01-9784431707936)
* [日本計算工学会流れの有限要素法研究委員会: 続・有限要素法による流れのシミュレーション](https://www.kinokuniya.co.jp/f/dsg-01-9784431100331)

# Methods
## Poisson方程式
以下のPoisson方程式を考える．
```math
\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}$は既知の関数とする．

## 弱形式化
関数空間
```math
\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$上で積分する．
```math
\begin{align}
    - \int_{\Omega} \nabla^2 u v \; dx
    &= \int_{\Omega} f v \; dx \\
\end{align}
```
Greenの第一恒等式（Gauss-Greenの式，付録参照）から，
```math
\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}
```
より，
```math
\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）と呼ぶ．なお，
```math
\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}
```
という記号を導入しておくと，弱形式は次のように書ける．
```math
\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追記)**
最近になって，やっと，試行関数を取り出す空間を
```math
\begin{align}
    V_{g_D} 
    &:= \left\{ 
        w \in H^1 \left( \Omega \right) 
        \mid 
        \left. w \right|_{\Gamma_D} = g_D
    \right\} \\
\end{align}
```
と書いておいた方が，ずっとすっきりすることが分かった．こうしておくと，試験関数を取り出す空間が
```math
\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}$を用意する．これには，
```math
\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次独立），
```math
\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}
```
を定める[^1]．ここから取り出す試行関数$\hat{u} \in \hat{U}$と試験関数$\hat{v} \in \hat{V}$を用いて弱形式を近似する．
[^1]: [Karniadakis&Sherwin2005](https://doi.org/10.1093/acprof:oso/9780198528692.001.0001)（の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$を用意する，という意味ではこの記述の方が明解だと思う．
```math
\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$に対しても成り立つ．したがって，
```math
\begin{align}
    \langle \hat{u}, \psi_i \rangle 
    &= \left( f, \psi_i \right) 
        + \left[ g_N, \psi_i \right] 
\end{align}
```
を得る．ここで，
```math
\begin{align}
    \hat{u} = \hat{g}_D + \sum_{i=1}^{N} \beta_i \psi_i
\end{align}
```
を代入すると，
```math
\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$元の連立方程式である．即ち，
```math
\begin{align}
    \boldsymbol{A} \boldsymbol{\beta} &= \boldsymbol{b}
\end{align}
```
と書ける．ここで，
```math
\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](https://www.saiensu.co.jp/search/?isbn=978-4-7819-0911-0&y=1999)より．

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

例えば，[菊池1999](https://www.saiensu.co.jp/search/?isbn=978-4-7819-0911-0&y=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次元の問題
```math
\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$を
```math
\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}
```
と定め，
```math
\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}
```
とおくと，弱形式
```math
\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$[^2]），各要素上でGalerkin法を適用する．1次元の問題では，$\bar{\Omega} = [0, 1]$を$K$個の要素（finite element，或いは単にelement）に分割する．各要素$e_k$は，
```math
\begin{align}
    e_k 
    = \left[ x_{k-1}, x_{k} \right] \quad (k = 1, 2, \ldots, K)
\end{align}
```
とする．即ち，$K$個の要素（element）と$K+1$個の節点（node）がある．
[^2]: $\Omega$が曲線や曲面で閉じられている場合，これを多角形や多面体で表現することになる．したがって，形状を全く厳密に表現することは難しく，区分的な直線や平面で近似している，という意味で近似の記号を用いた．
 
各要素$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)$．

具体的な形を書けば，
```math
\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}$を以下のように定める．
```math
\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$を用いて書くと，
```math
\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}$は，
```math
\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要素[^3]）を用いた有限要素解と呼ぶ．
[^3]: 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](https://doi.org/10.1115/1.3424474)らしい．詳しくは，例えば，[中山司: 流れ解析のための有限要素法入門](https://www.utp.or.jp/book/b305778.html)や[日本計算工学会流れの有限要素法研究委員会: 続・有限要素法による流れのシミュレーション](https://www.kinokuniya.co.jp/f/dsg-01-9784431100331)を参照．

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

### 長さ座標
要素$e_k = \[ x_{k-1}, x_k \]$について，以下の1次関数$L_0, L_1$を用意する．
```math
\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$の長さ座標，或いは形状関数と呼ぶ．これは，
```math
\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$ |
|:-:|:-:|
| ![phi.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/cc962b33-46fd-1028-637b-c1eb1c79f930.png) | ![L.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/09eeeaf5-9cf6-9119-65a9-9890299cc2a2.png) |

$\hat{w} \in \hat{W}$を$L_0, L_1$を用いて次のように展開する．
```math
\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)$とした．

弱形式に登場する
```math
\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$で
```math
\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}
```
という記号を導入すると，先に記した弱形式
```math
\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}
```
は，
```math
\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}$）．

さて，
```math
\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}$を具体的に書くと，
```math
\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}$は，
```math
\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$より，
```math
\begin{align}
    g_N \hat{v} (1)
    &= g_N v_K
\end{align}
```
を得る．

### 要素係数行列，要素自由項ベクトル
要素$e_k = \[ x_{k-1}, x_k \]$において，
```math
\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$について，
```math
\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$を計算してみる．
```math
\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}
```
即ち，
```math
\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$も計算してみる．
```math
\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}
```
こちらは，何らかの求積法を用いて計算しておく．例えば，
```math
\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}
```
という線形補間を利用し，
```math
\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)$と変数変換すれば，
```math
\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}
```
故に，
```math
\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]$は，以下のように改められる．
```math
\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}
```
或いは，
```math
\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}$を
```math
\begin{align}
    \boldsymbol{u} &= (u_0, u_1, \ldots, u_K)^{\top} \\
    \boldsymbol{v} &= (v_0, v_1, \ldots, v_K)^{\top}
\end{align}
```
と定める．次に，
```math
\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}
```
とすると，
```math
\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}
```
を得るから，弱形式は
```math
\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}
```
と書ける．ここで，
```math
\begin{align}
    \check{V} := \left\{ \boldsymbol{v} \in \mathbb{R}^{K+1} \mid v_0 = 0 \right\}
\end{align}
```
である．この式は，任意の$\boldsymbol{v} \in \check{V}$に対して
```math
\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}
```
と書ける．ここで，
```math
\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}
```
と書くと，
```math
\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$より，
```math
\begin{align}
    \left(
        \boldsymbol{A}^{*} \boldsymbol{u} - \boldsymbol{b}^{*}\text{の最初の成分以外}
    \right)
    &= 0
\end{align}
```
が要求される．即ち，
```math
\begin{align}
    \boldsymbol{A}^{**} \boldsymbol{u} 
    &= \boldsymbol{b}^{**}
\end{align}
```
ここで，
```math
\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$の正方行列になる．
```math
\begin{align}
    \boldsymbol{u}^{*}
    &:= \boldsymbol{u} \text{の第0成分を除いた} K \text{次元縦ベクトル} \\
    \boldsymbol{A}
    &:= \boldsymbol{A}^{**} \text{の第0列を除いた} K \text{次正方行列} \\
\end{align}
```
を定義する．即ち，
```math
\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}
```
上式を改めて見ると，
```math
\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}
```
より，
```math
\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方程式を考える．
```math
\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$）ならば，
```math
\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}
```
境界条件から，
```math
\begin{align}
    c_1
    &= g_N + 1 \\
    c_2
    &= g_D \\
\end{align}
```
よって，
```math
\begin{align}
    u 
    &= - \frac{1}{2} x^2 + (g_N + 1) x + g_D \\
\end{align}
```
$(g_D, g_N) = (0, 0)$の場合を，要素数$K$を変化させながら解く．ソルバーには，[昔自作した](https://qiita.com/ShotaDeguchi/items/bd10ad3e65922b797d77)CG法を使った（`numpy.linalg.solve`でも同じ結果であった）．

| $K=2$ | $K=4$ | $K=8$ |
|:-:|:-:|:-:|
| ![solution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/faf5cbef-e0c0-aa4c-3ae2-8a038e29eec7.png) | ![solution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/da614c2d-db60-5243-1226-8b385f5893ce.png) | ![solution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/946df81a-95d4-ff82-8108-aadd99c6c59b.png) |

$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)$ |
|:-:|:-:|:-:|
| ![solution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/f3e93943-c7da-a87e-b97f-5f643a0160ff.png) | ![solution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/6f7ed43b-26c4-b099-298f-07e2346e47d6.png) | ![solution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/f498ff59-5479-bdf0-f273-91c8a67737b6.png) |

正しそうだ．

[Galerkin法を学んだとき](https://qiita.com/ShotaDeguchi/items/4bfd94b11adea89d8a85)に見たように，Neumann条件は厳密に満たされていない．が，"柔らかい"意味で満たされている．また，Neumann条件を厳密に満たしていなくても，近似解は良質と言えるだろう．

## （比較的？）簡単でない問題
同じ問題を考える．
```math
\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](https://doi.org/10.1115/1.1436090)）を用いて，解$u$を用意する．以下のように定めておく．
```math
\begin{equation}
    u = \pi + x \sin(2 \pi x)
\end{equation}
```
よって，
```math
\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$ |
|:-:|:-:|
| ![analytic_solution.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/d089bead-ef85-fd88-34cb-ae715080b779.png) | ![source.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/28d8fdf8-0d4c-3b99-743b-674a39a705ea.png) |

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

| $K=8$ | $K=16$ | $K=32$ | Error |
|:-:|:-:|:-:|:-:|
| ![n_elem8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/3c2ed2cd-9729-81ca-92ed-5bd837ae2677.png) | ![n_elem16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/9592959b-f453-80d1-b6c3-de747b5184d4.png) | ![n_elem32.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/eb835e32-459f-ad1a-ed23-0fef0cf5e5fc.png) | ![convergence.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/ef5ea91d-6f7d-3e78-575f-58ef724fd844.png) |

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

因みに，同書の問題（の一つ）は，
```math
\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 |
|:-:|:-:|:-:|:-:|
| ![n_elem4.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/df16aefc-a4c8-6380-8963-dce87d24e814.png) | ![n_elem8.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/5b07b9d5-d58f-d6b9-2014-d0302c9bcb18.png) | ![n_elem16.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/fa96a67e-2dac-6950-c5d5-445b2d2cd689.png) | ![convergence.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/849ce24a-0435-2ca2-1692-421542380488.png) |

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

ところで，要素長を小さくすることを，$h$-refinementと呼ぶ．基底の次数を大きくすることを，$p$-refinementと呼ぶ．$h$-refinementと$p$-refinementを同時に使う場合，$hp$-refinementという呼び方をする（例えば，[松尾+2017](https://cir.nii.ac.jp/crid/1050574036225158272)）．

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

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

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

# Appendix
## Greenの恒等式（Green's identities）
### Greenの第一恒等式
弱形式を得る途中に用いたGreenの第一恒等式を記しておく．

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

このとき，
```math
\begin{align}
    \nabla \cdot \boldsymbol{A} 
    &= \nabla \cdot (v \nabla u) \\
    &= \nabla v \cdot \nabla u + v \nabla^2 u
\end{align}
```
および
```math
\begin{align}
    \boldsymbol{A} \cdot \boldsymbol{n}
    &= (v \nabla u) \cdot \boldsymbol{n} \\
    &= v \nabla u \cdot \boldsymbol{n}
\end{align}
```
より，元のGaussの発散の定理から，
```math
\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$で与える．

このとき，
```math
\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}
```
および
```math
\begin{align}
    \boldsymbol{A} \cdot \boldsymbol{n}
    &= (u \nabla v - v \nabla u) \cdot \boldsymbol{n} \\
\end{align}
```
より，元のGaussの発散の定理から，
```math
\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の第三恒等式
こちらも，本稿には登場しない．記そうか迷ったが，勉強している内に大変面白そうで，かつ大きい内容でありそうだと思ったので，議論の発散を避けて別に纏めることにする．[こちら](https://hooktail.sub.jp/vectoranalysis/GreensSecondTheorem/)と[こちら](https://eman-physics.net/math/differential24.html)が詳しい．

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

Poisson方程式の弱解$u \in U$は，
```math
\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}$は，
```math
\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$より，
```math
\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式の差分を取り，
```math
\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}$なので，
```math
\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}
```
以下のノルムを定義し，
```math
\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||$）で押さえつける．
```math
\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$に近い．
```math
\begin{align}
    ||| \hat{u} - u |||
    &= \min_{\hat{w} \in \hat{U}} ||| \hat{w} - u ||| \\
\end{align}
```
とも書ける．

## 関数空間
[Jiao+2024](https://arxiv.org/abs/2403.19090)がすっきり纏まっている．

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/f0d47101-633d-8072-8839-03e7d1b31057.png)

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

**(2024/12追記)**
[Evans2010](https://bookstore.ams.org/view?ProductCode=GSM/19.R)を読め（somehow, it is available [here](http://home.ustc.edu.cn/~huangty0407/[Evans,%20Lawrence%20Craig]%20Partial%20Differential%20Equation.pdf)）．

### Lebesgue空間
Lebesgue空間（$L^p (\Omega)$空間）を次のように定める．
```math
\begin{align}
    L^p (\Omega)
    := \lbrace
        u: \Omega \rightarrow \mathbb{R}
        \mid || u ||_{p} < \infty
    \rbrace
\end{align}
```
但し，
```math
\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）とかの概念を理解していないといけない．昔，大学院で受講した有限要素法の講義（数学科の講義を受講した）では，
```math
\begin{align}
    L^p (\Omega)
    := \lbrace
        u: \Omega \rightarrow \mathbb{R}
        \mid \text{measurable, } || u ||_{p} < \infty
    \rbrace
\end{align}
```
但し，
```math
\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空間を次のように定める．
```math
\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)$と書く．
```math
\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}$は以下で定める偏微分を表す．
```math
\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$の長さと呼ぶ．[こちら](https://www.youtube.com/watch?v=Hiwth6HsUq0)が分かり易い．

例えば，$|\alpha| = 1$のとき，これを満たすのは$\alpha = (1, 0, 0), (0, 1, 0), (0, 0, 1)$であるから，
```math
\begin{align}
    D^{\alpha}
    &= \frac{\partial}{\partial x_1}, 
        \frac{\partial}{\partial x_2}, 
        \frac{\partial}{\partial x_3} \\
\end{align}
```
$|\alpha| = 2$のとき，
```math
\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 \)$）を，
```math
\begin{align}
    || u ||_{L^{p}(\Omega)}
    := \left(
        \int_{\Omega}
            | u |^{p}
        \; dx
    \right)^{1/p}
\end{align}
```
と書くことにする（記号を改めただけ）．このとき，$H^1$ノルムは，
```math
\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}
```
と書ける．冗長だろうが，人生で一度は具体的に書き出してみたかった．一般的には，
```math
\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}
```
や
```math
\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}
```
などと書かれる印象がある．ここで，
```math
\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乗である．