Symplecticな時間積分法（+ Lotka-Volterra）

# 初めに
[Explicit RK familyを学んだことがある](https://qiita.com/ShotaDeguchi/items/639657f65ee136d17add)．

昔，シンプレクティックな数値積分法（symplectic integrator）のことは全く知らなくて（聞いたことすらなくて），初めて名前を知ったのが2023年だったと思う．その後，情報をつまみ食いしていたが，きちんと勉強したことが無い．

現状，理解が乏しく，かなり初歩的なところにいる．したがって，かなり冗長である．

参考:

* [Hairer2006](https://doi.org/10.1007/3-540-30666-8): Geometric Numerical Integration（本稿は，本書のch1程度）

# MAIN
## Harmonic oscillation
### 数値実験
ヒューリスティックに進む．

ばね定数$k$のばねに繋がれた質量$m$の質点の1次元的な運動を考える．Newton's 2nd lawとHooke's lawは，
```math
\begin{align}
    m \ddot{q}
    = f
    = - k q
\end{align}
```
を主張している（位置を$q$とした）．運動量$p$を，
```math
\begin{align}
    p
    = m v
    = m \dot{q}
\end{align}
```
と定めると，
```math
\begin{align}
    \dot{p}
    = m \ddot{q}
    = - k q
\end{align}
```
であるから，
```math
\begin{align}
    \begin{cases}
        \dot{q} = p / m \\
        \dot{p} = - k q
    \end{cases}
\end{align}
```
を得る．$\omega^2 = k / m$としておく．

落ち着いて解く．
```math
\begin{align}
    \frac{d}{d t}
    \begin{pmatrix}
        q \\
        p
    \end{pmatrix}
    &= \begin{pmatrix}
        0   & 1 / m \\
        - k & 0
    \end{pmatrix}
    \begin{pmatrix}
        q \\
        p
    \end{pmatrix}
    \\
    &= M
    \begin{pmatrix}
        q \\
        p
    \end{pmatrix}
    \\
    \therefore
    \begin{pmatrix}
        q \\
        p
    \end{pmatrix}
    &= \exp \left( M t \right)
    \begin{pmatrix}
        q_0 \\
        p_0
    \end{pmatrix}
    \\
    &= \left(
        I + M t + \frac{M^2}{2!} t^2 + \cdots
    \right)
    \begin{pmatrix}
        q_0 \\
        p_0
    \end{pmatrix}
\end{align}
```
ところで，
```math
\begin{align}
    M
    &= \begin{pmatrix}
        0   & 1 / m \\
        - k & 0
    \end{pmatrix}
    \\
    M^2
    &= \begin{pmatrix}
        0   & 1 / m \\
        - k & 0
    \end{pmatrix}
    \begin{pmatrix}
        0   & 1 / m \\
        - k & 0
    \end{pmatrix} = - \omega^2 I
    \\
    M^3
    &= - \omega^2 M
    \\
    M^4
    &= \omega^4 I
    \\
    \vdots
\end{align}
```
且つ，
```math
\begin{align}
    \cos \left( \omega t \right)
    &= 1
        - \frac{\omega^2}{2!} t^2
        + \frac{\omega^4}{4!} t^4
        + \cdots
    \\
    \sin \left( \omega t \right)
    &= \omega t
        - \frac{\omega^3}{3!} t^3
        + \frac{\omega^5}{5!} t^5
        + \cdots
\end{align}
```
より，
```math
\begin{align}
    \exp \left( M t \right)
    &= \left(
        I 
        + M t
        + \frac{M^2}{2!} t^2
        + \frac{M^3}{3!} t^3 
        + \frac{M^4}{4!} t^4 
        + \cdots
    \right)
    \\
    &= \left(
        I 
        + M t
        - \frac{\omega^2}{2!} t^2 I
        - \frac{\omega^2}{3!} t^2 M
        + \frac{\omega^4}{4!} t^4 I
        + \cdots
    \right)
    \\
    &= \cos \left( \omega t \right) I
        + \sin \left( \omega t \right) \frac{M}{\omega}
    \\
    &= \cos \left( \omega t \right) \begin{pmatrix}
        1 & 0 \\
        0 & 1
    \end{pmatrix}
        + \sin \left( \omega t \right) \begin{pmatrix}
        0 & 1 / m \omega \\
        - k / \omega & 0
    \end{pmatrix}
    \\
\end{align}
```
$\omega^2 = k / m$であるから，
```math
\begin{align}
    \exp \left( M t \right)
    &= \cos \left( \omega t \right) \begin{pmatrix}
        1 & 0 \\
        0 & 1
    \end{pmatrix}
        + \sin \left( \omega t \right) \begin{pmatrix}
        0 & 1 / m \omega \\
        - m \omega & 0
    \end{pmatrix}
    \\
    &= \begin{pmatrix}
        \cos \left( \omega t \right) & \sin \left( \omega t \right) / m \omega \\
        - m \omega \sin \left( \omega t \right) & \cos \left( \omega t \right)
    \end{pmatrix}
\end{align}
```
したがって，解は，
```math
\begin{align}
    \begin{pmatrix}
        q \\
        p
    \end{pmatrix}
    &= \exp \left( M t \right)
    \begin{pmatrix}
        q_0 \\
        p_0
    \end{pmatrix}
    \\
    &= \begin{pmatrix}
        \cos \left( \omega t \right) & \sin \left( \omega t \right) / m \omega \\
        - m \omega \sin \left( \omega t \right) & \cos \left( \omega t \right)
    \end{pmatrix}
    \begin{pmatrix}
        q_0 \\
        p_0
    \end{pmatrix}
\end{align}
```
より，$(q_0, p_0)^{\top}$から時計回りに回転する．$p_0 = 0$（静止した初期状態）を与えると，
```math
\begin{align}
    \begin{pmatrix}
        q \\
        p
    \end{pmatrix}
    &= \begin{pmatrix}
        \cos \left( \omega t \right) & \sin \left( \omega t \right) / m \omega \\
        - m \omega \sin \left( \omega t \right) & \cos \left( \omega t \right)
    \end{pmatrix}
    \begin{pmatrix}
        q_0 \\
        0
    \end{pmatrix}
    \\
    &= \begin{pmatrix}
        q_0 \cos \left( \omega t \right) \\
        - m \omega x_0 \sin \left( \omega t \right)
    \end{pmatrix}
\end{align}
```
である．

この系の全エネルギーは，運動エネルギーとポテンシャルの和であり，
```math
\begin{align}
    E
    &= K + U
    \\
    &= \frac{1}{2} m v^2 + \frac{1}{2} k q^2
    \\
    &= \frac{1}{2 m} p^2 + \frac{1}{2} k q^2
     =: H (q, p)
\end{align}
```
これをHamiltonian$H$と呼ぶ．

$m = 1, \ k = 1, \ \omega = 1, \ (q_0, p_0)^{\top} = (1, 0)^{\top}$と定め，解の動きを見る．

| position | Hamiltonian | phase |
|:-:|:-:|:-:|
| ![time_series.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/af5ecc58-a35a-4220-aa39-6e52d4b176d9.png) | ![energy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/1cacc525-f320-9d42-c055-4e3bbfbca43a.png) | ![phase_space.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/70f4c10c-4565-bc4c-f52f-b64fd12fd724.png) |

これを数値的にシミュレートする．陽的Euler法を採用すると，
```math
\begin{alignat}{2}
    \frac{d q}{d t}
    &= \frac{p}{m}
    \quad
    &&\rightarrow
    \quad
    q^{(n+1)}
    = q^{(n)} + \tau \frac{p^{(n)}}{m}
    \\
    \frac{d p}{d t}
    &= - k q
    \quad
    &&\rightarrow
    \quad
    p^{(n+1)}
    = p^{(n)} - \tau k q^{(n)}
\end{alignat}
```
である．$\tau = 0.1$として数値解を作ると，

| position | Hamiltonian | phase |
|:-:|:-:|:-:|
| ![time_series.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/3ae6f817-6cae-02a9-6162-b2062186f84e.png) | ![energy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/c56bb22f-1976-c2e1-952c-a428a6aefc83.png) | ![phase_space.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/3e1b9b81-0eaa-5c55-7301-acddbdd3ce9b.png) |

エネルギーが時間的に増大し，phaseが閉じない．一方，
```math
\begin{alignat}{2}
    q^{(n+1)}
    &= q^{(n)} + \tau \frac{p^{(n)}}{m}
    \quad
    &&\rightarrow
    \quad
    q^{(n+1)}
    = q^{(n)} + \tau \frac{p^{(n)}}{m}
    \\
    p^{(n+1)}
    &= p^{(n)} - \tau k q^{(n)}
    \quad
    &&\rightarrow
    \quad
    p^{(n+1)}
    = p^{(n)} - \tau k q^{(n+1)}
\end{alignat}
```
と改めると（Leap Frogとはいかないが，$q$が$p$を追いかけるものの，ずっと追いつかない），

| position | Hamiltonian | phase |
|:-:|:-:|:-:|
| ![time_series.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/fea15e8e-afa6-b027-50b1-51e49a3e67df.png) | ![energy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/0a1f0ce2-2a0a-bbd9-37b6-0c8943e94c20.png) | ![phase_space.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/3e7096a4-0703-9f41-4877-8d6b7b68a2ef.png) |

きちんと閉じる（但し，左斜め45度方向に長軸を持つ楕円）．どちらも1次精度なのに，大きく異なる．

速度Verlet法は，さらに精度が良い（大雑把に計算したところ，局所的に3次精度程度ありそうなので，大域的に2次精度だと思う）．

| position | Hamiltonian | phase |
|:-:|:-:|:-:|
| ![time_series.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/91376f4b-7180-5865-dd49-9d162a0b116e.png) | ![energy.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/b4c2bc0f-1386-4de5-6dc8-80ccfa04f7ee.png) | ![phase_space.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/0ac826d1-fce5-3b70-27cd-31291bc509c0.png) |

symplecticな方法を採用したとき，エネルギーは一定の範囲内で振動している．また，数値的にも，各手法が1次，1次，2次の速度であることが確認できる．

| Exp Euler | Symp Euler | Vel Verlet |
|:-:|:-:|:-:|
| ![energy_change_Euler.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/026aca68-1148-c417-374a-416c9ec1e52f.png) | ![energy_change_Symplectic.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/be98b505-53c5-a824-aa94-839790490859.png) | ![energy_change_Verlet.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/08e0dee0-6391-ab8a-d594-a8cbe57c38f4.png) |

### 何故か？
運動方程式を再掲する．
```math
\begin{align}
    \frac{d}{d t}
    \begin{pmatrix}
        q \\
        p
    \end{pmatrix}
    &= \begin{pmatrix}
        0   & 1 / m \\
        - k & 0
    \end{pmatrix}
    \begin{pmatrix}
        q \\
        p
    \end{pmatrix}
\end{align}
```
解は，
```math
\begin{align}
    \begin{pmatrix}
        q \\
        p
    \end{pmatrix}
    &= \begin{pmatrix}
        \cos \left( \omega t \right) & \sin \left( \omega t \right) / m \omega \\
        - m \omega \sin \left( \omega t \right) & \cos \left( \omega t \right)
    \end{pmatrix}
    \begin{pmatrix}
        q_0 \\
        p_0
    \end{pmatrix}
     = M_{\text{Solution}}
    \begin{pmatrix}
        q_0 \\
        p_0
    \end{pmatrix}
\end{align}
```
陽的Euler法は，
```math
\begin{align}
    \begin{pmatrix}
        q^{(n+1)} \\
        p^{(n+1)}
    \end{pmatrix}
    &= \begin{pmatrix}
        1   & \tau / m \\
        - \tau k & 1
    \end{pmatrix}
    \begin{pmatrix}
        q^{(n)} \\
        p^{(n)}
    \end{pmatrix}
     = M_{\text{ExpEuler}}
    \begin{pmatrix}
        q^{(n)} \\
        p^{(n)}
    \end{pmatrix}
\end{align}
```
Symplectic Euler法は，
```math
\begin{align}
    \begin{pmatrix}
        q^{(n+1)} \\
        p^{(n+1)}
    \end{pmatrix}
    &= \begin{pmatrix}
        1   & \tau / m \\
        - \tau k & 1 - \tau^2 k / m
    \end{pmatrix}
    \begin{pmatrix}
        q^{(n)} \\
        p^{(n)}
    \end{pmatrix}
     = M_{\text{SympEuler}}
    \begin{pmatrix}
        q^{(n)} \\
        p^{(n)}
    \end{pmatrix}
\end{align}
```
1ステップの間にどれだけの引き伸ばしがあるかを見る．
```math
\begin{align}
    | M_{\text{Solution}} |
    &= \begin{vmatrix}
        \cos \left( \omega t \right) & \sin \left( \omega t \right) / m \omega \\
        - m \omega \sin \left( \omega t \right) & \cos \left( \omega t \right)
    \end{vmatrix}
    = 1
    \\
    | M_{\text{ExpEuler}} |
    &= \begin{vmatrix}
        1   & \tau / m \\
        - \tau k & 1
    \end{vmatrix}
    = 1 + \tau^2 \frac{k}{m}
    = 1 + \tau^2 \omega^2
    \\
    | M_{\text{SympEuler}} |
    &= \begin{vmatrix}
        1   & \tau / m \\
        - \tau k & 1 - \tau^2 k / m
    \end{vmatrix}
    = 1
\end{align}
```
$M_{\text{Solution}}$はarea-preserving map．陽的Euler法は1ステップ当たり$\tau^2 \omega^2$だけのエネルギーの増加を生ずる．反復すると指数的な増加．$(1 + \tau^2 \omega^2)^n$を計算して初期値に乗ずると，既出の図と合致する．Symplectic Eulerは引き伸ばしが無い．

[Hairer2006](https://doi.org/10.1007/3-540-30666-8)を参照すれば，陰的Euler法はエネルギーの減少を生ずる．

## Lotka–Volterra
phaseが閉じない様子は，[過去のLotka–Volterraの結果](https://qiita.com/ShotaDeguchi/items/639657f65ee136d17add#lotkavolterra-system)を思い出させる．
```math
\begin{align}
    \frac{d q}{d t}
    &= \alpha q - \beta q p
     = \left( \alpha - \beta p \right) q
    \\
    \frac{d p}{d t}
    &= - \gamma p + \delta q p
     = - \left( \gamma - \delta q \right) p
\end{align}
```
$q$はprey population（例えば，ウサギ），$p$はpredator population（例えば，キツネ）．

1本目の式を2本目の式で除すと，変数が分離できる．
```math
\begin{align}
    \frac{\gamma - \delta q}{q} \frac{d q}{dt}
        + \frac{\alpha - \beta p}{p} \frac{d p}{d t}
    &= 0
    \\
\end{align}
```
このとき，左辺は，$I = I (q, p)$の全微分になっている．
```math
\begin{align}
    0 =
    \frac{d I}{d t}
    &= \frac{\partial I}{\partial q} \frac{d q}{d t}
        + \frac{\partial I}{\partial p} \frac{d p}{d t}
    \\
    \frac{\partial I}{\partial q}
    &= \frac{\gamma - \delta q}{q} 
    \\
    \frac{\partial I}{\partial p}
    &= \frac{\alpha - \beta p}{p} 
\end{align}
```
積分して，
```math
\begin{align}
    (\text{const.})
    &= I (q, p)
    \\
    &= \int d I
    \\
    &= \int \left(
        \frac{\partial I}{\partial q} d q
            + \frac{\partial I}{\partial p} d p
    \right)
    \\
    &= \gamma \ln (q) - \delta q
        + \alpha \ln (p) - \beta p
    \\
\end{align}
```
$I (q, p)$を系の不変量（invariant）と呼ぶ．

Explicit Eulerでは，
```math
\begin{align}
    q^{(n+1)}
    &= q^{(n)} + \tau \left( \alpha q^{(n)} - \beta q^{(n)} p^{(n)} \right)
    \\
    p^{(n+1)}
    &= p^{(n)} + \tau \left( - \gamma p^{(n)} + \delta q^{(n)} p^{(n)} \right)
\end{align}
```
Symplectic Eulerでは，
```math
\begin{align}
    q^{(n+1)}
    &= q^{(n)} + \tau \left( \alpha q^{(n)} - \beta q^{(n)} p^{(n)} \right)
    \\
    p^{(n+1)}
    &= p^{(n)} + \tau \left( - \gamma p^{(n)} + \delta q^{(n+1)} p^{(n)} \right)
\end{align}
```
パラメータに$(a, b, c, d) = (2/3, 4/3, 1, 1)$，初期値に$(q^{(0)}, p^{(0)}) = (1, 1)$を選び，$\tau = 0.1$でシミュレートした．

| prey | predator | invariant | phase |
|:-:|:-:|:-:|:-:|
| ![population_prey.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/60082669-d38d-9612-5161-ee09bbee800a.png) | ![population_predator.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/25971b4c-4693-c27b-531a-8597998b3a91.png) | ![invariant.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/02eb5bf2-992c-8f2a-71de-b162b41c4c03.png) | ![population_phase.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2529612/4b30929a-f139-3ec7-b258-ea0827328294.png) |

きちんと閉じた．しかし，[過去の実験](https://qiita.com/ShotaDeguchi/items/639657f65ee136d17add)で，HeunやRKが閉じた理由まではきちんと理解できていない．意図しない内に何か都合の良いコードになったか？或いは見え辛いだけかもしれない．

# 終わりに
謎が深まった．喜ばしい．

# Appendix
## Verlet / Velocity Verlet method
$q \ (t \pm \tau)$を$q \ (t)$周りでTaylor展開する．
```math
\begin{align}
    q \ (t \pm \tau)
    &= q \ (t)
        \pm              \tau     \dot{q}   \ (t)
        +   \frac{1}{2!} \tau^2 \ \ddot{q}  \ (t)
        \pm \frac{1}{3!} \tau^3 \ \dddot{q} \ (t)
        +   \mathcal{O} \ (\tau^4)
\end{align}
```
和から，
```math
\begin{align}
    q \ (t + \tau) + q \ (t - \tau)
    &= 2 q \ (t)
        + \tau^2 \ \ddot{q}  \ (t)
        + \mathcal{O} \ (\tau^4)
    \\
    \therefore
    q \ (t + \tau)
    &= 2 q \ (t)
        - q \ (t - \tau)
        + \tau^2 \ \ddot{q}  \ (t)
        + \mathcal{O} \ (\tau^4)
\end{align}
```
をVerlet法と呼ぶ．

ところで，素朴だが，
```math
\begin{align}
    q \ (t + \tau)
    &= q \ (t)
        +              \tau     \dot{q}   \ (t)
        + \frac{1}{2!} \tau^2 \ \ddot{q}  \ (t)
        + \mathcal{O} \ (\tau^3)
\end{align}
```
は位置の更新に使えそうである．また，
```math
\begin{align}
    d p \ (t)
    &= m \ddot{q} \ (t) \; d t
    \\
    p \ (t + \tau) - p \ (t)
    &= \int_{t}^{t + \tau} m \ddot{q} \ (s) \; d s
    \\
\end{align}
```
に対して，2次のHeun法と同じ考え方で台形近似を行い，
```math
\begin{align}
    p \ (t + \tau)
    &= p \ (t) + \int_{t}^{t + \tau} m \ddot{q} \ (s) \; d s
    \\
    &\simeq p \ (t) + \tau \frac{m \ddot{q} \ (t + \tau) + m \ddot{q} \ (t)}{2}
    \\
    &\simeq p \ (t) + \tau \frac{f \left( q \ (t + \tau) \right) + f \left( q \ (t) \right)}{2}
    \\
\end{align}
```
を運動量の更新に採用する．即ち，
```math
\begin{align}
    q \ (t + \tau)
    &\leftarrow q \ (t)
        + \frac{\tau}{m} p \ (t)
        + \frac{\tau^2}{2 m} f \left( q \ (t) \right)
    \\
    p \ (t + \tau)
    &\leftarrow p \ (t) + \tau \frac{f \left( q \ (t + \tau) \right) + f \left( q \ (t) \right)}{2}
    \\
\end{align}
```
を速度Verlet法と呼ぶ．

少し並び替えると，
```math
\begin{align}
    \begin{pmatrix}
        q^{(n+1)} \\
        p^{(n+1)}
    \end{pmatrix}
    &= \begin{pmatrix}
        1 - \tau^2 k / 2 m        & \tau / m \\
        \tau^3 k^2 / 4 m - \tau k & 1 - \tau^2 k / 2 m
    \end{pmatrix}
    \begin{pmatrix}
        q^{(n)} \\
        p^{(n)}
    \end{pmatrix}
     = M_{\text{VelVerlet}}
    \begin{pmatrix}
        q^{(n)} \\
        p^{(n)}
    \end{pmatrix}
\end{align}
```
と書ける．引き伸ばしを計算すると，$| M_{\text{VelVerlet}} | = 1$であり，area-preserving．