「下ごしらえした万有引力の2体問題の運動方程式を Maxima で数値的に解く」の Python & gnuplot 版。
無次元化された運動方程式
「万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ」にまとめたように，系に特徴的な長さと時間スケールで無次元化した変数に対する運動方程式は以下のようになるのであった。
\begin{eqnarray}
\frac{d^2 X}{dT^2} = – \frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{d^2 Y}{dT^2} = – \frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}
初期条件
無次元化された変数に対する初期条件は $T = 0$ で
\begin{eqnarray}
X &=& 1 \\
Y &=& 0\\
V_x &=& \frac{dX}{dT} = 0 \\
V_y &=& \frac{dY}{dT} = \frac{v_{y0}}{v_0} = k \ \ ( 0 < k < \sqrt{2})
\end{eqnarray}
楕円の軌道要素（長半径，離心率，周期）
このとき，数値計算する前にもう，軌道は規格化された長半径が
$$ A \equiv \frac{a}{x_0} = \frac{1}{2 – k^2} $$
離心率が
$$e = |k^2 – 1|$$
の楕円になることがわかってしまう。また，規格化された周期が
$$P \equiv \frac{p}{\tau} = 2 \pi A^{\frac{3}{2}} = \frac{2 \pi}{\left(2 – k^2\right)^{\frac{3}{2}}}$$
となることもわかってしまっているのであった。
じゃあ，そこまでわかっているのなら，なぜわざわざ数値計算するのかというと，それは楕円軌道の解が時間 $t$ の陽関数として与えられているわけではないから。時刻 $t$ のとき，どこにいるかが解析的にわかっていないので，それを知りたいから数値計算することになる。
数値解法用に連立1階微分方程式系に
\begin{eqnarray}
\frac{dX}{dt} &=& V \\
\frac{dY}{dt} &=& W \\
\frac{dV}{dt} &=&
-\frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{dW}{dt} &=&
-\frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}
scipy.integrate.solve_ivp() で解く
初期条件 $k = 1$ で円軌道となることを確認
まずは，$k = 1$ という初期条件で円軌道になることを確認。1周期を $N$ 等分し，$t = P(1) = 2 \pi$ まで計算したときの答えが $x_0 = 1, y_0 = 0$ になっているか。
import scipy.integrate
import numpy as np
import matplotlib.pyplot as plt
from IPython.display import set_matplotlib_formats
%matplotlib inline
set_matplotlib_formats('svg')
# scipy.integrate.solve_ivp() は
# dy/dt = F(t, y) の形を解く
def F(t, y):
X = y[0]
Y = y[1]
V = y[2]
W = y[3]
return [V, W, -X/np.sqrt(X**2+Y**2)**3, -Y/np.sqrt(X**2+Y**2)**3]
# 周期
def P(k):
return 2*np.pi/(np.sqrt(2 - k**2)**3)
N = 36
y_ini = [1, 0, 0, 1]
t_span = [0, P(1)]
t_list = np.linspace(0, P(1), N+1)
# scipy.integrate.solve_ivp() は
# dy/dt = F(t, y) の形を解く
ivp_c = scipy.integrate.solve_ivp(F, t_span, y_ini,
t_eval = t_list,
rtol = 1.e-12,
atol = 1.e-15)
print("X の誤差 ", 1 - ivp_c.y[0, -1])
print("Y の誤差 ", 0 - ivp_c.y[1, -1])
$e = 0.6$ となる初期条件で数値計算
初期条件 $k$ と楕円の離心率 $e$ との関係は $e = |k^2 – 1|$ であったから…
k = np.sqrt(1 + 0.6)
y_ini = [1, 0, 0, k]
t_span = [0, P(k)]
t_list = np.linspace(0, P(k), N+1)
# scipy.integrate.solve_ivp() は
# dy/dt = Fun(t, y) の形を解く
ivp_e = scipy.integrate.solve_ivp(F, t_span, y_ini,
t_eval = t_list,
rtol = 1.e-12,
atol = 1.e-15)
print("X の誤差 ", 1 - ivp_e.y[0, -1])
print("Y の誤差 ", 0 - ivp_e.y[1, -1])
一定時間間隔ごとの位置をグラフに
plt.figure(figsize=(6,6))
g = plt.subplot()
g.set_xlim (-4.5, 1.5)
g.set_ylim (-3, 3)
g.set_aspect('equal')
g.plot(ivp_e.y[0], ivp_e.y[1], 'o', color = "blue", markersize = 4)
g.plot(0, 0, 'o', color = "black", markersize = 6)
# 軸目盛を非表示
plt.xticks([])
plt.yticks([])
plt.savefig("./ivp_e_fig1.svg")
# グラフを表示
plt.show()
数値計算の結果をファイルへ書き出す
matplotlib に不慣れな場合（私です），数値計算の結果をファイルに書き出して，gnuplot などを使ってグラフにすることができる。
まず，$e=0.6$ となる初期条件で数値計算した結果の $x, y$ をファイル xy_e.txt に書き出す。
# ファイルへの書き出し
xy_e = []
for i in range(0, N+1):
xy_e.append([ivp_e.y[0, i], ivp_e.y[1, i]])
# 配列 xy_e をファイル xy_e.txt に書き込みます。
np.savetxt('xy_e.txt', xy_e)
# ファイルの中身を確認するときは...
# ! cat xy_e.txt
Gnuplot を Python Notebook で使う
弘大 JupyterHub には gnuplot kernel もインストールされている。別途 gnuplot kernel の Notebook を開くまでもなく，この Python の Notebook の中で gnuplot kernel を使うことができる。
# Python Notebook の中から gnuplot_kernel を使う場合，
# まず最初に1回，以下のように %load_ext します。
%load_ext gnuplot_kernel
一旦 %load_ext すれば，%%gnuplot から始まるセルの中では gnuplot のコマンドがそのまま使える。
inline の plot オプション（画像フォーマット，画像サイズ，フォントサイズ等）を変更するには，以下のように %gnuplot line magic を使う。
%gnuplot inline svg size 600,600 fixed enhanced font ',16'
計算結果の $x, y$ 座標値はファイル xy_e.txt に保存されている。
gnuplot でこのファイル xy_e.txt を読み込み，グラフにしてみる。
%%gnuplot
set xrange [-4.5:1.5]
set yrange [-3:3]
set zeroaxis
set size ratio 1
unset xtics
unset ytics
set title "一定時間間隔ごとの位置"
plot "xy_e.txt" pt 7 ps 0.6 lc "blue" notitle
%%gnuplot
set output "./ivp_e_fig01.svg"
replot
set output
Python による数値計算の結果 xy_e.txt だけでなく，楕円や原点を表す点も重ねて描いてみる。
焦点の1つを原点とした極座標では，楕円は
$$r(\phi) = \frac{a (1-e^2)}{1 + e \cos\phi} $$
（このままでも set polar で描けそうだが）以下のように媒介変数表示にして，set parametric としてグラフを描いてみる。
\begin{eqnarray}
x(\phi) &=& r(\phi) \cos\phi \\
y(\phi) &=& r(\phi) \sin\phi
\end{eqnarray}
%%gnuplot
# 極座標表示の楕円
r(a, e, phi) = a*(1-e**2)/(1+e*cos(phi))
x(a, e, phi) = r(a, e, phi) * cos(phi)
y(a, e, phi) = r(a, e, phi) * sin(phi)
# Python による数値計算の結果は以下のような楕円になるハズ
e = 0.6
a = 1/(1-e)
# 原点
set label 1 point pt 7 ps 1 lc "black" at 0, 0
set parametric
plot [phi=0:2*pi] x(a, e, phi), y(a, e, phi) lc "red" notitle, \
"xy_e.txt" pt 7 ps 0.6 lc "blue" notitle
%%gnuplot
set output "./ivp_e_fig02.svg"
replot
set output