夜空を見上げると、月は地球の周りを回り、地球は太陽の周りを回っています。では、なぜ月は地球に落ちてこないのでしょうか？ なぜ惑星は太陽の周りを楕円軌道で回るのでしょうか？ これらの問いに答えるのが二体問題（Two-Body Problem）です。
二体問題は、万有引力で相互作用する2つの天体の運動を数学的に記述する問題です。ニュートンが万有引力の法則を提唱した17世紀以来、天体力学の最も基本的な問題として研究されてきました。二体問題を理解すると、以下のような幅広い応用が見えてきます。
- 人工衛星の軌道設計: 地球周回軌道の計算、静止軌道への投入計画
- 惑星間航行: 火星探査機の遷移軌道の設計
- 天体物理学: 連星系の運動、太陽系外惑星の検出
本記事の内容
- 二体問題の定式化と換算質量による1体問題への帰着
- 運動方程式の導出と角運動量保存則
- ビネの公式を用いた軌道方程式の導出
- 円錐曲線解とエネルギー・軌道形状の関係
- Pythonによる楕円・放物線・双曲線軌道の数値積分と可視化
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
二体問題とは — 2天体の運動を1つにまとめる
問題の設定
二体問題を日常の言葉で言い換えると、「お互いに引力で引き合う2つの物体が、時間とともにどう動くか？」という問題です。たとえば、キャッチボールで投げたボールが放物線を描くのは、ボールと地球の二体問題を解いた結果です。ただし地球はボールに比べて圧倒的に重いため、地球の運動は無視できます。
しかし、月と地球のように質量が近い場合はどうでしょうか。月の質量は地球の約 $1/81$ あり、地球も実は月の引力でわずかに動いています。地球と月は共通の重心（地球中心から約4,670 km の位置）の周りを互いに回っているのです。
このように、2天体が互いに影響し合う運動を厳密に扱うのが二体問題です。質量 $m_1$ と $m_2$ の2天体が万有引力のみで相互作用する場合、それぞれの位置ベクトルを $\bm{r}_1$, $\bm{r}_2$ として運動方程式を書くと次のようになります。
$$ m_1 \ddot{\bm{r}}_1 = -\frac{G m_1 m_2}{|\bm{r}_2 – \bm{r}_1|^3}(\bm{r}_1 – \bm{r}_2) $$
$$ m_2 \ddot{\bm{r}}_2 = -\frac{G m_1 m_2}{|\bm{r}_2 – \bm{r}_1|^3}(\bm{r}_2 – \bm{r}_1) $$
ここで $G$ は万有引力定数（$G = 6.674 \times 10^{-11} \, \text{m}^3 \text{kg}^{-1} \text{s}^{-2}$）です。2つの連立微分方程式がありますが、これを直接解くのは大変です。そこで、巧妙な変数変換によって、この問題を「1つの天体の運動」に帰着させます。
重心運動と相対運動への分離
2天体の全運動は、重心の運動と相対運動に分離できます。これは日常の感覚でも理解できます。2人がスケートリンクでロープを引き合っている場面を想像してください。2人の重心は等速直線運動（外力がなければ）し、重心から見ると2人はそれぞれ重心の周りを回転しています。
まず、重心位置 $\bm{R}$ を定義します。
$$ \bm{R} = \frac{m_1 \bm{r}_1 + m_2 \bm{r}_2}{m_1 + m_2} $$
2つの運動方程式を足し合わせると、右辺は作用・反作用の法則（ニュートンの第3法則）により打ち消し合います。
$$ m_1 \ddot{\bm{r}}_1 + m_2 \ddot{\bm{r}}_2 = \bm{0} $$
すなわち
$$ (m_1 + m_2)\ddot{\bm{R}} = \bm{0} $$
これは重心が等速直線運動することを意味します。外力がなければ、2天体の重心は一定速度で直進します。
次に、相対位置ベクトル $\bm{r}$ を定義します。
$$ \bm{r} = \bm{r}_2 – \bm{r}_1 $$
$m_2$ の運動方程式を $m_2$ で割り、$m_1$ の運動方程式を $m_1$ で割って引くと、相対運動の方程式が得られます。$m_1$ の式から
$$ \ddot{\bm{r}}_1 = -\frac{G m_2}{|\bm{r}|^3}\bm{r} $$
$m_2$ の式から
$$ \ddot{\bm{r}}_2 = +\frac{G m_1}{|\bm{r}|^3}\bm{r} \cdot (-1) = -\frac{G m_1}{|\bm{r}|^3}(-\bm{r})= \frac{G m_1}{|\bm{r}|^3}\bm{r} \cdot (-1) $$
整理すると $\ddot{\bm{r}}_2 = \frac{Gm_1}{|\bm{r}|^3}(-\bm{r})$、つまり $\ddot{\bm{r}}_2 = -\frac{Gm_1}{|\bm{r}|^3}\bm{r}$ ではなく、符号に注意が必要です。丁寧にやり直しましょう。
$\bm{r}_2 – \bm{r}_1 = \bm{r}$ と定義したので、$m_2$ の運動方程式は
$$ m_2 \ddot{\bm{r}}_2 = -\frac{G m_1 m_2}{|\bm{r}|^3}\bm{r} $$
$m_1$ の運動方程式は
$$ m_1 \ddot{\bm{r}}_1 = +\frac{G m_1 m_2}{|\bm{r}|^3}\bm{r} $$
それぞれ質量で割ります。
$$ \ddot{\bm{r}}_2 = -\frac{G m_1}{|\bm{r}|^3}\bm{r}, \quad \ddot{\bm{r}}_1 = +\frac{G m_2}{|\bm{r}|^3}\bm{r} $$
$\ddot{\bm{r}} = \ddot{\bm{r}}_2 – \ddot{\bm{r}}_1$ を計算すると
$$ \ddot{\bm{r}} = -\frac{G m_1}{|\bm{r}|^3}\bm{r} – \frac{G m_2}{|\bm{r}|^3}\bm{r} = -\frac{G(m_1 + m_2)}{|\bm{r}|^3}\bm{r} $$
$\mu = G(m_1 + m_2)$ を重力パラメータと定義すると、相対運動の方程式は驚くほどシンプルになります。
$$ \boxed{\ddot{\bm{r}} = -\frac{\mu}{r^3}\bm{r}} $$
ここで $r = |\bm{r}|$ です。この式は、質量 $\mu$（重力パラメータ）の中心天体の周りを質点が運動する1体問題の運動方程式と全く同じ形です。つまり、二体問題は適切な変数変換により、1体問題に帰着できるのです。
換算質量
二体問題を1体問題として解釈する際に登場するもう一つの重要な概念が換算質量（Reduced Mass）です。換算質量 $m_r$ は次のように定義されます。
$$ m_r = \frac{m_1 m_2}{m_1 + m_2} $$
換算質量を使うと、相対運動の運動方程式は「質量 $m_r$ の粒子が、原点に固定された中心力場内を運動する」と解釈できます。
$$ m_r \ddot{\bm{r}} = -\frac{G m_1 m_2}{r^3}\bm{r} $$
人工衛星（質量 $m$）と地球（質量 $M$）の場合、$m \ll M$ なので
$$ m_r = \frac{mM}{m + M} \approx m $$
$$ \mu = G(m + M) \approx GM $$
となり、「地球を固定して衛星だけが運動する」という近似が正当化されます。地球の重力パラメータは $\mu_{\oplus} = GM_{\oplus} = 3.986 \times 10^{14} \, \text{m}^3/\text{s}^2$ です。
ここまでで、二体問題を1体問題に帰着させることができました。次に、この運動方程式を解いて軌道の形状を求めていきましょう。
角運動量の保存と平面運動
角運動量ベクトルの保存
二体問題において、角運動量が保存することは軌道の形状を決定する上で極めて重要です。日常の例でいえば、フィギュアスケーターが腕を縮めると回転速度が上がる現象が角運動量保存の表れです。天体力学でも同様に、衛星が地球に近づくと速度が上がり、遠ざかると速度が下がります。
角運動量ベクトルは次のように定義されます。
$$ \bm{h} = \bm{r} \times \dot{\bm{r}} $$
$\bm{h}$ の時間微分を計算してみましょう。
$$ \dot{\bm{h}} = \dot{\bm{r}} \times \dot{\bm{r}} + \bm{r} \times \ddot{\bm{r}} $$
第1項は同じベクトルの外積なのでゼロです。第2項について、$\ddot{\bm{r}} = -\frac{\mu}{r^3}\bm{r}$ を代入すると
$$ \bm{r} \times \ddot{\bm{r}} = \bm{r} \times \left(-\frac{\mu}{r^3}\bm{r}\right) = -\frac{\mu}{r^3}(\bm{r} \times \bm{r}) = \bm{0} $$
$\bm{r}$ と $\bm{r}$ は平行なので外積はゼロです。したがって
$$ \dot{\bm{h}} = \bm{0} $$
すなわち、角運動量ベクトル $\bm{h}$ は時間によらず一定です。
軌道は平面運動
角運動量が保存することには深い幾何学的意味があります。$\bm{h} = \bm{r} \times \dot{\bm{r}}$ は $\bm{r}$ と $\dot{\bm{r}}$ の両方に垂直なベクトルです。$\bm{h}$ が一定ということは、$\bm{r}$ と $\dot{\bm{r}}$ は常に $\bm{h}$ に垂直な平面内にあるということです。
つまり、二体問題における運動は常に1つの平面内で起こります。3次元の問題が2次元に帰着するのです。これは解析を大幅に簡単にしてくれます。
この軌道面内で極座標 $(r, \theta)$ を導入しましょう。位置ベクトルを極座標で表すと
$$ \bm{r} = r\bm{e}_r $$
速度ベクトルは
$$ \dot{\bm{r}} = \dot{r}\bm{e}_r + r\dot{\theta}\bm{e}_\theta $$
角運動量の大きさ $h$ は
$$ h = |\bm{r} \times \dot{\bm{r}}| = r \cdot r\dot{\theta} = r^2 \dot{\theta} $$
この $h = r^2 \dot{\theta} = \text{const.}$ がケプラーの第2法則（面積速度一定）に対応します。衛星が近地点（$r$ が小さい）にいるとき $\dot{\theta}$ が大きくなり、遠地点（$r$ が大きい）にいるとき $\dot{\theta}$ が小さくなるのです。
角運動量保存により運動が平面内に制約されることがわかりました。次に、この平面内での運動方程式を極座標で書き下し、軌道の形を求めていきます。
極座標における運動方程式
動径方向と接線方向の分解
軌道面内に極座標 $(r, \theta)$ を設定します。中心力場における加速度を極座標で表すと、動径方向（$\bm{e}_r$ 方向）と接線方向（$\bm{e}_\theta$ 方向）にそれぞれ次の成分を持ちます。
$$ \ddot{\bm{r}} = (\ddot{r} – r\dot{\theta}^2)\bm{e}_r + (r\ddot{\theta} + 2\dot{r}\dot{\theta})\bm{e}_\theta $$
万有引力は中心力なので $\bm{e}_r$ 方向のみに作用します。したがって、2つの成分方程式は
動径方向:
$$ \ddot{r} – r\dot{\theta}^2 = -\frac{\mu}{r^2} $$
接線方向:
$$ r\ddot{\theta} + 2\dot{r}\dot{\theta} = 0 $$
接線方向の方程式は角運動量の保存 $h = r^2\dot{\theta} = \text{const.}$ と等価です。実際、$h = r^2\dot{\theta}$ の時間微分を取ると
$$ \dot{h} = 2r\dot{r}\dot{\theta} + r^2\ddot{\theta} = r(2\dot{r}\dot{\theta} + r\ddot{\theta}) = 0 $$
これは接線方向の方程式そのものです。
$\dot{\theta}$ の消去
動径方向の方程式から軌道形状を求めるために、$\dot{\theta} = h/r^2$ を代入して $\dot{\theta}$ を消去します。
$$ \ddot{r} – r \cdot \frac{h^2}{r^4} = -\frac{\mu}{r^2} $$
整理すると
$$ \ddot{r} – \frac{h^2}{r^3} = -\frac{\mu}{r^2} $$
この方程式は $r$ と $t$ の関係を定めますが、このままでは解くのが困難です。そこで、天体力学で伝統的に用いられるビネの変数変換を導入します。
ここまでの議論で、運動方程式を極座標に変換し、角運動量保存を使って変数を減らしました。しかし、$r(t)$ を直接求めるのは難しいので、次に $r(\theta)$ つまり「軌道の形」を直接求める方法を考えます。
ビネの公式と軌道方程式の導出
ビネの変数変換
ビネの公式は、時間変数 $t$ を角度変数 $\theta$ に置き換えることで、軌道の形状 $r(\theta)$ を直接求める手法です。イメージとしては、「いつどこにいるか」ではなく「角度がいくつのときどのくらいの距離にいるか」を問う方向転換です。
新しい変数として $u = 1/r$ を導入します。すると $r = 1/u$ であり、$u$ は $\theta$ の関数として考えます。
まず、$\dot{r}$ を $u$ と $\theta$ で表します。$r = 1/u$ を時間で微分すると
$$ \dot{r} = -\frac{1}{u^2}\dot{u} = -\frac{1}{u^2}\frac{du}{d\theta}\dot{\theta} $$
$\dot{\theta} = hu^2$（$h = r^2\dot{\theta} = (1/u^2) \cdot \dot{\theta}$ より）を代入すると
$$ \dot{r} = -\frac{1}{u^2} \cdot \frac{du}{d\theta} \cdot hu^2 = -h\frac{du}{d\theta} $$
次に $\ddot{r}$ を計算します。$\dot{r} = -h\frac{du}{d\theta}$ を時間で微分すると
$$ \ddot{r} = -h\frac{d^2u}{d\theta^2}\dot{\theta} = -h\frac{d^2u}{d\theta^2} \cdot hu^2 = -h^2u^2\frac{d^2u}{d\theta^2} $$
ビネの運動方程式
これらを動径方向の運動方程式 $\ddot{r} – h^2/r^3 = -\mu/r^2$ に代入します。
$\ddot{r} = -h^2u^2\frac{d^2u}{d\theta^2}$ と $r = 1/u$ を使って
$$ -h^2u^2\frac{d^2u}{d\theta^2} – h^2u^3 = -\mu u^2 $$
左辺を $-h^2u^2$ で括ると
$$ -h^2u^2\left(\frac{d^2u}{d\theta^2} + u\right) = -\mu u^2 $$
両辺を $-h^2u^2$ で割ると（$u \neq 0$、つまり $r \neq \infty$ を仮定）
$$ \boxed{\frac{d^2u}{d\theta^2} + u = \frac{\mu}{h^2}} $$
これがビネの公式（Binet’s equation）です。驚くべきことに、右辺は定数です。これは単振動の方程式に定数項が加わった形であり、解析的に解くことができます。
軌道方程式の一般解
ビネの公式は、自由項 $\mu/h^2$ を持つ2階線形常微分方程式です。一般解は、斉次方程式 $d^2u/d\theta^2 + u = 0$ の一般解（正弦関数）と特殊解（定数 $\mu/h^2$）の和で与えられます。
$$ u(\theta) = \frac{\mu}{h^2} + C\cos(\theta – \omega) $$
ここで $C$ と $\omega$ は初期条件で決まる積分定数です。$\omega$ は近点方向の角度（近点引数）に対応します。$\theta$ の原点を近点方向に取ると $\omega = 0$ とでき
$$ u(\theta) = \frac{\mu}{h^2}(1 + e\cos\theta) $$
ここで $e = Ch^2/\mu$ を離心率（eccentricity）と定義しました。$r = 1/u$ に戻すと
$$ \boxed{r(\theta) = \frac{h^2/\mu}{1 + e\cos\theta} = \frac{p}{1 + e\cos\theta}} $$
ここで $p = h^2/\mu$ を半通径（semi-latus rectum）と呼びます。
この式は円錐曲線の極座標方程式に他なりません。つまり、万有引力の下での2天体の相対軌道は、必ず円錐曲線（円、楕円、放物線、双曲線のいずれか）になることが証明されました。これがケプラーの第1法則の数学的証明です。
軌道の形が円錐曲線であることがわかりました。では、どのような条件で楕円になり、どのような条件で双曲線になるのでしょうか？ それを決めるのが軌道のエネルギーです。
軌道エネルギーと軌道形状の関係
力学的エネルギーの保存
中心力場では、力学的エネルギーが保存します。運動エネルギーとポテンシャルエネルギーの和を計算しましょう。
極座標での速度の2乗は
$$ v^2 = \dot{r}^2 + (r\dot{\theta})^2 $$
力学的エネルギー（単位質量あたり）を $\mathcal{E}$ と書くと
$$ \mathcal{E} = \frac{1}{2}v^2 – \frac{\mu}{r} = \frac{1}{2}(\dot{r}^2 + r^2\dot{\theta}^2) – \frac{\mu}{r} $$
このエネルギーが保存することを示すために、時間微分を取ってゼロになることを確認します。
$$ \dot{\mathcal{E}} = \dot{r}\ddot{r} + r\dot{r}\dot{\theta}^2 + r^2\dot{\theta}\ddot{\theta} + \frac{\mu}{r^2}\dot{r} $$
動径方向の運動方程式 $\ddot{r} = r\dot{\theta}^2 – \mu/r^2$ を代入し、角運動量保存 $r\ddot{\theta} + 2\dot{r}\dot{\theta} = 0$（つまり $r^2\ddot{\theta} = -2r\dot{r}\dot{\theta}$）を用いると、全ての項が打ち消し合い $\dot{\mathcal{E}} = 0$ が得られます。
近点でのエネルギー
エネルギーと離心率の関係を導くため、近点（$\theta = 0$, $r = r_p$）での条件を使います。近点では $\dot{r} = 0$（動径方向速度がゼロ）なので
$$ \mathcal{E} = \frac{1}{2}\left(\frac{h}{r_p}\right)^2 – \frac{\mu}{r_p} $$
軌道方程式から $r_p = p/(1+e) = h^2/[\mu(1+e)]$ なので
$$ \frac{h}{r_p} = \frac{\mu(1+e)}{h} $$
代入して整理すると
$$ \mathcal{E} = \frac{\mu^2}{2h^2}(1+e)^2 – \frac{\mu^2(1+e)}{h^2} $$
$\mu^2/h^2$ を括り出すと
$$ \mathcal{E} = \frac{\mu^2}{2h^2}\left[(1+e)^2 – 2(1+e)\right] $$
$(1+e)^2 – 2(1+e) = (1+e)(1+e-2) = (1+e)(e-1) = e^2 – 1$ を使うと
$$ \boxed{\mathcal{E} = \frac{\mu^2}{2h^2}(e^2 – 1)} $$
エネルギーと軌道の分類
上の関係から、エネルギーの符号が離心率（つまり軌道の形状）を直接決定することがわかります。
| 離心率 $e$ | エネルギー $\mathcal{E}$ | 軌道の形状 | 物理的意味 |
|---|---|---|---|
| $e = 0$ | $\mathcal{E} < 0$ | 円 | 完全な円軌道 |
| $0 < e < 1$ | $\mathcal{E} < 0$ | 楕円 | 束縛軌道（周期運動） |
| $e = 1$ | $\mathcal{E} = 0$ | 放物線 | 脱出の臨界条件 |
| $e > 1$ | $\mathcal{E} > 0$ | 双曲線 | 非束縛軌道（脱出） |
エネルギーが負のとき、天体は中心天体に束縛されて周期的な楕円運動をします。エネルギーがちょうどゼロのとき、天体は無限遠で速度ゼロになる放物線軌道を描きます。エネルギーが正のとき、天体は中心天体の重力圏から脱出する双曲線軌道を描きます。
楕円軌道の場合のビザビバ方程式
楕円軌道（$0 \leq e < 1$）の場合、半長軸 $a$ との関係を導きます。楕円の幾何学から
$$ r_p = a(1-e), \quad r_a = a(1+e) $$
$$ p = a(1-e^2) $$
したがって $h^2 = \mu p = \mu a(1-e^2)$ であり、エネルギーの式に代入すると
$$ \mathcal{E} = \frac{\mu^2(e^2-1)}{2\mu a(1-e^2)} = \frac{\mu(e^2-1)}{2a(1-e^2)} = -\frac{\mu}{2a} $$
$$ \boxed{\mathcal{E} = -\frac{\mu}{2a}} $$
この式を速度の形で書き直すと、有名なビザビバ方程式（vis-viva equation）が得られます。
$$ \boxed{v^2 = \mu\left(\frac{2}{r} – \frac{1}{a}\right)} $$
ビザビバ方程式は、軌道上の任意の位置での速度を、その位置 $r$ と半長軸 $a$ だけから求められる非常に便利な式です。
ここまでの理論的な導出を踏まえて、次にPythonを使って楕円・放物線・双曲線の3種類の軌道を数値的にシミュレーションし、理論の正しさを確認しましょう。
Pythonで軌道方程式を可視化する
円錐曲線軌道の描画
まず、軌道方程式 $r(\theta) = p/(1+e\cos\theta)$ を使って、離心率を変えた場合の軌道形状を可視化します。
import numpy as np
import matplotlib.pyplot as plt
# 半通径を固定して離心率を変える
p = 1.0 # 半通径（正規化）
# 各離心率に対する軌道
eccentricities = [0.0, 0.3, 0.6, 0.9, 1.0, 1.5]
labels = [
f"e = {e} ({'円' if e == 0 else '楕円' if e < 1 else '放物線' if e == 1 else '双曲線'})"
for e in eccentricities
]
fig, ax = plt.subplots(figsize=(10, 10))
for e, label in zip(eccentricities, labels):
if e < 1:
theta = np.linspace(0, 2 * np.pi, 500)
elif e == 1:
theta = np.linspace(-2.8, 2.8, 500)
else:
# 双曲線: cos(theta) > -1/e の範囲
theta_max = np.arccos(-1 / e) - 0.05
theta = np.linspace(-theta_max, theta_max, 500)
r = p / (1 + e * np.cos(theta))
# 近すぎる点や負のrを除外
mask = r > 0
r = r[mask]
theta = theta[mask]
x = r * np.cos(theta)
y = r * np.sin(theta)
# 表示範囲の制限
valid = (np.abs(x) < 5) & (np.abs(y) < 5)
ax.plot(x[valid], y[valid], linewidth=2, label=label)
# 焦点（原点）をマーク
ax.plot(0, 0, 'ko', markersize=8, label='焦点（中心天体）')
ax.set_xlim(-4, 3)
ax.set_ylim(-3.5, 3.5)
ax.set_xlabel('x [正規化単位]', fontsize=12)
ax.set_ylabel('y [正規化単位]', fontsize=12)
ax.set_title('離心率による軌道形状の変化', fontsize=14)
ax.legend(fontsize=10, loc='upper right')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このグラフからいくつかの重要なことが読み取れます。
- 離心率 $e=0$ の円軌道は焦点を中心とする真円になっています。これは軌道方程式で $r = p/(1+0) = p = \text{const.}$ となることに対応します。
- 離心率が大きくなるにつれて、軌道は徐々に引き伸ばされることがわかります。$e=0.9$ の楕円は非常に細長い形をしており、彗星の軌道を連想させます。
- $e=1$ の放物線と $e=1.5$ の双曲線は開いた軌道であり、焦点に最も近づいた後は無限遠に飛び去ります。恒星間天体オウムアムアはこのような双曲線軌道で太陽系を通過しました。
- 全ての軌道は焦点（中心天体）の位置で共通の近点距離を持つことが確認できます。これは半通径 $p$ を固定したためです。
数値積分による軌道シミュレーション
次に、運動方程式を数値積分して軌道を計算します。理論的に得られた円錐曲線と一致するか確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def two_body_eom(t, state, mu):
"""二体問題の運動方程式（2次元）"""
x, y, vx, vy = state
r = np.sqrt(x**2 + y**2)
ax = -mu * x / r**3
ay = -mu * y / r**3
return [vx, vy, ax, ay]
# 重力パラメータ（正規化）
mu = 1.0
# 3種類の初期条件を設定（近点から出発）
# 近点距離 r_p = 1.0, 近点で動径方向速度ゼロ
r_p = 1.0
cases = {}
# 楕円軌道: e = 0.5
e_ellipse = 0.5
p_ellipse = r_p * (1 + e_ellipse)
h_ellipse = np.sqrt(mu * p_ellipse)
v_ellipse = h_ellipse / r_p
cases['楕円 (e=0.5)'] = {
'state0': [r_p, 0, 0, v_ellipse],
'T': 2 * np.pi * (p_ellipse / (1 - e_ellipse**2))**1.5 / np.sqrt(mu),
'color': '#2196F3',
'e': e_ellipse
}
# 放物線軌道: e = 1.0
v_parabola = np.sqrt(2 * mu / r_p)
cases['放物線 (e=1.0)'] = {
'state0': [r_p, 0, 0, v_parabola],
'T': 15.0,
'color': '#4CAF50',
'e': 1.0
}
# 双曲線軌道: e = 1.5
e_hyp = 1.5
p_hyp = r_p * (1 + e_hyp)
h_hyp = np.sqrt(mu * p_hyp)
v_hyp = h_hyp / r_p
cases['双曲線 (e=1.5)'] = {
'state0': [r_p, 0, 0, v_hyp],
'T': 8.0,
'color': '#F44336',
'e': e_hyp
}
fig, ax = plt.subplots(figsize=(10, 10))
for name, case in cases.items():
# 数値積分
sol = solve_ivp(
two_body_eom, [0, case['T']], case['state0'],
args=(mu,), method='RK45', max_step=0.01,
rtol=1e-10, atol=1e-12
)
x_num = sol.y[0]
y_num = sol.y[1]
# 理論曲線
e = case['e']
p = r_p * (1 + e)
if e < 1:
theta_th = np.linspace(0, 2 * np.pi, 500)
elif e == 1:
theta_th = np.linspace(-2.8, 2.8, 500)
else:
theta_max = np.arccos(-1 / e) - 0.05
theta_th = np.linspace(-theta_max, theta_max, 500)
r_th = p / (1 + e * np.cos(theta_th))
mask = r_th > 0
x_th = r_th[mask] * np.cos(theta_th[mask])
y_th = r_th[mask] * np.sin(theta_th[mask])
# プロット（数値解は実線、理論は破線）
valid = (np.abs(x_th) < 8) & (np.abs(y_th) < 8)
ax.plot(x_th[valid], y_th[valid], '--', color=case['color'],
linewidth=1.5, alpha=0.7)
ax.plot(x_num, y_num, '-', color=case['color'],
linewidth=2.5, label=f'{name}')
# 中心天体
ax.plot(0, 0, 'o', color='#FFC107', markersize=15,
markeredgecolor='black', markeredgewidth=1, label='中心天体')
ax.set_xlim(-6, 4)
ax.set_ylim(-5, 5)
ax.set_xlabel('x [正規化単位]', fontsize=12)
ax.set_ylabel('y [正規化単位]', fontsize=12)
ax.set_title('数値積分による軌道（実線）と理論曲線（破線）の比較',
fontsize=14)
ax.legend(fontsize=11, loc='upper right')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このシミュレーション結果から、以下のことが確認できます。
- 数値積分の結果（実線）と理論的な円錐曲線（破線）が完全に重なっていることから、ビネの公式から得た軌道方程式の正しさが数値的に検証されました。
- 楕円軌道は閉じた軌道であり、衛星は同じ経路を繰り返し周回します。初期条件から1周期分の積分で軌道が閉じることも確認できます。
- 放物線・双曲線軌道は開いた軌道であり、中心天体に最も近づいた後は永遠に遠ざかっていきます。放物線は双曲線よりも中心天体に長くとどまる（曲がりがゆるやか）ことが見て取れます。
エネルギーと角運動量の保存確認
数値積分の精度を検証するために、エネルギーと角運動量が保存されているか確認しましょう。
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def two_body_eom(t, state, mu):
"""二体問題の運動方程式（2次元）"""
x, y, vx, vy = state
r = np.sqrt(x**2 + y**2)
ax = -mu * x / r**3
ay = -mu * y / r**3
return [vx, vy, ax, ay]
mu = 1.0
r_p = 1.0
e = 0.5
p = r_p * (1 + e)
h0 = np.sqrt(mu * p)
v0 = h0 / r_p
# 楕円軌道の数値積分（5周期分）
a = p / (1 - e**2)
T = 2 * np.pi * a**1.5 / np.sqrt(mu)
sol = solve_ivp(
two_body_eom, [0, 5 * T], [r_p, 0, 0, v0],
args=(mu,), method='RK45', max_step=0.01,
rtol=1e-12, atol=1e-14, dense_output=True
)
t_eval = np.linspace(0, 5 * T, 5000)
states = sol.sol(t_eval)
x, y, vx, vy = states
# エネルギーと角運動量の計算
r = np.sqrt(x**2 + y**2)
v2 = vx**2 + vy**2
energy = 0.5 * v2 - mu / r
angular_mom = x * vy - y * vx
# 初期値からの相対誤差
E0 = energy[0]
h_init = angular_mom[0]
energy_err = np.abs((energy - E0) / E0)
h_err = np.abs((angular_mom - h_init) / h_init)
fig, axes = plt.subplots(2, 1, figsize=(10, 6))
axes[0].semilogy(t_eval / T, energy_err, color='#E91E63', linewidth=1)
axes[0].set_ylabel('|ΔE/E₀|', fontsize=12)
axes[0].set_title('エネルギー保存の数値誤差', fontsize=13)
axes[0].grid(True, alpha=0.3)
axes[1].semilogy(t_eval / T, h_err, color='#3F51B5', linewidth=1)
axes[1].set_xlabel('時間 [周期数]', fontsize=12)
axes[1].set_ylabel('|Δh/h₀|', fontsize=12)
axes[1].set_title('角運動量保存の数値誤差', fontsize=13)
axes[1].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
このグラフからは、RK45法（4-5次のルンゲ・クッタ法）による数値積分が高い精度を達成していることが読み取れます。
- エネルギーの相対誤差は $10^{-10}$ 程度に抑えられており、5周期分の積分でも良好な保存性を示しています。長時間の積分ではシンプレクティック積分法（Verlet法など）を使うと、エネルギー誤差がさらに小さく保たれます。
- 角運動量の保存精度もエネルギーと同程度であり、数値誤差が物理量の保存則を損なっていないことが確認できます。
速度ベクトルの可視化
楕円軌道上の各点における速度ベクトルを描画して、ケプラーの第2法則（面積速度一定）を直感的に理解しましょう。
import numpy as np
import matplotlib.pyplot as plt
mu = 1.0
e = 0.6
p = 1.5
a = p / (1 - e**2)
h = np.sqrt(mu * p)
# 軌道全体
theta_full = np.linspace(0, 2 * np.pi, 500)
r_full = p / (1 + e * np.cos(theta_full))
x_full = r_full * np.cos(theta_full)
y_full = r_full * np.sin(theta_full)
# 速度ベクトルを表示する点
n_arrows = 12
theta_arrows = np.linspace(0, 2 * np.pi, n_arrows, endpoint=False)
r_arrows = p / (1 + e * np.cos(theta_arrows))
x_arr = r_arrows * np.cos(theta_arrows)
y_arr = r_arrows * np.sin(theta_arrows)
# 極座標での速度成分
vr = mu / h * e * np.sin(theta_arrows) # 動径方向速度
vt = mu / h * (1 + e * np.cos(theta_arrows)) # 接線方向速度
# 直交座標に変換
vx = vr * np.cos(theta_arrows) - vt * np.sin(theta_arrows)
vy = vr * np.sin(theta_arrows) + vt * np.cos(theta_arrows)
fig, ax = plt.subplots(figsize=(10, 10))
# 軌道
ax.plot(x_full, y_full, 'b-', linewidth=2, label='楕円軌道 (e=0.6)')
# 速度ベクトル
scale = 0.7 # 矢印のスケーリング
for i in range(n_arrows):
ax.annotate('', xy=(x_arr[i] + scale * vx[i], y_arr[i] + scale * vy[i]),
xytext=(x_arr[i], y_arr[i]),
arrowprops=dict(arrowstyle='->', color='red', lw=2))
ax.plot(x_arr[i], y_arr[i], 'ro', markersize=5)
# 中心天体と楕円の中心
ax.plot(0, 0, 'o', color='#FFC107', markersize=15,
markeredgecolor='black', markeredgewidth=1, label='焦点')
# 近点・遠点ラベル
ax.annotate('近点', xy=(p / (1 + e), 0), fontsize=11,
textcoords="offset points", xytext=(10, -15))
ax.annotate('遠点', xy=(-p / (1 - e), 0), fontsize=11,
textcoords="offset points", xytext=(-40, -15))
ax.set_xlabel('x [正規化単位]', fontsize=12)
ax.set_ylabel('y [正規化単位]', fontsize=12)
ax.set_title('楕円軌道上の速度ベクトル（ケプラーの第2法則）',
fontsize=14)
ax.legend(fontsize=11, loc='upper right')
ax.set_aspect('equal')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
この図から、ケプラーの第2法則の直感的な意味が視覚的に理解できます。
- 近点付近では速度ベクトルが長い（速度が大きい）のに対し、遠点付近では速度ベクトルが短い（速度が小さい）ことがはっきり見えます。これは面積速度一定の法則 $h = r^2\dot{\theta} = \text{const.}$ の直接的な帰結です。
- 速度ベクトルの方向も変化しています。近点と遠点では速度は純粋に接線方向（動径方向成分ゼロ）ですが、中間地点では動径方向にも成分を持ちます。これは衛星が楕円軌道上で中心天体に近づいたり遠ざかったりすることに対応します。
- 全ての矢印が等間隔の真近点角 $\theta$ で配置されているにもかかわらず、近点付近では矢印同士の間隔が密であり、遠点付近では疎です。同じ角度だけ掃くのに近点では短い時間、遠点では長い時間がかかることを反映しています。
ケプラーの法則との対応
ここまでの導出結果を、ケプラーの3法則と対応付けて整理しましょう。
第1法則（楕円軌道の法則）
「惑星は太陽を1つの焦点とする楕円軌道上を運動する」
ビネの公式を解いて得た軌道方程式 $r = p/(1+e\cos\theta)$ は、焦点を原点とする楕円（$0 \leq e < 1$）の極座標方程式そのものです。第1法則はニュートンの万有引力の法則（逆2乗則）から数学的に証明されました。
第2法則（面積速度一定の法則）
「惑星と太陽を結ぶ線分が単位時間に掃く面積は一定」
角運動量 $h = r^2\dot{\theta}$ の保存から直接導かれます。面積速度 $dA/dt = h/2 = \text{const.}$ です。第2法則は中心力（万有引力に限らない）の一般的な性質です。
第3法則（調和の法則）
「軌道周期の2乗は半長軸の3乗に比例する」
1周期で掃く面積は楕円の面積 $\pi ab$ に等しいので
$$ \frac{h}{2}T = \pi ab $$
$b = a\sqrt{1-e^2}$ と $h = \sqrt{\mu p} = \sqrt{\mu a(1-e^2)}$ を代入して $T$ について解くと
$$ T = \frac{2\pi a \cdot a\sqrt{1-e^2}}{\sqrt{\mu a(1-e^2)}} = \frac{2\pi a^{3/2}}{\sqrt{\mu}} $$
したがって
$$ \boxed{T^2 = \frac{4\pi^2}{\mu}a^3} $$
周期は離心率 $e$ によらず、半長軸 $a$ のみで決まるという驚くべき結果です。
まとめ
本記事では、二体問題の定式化からケプラー軌道の数学的導出までを解説しました。
- 二体問題は1体問題に帰着できる — 重心運動と相対運動に分離し、相対運動は換算質量を用いた中心力問題になる
- 角運動量保存により運動は平面内に限られる — 3次元問題が2次元に帰着し、ケプラーの第2法則が導かれる
- ビネの公式により軌道は円錐曲線になることが示される — ケプラーの第1法則の数学的証明
- 軌道エネルギーの符号が軌道の種類を決定する — 負なら楕円（束縛軌道）、ゼロなら放物線、正なら双曲線（脱出軌道）
- ビザビバ方程式は軌道上の任意の点での速度を与える — 軌道遷移の計算に不可欠
二体問題の理解は、軌道力学の全ての議論の出発点です。現実の宇宙では、第三の天体の重力摂動、大気抵抗、太陽輻射圧などの効果が加わりますが、それらは全て二体問題の解を基準として摂動として扱われます。
次のステップとして、以下の記事も参考にしてください。