ケプラーの法則は、惑星の運動を記述する3つの経験法則です。ニュートンの万有引力の法則と運動の法則から、これらすべてを数学的に導出できます。この導出は、古典力学の最も美しい成果のひとつであり、軌道力学の基礎を成しています。
本記事では、万有引力の法則から出発し、ケプラーの3法則を1つずつ数学的に証明します。
本記事の内容
- 万有引力の法則
- 第2法則（面積速度一定）の証明
- 第1法則（楕円軌道）の導出
- 第3法則（周期と半長径の関係）の証明
- Pythonでの楕円軌道の描画
万有引力の法則
$$ \begin{equation} \bm{F} = -\frac{GMm}{r^2}\hat{\bm{r}} \end{equation} $$
$G$: 万有引力定数、$M$: 中心天体の質量、$m$: 軌道上の物体の質量、$r$: 距離。
重力パラメータ $\mu = GM$ を導入すると、運動方程式は：
$$ \ddot{\bm{r}} = -\frac{\mu}{r^2}\hat{\bm{r}} = -\frac{\mu}{r^3}\bm{r} $$
第2法則: 面積速度一定
定理
中心力のもとで運動する物体の 面積速度 $dA/dt$ は一定である。
証明
角運動量 $\bm{L} = m\bm{r} \times \bm{v}$ の時間微分：
$$ \dot{\bm{L}} = m(\dot{\bm{r}} \times \bm{v} + \bm{r} \times \dot{\bm{v}}) = m(\bm{v} \times \bm{v} + \bm{r} \times \bm{a}) $$
$\bm{v} \times \bm{v} = \bm{0}$ であり、中心力では $\bm{a} \parallel \bm{r}$ なので $\bm{r} \times \bm{a} = \bm{0}$。よって：
$$ \dot{\bm{L}} = \bm{0} \implies \bm{L} = \text{const} $$
面積速度は $\frac{dA}{dt} = \frac{1}{2}|\bm{r} \times \bm{v}| = \frac{L}{2m} = \text{const}$。$\square$
第1法則: 楕円軌道の導出
極座標 $(r, \theta)$ での運動方程式から導出します。
極座標での運動方程式
$r$ 方向: $\ddot{r} – r\dot{\theta}^2 = -\frac{\mu}{r^2}$
$\theta$ 方向: $r\ddot{\theta} + 2\dot{r}\dot{\theta} = 0$（角運動量保存 $h = r^2\dot{\theta} = \text{const}$）
ビネの変数変換
$u = 1/r$ とおくと、$\dot{\theta} = hu^2$ より：
$$ \dot{r} = -\frac{h \, du}{d\theta}, \quad \ddot{r} = -h^2u^2\frac{d^2u}{d\theta^2} $$
$r$ 方向の方程式に代入すると ビネの方程式：
$$ \begin{equation} \frac{d^2u}{d\theta^2} + u = \frac{\mu}{h^2} \end{equation} $$
解
これは非斉次の調和振動子方程式で、一般解は：
$$ u = \frac{\mu}{h^2} + C\cos(\theta – \theta_0) $$
$r = 1/u$ に戻すと：
$$ \begin{equation} r = \frac{h^2/\mu}{1 + e\cos(\theta – \theta_0)} = \frac{p}{1 + e\cos\nu} \end{equation} $$
ここで $p = h^2/\mu$ は半通径、$e = Ch^2/\mu$ は離心率、$\nu = \theta – \theta_0$ は真近点角です。
これは 円錐曲線の極座標表現 であり、$0 \leq e < 1$ のとき楕円です。$\square$
楕円の幾何学
$$ a = \frac{p}{1 – e^2}, \quad b = \frac{p}{\sqrt{1 – e^2}}, \quad r_{\text{min}} = a(1-e), \quad r_{\text{max}} = a(1+e) $$
第3法則: 周期と半長径の関係
楕円の面積 $A = \pi ab$ を面積速度 $dA/dt = h/2$ で割ると周期が求まります：
$$ T = \frac{A}{dA/dt} = \frac{\pi ab}{h/2} = \frac{2\pi ab}{h} $$
$b = a\sqrt{1-e^2}$、$p = a(1-e^2) = h^2/\mu$ より $h = \sqrt{\mu a(1-e^2)}$：
$$ T = \frac{2\pi a \cdot a\sqrt{1-e^2}}{\sqrt{\mu a(1-e^2)}} = 2\pi\sqrt{\frac{a^3}{\mu}} $$
$$ \begin{equation} T^2 = \frac{4\pi^2}{\mu} a^3 \end{equation} $$
周期の2乗は半長径の3乗に比例します。$\square$
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
mu = 3.986e14 # 地球 [m^3/s^2]
# 楕円軌道のパラメータ
a = 20000e3 # 半長径 [m]
e = 0.6 # 離心率
p = a * (1 - e**2)
T = 2 * np.pi * np.sqrt(a**3 / mu)
# 真近点角に対する r
nu = np.linspace(0, 2*np.pi, 1000)
r = p / (1 + e * np.cos(nu))
x = r * np.cos(nu)
y = r * np.sin(nu)
fig, axes = plt.subplots(1, 2, figsize=(14, 6))
# 左: 楕円軌道
axes[0].plot(x/1e6, y/1e6, 'b-', linewidth=2)
axes[0].plot(0, 0, 'yo', markersize=15, label='焦点（地球）')
axes[0].plot(a*(1-e)/1e6, 0, 'r^', markersize=10, label='近地点')
axes[0].plot(-a*(1+e)/1e6, 0, 'gv', markersize=10, label='遠地点')
axes[0].set_xlabel('x [×1000 km]', fontsize=12)
axes[0].set_ylabel('y [×1000 km]', fontsize=12)
axes[0].set_title(f'楕円軌道 (a={a/1e3:.0f}km, e={e})', fontsize=14)
axes[0].set_aspect('equal')
axes[0].legend(fontsize=10)
axes[0].grid(True, alpha=0.3)
# 右: 第3法則の検証
a_range = np.linspace(7000e3, 50000e3, 100)
T_range = 2 * np.pi * np.sqrt(a_range**3 / mu)
axes[1].plot(a_range/1e6, T_range/3600, 'b-', linewidth=2)
axes[1].set_xlabel('半長径 a [×1000 km]', fontsize=12)
axes[1].set_ylabel('周期 T [時間]', fontsize=12)
axes[1].set_title(r'ケプラーの第3法則: $T^2 \propto a^3$', fontsize=14)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"周期: {T/3600:.2f} 時間")
print(f"近地点: {a*(1-e)/1e3:.0f} km, 遠地点: {a*(1+e)/1e3:.0f} km")
まとめ
- 万有引力: $F = -GMm/r^2$
- 第2法則: 角運動量保存 → 面積速度一定
- 第1法則: ビネの方程式 → 円錐曲線 → 楕円軌道 $r = p/(1+e\cos\nu)$
- 第3法則: $T^2 = 4\pi^2 a^3/\mu$
次のステップとして、以下の記事も参考にしてください。