最小二乗法とGalerkin法

# 初めに
重み付き残差法（特に，最小二乗法），およびGalerkin法について，自分なりに纏めてみる．

以下を参考にしている:
* [明治大学 桂田先生: 応用数値解析特論 講義ノート](http://nalab.mind.meiji.ac.jp/~mk/lecture/ouyousuuchikaisekitokuron-2023/)
* 竹内 則雄，樫山 和男，寺田 賢二郎: 計算力学（第2版）有限要素法の基礎
* 地盤工学会: 地盤技術者のためのFEMシリーズ 有限要素法を学ぶ

# モデル問題
## 問題設定
厚さ$h$の均一な地盤（単位体積重量$\gamma$，ヤング率$E$）を考える．

地盤は下端境界で固定支持されており，上端境界で分布荷重$q$が作用しているものとする．

このとき，地盤の深さ方向の変位$u(z)$を求める問題を考える．

微小要素における釣り合いを考えると，
```math
\begin{align}
    (- \sigma_{zz}(z) + \sigma_{zz}(z+\Delta z)) \Delta x \Delta y - \gamma \Delta x \Delta y \Delta z &= 0 \\
    \frac{d \sigma_{zz}}{d z} \Delta z - \gamma \Delta z &= 0 \\
    \therefore \frac{d \sigma_{zz}}{d z} &= \gamma
\end{align}
```
を得る．応力$\sigma_{zz}$と変位$u$の関係
```math
\begin{equation}
    \sigma_{zz} = E \varepsilon_{zz} = E \frac{du}{dz}
\end{equation}
```
より，以下の境界値問題を得る．
```math
\begin{alignat}{2}
    E \frac{d^2 u}{d z^2}
    &= \gamma
    \quad &&(0 < z < h) \\
    u &= 0
    \quad &&(z = 0) \\
    E \frac{d u}{d z} &= - q
    \quad &&(z = h) \\
\end{alignat}
```

## 厳密解
まずは，厳密解を求めておく．

元の微分方程式を積分すると，
```math
\begin{align}
    \frac{d^2 u}{d z^2} &= \frac{\gamma}{E} \\
    \frac{du}{dz} &= c_1 + \frac{\gamma}{E} z \\
    u &= c_0 + c_1 z + \frac{\gamma}{2E} z^2
\end{align}
```
として，一般解を得る．境界条件を考慮すると，
```math
\begin{align}
    0 &= u(z=0) = c_0 \\
    - \frac{q}{E} &= \frac{du}{dz}(z=h) = c_1 + \frac{\gamma}{E} h \\
\end{align}
```
より，
```math
\begin{align}
    c_0 &= 0 \\
    c_1 &= - \frac{q + \gamma h}{E}
\end{align}
```
であるから，厳密解は
```math
\begin{equation}
    u = - \frac{q + \gamma h}{E} z + \frac{\gamma}{2E} z^2 
\end{equation}
```
である．

例えば，$z=h$における変位は
```math
\begin{align}
    u(h)
    &= \frac{\gamma}{2E} h^2 - \frac{q + \gamma h}{E} h \\
    &= - \frac{2 q h + \gamma h^2}{2E}
\end{align}
```
であり，負の値であることは，沈下していることを示している．

## 重み付き残差法（最小二乗法）
### 近似解の用意
厳密解$u$の代わりに，近似解$\hat{u}$を用意する．
```math
\begin{align}
    \hat{u}(z; \boldsymbol{\beta})
    &= \sum_{i=0}^{N} \beta_i \psi_i(z) \\
    &= \sum_{i=0}^{N} \beta_i z^i \\
\end{align}
```
ここで，$\boldsymbol{\beta} = \lbrace \beta_i \rbrace_{i=0}^{N}$はパラメータ（未知数）である[^1]．$\psi_i(z)$は既知の基底関数であり，今回は$\psi_i(z) = z^i$としてみる．近似関数$\hat{u}$を決定することは，$N+1$個の未知数$\lbrace \beta_i \rbrace_{i=0}^{N}$を決定することである．
[^1]: これをパラメータと呼ぶのは嫌われる気がするが，本稿は論文でもないので，少し大らかに考え，私が解釈し易い言葉を使うことにした．

まずは，境界条件を満たすために，$\hat{u}$の動ける範囲を絞る．$N \ge 1$とし，$\hat{u}$は少なくとも1次以上の多項式とする．

Dirichlet境界条件$u(z=0) = 0$を$\hat{u}$に課す．
```math
\begin{align}
    0 = \hat{u} 
    &= \sum_{i=0}^{N} \beta_i z^i \\
    &= \beta_0 + \beta_1 \cdot 0 + \beta_2 \cdot 0^2 + \cdots + \beta_N \cdot 0^N \\
    &= \beta_0
\end{align}
```
すなわち，$\beta_0 = 0$である．

続いて，$\hat{u}$の1階の導関数を求めておく．
```math
\begin{align}
    \frac{d \hat{u}}{d z}(z; \boldsymbol{\beta})
    &= \beta_1 + 2 \beta_2 z + 3 \beta_3 z^2 + \cdots + N \beta_N z^{N-1} \\
    &= \sum_{i=1}^{N} i \beta_i z^{i-1}
\end{align}
```

これを以て，Neumann境界条件$d u/d z = -q/E$を$\hat{u}$に課す．
```math
\begin{align}
    - \frac{q}{E} = \frac{d \hat{u}}{d z}
    &= \sum_{i=1}^{N} i \beta_i h^{i-1} \\
    &= \beta_1 + \sum_{i=2}^{N} i \beta_i h^{i-1}
\end{align}
```
すなわち，$\beta_1 = -q/E - \sum_{i=2}^{N} i \beta_i h^{i-1}$である．

以上より，近似解$\hat{u}$は，
```math
\begin{align}
    \hat{u}(z; \boldsymbol{\beta})
    &= \sum_{i=0}^{N} \beta_i z^i \\
    &= \beta_0 + \beta_1 z + \sum_{i=2}^{N} \beta_i z^i \\
    &= \beta_1 z + \sum_{i=2}^{N} \beta_i z^i
\end{align}
```
に改められる．

$N+1$個の未知係数$\lbrace \beta_i \rbrace_{i=0}^{N}$を定める問題から，$N-1$個の未知係数$\lbrace \beta_i \rbrace_{i=2}^{N}$を定める問題に置き換わった．

### 残差
近似解$\hat{u}$を元の微分方程式に代入し，以下の残差$r$を定義する．
```math
\begin{align}
    r := E \frac{d^2 \hat{u}}{d z^2} - \gamma
\end{align}
```
ここで，$\hat{u} = \hat{u}(z; \boldsymbol{\beta})$より，$r = r(z; \boldsymbol{\beta})$である．

また，$\hat{u}$は，真なる解$u$とは異なると考えられるため，残差$r$は非ゼロであろうと予想する．ただし，真なる解$u$が残差ゼロを達成するのだから，近似解$\hat{u}$に対する要求は，残差$r$をなるべく小さくせよ，ということである．

すなわち，残差$r$を最小化するような$\boldsymbol{\beta}$を求めよ，という問題に置き換わった．

### 最小二乗法（2次多項式）
ところで，残差$r$を最小化せよ，という問題は少し都合が悪い．

なぜなら，$r$は$z$に依存する関数であるため，例えば，

* $r(z=0.2; \boldsymbol{\beta}_1) = 0.1$，かつ $r(z=0.8; \boldsymbol{\beta}_1) = 0.2$，
* $r(z=0.2; \boldsymbol{\beta}_2) = 0.2$，かつ $r(z=0.8; \boldsymbol{\beta}_2) = 0.1$，

という，2つの解の候補があるとき，$\boldsymbol{\beta}_1$と$\boldsymbol{\beta}_2$のどちらがを採用すべきか難しい．

そこで，$r$の大きさを測る尺度を用意する．具体的には，以下の二乗和を考える．
```math
\begin{align}
    R(\boldsymbol{\beta}) 
    &:= \int_{\Omega} r(z; \boldsymbol{\beta})^2 dz
\end{align}
```
ここで，$\Omega = (0, h)$である．$z$について積分することで，$z$の次元を押し潰している．得られた$R(\boldsymbol{\beta})$のことを，**飽くまでも本稿では**「積分残差」と呼ぶことにする[^2]．
[^2]: 恐らく，このような名前は存在しないが，区別のためにこう呼ぶことにする．

先ほどは，残差$r$を最小化するような$\boldsymbol{\beta}$を探そう，と記したが，上述の不都合などを解消するために，積分残差$R$を最小化するような$\boldsymbol{\beta}$を探せ，という問題に置き換えた．

さて，近似解の形を
```math
\begin{equation}
    \hat{u}(z; \boldsymbol{\beta})
    = \sum_{i=0}^{N} \beta_i z^i
    = \beta_1 z + \sum_{i=2}^{N} \beta_i z^i
\end{equation}
```
などとぼかしていたが，そろそろ具体的な形（次数）を定めて話を進めることにしよう．

まずは，$N=2$としてみる．すなわち，
```math
\begin{equation}
    \hat{u}(z; \boldsymbol{\beta}) = \beta_1 z + \beta_2 z^2
\end{equation}
```
である（既に$\beta_1$は既知であるため，$\beta_2$のみ定めれば良い）．このとき，残差$r$は
```math
\begin{align}
    r(z; \boldsymbol{\beta})
    &= E \frac{d^2 \hat{u}}{d z^2} - \gamma \\
    &= E \frac{d^2}{d z^2} (\beta_1 z + \beta_2 z^2) - \gamma \\
    &= 2 E \beta_2 - \gamma
\end{align}
```
である．よって，積分残差$R$は
```math
\begin{align}
    R(\boldsymbol{\beta})
    &= \int_{\Omega} r(z; \boldsymbol{\beta})^2 dz \\
    &= \int_{\Omega} (2 E \beta_2 - \gamma)^2 dz \\
    &= (2 E \beta_2 - \gamma)^2 h
\end{align}
```
である．

強調すべきは，$R$は$\beta_2$の2次式となっている点である．すなわち，$R$はパラメータ空間で（$\beta_2$について）綺麗なお椀型の（単峰の）下に凸な関数となっており，$R$を最小化する$\beta_2$は唯一に定まる（「$R$を最小化する$\beta_2$を探す」と「$\partial_{\beta_2} R = 0$なる$\beta_2$を探す」は同値である）．

ヤング率$E$や厚さ$h$は明らかに非ゼロなので，
```math
\begin{align}
    \frac{d R}{d \beta_2} = 0
    \Rightarrow 2 E (2 E \beta_2 - \gamma) h &= 0 \\
    \therefore \beta_2 &= \frac{\gamma}{2E}
\end{align}
```
より，$\beta_1$の具体的な形（$\beta_1 = -q/E - \sum_{i=2}^{N} i \beta_i h^{i-1}$）を思い出し，以下の近似解を得る．
```math
\begin{align}
    \hat{u}(z; \boldsymbol{\beta})
    &= \beta_1 z + \beta_2 z^2 \\
    &= \left( - \frac{q}{E} - 2 \beta_2 h \right) z + \beta_2 z^2 \\
    &= - \frac{q + \gamma h}{E} z + \frac{\gamma}{2E} z^2
\end{align}
```
これは，厳密解$u$と一致する（何ということだろう！）．

### 最小二乗法（3次多項式）
続いて，$N=3$としてみる．すなわち，近似解を以下の3次多項式とする[^3]．
[^3]: 我々は既に厳密解が2次式であることを知っているが，これは偶々であって，現実には事前に厳密解を知るどころか，次数を推定することもできない，と考えると，えいやっと3次式を試すのも自然と考える．

```math
\begin{equation}
    \hat{u}(z; \boldsymbol{\beta}) = \beta_1 z + \beta_2 z^2 + \beta_3 z^3
\end{equation}
```
このとき，残差$r$は
```math
\begin{align}
    r(z; \boldsymbol{\beta})
    &= E \frac{d^2 \hat{u}}{d z^2} - \gamma \\
    &= E \frac{d^2}{d z^2} (\beta_1 z + \beta_2 z^2 + \beta_3 z^3) - \gamma \\
    &= 2 E \beta_2 - \gamma + 6 E \beta_3 z
\end{align}
```
である．積分残差$R$は
```math
\begin{align}
    R(\boldsymbol{\beta})
    &= \int_{\Omega} r(z; \boldsymbol{\beta})^2 dz \\
    &= \int_{\Omega} (2 E \beta_2 - \gamma + 6 E \beta_3 z)^2 dz \\
    &= (2 E \beta_2 - \gamma)^2 h
        + 2 (2 E \beta_2 - \gamma) 6 E \beta_3 \frac{h^2}{2}
        + (6 E \beta_3)^2 \frac{h^3}{3} \\
    &= (2 E \beta_2 - \gamma)^2 h
        + 6 (2 E \beta_2 - \gamma) E \beta_3 h^2
        + 12 E^2 \beta_3^2 h^3
\end{align}
```
であり，やはりパラメータについて2次形式となっているから，$\partial_{\boldsymbol{\beta}} R = 0$を解けば良い．

実際に計算すれば，
```math
\begin{align}
    \frac{d R}{d \beta_2} = 0
    &\Rightarrow
    (2 E \beta_2 - \gamma) + 3 E \beta_3 h = 0 \\
    \frac{d R}{d \beta_3} = 0
    &\Rightarrow
    (2 E \beta_2 - \gamma) + 4 E \beta_3 h = 0
\end{align}
```
より，
```math
\begin{align}
    E
    \begin{bmatrix}
        2 & 3h \\
        2 & 4h
    \end{bmatrix}
    \begin{bmatrix}
        \beta_2 \\
        \beta_3
    \end{bmatrix}
    &=
    \begin{bmatrix}
        \gamma \\
        \gamma
    \end{bmatrix} \\
    \therefore
    \begin{bmatrix}
        \beta_2 \\
        \beta_3
    \end{bmatrix}
    &=
    E^{-1}
    \begin{bmatrix}
        2 & 3h \\
        2 & 4h
    \end{bmatrix}^{-1}
    \begin{bmatrix}
        \gamma \\
        \gamma
    \end{bmatrix} \\
    &=
    \frac{1}{2Eh}
    \begin{bmatrix}
        4h & -3h \\
        -2 & 2
    \end{bmatrix}
    \begin{bmatrix}
        \gamma \\
        \gamma
    \end{bmatrix} \\
    &=
    \frac{1}{2Eh}
    \begin{bmatrix}
        \gamma h \\
        0
    \end{bmatrix} \\
\end{align}
```
である．これは，先ほどの2次多項式の場合と全く同じである（これは統計学の最小二乗法にも登場するNormal Equationではないか，名前が偶然一致しているわけではないのか，うおお）．

### 最小二乗法（1次多項式）
先ほどは，特に理由を付けずに$N=2$や$N=3$としてみたが，$N=1$とするとどうなるだろうか（試行関数の定義$\hat{u} = \beta_1 z + \sum_{i=2}^{N} \beta_i z^i$上，こんなことをして良いのか？という疑問はあるがえいやっと進む）．近似解は，
```math
\begin{equation}
    \hat{u}(z; \boldsymbol{\beta}) = \beta_1 z
\end{equation}
```
である．$\beta_1$は既に定まっており，$\beta_1 = -q/E - \sum_{i=2}^{N} i \beta_i h^{i-1}$であるから，結局，
```math
\begin{equation}
    \hat{u}(z; \boldsymbol{\beta}) = - \frac{q}{E} z
\end{equation}
```
が近似解となろう．これは，$z=0$で$u=0$，$z=h$で$du/dz = -q/E$であり，境界条件を満たす．

ところで，元の微分方程式は2階の微分方程式であることを思い出そう．$d^2u / dz^2 = \gamma / E$において，例えば，右辺が$z$に依存する滑らかな関数のとき，$u$は2回微分して，まだ滑らかな（かつ$z$に依存する）関数でなければ等式が成り立たない．しかし，いま用意した近似解は1次多項式であり，そのラプラシアンは至るところでゼロである．

こうしたことを考えると，今回の近似解は境界条件を満たすものの，何だか違和感が残る（かつ，関数自体の定義の怪しさがある）．

## 弱形式
### 弱形式化
元の境界値問題を再掲する．
```math
\begin{alignat}{2}
    E \frac{d^2 u}{d z^2}
    &= \gamma
    \quad &&(0 < z < h) \\
    u &= 0
    \quad &&(z = 0) \\
    E \frac{d u}{d z} &= - q
    \quad &&(z = h) \\
\end{alignat}
```
ここで，関数空間$V$を
```math
\begin{align}
    V = \{ 
        v \in H^1 (\Omega) \mid 
        v: \bar{\Omega} \to \mathbb{R}, 
        v = 0 \text{ on } \Gamma_D
    \}
\end{align}
```
と定めておく（$\Gamma_D$はDirichlet境界であり，今回の例では$z=0$がこれに当たる）．

微分方程式の両辺に任意の$v \in V$を掛けて$\Omega$上で積分すると，
```math
\begin{align}
    \int_{\Omega} E \frac{d^2 u}{d z^2} v dz
    &= \int_{\Omega} \gamma v dz
\end{align}
```
である．ここで，部分積分法から，
```math
\begin{align}
    \int_{\Omega} E \frac{d^2 u}{d z^2} v dz
    &= E \left. \frac{d u}{d z} v \right|_{\Gamma}
        - \int_{\Omega} E \frac{d u}{d z} \frac{d v}{d z} dz \\
    &= - q v(h)
        - \int_{\Omega} E \frac{d u}{d z} \frac{d v}{d z} dz
\end{align}
```
より，先の積分の式は
```math
\begin{align}
    - \int_{\Omega} E \frac{d u}{d z} \frac{d v}{d z} dz
    &= \int_{\Omega} \gamma v dz + q v(h)
\end{align}
```
となる[^4]．
[^4]: 符号がヘンテコな気がするが，鉛直上向きが正の形式で，上端から鉛直下向きの$q$が作用するためである．今の問題設定ではこれが正しい．

強調すべきは，解$u$の微分の階数が1階減っていることである．すなわち，$u$に課される滑らかさの条件が緩和されている．また，Neumann境界条件は，境界条件として陽に課されるというよりも，方程式の中に取り込まれた形になっている．

このようにして得られた方程式を，元の方程式の弱形式（weak form）と呼ぶ．これに対して，元の方程式は強形式（strong form）と呼ばれる．弱形式の解を，元の方程式の弱解（weak solution）と呼ぶ．

また，この分野では，

* $\hat{u}$のことを試行関数（trial function），あるいは近似関数（approximate solution）
* $v$のことを試験関数（test function）[^5]，あるいは隆起関数（bump function），あるいは重み関数（weight function）

[^5]: [趣味の大学数学: テスト関数とは](https://math-fun.net/20220821/26248/)

と呼ぶ．

### Galerkin法（1次多項式）
弱解の有限次元近似を求める方法に，Galerkin法がある[^6]．
[^6]: Galerkinの読み方は，日本語では，「ガラーキン」とすることが多いが，英語の発音を真似るならば「ガレルキン」，「ガローキン」，「ガロールキン」辺りが近いだろう．

Galerkin法では，幾つかの関数の1次結合で試行関数や試験関数の近似関数を用意する（この意味で，先ほどまで使っていた近似解と似ている）．具体的には，関数空間
```math
\begin{align}
    U_h &= \left\{ \hat{g}_D + \sum_{i=1}^{N} \beta_i \psi_i \mid (\beta_1, \ldots, \beta_N) \in \mathbb{R}^N \right\} \\
    V_h &= \left\{ \sum_{i=1}^{N} \alpha_i \psi_i \mid (\alpha_1, \ldots, \alpha_N) \in \mathbb{R}^N \right\}
\end{align}
```
から，それぞれ試行関数$\hat{u}_h$と試験関数$v_h$を選ぶ．ただし，$\hat{g}_D$はDirichlet境界条件を満たす関数である（今回の例では，$\hat{g}_D = 0$が対応する）．

ここで，$U_h$と$V_h$はどちらも基底関数$\psi_i$で張られている．これがGalerkin法の特徴であり，別々のものを用いても良い（$U_h$に$\psi_i$，$V_h$に$\phi_i$を用いるなど）．実際，別の基底で張る方法が重み付き残差法と呼ばれるものであり，先に記した最小二乗法はその一つである．

さて，いよいよGalerkin法で1次多項式による近似を試みる．基底関数$\psi_i$には，先ほどまでと同様に，$\psi_i = z^i$を用いる．試行関数$\hat{u}_h$と試験関数$v_h$を以下のように選ぶ．
```math
\begin{align}
    \hat{u}_h &= \beta_1 z \\
    v_h &= \alpha_1 z
\end{align}
```
隠れているが，$\beta_0 = 0$がDirichlet条件を満たす関数$\hat{g}_D$である．また，$v_h$は$v_h = 0$ on $\Gamma_D$を満たす．

$\hat{u}_h$と$v_h$を弱形式に代入すると，
```math
\begin{align}
    - \int_{\Omega} E \frac{d \hat{u}_h}{d z} \frac{d v_h}{d z} dz
    &= \int_{\Omega} \gamma v_h dz + q v_h(h) \\
    - \int_{\Omega} E \beta_1 \alpha_1 dz
    &= \int_{\Omega} \gamma \alpha_1 z dz + q \alpha_1 h \\
    - E \beta_1 \alpha_1 h
    &= \frac{1}{2} \gamma \alpha_1 h^2 + q \alpha_1 h \\
    \alpha_1 \left(
        E \beta_1 h + \frac{1}{2} \gamma h^2 + q h
     \right)
    &= 0
\end{align}
```
ここで，$v_h$は任意の関数であるから，上式は$\alpha_1$に依らずに成り立つ．すなわち，
```math
\begin{align}
    E \beta_1 h + \frac{1}{2} \gamma h^2 + q h
    &= 0 \\
    \therefore \beta_1
    &= - \frac{2 q + \gamma h}{2 E}
\end{align}
```
を得る．よって，Galerkin法（1次近似）による近似解は，
```math
\begin{align}
    \hat{u}_h
    &= \beta_1 z \\
    &= - \frac{2 q + \gamma h}{2 E} z
\end{align}
```
である．これは，厳密解とは異なるものの，良く似ている．

実際，$z=h$における変位を計算すると，
```math
\begin{align}
    \hat{u}_h(h)
    &= - \frac{2 q + \gamma h}{2 E} h \\
    &= - \frac{2 q h + \gamma h^2}{2 E}
\end{align}
```
であり，これは厳密解と一致する．

### Galerkin法（2次多項式）
続いて，Galerkin法で2次多項式による近似を考える．
試行関数$\hat{u}_h$と試験関数$v_h$を以下のように定める．
```math
\begin{align}
    \hat{u}_h &= \beta_1 z + \beta_2 z^2 \\
    v_h &= \alpha_1 z + \alpha_2 z^2
\end{align}
```
Dirichlet条件を満たす関数として，$\beta_0 = 0$を選んでいる．

$\hat{u}_h$と$v_h$を弱形式に代入する.
```math
\begin{align}
    - \int_{\Omega} E \frac{d \hat{u}_h}{d z} \frac{d v_h}{d z} dz
    &= \int_{\Omega} \gamma v_h dz + q v_h(h) \\
    - \int_{\Omega} E (\beta_1 + 2 \beta_2 z) (\alpha_1 + 2 \alpha_2 z) dz
    &= \int_{\Omega} \gamma (\alpha_1 z + \alpha_2 z^2) dz + q (\alpha_1 h + \alpha_2 h^2) \\
\end{align}
```
分けて計算する．
```math
\begin{align}
    (\mathrm{LHS})
    &= - \int_{\Omega} E (\beta_1 + 2 \beta_2 z) (\alpha_1 + 2 \alpha_2 z) dz \\
    &= - \int_{\Omega} E 
        \left(
            \beta_1 \alpha_1 
            + 2 (\beta_1 \alpha_2 + \beta_2 \alpha_1) z
            + 4 \beta_2 \alpha_2 z^2
        \right) dz \\
    &= - E 
        \left(
            \beta_1 \alpha_1 h
            + (\beta_1 \alpha_2 + \beta_2 \alpha_1) h^2
            + \frac{4}{3} \beta_2 \alpha_2 h^3
        \right) \\
    (\mathrm{RHS})
    &= \int_{\Omega} \gamma (\alpha_1 z + \alpha_2 z^2) dz 
        + q (\alpha_1 h + \alpha_2 h^2) \\
    &= \gamma \left( \frac{1}{2} \alpha_1 h^2 + \frac{1}{3} \alpha_2 h^3 \right)
        + q (\alpha_1 h + \alpha_2 h^2)
\end{align}
```
よって，$\alpha_1$と$\alpha_2$について整理すると，
```math
\begin{align}
    \begin{bmatrix}
        \alpha_1 \\
        \alpha_2
    \end{bmatrix}^{\top}
    \begin{bmatrix}
        E \beta_1 h
        + E \beta_2 h^2
        + \frac{1}{2} \gamma h^2
        + q h
        \\
        E \beta_1 h^2
        + \frac{4}{3} E \beta_2 h^3
        + \frac{1}{3} \gamma h^3
        + q h^2
    \end{bmatrix}
    = 0
\end{align}
```
である．任意の$\alpha_1$と$\alpha_2$について成り立つから，
```math
\begin{align}
    E h
    \begin{bmatrix}
        1 & h \\
        h & \frac{4}{3} h^2
    \end{bmatrix}
    \begin{bmatrix}
        \beta_1 \\
        \beta_2
    \end{bmatrix}
    &=
    \begin{bmatrix}
        - \frac{1}{2} \gamma h^2 - q h \\
        - \frac{1}{3} \gamma h^3 - q h^2
    \end{bmatrix} \\
\end{align}
```
を解けば，
```math
\begin{align}
    \begin{bmatrix}
        \beta_1 \\
        \beta_2
    \end{bmatrix}
    &=
    \begin{bmatrix}
        - \frac{q + \gamma h}{E} \\
        \frac{\gamma}{2 E}
    \end{bmatrix}
\end{align}
```
より，Galerkin法（2次近似）による近似解は，
```math
\begin{align}
    \hat{u}_h
    &= \beta_1 z + \beta_2 z^2 \\
    &= - \frac{q + \gamma h}{E} z + \frac{\gamma}{2 E} z^2
\end{align}
```
となり，厳密解と一致する．

## 簡単な比較
乱暴だが，パラメータを全て1としたとき（$h = \gamma = E = q = 1$），各種近似解は以下のようになる．

| $\mathrm{Least \ square}$ | $\mathrm{Galerkin}$ |
|:-:|:-:|
| ![LeastSquare.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/10e680ea-1ff6-5e29-58e6-d4e36a7472c4.png) | ![Galerkin.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/db44281a-661d-9294-8352-c31e57480dca.png) |

興味深いことは，最小二乗法の1次近似は，きちんとNeumann条件を満たしている，ということだ．にも関わらず，近似関数全体としてズレが大きい．一方，Galerkin法における1次近似は，Neumann条件をしっかりと満たしている訳ではないものの，これを弱形式の内側に取り込み，近似関数全体としてのズレは比較的小さい．ここから，Galerkin法における1次近似の強力さを感じる（最小二乗法の1次近似のような，関数の定義自体の怪しさも無い）．

# 終わりに
昔から気になっていたが，何となく棲み分けがはっきりしなかった事項を，自分なりに整理した．

今回は（所謂）重み付き残差法の例として，最小二乗法のみを取り上げたが，他にも様々な方法があるため，それらについても勉強して纏めてみたい．また，空間次元をもう少し大きくしたり，Galerkin法とRitz法の関係[^7]，有限要素法の勉強などもしたい．

[^7]: [桂田先生の資料](http://nalab.mind.meiji.ac.jp/~mk/lecture/applied-complex-function-2019/members/fem.pdf)が参考になる，とだけメモしておく．