万有引力の法則は、ニュートンが1687年に発表した物理学の基本法則の一つです。この法則は、あらゆる物体間に作用する引力を記述し、リンゴの落下から惑星の公転運動まで統一的に説明します。
万有引力を理解することは、軌道力学、天体力学、衛星工学の基盤となります。本記事では、万有引力の法則からケプラーの3法則を導出し、Pythonで惑星軌道をシミュレーションします。
本記事の内容
- 万有引力の法則の数学的表現
- ケプラーの3法則の導出
- Pythonによる惑星軌道のシミュレーション
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
万有引力の法則
質量 $M$ と $m$ の2つの物体が距離 $r$ だけ離れているとき、両者間に働く引力の大きさは、
$$ F = G\frac{Mm}{r^2} $$
ここで $G = 6.674 \times 10^{-11} \, \text{N}\cdot\text{m}^2/\text{kg}^2$ は万有引力定数です。
ベクトル形式では、物体 $M$ から物体 $m$ への引力は、
$$ \bm{F} = -G\frac{Mm}{r^2}\hat{\bm{r}} $$
$\hat{\bm{r}}$ は $M$ から $m$ へ向かう単位ベクトルで、マイナス符号は引力（引き合う方向）を示します。
ケプラーの第1法則（楕円軌道の法則）の導出
運動方程式の設定
中心天体（質量 $M$）のまわりの天体（質量 $m$）の運動を考えます。極座標 $(r, \theta)$ での運動方程式は、
$$ m\ddot{r} – mr\dot{\theta}^2 = -\frac{GMm}{r^2} $$
$$ m(r\ddot{\theta} + 2\dot{r}\dot{\theta}) = 0 $$
角運動量の保存
第2式より、
$$ \frac{1}{r}\frac{d}{dt}(mr^2\dot{\theta}) = 0 $$
角運動量 $L = mr^2\dot{\theta}$ が保存されます。
$$ L = mr^2\dot{\theta} = \text{const.} $$
軌道方程式の導出
$u = 1/r$ と置き、$\theta$ を独立変数とします。$\dot{r}$ を変換すると、
$$ \dot{r} = \frac{dr}{dt} = \frac{dr}{d\theta}\dot{\theta} = -\frac{L}{m}\frac{du}{d\theta} $$
$$ \ddot{r} = -\frac{L^2u^2}{m^2}\frac{d^2u}{d\theta^2} $$
第1式の運動方程式に代入して整理すると、ビネの方程式が得られます。
$$ \frac{d^2u}{d\theta^2} + u = \frac{GMm^2}{L^2} $$
この微分方程式の一般解は、
$$ u = \frac{GMm^2}{L^2}(1 + e\cos(\theta – \theta_0)) $$
$r = 1/u$ に戻すと、
$$ r = \frac{L^2/(GMm^2)}{1 + e\cos(\theta – \theta_0)} $$
$p = L^2/(GMm^2)$ を半直弦とおくと、
$$ \boxed{r = \frac{p}{1 + e\cos\theta}} $$
これは 離心率 $e$ の円錐曲線の方程式です。$0 \le e < 1$ のとき楕円となり、ケプラーの第1法則が導出されました。
ケプラーの第2法則（面積速度一定の法則）
面積速度 $\dot{S}$ は、微小時間 $dt$ の間に掃く面積 $dS$ の時間変化率です。
$$ dS = \frac{1}{2}r^2 d\theta $$
$$ \dot{S} = \frac{dS}{dt} = \frac{1}{2}r^2\dot{\theta} = \frac{L}{2m} = \text{const.} $$
角運動量 $L$ が保存されるため、面積速度は一定です。これがケプラーの第2法則です。
ケプラーの第3法則（調和の法則）
楕円の面積は $S = \pi ab$（$a$: 長半径、$b$: 短半径）です。周期 $T$ で一周するので、
$$ \pi ab = \dot{S} \cdot T = \frac{L}{2m}T $$
$$ T = \frac{2\pi m ab}{L} $$
楕円の関係式 $b = a\sqrt{1-e^2}$、$p = a(1-e^2) = L^2/(GMm^2)$ を用いると、
$$ \begin{align} T^2 &= \frac{4\pi^2 m^2 a^2 b^2}{L^2} \\ &= \frac{4\pi^2 m^2 a^2 \cdot a^2(1-e^2)}{L^2} \\ &= \frac{4\pi^2 m^2 a^4(1-e^2)}{L^2} \end{align} $$
$a(1-e^2) = L^2/(GMm^2)$ より $L^2 = GMm^2 a(1-e^2)$ を代入すると、
$$ \boxed{T^2 = \frac{4\pi^2}{GM}a^3} $$
周期の2乗は長半径の3乗に比例します。これがケプラーの第3法則です。
Pythonでの実装
万有引力による軌道運動を数値シミュレーションし、楕円軌道とケプラーの第2法則を可視化します。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# 定数（正規化: G*M = 1）
GM = 1.0
# 初期条件（楕円軌道: e=0.6）
r0 = 1.0 # 近点距離
e = 0.6 # 離心率
a = r0 / (1 - e) # 長半径
v0 = np.sqrt(GM * (2 / r0 - 1 / a)) # ビスビバ方程式
# 運動方程式
def equations(t, state):
x, y, vx, vy = state
r = np.sqrt(x**2 + y**2)
ax = -GM * x / r**3
ay = -GM * y / r**3
return [vx, vy, ax, ay]
# 数値積分
T_orbit = 2 * np.pi * a**1.5 / np.sqrt(GM)
state0 = [r0, 0, 0, v0]
sol = solve_ivp(equations, [0, T_orbit], state0, max_step=0.01, rtol=1e-10, atol=1e-12)
# 可視化
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 軌道
axes[0].plot(sol.y[0], sol.y[1], 'b-', linewidth=1.5)
axes[0].plot(0, 0, 'yo', markersize=12, label='Central body')
axes[0].set_xlabel('x')
axes[0].set_ylabel('y')
axes[0].set_title('Elliptical Orbit (e=0.6)')
axes[0].set_aspect('equal')
axes[0].legend()
axes[0].grid(True)
# ケプラーの第2法則: 面積速度の時間変化
x, y = sol.y[0], sol.y[1]
vx, vy = sol.y[2], sol.y[3]
dS_dt = 0.5 * np.abs(x * vy - y * vx) # 面積速度
axes[1].plot(sol.t, dS_dt, 'r-', linewidth=1.5)
axes[1].set_xlabel('Time')
axes[1].set_ylabel('Areal velocity $dS/dt$')
axes[1].set_title("Kepler's 2nd Law: Constant Areal Velocity")
axes[1].grid(True)
plt.tight_layout()
plt.show()
左図では楕円軌道が描かれ、右図では面積速度が時間によらずほぼ一定であることが確認できます。
まとめ
本記事では、万有引力の法則とケプラーの3法則について解説しました。
- 万有引力 $F = GMm/r^2$ から出発し、極座標での運動方程式を立てた
- 第1法則: 角運動量保存とビネの方程式から楕円軌道 $r = p/(1+e\cos\theta)$ を導出した
- 第2法則: 角運動量保存から面積速度一定 $\dot{S} = L/(2m)$ が導かれた
- 第3法則: 楕円面積と面積速度の関係から $T^2 \propto a^3$ を導出した
次のステップとして、以下の記事も参考にしてください。