1つの質点の力学を理解した次のステップは、複数の質点からなる系（質点系） の力学です。質点系では「重心」が系全体の運動を代表し、内力は重心の運動に影響しないという重要な性質があります。
天体の運動、分子間の相互作用、衝突問題など、物理学のあらゆる場面で質点系の概念は登場します。本記事では質点系の基本法則を導出し、2体問題を例にPythonでシミュレーションを行います。
本記事の内容
- 質点系の重心の定義と運動方程式
- 内力と外力の分離
- 運動量・角運動量・エネルギーの保存則
- 2体問題と換算質量
- Pythonで2体運動のシミュレーション
前提知識
質点系の重心
$N$ 個の質点からなる系を考えます。各質点の質量を $m_i$、位置を $\bm{r}_i$（$i = 1, 2, \dots, N$）とします。
重心（質量中心）は次のように定義されます：
$$ \boxed{\bm{R} = \frac{\sum_{i=1}^{N} m_i \bm{r}_i}{\sum_{i=1}^{N} m_i} = \frac{1}{M} \sum_{i=1}^{N} m_i \bm{r}_i} $$
ここで $M = \sum_{i} m_i$ は系の全質量です。
重心の運動方程式
各質点に作用する力を考えます。質点 $i$ に作用する力は：
- 外力 $\bm{F}_i^{(\text{ext})}$：系の外部から作用する力
- 内力 $\bm{f}_{ij}$：質点 $j$ から質点 $i$ への力
ニュートンの第3法則（作用・反作用の法則）より：
$$ \bm{f}_{ij} = -\bm{f}_{ji} $$
質点 $i$ の運動方程式は：
$$ m_i \ddot{\bm{r}}_i = \bm{F}_i^{(\text{ext})} + \sum_{j \neq i} \bm{f}_{ij} $$
すべての質点について和をとると：
$$ \sum_{i=1}^{N} m_i \ddot{\bm{r}}_i = \sum_{i=1}^{N} \bm{F}_i^{(\text{ext})} + \sum_{i=1}^{N}\sum_{j \neq i} \bm{f}_{ij} $$
作用・反作用の法則から、内力の和は：
$$ \sum_{i=1}^{N}\sum_{j \neq i} \bm{f}_{ij} = \sum_{i < j}(\bm{f}_{ij} + \bm{f}_{ji}) = \bm{0} $$
左辺は $M\ddot{\bm{R}}$ に等しいので：
$$ \boxed{M\ddot{\bm{R}} = \bm{F}^{(\text{ext})} \equiv \sum_{i=1}^{N} \bm{F}_i^{(\text{ext})}} $$
これが重心の運動方程式です。重心は全質量 $M$ が集中した1つの質点のように、外力の合力 $\bm{F}^{(\text{ext})}$ によって運動します。内力は重心の運動に一切影響しません。
運動量の保存則
系の全運動量は：
$$ \bm{P} = \sum_{i=1}^{N} m_i \dot{\bm{r}}_i = M\dot{\bm{R}} $$
重心の運動方程式から：
$$ \frac{d\bm{P}}{dt} = \bm{F}^{(\text{ext})} $$
外力の合力がゼロのとき：
$$ \boxed{\bm{F}^{(\text{ext})} = \bm{0} \quad \Longrightarrow \quad \bm{P} = \text{const.}} $$
これが運動量保存則です。
角運動量
全角運動量
原点まわりの全角運動量は：
$$ \bm{L} = \sum_{i=1}^{N} \bm{r}_i \times m_i \dot{\bm{r}}_i $$
この時間微分は：
$$ \frac{d\bm{L}}{dt} = \sum_{i=1}^{N} \bm{r}_i \times m_i \ddot{\bm{r}}_i = \sum_{i=1}^{N} \bm{r}_i \times \bm{F}_i^{(\text{ext})} + \sum_{i=1}^{N}\sum_{j \neq i} \bm{r}_i \times \bm{f}_{ij} $$
内力が中心力（2質点を結ぶ線分に沿う力）のとき、内力のトルクの和はゼロになります。これを示します。
$$ \sum_{i}\sum_{j \neq i} \bm{r}_i \times \bm{f}_{ij} = \sum_{i < j}(\bm{r}_i \times \bm{f}_{ij} + \bm{r}_j \times \bm{f}_{ji}) $$
$$ = \sum_{i < j}(\bm{r}_i - \bm{r}_j) \times \bm{f}_{ij} $$
$\bm{f}_{ij}$ が $\bm{r}_i – \bm{r}_j$ に平行なとき、この外積はゼロになります。したがって：
$$ \boxed{\frac{d\bm{L}}{dt} = \bm{N}^{(\text{ext})} \equiv \sum_{i=1}^{N} \bm{r}_i \times \bm{F}_i^{(\text{ext})}} $$
外力のトルクがゼロなら角運動量は保存されます。
ケーニッヒの定理（角運動量の分解）
重心まわりの相対位置 $\bm{r}_i’ = \bm{r}_i – \bm{R}$ を用いると、全角運動量は：
$$ \boxed{\bm{L} = \bm{R} \times M\dot{\bm{R}} + \sum_{i=1}^{N} \bm{r}_i’ \times m_i \dot{\bm{r}}_i’} $$
第1項は「重心の公転」の角運動量、第2項は「重心まわりの回転（自転）」の角運動量です。
エネルギー
運動エネルギーの分解
系の全運動エネルギーも同様に分解できます：
$$ T = \sum_{i=1}^{N} \frac{1}{2}m_i |\dot{\bm{r}}_i|^2 $$
$\dot{\bm{r}}_i = \dot{\bm{R}} + \dot{\bm{r}}_i’$ を代入すると：
$$ T = \frac{1}{2}M|\dot{\bm{R}}|^2 + \sum_{i=1}^{N}\frac{1}{2}m_i|\dot{\bm{r}}_i’|^2 + \dot{\bm{R}} \cdot \sum_{i} m_i \dot{\bm{r}}_i’ $$
重心の定義から $\sum_i m_i \bm{r}_i’ = \bm{0}$ であり、微分して $\sum_i m_i \dot{\bm{r}}_i’ = \bm{0}$ です。したがって第3項は消え：
$$ \boxed{T = \frac{1}{2}M|\dot{\bm{R}}|^2 + \sum_{i=1}^{N}\frac{1}{2}m_i|\dot{\bm{r}}_i’|^2} $$
第1項は重心の運動エネルギー、第2項は重心まわりの運動エネルギーです。これはケーニッヒの定理と呼ばれます。
2体問題と換算質量
2体問題の定式化
質量 $m_1$, $m_2$ の2質点が相互作用のみで運動する場合を考えます。外力はないものとします。
$$ m_1 \ddot{\bm{r}}_1 = \bm{f}_{12}, \quad m_2 \ddot{\bm{r}}_2 = \bm{f}_{21} = -\bm{f}_{12} $$
相対座標を導入します：
$$ \bm{r} = \bm{r}_1 – \bm{r}_2 $$
第1式を $m_1$ で、第2式を $m_2$ で割って差をとると：
$$ \ddot{\bm{r}}_1 – \ddot{\bm{r}}_2 = \bm{f}_{12}\left(\frac{1}{m_1} + \frac{1}{m_2}\right) $$
$$ \ddot{\bm{r}} = \frac{\bm{f}_{12}}{\mu} $$
ここで $\mu$ は換算質量（reduced mass）：
$$ \boxed{\frac{1}{\mu} = \frac{1}{m_1} + \frac{1}{m_2} \quad \Longleftrightarrow \quad \mu = \frac{m_1 m_2}{m_1 + m_2}} $$
1体問題への帰着
相対運動の方程式は：
$$ \boxed{\mu \ddot{\bm{r}} = \bm{f}_{12}(\bm{r})} $$
これは質量 $\mu$ の1質点が力 $\bm{f}_{12}$ のもとで運動する方程式と同じ形です。つまり2体問題は、重心の等速直線運動と、換算質量 $\mu$ の1体問題に分離できます。
万有引力の場合
$\bm{f}_{12} = -\frac{Gm_1 m_2}{r^2}\hat{r}$ のとき、相対運動の方程式は：
$$ \mu \ddot{\bm{r}} = -\frac{Gm_1 m_2}{r^2}\hat{r} $$
両辺を $\mu$ で割ると：
$$ \ddot{\bm{r}} = -\frac{G(m_1 + m_2)}{r^2}\hat{r} $$
これはケプラー問題と同じ形であり、相対座標 $\bm{r}$ は楕円軌道を描きます。
Pythonでの可視化
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
# ============================
# 万有引力の2体問題シミュレーション
# ============================
G = 1.0 # 万有引力定数（無次元化）
def two_body_eom(t, state, m1, m2):
"""2体問題の運動方程式"""
x1, y1, vx1, vy1, x2, y2, vx2, vy2 = state
dx = x2 - x1
dy = y2 - y1
r = np.sqrt(dx**2 + dy**2)
r3 = r**3
# 万有引力
fx = G * m1 * m2 * dx / r3
fy = G * m1 * m2 * dy / r3
return [vx1, vy1, fx/m1, fy/m1,
vx2, vy2, -fx/m2, -fy/m2]
fig, axes = plt.subplots(1, 3, figsize=(18, 6))
# === (1) 等質量の2体問題 ===
m1, m2 = 1.0, 1.0
# 重心が原点になるように初期条件を設定
r0 = 2.0 # 初期距離
v0 = 0.5 # 初期速度（円軌道に近い値）
y0_case1 = [r0/2, 0, 0, v0/2,
-r0/2, 0, 0, -v0/2]
sol1 = solve_ivp(two_body_eom, (0, 40), y0_case1,
args=(m1, m2), method='RK45',
rtol=1e-10, atol=1e-12, max_step=0.01)
axes[0].plot(sol1.y[0], sol1.y[1], 'b-', lw=1.2, label=f'$m_1={m1}$')
axes[0].plot(sol1.y[4], sol1.y[5], 'r-', lw=1.2, label=f'$m_2={m2}$')
# 重心
Rx = (m1*sol1.y[0] + m2*sol1.y[4]) / (m1+m2)
Ry = (m1*sol1.y[1] + m2*sol1.y[5]) / (m1+m2)
axes[0].plot(Rx[0], Ry[0], 'k+', ms=15, mew=2, label='CM')
axes[0].set_xlabel('x'); axes[0].set_ylabel('y')
axes[0].set_title(f'Equal mass ($m_1={m1}$, $m_2={m2}$)', fontsize=12)
axes[0].legend(fontsize=10); axes[0].set_aspect('equal')
axes[0].grid(True, alpha=0.3)
# === (2) 質量比 1:3 ===
m1, m2 = 1.0, 3.0
mu = m1*m2/(m1+m2)
# 重心系で初期条件設定
v_rel = 0.7
y0_case2 = [m2/(m1+m2)*r0, 0, 0, m2/(m1+m2)*v_rel,
-m1/(m1+m2)*r0, 0, 0, -m1/(m1+m2)*v_rel]
sol2 = solve_ivp(two_body_eom, (0, 40), y0_case2,
args=(m1, m2), method='RK45',
rtol=1e-10, atol=1e-12, max_step=0.01)
axes[1].plot(sol2.y[0], sol2.y[1], 'b-', lw=1.2, label=f'$m_1={m1}$')
axes[1].plot(sol2.y[4], sol2.y[5], 'r-', lw=1.2, label=f'$m_2={m2}$')
Rx2 = (m1*sol2.y[0] + m2*sol2.y[4]) / (m1+m2)
Ry2 = (m1*sol2.y[1] + m2*sol2.y[5]) / (m1+m2)
axes[1].plot(Rx2[0], Ry2[0], 'k+', ms=15, mew=2, label='CM')
axes[1].set_xlabel('x'); axes[1].set_ylabel('y')
axes[1].set_title(f'Unequal mass ($m_1={m1}$, $m_2={m2}$)', fontsize=12)
axes[1].legend(fontsize=10); axes[1].set_aspect('equal')
axes[1].grid(True, alpha=0.3)
# === (3) 換算質量 vs 質量比 ===
m1_fixed = 1.0
m2_range = np.linspace(0.01, 10, 500)
mu_range = m1_fixed * m2_range / (m1_fixed + m2_range)
axes[2].plot(m2_range/m1_fixed, mu_range/m1_fixed, 'b-', lw=2)
axes[2].axhline(1, color='gray', ls='--', alpha=0.5, label='$\\mu/m_1 = 1$')
axes[2].set_xlabel('$m_2 / m_1$', fontsize=12)
axes[2].set_ylabel('$\\mu / m_1$', fontsize=12)
axes[2].set_title('Reduced mass ratio', fontsize=12)
axes[2].legend(fontsize=10)
axes[2].grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('two_body_problem.png', dpi=150, bbox_inches='tight')
plt.show()
# === エネルギー保存の確認 ===
fig2, ax = plt.subplots(figsize=(10, 5))
m1, m2 = 1.0, 3.0
dx = sol2.y[4] - sol2.y[0]
dy = sol2.y[5] - sol2.y[1]
r = np.sqrt(dx**2 + dy**2)
T1 = 0.5 * m1 * (sol2.y[2]**2 + sol2.y[3]**2)
T2 = 0.5 * m2 * (sol2.y[6]**2 + sol2.y[7]**2)
T_total = T1 + T2
V = -G * m1 * m2 / r
E_total = T_total + V
ax.plot(sol2.t, T_total, 'r-', lw=1.5, label='Kinetic $T$')
ax.plot(sol2.t, V, 'b-', lw=1.5, label='Potential $V$')
ax.plot(sol2.t, E_total, 'k--', lw=2, label='Total $E = T + V$')
ax.set_xlabel('Time', fontsize=12)
ax.set_ylabel('Energy', fontsize=12)
ax.set_title('Energy conservation in two-body problem', fontsize=13)
ax.legend(fontsize=11)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig('two_body_energy.png', dpi=150, bbox_inches='tight')
plt.show()
上段の左と中央の図は、等質量と質量比 1:3 の2体問題の軌道を示しています。等質量の場合、2つの質点は重心を中心として対称に周回します。質量比が異なる場合、重い質点の軌道は小さく、軽い質点の軌道は大きくなります。右の図は、換算質量 $\mu$ が質量比とともに単調増加し、$m_2 \gg m_1$ では $\mu \to m_1$ に近づくことを示しています。
下の図では、運動エネルギー・ポテンシャルエネルギーが個別に変動しても全エネルギーが一定に保たれることを確認できます。
まとめ
本記事では、質点系の力学について解説しました。
- 重心の運動方程式: $M\ddot{\bm{R}} = \bm{F}^{(\text{ext})}$ — 内力は重心に影響しない
- 運動量保存則: 外力がゼロなら全運動量は保存
- 角運動量: 外力のトルクがゼロなら角運動量は保存
- ケーニッヒの定理: 運動エネルギーは重心の運動と重心まわりの運動に分解できる
- 2体問題: 換算質量 $\mu = m_1 m_2/(m_1 + m_2)$ を用いて1体問題に帰着できる
次のステップとして、以下の記事も参考にしてください。