はじめに
前回の記事では、万有引力を仮定した数値計算によって、
惑星の公転周期の2乗と楕円軌道の長半径の3乗が比例するというケプラーの第3法則を視覚的に確認した。
しかし、そこで得られた楕円軌道は、あくまで数値計算の結果として「そうなった」ものであり、なぜ惑星が恒星から一定の範囲内を往復し続けるのか、すなわち、なぜ軌道が閉じ込められるのかという問いには答えていない。
2次元の惑星の運動は、中心力のもとでは角運動量が保存されるため、動径 $r$ のみを変数とする1次元の運動に帰着させることができる。
このとき、動径方向の運動エネルギーと足し合わせると力学的エネルギーとして保存されるポテンシャル、すなわち有効ポテンシャル $W(r)$ が現れる。この $W(r)$ の形を見れば、惑星がどの範囲の $r$ を運動するのかが一目でわかる。
そこで本記事では、有効ポテンシャル $W(r)$ を理論的に導出し、それと前回の記事と同じ数値計算の結果とを比較することで、両者が一致することを確認する。
理論
運動方程式と単位系
前回の記事と同様に、原点に質量 $M$ の恒星を固定し、その周りを質量 $m$ $(m \ll M)$ の惑星が運動するものとする。惑星に働く万有引力は、
\vec{F}(\vec{r})=-G\frac{Mm}{r^2}\frac{\vec{r}}{r}
である。以下では簡単のため、前回の記事のプログラムと同じく $GM=1$、$m=1$ とする単位系を用いる。すなわち、エネルギーや角運動量はすべて単位質量あたりの量として扱う。このとき、万有引力のポテンシャルエネルギーは、
U(r)=-\frac{1}{r}
となる。
角運動量の保存
極座標 $(r,\theta)$ を用いると、惑星の速度の2乗は動径方向成分と角度方向成分に分解でき、
v^2=\dot{r}^2+r^2\dot{\theta}^2
と書ける。万有引力は常に原点を向く中心力であるため、原点まわりのトルクはゼロとなり、角運動量
L=r^2\dot{\theta}
は時間によらず一定である。
力学的エネルギーの保存と有効ポテンシャル
力学的エネルギー $E$ は運動エネルギーとポテンシャルエネルギーの和であり、保存される。
E=\frac{1}{2}\left(\dot{r}^2+r^2\dot{\theta}^2\right)-\frac{1}{r}
ここで、角運動量の保存則 $\dot{\theta}=L/r^2$ を代入して $\dot{\theta}$ を消去すると、
E=\frac{1}{2}\dot{r}^2+\frac{L^2}{2r^2}-\frac{1}{r}
を得る。そこで、
W(r)=-\frac{1}{r}+\frac{L^2}{2r^2}
と定義すれば、力学的エネルギーは
E=\frac{1}{2}\dot{r}^2+W(r)
と表される。これは、ポテンシャル $W(r)$ の中を1次元的に運動する質点のエネルギー保存則そのものである。つまり、動径方向の運動エネルギー $\frac{1}{2}\dot{r}^2$ と $W(r)$ の和が保存される。この $W(r)$ を有効ポテンシャルと呼ぶ。
$W(r)$ の第1項は万有引力による引力のポテンシャル、第2項は角度方向の運動に由来する遠心力のポテンシャルである。$r$ が小さい領域では第2項が支配的となって $W(r)\to+\infty$ となり、$r$ が大きい領域では第1項が支配的となって $W(r)\to 0^-$ となる。したがって、$W(r)$ は井戸型の形状をもつ。
井戸の底
$W(r)$ を $r$ で微分してゼロとおくと、
\frac{dW}{dr}=\frac{1}{r^2}-\frac{L^2}{r^3}=0
より、$r=L^2$ で最小値
W_{\min}=-\frac{1}{2L^2}
をとる。$E=W_{\min}$ のときは $\dot{r}=0$ が常に成り立つため、惑星は半径 $L^2$ の円軌道を描く。
転回点
$\frac{1}{2}\dot{r}^2\ge 0$ であるから、惑星が運動できるのは $W(r)\le E$ を満たす領域に限られる。その境界、すなわち $\dot{r}=0$ となる点を転回点と呼ぶ。$W(r)=E$ を整理すると、
Er^2+r-\frac{L^2}{2}=0
という $r$ についての2次方程式になる。$E<0$ のとき、この方程式は2つの正の実数解
r_{\min}=\frac{1-\sqrt{1+2EL^2}}{-2E},\qquad r_{\max}=\frac{1+\sqrt{1+2EL^2}}{-2E}
をもつ。惑星は $r_{\min}\le r\le r_{\max}$ の範囲に閉じ込められ、その間を往復する。これが束縛運動、すなわち楕円軌道に対応する。
なお、楕円の長半径 $a$ は $r_{\min}$ と $r_{\max}$ の平均であるから、
a=\frac{r_{\min}+r_{\max}}{2}=-\frac{1}{2E}
となる。すなわち、長半径は力学的エネルギー $E$ のみで決まり、角運動量 $L$ にはよらない。前回の記事で扱ったケプラーの第3法則において、周期が長半径のみで決まることと対応している。
また、離心率 $e$ は
e=\sqrt{1+2EL^2}
であり、$r_{\min}=a(1-e)$、$r_{\max}=a(1+e)$ と書ける。
数値計算
計算条件
前回の記事の初期条件の範囲 $1\le x\le 2$、$0.5\le v_y\le 0.8$ から、以下の1例を選んだ。
| 物理量 | 値 |
|---|---|
| 初期位置 $(x_0,y_0)$ | $(1.5,\ 0)$ |
| 初速度 $(v_{x0},v_{y0})$ | $(0,\ 0.65)$ |
| 時間刻み $h$ | $1.0\times10^{-3}$ |
| 計算時間 | $0\le t\le 50$ |
このとき、保存量は
L=x_0v_{y0}-y_0v_{x0}=0.975
E=\frac{1}{2}\left(v_{x0}^2+v_{y0}^2\right)-\frac{1}{r_0}\approx-0.4554
である。$E<0$ であるから、惑星は束縛運動を行う。
時間発展には前回の記事と同じく、速度を先に更新し、更新後の速度で位置を更新する半陰的オイラー法（シンプレクティック・オイラー法）を用いた。
数値解との比較には、動径方向の速度
\dot{r}=\frac{x\dot{x}+y\dot{y}}{r}
を各時刻で求め、エネルギー保存則から逆算した量 $E-\frac{1}{2}\dot{r}^2$ をその時刻の $r$ に対してプロットする。これが理論式 $W(r)$ の曲線上に乗れば、理論と数値計算が一致していることになる。
プログラム
import numpy as np
import matplotlib.pyplot as plt
import japanize_matplotlib
# 前回の記事と同じ単位系: GM = 1, m = 1 (加速度 a = -r^alpha * r_vec/r, alpha = -2)
h = 1.0e-3
alpha = -2.0
# 初期条件（前回の記事の範囲 1<=x<=2, 0.5<=v_y<=0.8 から1つ選ぶ）
x, y = 1.5, 0.0
v_x, v_y = 0.0, 0.65
# 保存量
L = x * v_y - y * v_x # 角運動量（単位質量あたり）
E = 0.5 * (v_x**2 + v_y**2) - 1.0 / np.hypot(x, y) # 力学的エネルギー
# 有効ポテンシャル W(r) = -1/r + L^2/(2 r^2)
def W(r):
return -1.0 / r + L**2 / (2.0 * r**2)
# 半陰的オイラー法で軌道計算
r_ary, rdot_ary, t_ary = [], [], []
t = 0.0
while t <= 50:
r = np.hypot(x, y)
a_x = -r**alpha * x / r
a_y = -r**alpha * y / r
v_x += a_x * h
v_y += a_y * h
x += v_x * h
y += v_y * h
t += h
r = np.hypot(x, y)
r_ary.append(r)
rdot_ary.append((x * v_x + y * v_y) / r) # 動径速度 dr/dt
t_ary.append(t)
r_ary = np.array(r_ary)
rdot_ary = np.array(rdot_ary)
t_ary = np.array(t_ary)
K_r = 0.5 * rdot_ary**2 # 動径方向の運動エネルギー
W_num = E - K_r # エネルギー保存から逆算した W
# 転回点（W(r) = E の解）
r_min = (1 - np.sqrt(1 + 2 * E * L**2)) / (-2 * E)
r_max = (1 + np.sqrt(1 + 2 * E * L**2)) / (-2 * E)
# ---- 描画 ----
fig, axes = plt.subplots(1, 2, figsize=(13, 5.2))
ax = axes[0]
rr = np.linspace(0.25, 4.0, 1000)
ax.plot(rr, W(rr), color='black', lw=2, label=r'$W(r)=-\frac{1}{r}+\frac{L^2}{2r^2}$')
ax.plot(rr, -1 / rr, '--', color='gray', lw=1, label=r'重力ポテンシャル $-1/r$')
ax.plot(rr, L**2 / (2 * rr**2), ':', color='gray', lw=1, label=r'遠心ポテンシャル $L^2/2r^2$')
ax.axhline(E, color='red', lw=1.5, label=f'力学的エネルギー E = {E:.4f}')
ax.scatter(r_ary[::200], W_num[::200], s=12, color='blue', zorder=5,
label=r'数値解: $E-\frac{1}{2}\dot{r}^2$')
for rt in (r_min, r_max):
ax.axvline(rt, color='red', ls=':', lw=1)
ax.fill_between(rr, W(rr), E, where=(W(rr) <= E), color='orange', alpha=0.25,
label=r'動径運動エネルギー $\frac{1}{2}\dot{r}^2$')
ax.annotate(f'$r_{{min}}$={r_min:.3f}', (r_min, E), xytext=(r_min + 0.05, E + 0.25))
ax.annotate(f'$r_{{max}}$={r_max:.3f}', (r_max, E), xytext=(r_max + 0.05, E + 0.25))
ax.set_xlim(0.25, 4.0)
ax.set_ylim(-1.2, 1.0)
ax.set_xlabel('動径 r')
ax.set_ylabel('エネルギー')
ax.set_title(f'有効ポテンシャル W(r)（L = {L:.3f}）')
ax.legend(fontsize=8, loc='upper right')
ax.grid(alpha=0.3)
ax = axes[1]
ax.plot(t_ary, K_r + W(r_ary) - E, color='blue', lw=1)
ax.set_xlabel('時間 t')
ax.set_ylabel(r'$\frac{1}{2}\dot{r}^2+W(r)-E$')
ax.set_title('エネルギー保存の確認（誤差）')
ax.grid(alpha=0.3)
plt.tight_layout()
plt.savefig('有効ポテンシャル.png', dpi=150)
plt.show()
print(f'L={L:.4f}, E={E:.4f}, r_min={r_min:.4f}, r_max={r_max:.4f}')
print(f'数値解 r範囲: {r_ary.min():.4f} ~ {r_ary.max():.4f}')
print(f'最大誤差: {np.abs(K_r + W(r_ary) - E).max():.2e}')
結果
プログラムを実行すると、以下のグラフが出力される。
左図は有効ポテンシャル $W(r)$（黒実線）と、それを構成する重力ポテンシャル $-1/r$（灰破線）および遠心ポテンシャル $L^2/2r^2$（灰点線）である。赤実線は力学的エネルギー $E$ を表し、赤実線と $W(r)$ に挟まれたオレンジ色の領域の高さが、各 $r$ における動径方向の運動エネルギー $\frac{1}{2}\dot{r}^2$ に対応する。青点は数値解から逆算した $E-\frac{1}{2}\dot{r}^2$ である。
右図は $\frac{1}{2}\dot{r}^2+W(r)-E$ の時間変化であり、エネルギー保存則からのずれ（数値誤差）を表す。
理論値と数値計算の結果を比較すると、以下の表のようになる。
| 物理量 | 理論値 | 数値計算 |
|---|---|---|
| $r_{\min}$ | 0.6958 | 0.6958 |
| $r_{\max}$ | 1.5000 | 1.5000 |
| $a=(r_{\min}+r_{\max})/2$ | 1.0979 | 1.0979 |
| $W_{\min}$（$r=L^2\approx0.951$） | $-0.526$ | − |
| エネルギー保存の最大誤差 | 0 | $2.5\times10^{-4}$ |
考察
理論と数値計算の一致
青点は $r_{\min}\le r\le r_{\max}$ の範囲で $W(r)$ の曲線上に正確に乗っている。また、数値計算で得られた $r$ の最小値・最大値は、$W(r)=E$ から求めた転回点と小数第4位まで一致した。これにより、2次元の惑星運動が、有効ポテンシャル $W(r)$ の中の1次元運動として正しく記述できることが確認できた。
なお、今回の初期条件では初期位置が $r_{\max}$ と一致している。これは、初速度が $y$ 方向、すなわち動径に垂直な方向に与えられているため、初期時刻で $\dot{r}=0$ となり、出発点そのものが転回点（遠日点）になっているためである。
束縛運動となる理由
惑星が恒星から逃げ去らずに一定の範囲を往復し続けるのは、$E<0$ であり、かつ $W(r)$ が井戸型をしているためである。$r$ が小さい側では遠心ポテンシャルが壁となって恒星への落下を防ぎ、$r$ が大きい側では $W(r)\to 0^-$ であるため、$E<0$ の惑星はそれを越えられない。この2つの壁の間に閉じ込められることが、楕円軌道の本質である。
一方、$E\ge0$ となる場合は $r_{\max}$ が存在せず、惑星は無限遠に飛び去る（放物線軌道または双曲線軌道）。また、$L=0$ の場合は遠心ポテンシャルが消えるため、惑星は恒星に向かってまっすぐ落下する。
数値誤差について
右図のように、エネルギー保存の誤差は $10^{-4}$ 程度の振幅で周期的に振動しているが、時間とともに増大していない。これは、用いた半陰的オイラー法がシンプレクティック積分法であり、長時間の計算においてもエネルギーが系統的にずれないという性質をもつためである。誤差が近日点付近（$r$ が小さく、速度変化が大きい領域）で大きくなるのは、1ステップあたりの局所誤差が加速度の大きさに依存するためと考えられる。
ケプラーの第3法則との関係
理論の節で示したように、長半径は $a=-1/(2E)$ となり、力学的エネルギーのみで決まる。今回の数値計算でも $a=1.0979$ となり、$-1/(2E)$ と一致した。前回の記事では、様々な初期条件の惑星について $T^2\propto a^3$ が成り立つことを確認したが、有効ポテンシャルの観点からは、周期と長半径がいずれも角運動量によらずエネルギーで決まる量であることがわかる。
まとめ
本記事では、万有引力のもとでの惑星の運動について、角運動量保存則を用いて動径方向の1次元運動に帰着させ、有効ポテンシャル
W(r)=-\frac{1}{r}+\frac{L^2}{2r^2}
を導出した。そして、動径方向の運動エネルギー $\frac{1}{2}\dot{r}^2$ と $W(r)$ の和が力学的エネルギーとして保存されることを利用して、数値計算の結果と理論とを比較した。
その結果、数値解から逆算した $E-\frac{1}{2}\dot{r}^2$ は $W(r)$ の曲線に一致し、転回点 $r_{\min}$、$r_{\max}$ も理論値と一致した。有効ポテンシャルを描くことで、惑星がなぜ楕円軌道という束縛運動を行うのかを視覚的に理解することができた。
有効ポテンシャルの考え方は、万有引力に限らず、任意の中心力のもとでの運動に適用できる。今後は、中心力の $r$ 依存性を変えた場合に、$W(r)$ の形状や軌道がどのように変化するかについても調べたい。
また、電磁気学や量子化学の分野でもクーロン力の有効ポテンシャルとして活用される。
参考文献