二体問題（two-body problem）は、万有引力で相互作用する2つの天体の運動を求める問題です。天体力学の最も基本的な問題であり、解析的に解ける数少ない多体問題のひとつです。ケプラー軌道はこの解として得られます。
本記事の内容
- 二体問題の定式化
- 換算質量による1体問題への帰着
- 軌道方程式 $r = p/(1 + e\cos\nu)$ の導出
- 比エネルギーと軌道形状の関係
- Pythonでの数値積分による検証
前提知識
この記事を読む前に、以下の記事を読んでおくと理解が深まります。
二体問題の定式化
質量 $m_1$, $m_2$ の2つの天体が万有引力のみで相互作用するとき、ニュートンの運動方程式は:
$$ m_1 \ddot{\bm{r}}_1 = \frac{Gm_1 m_2}{|\bm{r}_2 – \bm{r}_1|^3}(\bm{r}_2 – \bm{r}_1) $$
$$ m_2 \ddot{\bm{r}}_2 = -\frac{Gm_1 m_2}{|\bm{r}_2 – \bm{r}_1|^3}(\bm{r}_2 – \bm{r}_1) $$
重心座標と相対座標
重心 $\bm{R}$ と相対位置 $\bm{r}$ を定義します:
$$ \bm{R} = \frac{m_1 \bm{r}_1 + m_2 \bm{r}_2}{m_1 + m_2}, \quad \bm{r} = \bm{r}_2 – \bm{r}_1 $$
2式を足すと重心の運動方程式:
$$ (m_1 + m_2)\ddot{\bm{R}} = \bm{0} \implies \bm{R}(t) = \bm{R}_0 + \dot{\bm{R}}_0 t $$
重心は等速直線運動します。
換算質量と相対運動
2式を変形すると、相対運動の方程式が得られます:
$$ \ddot{\bm{r}} = \ddot{\bm{r}}_2 – \ddot{\bm{r}}_1 = -\frac{G(m_1 + m_2)}{r^3}\bm{r} $$
換算質量 $\mu_r$ と重力パラメータ $\mu$ を定義:
$$ \mu_r = \frac{m_1 m_2}{m_1 + m_2}, \quad \mu = G(m_1 + m_2) $$
相対運動の方程式:
$$ \boxed{\ddot{\bm{r}} = -\frac{\mu}{r^3}\bm{r}} $$
二体問題は、質量 $\mu_r$ の仮想的な1つの粒子が中心力場 $\mu$ の下で運動する1体問題に帰着されました。
保存量
角運動量
中心力のためトルクはゼロ:
$$ \bm{h} = \bm{r} \times \dot{\bm{r}} = \text{const} $$
$\bm{h}$ は一定ベクトルなので、運動は $\bm{h}$ に垂直な平面内に限られます。
比力学的エネルギー
$$ \boxed{\varepsilon = \frac{v^2}{2} – \frac{\mu}{r} = \text{const}} $$
軌道方程式の導出
極座標 $(r, \nu)$ を導入します。比角運動量は:
$$ h = r^2 \dot{\nu} $$
$u = 1/r$ と変換し、独立変数を $t$ から $\nu$ に変更します。
$$ \dot{r} = \frac{dr}{dt} = \frac{dr}{d\nu}\dot{\nu} = -h\frac{du}{d\nu} $$
$$ \ddot{r} = -h^2 u^2 \frac{d^2u}{d\nu^2} $$
動径方向の運動方程式 $\ddot{r} – r\dot{\nu}^2 = -\mu/r^2$ に代入:
$$ -h^2 u^2 \frac{d^2u}{d\nu^2} – h^2 u^3 = -\mu u^2 $$
$$ \frac{d^2 u}{d\nu^2} + u = \frac{\mu}{h^2} $$
この2階線形ODEの一般解は:
$$ u = \frac{\mu}{h^2} + A\cos(\nu – \nu_0) $$
$\nu_0 = 0$（近地点方向を基準）とし、$r = 1/u$ に戻すと:
$$ r = \frac{h^2/\mu}{1 + (Ah^2/\mu)\cos\nu} $$
半通径 $p = h^2/\mu$ と離心率 $e = Ah^2/\mu$ を定義して:
$$ \boxed{r = \frac{p}{1 + e\cos\nu} = \frac{a(1-e^2)}{1 + e\cos\nu}} $$
これが軌道方程式（円錐曲線の極座標表示）です。
比エネルギーと軌道形状
軌道方程式をエネルギーで表すと:
$$ \boxed{\varepsilon = -\frac{\mu}{2a}} $$
| 軌道形状 | 離心率 $e$ | 半長軸 $a$ | エネルギー $\varepsilon$ |
|---|---|---|---|
| 円 | $e = 0$ | $a > 0$ | $\varepsilon < 0$ |
| 楕円 | $0 < e < 1$ | $a > 0$ | $\varepsilon < 0$ |
| 放物線 | $e = 1$ | $a \to \infty$ | $\varepsilon = 0$ |
| 双曲線 | $e > 1$ | $a < 0$ | $\varepsilon > 0$ |
Pythonでの実装
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
mu = 3.986004418e14 # 地球 [m^3/s^2]
R_earth = 6371e3 # [m]
def two_body_ode(t, state, mu):
"""二体問題の運動方程式"""
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]
fig, axes = plt.subplots(2, 2, figsize=(13, 11))
# (a) 軌道方程式: 異なる離心率の軌道
nu = np.linspace(0, 2*np.pi, 1000)
p = 2 * R_earth # 半通径
for e_val, ls, label in [(0.0, '-', '円 (e=0)'),
(0.3, '-', '楕円 (e=0.3)'),
(0.7, '-', '楕円 (e=0.7)'),
(0.95, '--', '楕円 (e=0.95)')]:
r_orbit = p / (1 + e_val * np.cos(nu))
x_orbit = r_orbit * np.cos(nu) / R_earth
y_orbit = r_orbit * np.sin(nu) / R_earth
axes[0, 0].plot(x_orbit, y_orbit, ls, linewidth=2, label=label)
# 地球
theta_e = np.linspace(0, 2*np.pi, 100)
axes[0, 0].fill(np.cos(theta_e), np.sin(theta_e), color='cyan', alpha=0.3)
axes[0, 0].plot(0, 0, 'ko', markersize=5)
axes[0, 0].set_xlabel('x / $R_E$')
axes[0, 0].set_ylabel('y / $R_E$')
axes[0, 0].set_title('軌道方程式 $r = p/(1+e\\cos\\nu)$')
axes[0, 0].set_aspect('equal')
axes[0, 0].legend(fontsize=9)
axes[0, 0].grid(True, alpha=0.3)
axes[0, 0].set_xlim(-8, 5)
# (b) 数値解と解析解の比較
a = 2 * R_earth
e = 0.5
r_p = a * (1 - e)
v_p = np.sqrt(mu * (2/r_p - 1/a))
T = 2 * np.pi * np.sqrt(a**3 / mu)
sol = solve_ivp(two_body_ode, [0, T], [r_p, 0, 0, v_p],
args=(mu,), t_eval=np.linspace(0, T, 1000),
rtol=1e-12, atol=1e-14)
axes[0, 1].plot(sol.y[0]/R_earth, sol.y[1]/R_earth, 'b-', linewidth=2, label='数値解')
r_analytic = a * (1 - e**2) / (1 + e * np.cos(nu))
x_a = r_analytic * np.cos(nu) / R_earth
y_a = r_analytic * np.sin(nu) / R_earth
axes[0, 1].plot(x_a, y_a, 'r--', linewidth=1, label='解析解')
axes[0, 1].fill(np.cos(theta_e), np.sin(theta_e), color='cyan', alpha=0.3)
axes[0, 1].set_xlabel('x / $R_E$')
axes[0, 1].set_ylabel('y / $R_E$')
axes[0, 1].set_title('数値解 vs 解析解 (e=0.5)')
axes[0, 1].set_aspect('equal')
axes[0, 1].legend()
axes[0, 1].grid(True, alpha=0.3)
# (c) エネルギーと角運動量の保存
r_vals = np.sqrt(sol.y[0]**2 + sol.y[1]**2)
v_vals = np.sqrt(sol.y[2]**2 + sol.y[3]**2)
energy = 0.5 * v_vals**2 - mu / r_vals
h_vals = sol.y[0] * sol.y[3] - sol.y[1] * sol.y[2]
t_norm = sol.t / T
axes[1, 0].plot(t_norm, (energy - energy[0])/np.abs(energy[0]), 'b-', linewidth=1.5,
label='比エネルギー偏差')
ax_twin = axes[1, 0].twinx()
ax_twin.plot(t_norm, (h_vals - h_vals[0])/np.abs(h_vals[0]), 'r-', linewidth=1.5,
label='比角運動量偏差')
axes[1, 0].set_xlabel('$t / T$')
axes[1, 0].set_ylabel('エネルギー相対偏差', color='b')
ax_twin.set_ylabel('角運動量相対偏差', color='r')
axes[1, 0].set_title('保存量の数値精度')
axes[1, 0].grid(True, alpha=0.3)
# (d) 放物線・双曲線軌道
nu_open = np.linspace(-2.5, 2.5, 1000)
for e_val, label in [(1.0, '放物線 (e=1.0)'),
(1.5, '双曲線 (e=1.5)'),
(3.0, '双曲線 (e=3.0)')]:
r_open = p / (1 + e_val * np.cos(nu_open))
mask = r_open > 0
x_open = r_open[mask] * np.cos(nu_open[mask]) / R_earth
y_open = r_open[mask] * np.sin(nu_open[mask]) / R_earth
valid = np.abs(x_open) < 15
axes[1, 1].plot(x_open[valid], y_open[valid], linewidth=2, label=label)
axes[1, 1].fill(np.cos(theta_e), np.sin(theta_e), color='cyan', alpha=0.3)
axes[1, 1].plot(0, 0, 'ko', markersize=5)
axes[1, 1].set_xlabel('x / $R_E$')
axes[1, 1].set_ylabel('y / $R_E$')
axes[1, 1].set_title('放物線・双曲線軌道')
axes[1, 1].set_aspect('equal')
axes[1, 1].legend(fontsize=9)
axes[1, 1].grid(True, alpha=0.3)
axes[1, 1].set_xlim(-12, 8)
axes[1, 1].set_ylim(-10, 10)
plt.tight_layout()
plt.savefig('two_body_problem.png', dpi=150, bbox_inches='tight')
plt.show()
まとめ
本記事では、二体問題の解法について解説しました。
- 二体問題は換算質量を用いて1体の中心力問題に帰着される
- 相対運動方程式: $\ddot{\bm{r}} = -\mu\bm{r}/r^3$
- 軌道方程式: $r = p/(1 + e\cos\nu)$（円錐曲線）
- 比エネルギー $\varepsilon = -\mu/(2a)$ と角運動量 $h$ が保存量
- $e < 1$: 束縛軌道（楕円）、$e \geq 1$: 非束縛軌道（放物線・双曲線）
次のステップとして、以下の記事も参考にしてください。