「ケプラーの3法則」「万有引力 F = G m₁m₂/r²」、教科書では文字や式で覚えますが、Python で太陽系を実際に動かしてみると、なぜそうなるのかが直感的に見えてきます。
ケプラーの3法則（おさらい）
- 1第1法則：惑星は太陽を1焦点とする楕円軌道を描く
- 2第2法則：惑星と太陽を結ぶ線分が単位時間に描く面積は一定（面積速度一定）
- 3第3法則：公転周期² ∝ 軌道長半径³
Step 1: 万有引力で2体問題を解く
太陽（質量 M）から距離 r にある惑星（質量 m）には、太陽方向に引力 F = GMm/r² が働く。これを Python で時間刻み dt ごとに更新します。
地球の公転シミュレーション
Python
import numpy as npimport matplotlib.pyplot as pltfrom matplotlib.animation import FuncAnimationfrom IPython.display import HTML
# 単位を「天文単位」「年」「太陽質量」で表記GM = 4 * np.pi**2 # 太陽の引力定数（天文単位/年² で 4π²）
# 地球の初期状態x, y = 1.0, 0.0 # 太陽から1天文単位vx, vy = 0.0, 2 * np.pi # 公転速度（円軌道なら 2π/年）
dt = 0.01
fig, ax = plt.subplots(figsize=(6, 6))ax.set_xlim(-1.5, 1.5); ax.set_ylim(-1.5, 1.5)ax.scatter([0], [0], s=200, color="orange", label="太陽")trace, = ax.plot([], [], "b-", linewidth=0.5)planet, = ax.plot([x], [y], "o", color="blue", markersize=10)
xs, ys = [x], [y]
def step(frame): global x, y, vx, vy r2 = x**2 + y**2 r = np.sqrt(r2) # 万有引力 ax_g = -GM * x / r**3 ay_g = -GM * y / r**3 vx += ax_g * dt; vy += ay_g * dt x += vx * dt; y += vy * dt xs.append(x); ys.append(y) trace.set_data(xs, ys) planet.set_data([x], [y]) return trace, planet
ani = FuncAnimation(fig, step, frames=300, interval=30, blit=True)HTML(ani.to_jshtml())
実行すると
青い点（地球）が太陽の周りを楕円軌道で回ります。1周がだいたい1年（自然単位系で）になっているのが確認できます。
Step 2: 軌道を変えて楕円に
初期速度を変えると、軌道の形が変わります：
| 初期速度 vy | 軌道 | 解釈 |
|---|---|---|
| 2π × 0.7 | つぶれた楕円 | 近日点で速い |
| 2π（円速度） | 円 | 等速円運動 |
| 2π × 1.2 | 伸びた楕円 | 遠日点で遅い |
| 2π × √2 ≈ 8.89 | 放物線 | 脱出速度（戻ってこない） |
| 10以上 | 双曲線 | 完全に逃げる |
Step 3: ケプラーの第3法則を確認
複数惑星で T² ∝ a³ を検証
Python
import numpy as np
# 太陽系の主要惑星のデータplanets = { "水星": {"a": 0.387, "T": 0.241}, "金星": {"a": 0.723, "T": 0.615}, "地球": {"a": 1.000, "T": 1.000}, "火星": {"a": 1.524, "T": 1.881}, "木星": {"a": 5.203, "T": 11.86}, "土星": {"a": 9.537, "T": 29.46},}
import matplotlib.pyplot as pltxs, ys, names = [], [], []for name, p in planets.items(): xs.append(p["a"]**3); ys.append(p["T"]**2) names.append(name)
plt.figure(figsize=(8, 6))plt.loglog(xs, ys, "o-", markersize=10)for x, y, name in zip(xs, ys, names): plt.annotate(name, (x, y), fontsize=11)plt.xlabel("軌道長半径 a の3乗 (天文単位³)")plt.ylabel("公転周期 T の2乗 (年²)")plt.title("ケプラーの第3法則: T² ∝ a³")plt.grid(True)plt.show()
両対数プロットで直線になり、傾きが1（つまり T² = a³）であることが確認できます。これがケプラーの第3法則。
次回予告
EP.05 は波動の世界。音をして、ドミソの和音が3つの周波数に分解できることを実験します。
この記事の感想を教えてください
あなたの 1 クリックで、本当にこの記事は更新されます。「もっと詳しく」「続編希望」が一定数集まった記事は、 ふくふくが 実際に内容を拡充したり続編記事を公開 します。 送信したリアクションはお使いのブラウザに記録され、再カウントされません。