「万有引力の運動方程式を数値的に解く準備」で下ごしらえした式を Python で数値的に解き，楕円運動のアニメをつくる。
無次元化された運動方程式
「万有引力の運動方程式を数値的に解く準備」にまとめたように，軌道長半径 $a$ と周期 $P$ で無次元化した変数に対する運動方程式は以下のようになるのであった。
\begin{eqnarray}
\frac{d^2 X}{dT^2} = – 4\pi^2\frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{d^2 Y}{dT^2} = – 4\pi^2\frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}
無次元化された初期条件
無次元化された変数に対する初期条件は $T = 0$ で
\begin{eqnarray}
X(0) &=& 1-e \\
Y(0) &=& 0\\
V_x(0) &=& 0 \\
V_y(0) &=& 2\pi\sqrt{\frac{1+e}{1-e}}
\end{eqnarray}
この初期条件で数値計算すると，（規格化された）長半径 $1$，離心率 $e$ の楕円の軌道が得られるハズ。
じゅあ，そこまでわかっているのなら，なぜわざわざ数値計算するのかというと，それは楕円軌道の解が時間 $t$ の陽関数として与えられているわけではないから。時刻 $t$ のとき，どこにいるかが解析的にわかっていないので，それを知りたいから数値計算することになる。
必要なモジュールの import
In [1]:
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
from scipy.integrate import solve_ivp
# 以下はグラフを SVG で Notebook にインライン表示させる設定
%config InlineBackend.figure_formats = ['svg']
数値計算用に連立1階微分方程式系に
\begin{eqnarray}
\frac{dX}{dT} &=& F_1(V_x) = V_x \\
\frac{dY}{dT} &=& F_2(V_y) = V_y \\
\frac{dV_x}{dT} &=& F_3(X, Y)
= -4\pi^2 \frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{dV_y}{dT} &=& F_4(X, Y)
= -4\pi^2 \frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}
In [2]:
# scipy.integrate.solve_ivp() は
# dy/dt = Fun(t, y) の形を解く。変数名は決め打ち。
def Fun(t, y):
# solve.ivp() に渡す関数の引数の順番に注意。t が先
X = y[0]
Y = y[1]
Vx = y[2]
Vy = y[3]
F1 = Vx
F2 = Vy
F3 = -4*np.pi**2 * X/(X**2 + Y**2)**(3/2)
F4 = -4*np.pi**2 * Y/(X**2 + Y**2)**(3/2)
return [F1, F2, F3, F4]
# 離心率 e = 0.6
e = 0.6
t0 = 0
tend = 1
t_span = [t0, tend]
# 初期条件 X0, Y0, Vx0, Vy0
X0 = 1-e
Y0 = 0
Vx0 = 0
Vy0 = 2*np.pi*np.sqrt((1+e)/(1-e))
y_ini = [X0, Y0, Vx0, Vy0]
Ndiv = 36
t_list = np.linspace(t0, tend, Ndiv+1)
sol = solve_ivp(Fun, t_span, y_ini,
t_eval = t_list,
rtol = 1.e-12,
atol = 1.e-14)
ちょうど1周期 $T=t/P = 1$ までいくと，$T=0$ の時の初期値になっているはずだから，数値誤差を確認する。
数値解 sol から X すなわち y[0] を取り出すには，sol.y[0]。最後の項は sol.y[0][-1]。最初の項（初期条件）sol.y[0][0] との差を確認する。
In [3]:
X = sol.y[0]
Y = sol.y[1]
# X の誤差
print(X[-1]-X[0])
# Y の誤差
print(Y[-1]-Y[0])
一定時間間隔ごとの位置をグラフに
In [4]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
plt.scatter(X, Y, s=20, c='blue')
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);
楕円軌道と重ねてグラフに
得られた数値解が確かに楕円軌道になっていることを確認してみる。
軌道長半径 $a$，離心率 $e$ の楕円の式は
\begin{eqnarray}
r &=& \frac{a(1-e^2)}{1 + e \cos\phi} \\
x &=& r \cos\phi \\
y &=& r \sin\phi
\end{eqnarray}
数値解とあわせるため，軌道長半径 $a$ で無次元化した座標は
\begin{eqnarray}
X &=& \frac{r}{a} \cos\phi = \frac{1-e^2}{1 + e \cos\phi} \cos\phi\\
Y &=& \frac{r}{a} \sin\phi = \frac{1-e^2}{1 + e \cos\phi} \sin\phi
\end{eqnarray}
In [5]:
def Xdaen(e, phi):
R = (1-e**2)/(1+e*np.cos(phi))
return R * np.cos(phi)
def Ydaen(e, phi):
R = (1-e**2)/(1+e*np.cos(phi))
return R * np.sin(phi)
数値解と楕円のグラフを重ねて描いてみると…
In [6]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
# 縦横のアスペクト比を等しくして
plt.axes().set_aspect('equal')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
phi = np.linspace(0, 2*np.pi, 361)
plt.plot(Xdaen(e, phi), Ydaen(e, phi),
linewidth = 0.5, c='red', label="楕円")
plt.plot(X, Y, 'ob', markersize = 2, label='数値解')
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.legend();
問題
- 離心率 $e$ を変えて数値計算を行う。
- 一定時間間隔ごとの位置をアニメにする。
- シーン1：軌道が楕円になる様子のアニメ。
- シーン2：楕円上に一定時間間隔ごとを位置を示すアニメ。
- シーン3：一定時間間隔ごとに掃く扇形を塗りつぶすアニメ。
アニメーション作成
滑らかにするため Ndiv は大きめにとって数値積分
In [7]:
# 離心率 e = 0.6
e = 0.6
t0 = 0
tend = 1
t_span = [t0, tend]
# 初期条件 X0, Y0, Vx0, Vy0
X0 = 1-e
Y0 = 0
Vx0 = 0
Vy0 = 2*np.pi*np.sqrt((1+e)/(1-e))
y_ini = [X0, Y0, Vx0, Vy0]
# 滑らかにするため Ndiv は大きめに
Ndiv = 360
t_list = np.linspace(t0, tend, Ndiv+1)
sol = solve_ivp(Fun, t_span, y_ini,
t_eval = t_list,
rtol = 1.e-12,
atol = 1.e-14)
# 後々のことを考えてリストに変換しておく
X = sol.y[0].tolist()
Y = sol.y[1].tolist()
まずは途中までの軌道を描く
数値計算の結果は X と Y に格納されたので，適宜（間引きしたりして）グラフにする。
In [8]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()
# 外枠と目盛を非表示に
ax.axis('off')
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
i = 3
# i*10 番目までの軌道
xi = X[:i*10]
yi = Y[:i*10]
plt.plot(xi, yi)
# i*10 番目の位置に赤丸
plt.plot(X[i*10], Y[i*10], "or", markersize = 3)
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);
FuncAnimation() でアニメーション作成
特定の i の値のところのグラフができたら，当該コードを def func(i): の中にコピぺして FuncAnimation() を呼べば，パラパラアニメをつくることができるんだったよね。
In [9]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
def func(i):
# 前の frame を消す
plt.cla()
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
# i*10 番目までの軌道
xi = X[:i*10]
yi = Y[:i*10]
plt.plot(xi, yi, linewidth = 0.5)
# i*10 番目の位置に赤丸
plt.plot(X[i*10], Y[i*10], "or", markersize = 3)
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.title("ケプラーの第1法則：楕円軌道")
frames = 36
ani = FuncAnimation(fig, func,
interval = 100,
frames = range(frames))
ani.save("anim01.mp4", dpi=300)
ani.save("anim01.gif", dpi=150)
楕円軌道全体に一定時間ごとの位置を描く
In [10]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
# 楕円軌道は全体
plt.plot(X, Y, linewidth = 0.5)
# i*10 番目までの 10 ごとの位置に赤丸
xi = X[:i*10:10]
yi = Y[:i*10:10]
plt.plot(xi, yi, "or", markersize = 3)
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);
FuncAnimation() でアニメーション作成
特定の i の値のところのグラフができたら，当該コードを def func(i): の中にコピぺして FuncAnimation() を呼べば，パラパラアニメをつくることができるんだったよね。
In [11]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
def func(i):
# 前の frame を消す
plt.cla()
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
# 楕円軌道は全体
plt.plot(X, Y, linewidth = 0.5)
# i*10 番目までの 10 ごとの位置に赤丸
xi = X[:i*10+1:10]
yi = Y[:i*10+1:10]
plt.plot(xi, yi, "or", markersize=3)
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.title("ケプラーの第2法則：一定時間間隔ごとの位置")
frames = 36
ani = FuncAnimation(fig, func,
interval = 100,
frames = range(frames))
ani.save("anim02.mp4", dpi=300)
ani.save("anim02.gif", dpi=150)
Δt の間に掃く扇形を塗りつぶす
In [12]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
i = 0
# i 番目と i+1 番目でできる扇形を塗りつぶす
Xougi = [0] + X[i*10:(i+1)*10+1]
Yougi = [0] + Y[i*10:(i+1)*10+1]
plt.fill(Xougi, Yougi, fc = "yellow")
# 原点と i 番目を結ぶ直線
plt.plot([Xougi[0], Xougi[1]],
[Yougi[0], Yougi[1]], 'tab:blue', linewidth = 0.5)
# 原点と i+1 番目を結ぶ直線
plt.plot([Xougi[0], Xougi[-1]],
[Yougi[0], Yougi[-1]], 'tab:blue', linewidth = 0.5)
# 楕円軌道は全体
plt.plot(X, Y, linewidth = 0.5)
# 10 ごとの全ての位置に赤丸
plt.plot(X[::10], Y[::10], "or", markersize = 3)
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6);
FuncAnimation() でアニメーション作成
特定の i の値のところのグラフができたら，当該コードを def func(i): の中にコピぺして FuncAnimation() を呼べば，パラパラアニメをつくることができるんだったよね。
In [13]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
def func(i):
# 前の frame を消す
plt.cla()
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
# i 番目と i+1 番目でできる扇形を塗りつぶす
Xougi = [0] + X[i*10:(i+1)*10+1]
Yougi = [0] + Y[i*10:(i+1)*10+1]
plt.fill(Xougi, Yougi, fc = "yellow")
# 原点と i 番目を結ぶ直線
plt.plot([Xougi[0], Xougi[1]],
[Yougi[0], Yougi[1]], 'tab:blue', linewidth = 0.5)
# 原点と i+1 番目を結ぶ直線
plt.plot([Xougi[0], Xougi[-1]],
[Yougi[0], Yougi[-1]], 'tab:blue', linewidth = 0.5)
# 楕円軌道は全体
plt.plot(X, Y, linewidth = 0.5)
# 10 ごとの全ての位置に赤丸
plt.plot(X[::10], Y[::10], "or", markersize = 3)
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.title("ケプラーの第2法則：面積速度一定")
frames = 36
ani = FuncAnimation(fig, func,
interval = 100,
frames = range(frames))
ani.save("anim03.mp4", dpi=300)
ani.save("anim03.gif", dpi=150)
面積を比較する静止画パート
面積速度が一定であることをじっくり目視で確認するために，2つの面積を塗りつぶした静止画を最後に追加しておこう。
In [14]:
# グラフのサイズの設定
fig = plt.figure(figsize=[6,6])
ax = fig.add_subplot()
# グラフの縦横のアスペクト比を equal に。
ax.set_aspect('equal')
def func(i):
# 前の frame を消す
plt.cla()
# 外枠と目盛を非表示に
ax.axis('off')
# 軸の設定
plt.axis([-2, 1, -1.5, 1.5])
# 近点付近の扇形
Xougi = [0] + X[-11:-1] + X[0:11]
Yougi = [0] + Y[-11:-1] + Y[0:11]
plt.fill(Xougi, Yougi, fc = "yellow")
# 扇形の直線2本
plt.plot([Xougi[0], Xougi[1]],
[Yougi[0], Yougi[1]], 'tab:blue', linewidth = 0.5)
plt.plot([Xougi[0], Xougi[-1]],
[Yougi[0], Yougi[-1]], 'tab:blue', linewidth = 0.5)
# 遠点付近の扇形
Xougi = [0] + X[170:191]
Yougi = [0] + Y[170:191]
plt.fill(Xougi, Yougi, fc = "yellow")
# 扇形の直線2本
plt.plot([Xougi[0], Xougi[1]],
[Yougi[0], Yougi[1]], 'tab:blue', linewidth = 0.5)
plt.plot([Xougi[0], Xougi[-1]],
[Yougi[0], Yougi[-1]], 'tab:blue', linewidth = 0.5)
# 楕円軌道は全体
plt.plot(X, Y, linewidth = 0.5)
# 10 ごとの全ての位置に赤丸
plt.plot(X[::10], Y[::10], "or", markersize = 3)
# x軸 y軸
plt.axhline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.axvline(0, color='black', dashes=(3, 3), linewidth=0.6)
plt.title("ケプラーの第2法則：面積は等しい")
frames = 36
ani = FuncAnimation(fig, func,
interval = 100,
frames = range(frames))
ani.save("anim04.mp4", dpi=300)
ani.save("anim04.gif", dpi=150)
参考：ffmpeg で動画を連結する
gif ファイルを連結する方法がうまくいかず，とりあえずは mp4 ファイルを連結して gif に変換することに。
In [15]:
# mp4 ファイルを連結する
dat = ["file anim01.mp4",
"file anim02.mp4",
"file anim03.mp4",
"file anim04.mp4"]
np.savetxt('input.txt', dat, fmt='%s') # 文字列として書き込む
# すでに（古い）ファイルがある場合は削除
!rm -f outfile1234.mp4
# input.txt の内容にしたがって連結する
!ffmpeg -hide_banner -loglevel error -f concat -i input.txt -c copy outfile1234.mp4
# すでに（古い）ファイルがある場合は削除
!rm -f outfile1234.gif
# mp4 を gif に変換する
!ffmpeg -hide_banner -loglevel error -i outfile1234.mp4 -vf scale=1280:-1 outfile1234.gif