ニュートンの万有引力の法則は、りんごの落下から惑星の運動までを統一的に記述する物理学の偉大な成果です。ケプラーが観測データから経験的に見出した3法則が、万有引力から数学的に導出できることを示したのがニュートンの最大の功績でした。
本記事の内容
- 万有引力の法則
- ケプラーの第1法則（楕円軌道）の導出
- 第2法則（面積速度一定）の導出
- 第3法則（調和の法則）の導出
- Pythonでの軌道シミュレーション
前提知識
万有引力の法則
質量 $M$ と $m$ の2物体間の引力:
$$ \bm{F} = -\frac{GMm}{r^2}\hat{\bm{r}} $$
重力ポテンシャル:
$$ U(r) = -\frac{GMm}{r} $$
ケプラーの第2法則の導出
万有引力は中心力（$\bm{F} \parallel \hat{\bm{r}}$）なので、トルク $\bm{\tau} = \bm{r} \times \bm{F} = \bm{0}$。
$$ \frac{d\bm{L}}{dt} = \bm{0} \implies \bm{L} = m\bm{r} \times \bm{v} = \text{const} $$
面積速度:
$$ \frac{dA}{dt} = \frac{1}{2}r^2\dot{\theta} = \frac{|\bm{L}|}{2m} = \text{const} $$
これが面積速度一定（ケプラーの第2法則）です。
ケプラーの第1法則の導出
極座標での運動方程式:
$$ m(\ddot{r} – r\dot{\theta}^2) = -\frac{GMm}{r^2} $$
角運動量 $L = mr^2\dot{\theta}$ を用いて $\dot{\theta}$ を消去し、$u = 1/r$ と変数変換すると:
$$ \frac{d^2u}{d\theta^2} + u = \frac{GMm^2}{L^2} $$
この2階線形ODEの一般解は:
$$ u = \frac{1}{r} = \frac{GMm^2}{L^2}(1 + e\cos(\theta – \theta_0)) $$
$$ \boxed{r = \frac{L^2/(GMm^2)}{1 + e\cos\theta} = \frac{p}{1 + e\cos\theta}} $$
これは離心率 $e$ の円錐曲線の極座標表示です。$0 \leq e < 1$ で楕円、$e = 0$ で円、$e = 1$ で放物線、$e > 1$ で双曲線。
ケプラーの第3法則の導出
楕円の面積 $A = \pi ab$ と面積速度から周期を求めます。
$$ T = \frac{A}{dA/dt} = \frac{\pi ab}{L/(2m)} $$
$a$ と $b$ を軌道パラメータで表すと $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 a^3}{GM} \cdot \frac{m^2 a(1-e^2)}{L^2} = \frac{4\pi^2 a^3}{GM} \end{align} $$
$$ \boxed{T^2 = \frac{4\pi^2}{GM}a^3} $$
周期の2乗は半長軸の3乗に比例します。
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
G = 6.674e-11
M_sun = 1.989e30
AU = 1.496e11
def two_body(t, state):
x, y, vx, vy = state
r = np.sqrt(x**2 + y**2)
ax = -G*M_sun*x/r**3
ay = -G*M_sun*y/r**3
return [vx, vy, ax, ay]
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
# 異なる離心率の軌道
for e, color, name in [(0.0, 'blue', '円 e=0'),
(0.2, 'green', '楕円 e=0.2'),
(0.6, 'red', '楕円 e=0.6'),
(0.9, 'purple', '楕円 e=0.9')]:
a = AU
r_p = a*(1-e)
v_p = np.sqrt(G*M_sun*(1+e)/(a*(1-e)))
T = 2*np.pi*np.sqrt(a**3/(G*M_sun))
sol = solve_ivp(two_body, [0, T*1.05], [r_p, 0, 0, v_p],
t_eval=np.linspace(0, T, 500), rtol=1e-10, atol=1e-12)
axes[0].plot(sol.y[0]/AU, sol.y[1]/AU, color=color, label=name)
axes[0].plot(0, 0, 'yo', markersize=12)
axes[0].set_title('異なる離心率の軌道')
axes[0].set_xlabel('x [AU]')
axes[0].set_ylabel('y [AU]')
axes[0].legend()
axes[0].set_aspect('equal')
axes[0].grid(True, alpha=0.3)
# 面積速度一定の可視化
e = 0.6
a = AU
r_p = a*(1-e)
v_p = np.sqrt(G*M_sun*(1+e)/(a*(1-e)))
T = 2*np.pi*np.sqrt(a**3/(G*M_sun))
sol = solve_ivp(two_body, [0, T], [r_p, 0, 0, v_p],
t_eval=np.linspace(0, T, 100), rtol=1e-10, atol=1e-12)
for i in range(0, len(sol.t)-1, 10):
axes[1].fill([0, sol.y[0][i]/AU, sol.y[0][i+1]/AU],
[0, sol.y[1][i]/AU, sol.y[1][i+1]/AU], alpha=0.3)
axes[1].plot(sol.y[0]/AU, sol.y[1]/AU, 'k-', linewidth=1)
axes[1].plot(0, 0, 'yo', markersize=12)
axes[1].set_title('面積速度一定（等時間間隔の掃過面積）')
axes[1].set_aspect('equal')
axes[1].grid(True, alpha=0.3)
# ケプラーの第3法則
a_vals = np.array([0.387, 0.723, 1.0, 1.524, 5.203]) * AU # 水星〜木星
T_vals = 2*np.pi*np.sqrt(a_vals**3/(G*M_sun))
T_years = T_vals / (365.25*24*3600)
a_AU = a_vals / AU
names = ['水星', '金星', '地球', '火星', '木星']
axes[2].loglog(a_AU, T_years, 'ro', markersize=8)
a_range = np.logspace(-0.5, 1, 100)
axes[2].loglog(a_range, a_range**1.5, 'b-', label='$T \\propto a^{3/2}$')
for i, name in enumerate(names):
axes[2].annotate(name, (a_AU[i], T_years[i]), textcoords="offset points",
xytext=(8, 5), fontsize=9)
axes[2].set_title('ケプラーの第3法則')
axes[2].set_xlabel('半長軸 [AU]')
axes[2].set_ylabel('周期 [年]')
axes[2].legend()
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('kepler_laws.png', dpi=150, bbox_inches='tight')
plt.show()
まとめ
- 万有引力 $F = GMm/r^2$ から3法則が全て導出される
- 第1法則: 軌道は楕円（$r = p/(1+e\cos\theta)$）
- 第2法則: 面積速度一定（角運動量保存の帰結）
- 第3法則: $T^2 = 4\pi^2 a^3/(GM)$
次のステップとして、以下の記事も参考にしてください。