ゆるふわ宇宙・航空
二体問題をPythonでシミュレーション｜衛星の軌道を数値積分で描く
2026-09-20 ・ 実践
vis-vivaの式で軌道速度は求まりますが、「軌道の形そのもの」を描くには運動方程式を数値で解く必要があります。この記事ではPythonで、万有引力だけで動く衛星の軌道をゼロからシミュレーションし、楕円が浮かび上がる様子を確かめます。
発想：力から位置を積み上げる
やることはシンプル。「引力から加速度を求め → 速度を更新 → 位置を更新」を細かい時間刻みで繰り返すだけ。これで軌道が描けます。
時間刻みで積み上げる
引力→加速度
各瞬間に計算
速度を更新
加速度×dt
位置を更新
速度×dt
準備
pip install numpy
① 二体問題を数値積分
安定性の良い「シンプレクティック・オイラー法」で解きます。
import numpy as np
GM = 3.986e14
R = 6.371e6
# 初期条件：近地点400km、速度は円軌道より少し速く（楕円にする）
pos = np.array([R + 400e3, 0.0])
v_circ = np.sqrt(GM / pos[0])
vel = np.array([0.0, v_circ * 1.2]) # 円速度の1.2倍 → 楕円軌道
dt = 1.0 # 1秒刻み
trajectory = []
for step in range(12000): # 約3.3時間ぶん
r = np.linalg.norm(pos)
acc = -GM * pos / r**3 # 万有引力による加速度
vel = vel + acc * dt # 速度を先に更新（シンプレクティック）
pos = pos + vel * dt
trajectory.append(pos.copy())
trajectory = np.array(trajectory)
r_all = np.linalg.norm(trajectory, axis=1)
print(f"近地点高度: {(r_all.min()-R)/1000:6.0f} km")
print(f"遠地点高度: {(r_all.max()-R)/1000:6.0f} km")
出力例:
近地点高度: 400 km
遠地点高度: 11040 km
初速を円軌道の1.2倍にしただけで、遠地点が約11040kmまで伸びた楕円軌道が再現できました。式を解いたのではなく、力から一歩ずつ積み上げた結果です。
② エネルギー保存で答え合わせ
物理的に正しければ、全エネルギー（運動＋位置）は一定に保たれるはずです。
def energy(pos, vel):
r = np.linalg.norm(pos)
return 0.5 * np.dot(vel, vel) - GM / r # 単位質量あたり
# 最初と最後のエネルギーを比較
pos0, vel0 = np.array([R+400e3, 0.0]), np.array([0.0, v_circ*1.2])
print("初期エネルギー:", round(energy(pos0, vel0), 1))
print("最終エネルギー:", round(energy(pos, vel), 1))
両者はほぼ一致します。エネルギーが保たれている＝シミュレーションが正しく動いている証拠。数値計算の妥当性は、こうした保存量で必ずチェックします。
積分法の選び方が命
単純なオイラー法だと誤差が溜まって軌道が徐々にずれます（エネルギーが増えていく）。今回の「速度を先に更新する」シンプレクティック法は、長時間でもエネルギーが安定するので軌道計算向き。より高精度にはルンゲ・クッタ法などを使います。
まとめ
- 二体問題は「引力→加速度→速度→位置」を時間刻みで積み上げて解ける
- 初速を変えるだけで円・楕円の軌道が再現できる
- エネルギー保存で計算の正しさを検証する
- 積分法の選択が精度を左右（シンプレクティック法が軌道向き）
もう少し詳しく（背景と理論）
二体問題は、万有引力で引き合う2物体の運動で、Newton が解いてケプラーの法則を理論的に導いた古典力学の金字塔です1。解は円・楕円・放物線・双曲線の円錐曲線になり、角運動量とエネルギーが保存します。数値シミュレーションで注意すべきは積分法の選択——単純なオイラー法は誤差でエネルギーが増減し軌道が崩れますが、シンプレクティック積分法（Leapfrog/Verlet 等）は長時間でもエネルギーをよく保存します2。現実には第三の天体やJ2摂動、抵抗などが効くため三体以上は解析解を持たず、数値計算が必須になります3。
次の一歩 🌸
速度の式はvis-viva方程式、円軌道は軌道速度、軌道遷移はホーマン遷移へどうぞ。