万有引力による惑星の軌道(2次曲線)のシミュレーション

はじめに

恒星の周りを周回する惑星の動きは、万有引力を中心力として考えた場合、2次曲線の軌道を描く。これは、エネルギー保存則と角運動量保存則から示される。それによると、無限遠を基準とした惑星の力学的エネルギーが正となるときは、双曲線を描き、0のときは放物線を負のときは楕円を描く。そこで、惑星の速度によってどのような惑星軌道を描くのかをPythonを用いた数値計算でシミュレーションする。結果として、初速度を変えると以下のように軌道が変化した。
(1)双曲線

(2)放物線

(3)楕円

導入

以下の記事を参考にした。

定点$0$から距離$r$だけ離れた物体に働く中心力を$F(r)=-r^{-2}$とすると、$xy$平面における微分方程式は以下のように表される。ただし質量$m=1$とする。

\frac{d^2 x}{d t^2}= -r^{-2} \frac{x}{r},\frac{d^2 y}{d t^2}= -r^{-2} \frac{y}{r},r=\sqrt{x^2+y^2}

これにより、ポテンシャルエネルギー$U(r)$は以下のように表される。ただし、無限遠を基準0とする。

U(r)=-\frac{1}{r}

したがって、力学的エネルギー$E$は以下のように表すことができる。

E=\frac{v^2}{2}-\frac{1}{r}

これにより、

v=\sqrt{2(E+r^{-1})}

となるので、
(1)楕円($E<0$)

0<v<\sqrt{2r^{-1}}

(2)放物線($E=0$)

v=\sqrt{2r^{-1}}

(3)双曲線($E>0$)

v>\sqrt{2r^{-1}}

となる。

プログラム

python Euler_Planet_2.py

import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
import math
#微小時間
h=1.0*10**-5
# 中心力の定数
alpha=-2.0

# 初期値
## 時間
t=0
## 位置
x=2
y=0
## 速度
v_x=0
### 双曲線はv_y=1.2,放物線はv_y=1.0,楕円はv_y=0.6としてシミュレーションを行った。
v_y=1.0

# 計算結果を入れるためのリスト
t_ary=[]
a_x_ary=[]
a_y_ary=[]
v_x_ary=[]
v_y_ary=[]
x_ary=[]
y_ary=[]
V_ary=[]
E_ary=[]
E_0_ary=[]
L_ary=[]
V_0_ary=[]
V_1_ary=[]
V_x=0
V_y=0
V_0=0
while t<=50:
  r=math.sqrt(x**2+y**2)
  a_x=-r**alpha*x/r
  a_y=-r**alpha*y/r
  v_x=v_x+a_x*h
  v_y=v_y+a_y*h
  x=x+v_x*h
  y=y+v_y*h
  t=t+h
  # オイラー法で求めたポテンシャルエネルギー
  #V_0=-math.sqrt(V_x**2+V_y**2)
  #V_0=V_0-np.array([a_x,a_y])@(np.array([v_x,v_y])*h)
  V_0=V_0-(a_x*v_x+a_y*v_y)*h
  # 理論値のポテンシャルエネルギー
  if alpha != -1:
    V_1=r**(alpha+1)/(alpha+1)
  else:
    V_1=np.log(r)
  E=0.5*(v_x**2+v_y**2)+V_1
  E_0=0.5*(v_x**2+v_y**2)+V_0
  L=x*v_y-y*v_x
  t_ary.append(t)
  a_x_ary.append(a_x)
  a_y_ary.append(a_y)
  v_x_ary.append(v_x)
  v_y_ary.append(v_y)
  x_ary.append(x)
  y_ary.append(y)
  V_0_ary.append(V_0)
  V_1_ary.append(V_1)
  E_ary.append(E)
  E_0_ary.append(E_0)
  L_ary.append(L)

## 時間と変位の関係
# plt.plot(t_ary,x_ary,label="x")
# plt.plot(t_ary,y_ary,label="y")
# plt.legend()
# plt.xlabel("t")
# plt.ylabel("x,y")
# plt.savefig("x_t_y_t.png")
# plt.show()

## 位置関係
plt.plot(x_ary,y_ary)
plt.plot(0,0,marker="o",color="red")
plt.xlabel("x")
plt.ylabel("y")
plt.savefig("x_y5_para.png")
plt.show()

# 力学的エネルギーの推移
# plt.plot(t_ary,E_ary)
# plt.xlabel("t")
# plt.ylabel("E")
# plt.savefig("エネルギー2.png")
# plt.show()

## 角運動量の推移
# plt.plot(t_ary,L_ary)
# plt.xlabel("t")
# plt.ylabel("L")
# plt.savefig("角運動量.png")
# plt.show()

## ポテンシャルエネルギーの理論値とオイラー法で求めた値の比較
# plt.plot(V_1_ary,V_0_ary)
# plt.xlabel("V_1:理論値")
# plt.ylabel("V_0:オイラー法")
# plt.savefig("V_0_V_1.png")
# plt.show()

# # 力学的エネルギーの推移
# plt.plot(t_ary,E_0_ary)
# plt.xlabel("t")
# plt.ylabel("E_0")
# plt.savefig("エネルギー_推定値.png")
# plt.show()

結果

結果として初速度の大きさにより以下のような双曲線、放物線、楕円の3種類の2次曲線を描写することができた。ただし、赤点が恒星で青い軌道が惑星軌道である。

(1)双曲線

(2)放物線

(3)楕円

まとめ

今回は、オイラー法を用いた数値解析により万有引力による惑星の軌道をPythonで描写した。結果、初速度によって惑星の力学的エネルギーが変化するため惑星の軌道が変化するということが分かった。具体的には、双曲線、放物線、楕円のうちのどれかになる。

参考文献