初めに

Explicit RK familyを学んだことがある．

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

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

参考:

Hairer2006: Geometric Numerical Integration（本稿は，本書のch1程度）

MAIN

Harmonic oscillation

数値実験

ヒューリスティックに進む．

ばね定数$k$のばねに繋がれた質量$m$の質点の1次元的な運動を考える．Newton's 2nd lawとHooke's lawは，

\begin{align}
    m \ddot{q}
    = f
    = - k q
\end{align}

を主張している（位置を$q$とした）．運動量$p$を，

\begin{align}
    p
    = m v
    = m \dot{q}
\end{align}

と定めると，

\begin{align}
    \dot{p}
    = m \ddot{q}
    = - k q
\end{align}

であるから，

\begin{align}
    \begin{cases}
        \dot{q} = p / m \\
        \dot{p} = - k q
    \end{cases}
\end{align}

を得る．$\omega^2 = k / m$としておく．

落ち着いて解く．

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

ところで，

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

且つ，

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

より，

\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$であるから，

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

したがって，解は，

\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$（静止した初期状態）を与えると，

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

である．

この系の全エネルギーは，運動エネルギーとポテンシャルの和であり，

\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

これを数値的にシミュレートする．陽的Euler法を採用すると，

\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

エネルギーが時間的に増大し，phaseが閉じない．一方，

\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

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

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

position	Hamiltonian	phase

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

Exp Euler	Symp Euler	Vel Verlet

何故か？

運動方程式を再掲する．

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

解は，

\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法は，

\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法は，

\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ステップの間にどれだけの引き伸ばしがあるかを見る．

\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を参照すれば，陰的Euler法はエネルギーの減少を生ずる．

Lotka–Volterra

phaseが閉じない様子は，過去のLotka–Volterraの結果を思い出させる．

\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本目の式で除すと，変数が分離できる．

\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)$の全微分になっている．

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

積分して，

\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では，

\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では，

\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

きちんと閉じた．しかし，過去の実験で，HeunやRKが閉じた理由まではきちんと理解できていない．意図しない内に何か都合の良いコードになったか？或いは見え辛いだけかもしれない．

終わりに

謎が深まった．喜ばしい．

Appendix

Verlet / Velocity Verlet method

$q \ (t \pm \tau)$を$q \ (t)$周りでTaylor展開する．

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

和から，

\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法と呼ぶ．

ところで，素朴だが，

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

は位置の更新に使えそうである．また，

\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法と同じ考え方で台形近似を行い，

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

を運動量の更新に採用する．即ち，

\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法と呼ぶ．

少し並び替えると，

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