中心力場における単一質点の軌道の可視化

最近、V.I.アーノルド「古典力学の数学的方法」を読んでいるのだが、かっこいい図があったので自分で描いてみた。

古典力学の数学的方法

作者:V.I.アーノルド
発売日: 2003/05/28
メディア: 単行本

f:id:estshorter:20210202223752p:plain — 円環の中、至るところ稠密な軌道

エディタのAtomみたいな画像である。以下、どう上の画像を描いたかかの説明。

いま、平面上の中心力場（質点にはたらく力の向きが「原点から質点へ向かう単位ベクトルかその逆ベクトル」であり、力の大きさも原点からの距離のみに依存する場。例えば重力場。）における、質量1の質点の運動を考える。支配方程式は、

$\displaystyle{ \ddot{\boldsymbol{r}} = - \dfrac{\partial U}{\partial \boldsymbol{r}}, \quad U=U(r) }$

である。この場合、原点からの距離 $r$ は、次のポテンシャルエネルギー（有効ポテンシャルエネルギー）をもつ一次元の問題に従う。

$\displaystyle{ V(r) = U(r) + \dfrac{M^2}{2r^2} }$

すなわち、

$\displaystyle{ \ddot{r} = - \dfrac{\partial V}{\partial {r}} = - \dfrac{\partial U}{\partial r} + \dfrac{M^2}{r^3} }$

である。簡単にいうと、 $r$ は、有効ポテンシャルエネルギー $V$ がなす「坂」を滑らかに転がったときの、地面に平行な方向への変位と同様に変化する。

また、今回のケースにおいては、極座標における角度 $\varphi$ と $r$ は次の関係がある。

$\displaystyle{ \dfrac{d\varphi}{dr} = \dfrac{M/r^2}{\sqrt{2(E-V(r))}} }$

ここで、 $E$ は全エネルギー

$\displaystyle{ E=\dfrac{\dot{r}^2}{2} + V(r) = \dfrac{\dot{\boldsymbol{r}}^2}{2} + U(r) }$

である。したがって、ポテンシャル $U$ の形状、 $M$ および、 $E$ を決めれば数値的に軌道 $\boldsymbol{r}$ を求めることができる。なお、 $U$ の形状によっては、軌道を解析に求められる。例えば、重力場（ $U(r)=-Gm/r$ ）では円錐曲線が $\frac{d\varphi}{dr}$ の式より導かれる。

ここまで説明してようやく画像の話に戻れる。当該画像は、

$\displaystyle{ U(r) = \dfrac{r^3}{2}, M=1, E=2 }$

と設定したときの軌道である。以下のPython スクリプトにより、 $\frac{d\varphi}{dr}$ をRunge-Kutta45で積分して描画した。 $\dot{r}=0$ のときに $\frac{d\varphi}{dr}$ の分母が0になってしまいうまく計算できないので多少の誤差はある（ $U(r)=-1/r$ と設定し、楕円曲線からズレていくのを見ると明らかだ）が、解析解は（多分）ないので仕方ないだろう。

from functools import partial

import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp


def U(r):
    return 0.5 * r ** 3


def V(r):
    return U(r) + 0.5 / r ** 2


def dphi(r, E):
    return 1 / r ** 2 * (1 / np.sqrt(2 * (E - V(r))))


r_min = 0.50842208660526406332007805917088124001667877021803
r_max = 1.5286429152095940020185327716896748970412268021517
E = 2

fun = partial(dphi, E=E)
r_eval = np.arange(r_min, r_max, 0.001)
sol = solve_ivp(lambda r_, phi: fun(r_), [r_min, r_max], [0], t_eval=r_eval)
phi = sol.y.flatten()
r = sol.t

ax = plt.subplot(111, projection="polar")
for i in range(26):
    ax.plot(phi, r, color="k")
    r = r[::-1]
    phi = 2 * phi[-1] - phi[::-1]

circle_phi = np.arange(0, 2 * np.pi, 0.01)
ax.plot(circle_phi, r_min * np.ones(circle_phi.shape), color="k")
ax.plot(circle_phi, r_max * np.ones(circle_phi.shape), color="k")
ax.axis("off")
plt.savefig("orbit.png", bbox_inches="tight", dpi=200)
plt.show()

ちなみに、最初の画像の軌道は「ドーナツ」の部分を稠密に満たす。以下は十分に時間がたった時(range(26*6))の、初期時刻からの軌道である。

f:id:estshorter:20210202233919p:plain — 軌道の長期的な発展

われがわログ

最適化アルゴリズムとかプログラミングについて書きたい

中心力場における単一質点の軌道の可視化