二体問題は、重心運動と相対運動に完全に分離できます。相対運動は換算質量 μ = m 1 m 2 / ( m 1 + m 2 ) \mu = m_1 m_2/(m_1+m_2) μ = m 1 m 2 / ( m 1 + m 2 ) を持つ質点 1 個の運動と同じ方程式に従います。
力が中心力であること(力が二体を結ぶ直線に沿うこと)だけから、角運動量 L \boldsymbol{L} L が保存 します。その帰結として運動は平面内に限られ、面積速度が一定 になります。これがケプラーの第二法則で、力の大きさが距離のどんな関数であっても成り立ちます。
角運動量保存を使うと、二次元の運動は有効ポテンシャル U e f f ( r ) = U ( r ) + L 2 / ( 2 μ r 2 ) U_{\mathrm{eff}}(r) = U(r) + L^2/(2\mu r^2) U eff ( r ) = U ( r ) + L 2 / ( 2 μ r 2 ) の中の一次元運動に帰着します。軌道の定性的な分類(束縛・非束縛・円軌道)はこの図だけで読み取れます。
軌道の形はビネの方程式 から決まります。逆二乗力 f ( r ) = − k / r 2 f(r) = -k/r^2 f ( r ) = − k / r 2 のときだけ方程式が線形になり、解が円錐曲線 r = ℓ / ( 1 + e cos θ ) r = \ell/(1 + e\cos\theta) r = ℓ / ( 1 + e cos θ ) になります。これがケプラーの第一法則です。離心率はエネルギーと角運動量で e 2 = 1 + 2 E L 2 / ( μ k 2 ) e^2 = 1 + 2EL^2/(\mu k^2) e 2 = 1 + 2 E L 2 / ( μ k 2 ) と表されます。
面積速度一定則と楕円の面積を組み合わせると T 2 = 4 π 2 a 3 / ( G ( m 1 + m 2 ) ) T^2 = 4\pi^2 a^3 / \bigl(G(m_1+m_2)\bigr) T 2 = 4 π 2 a 3 / ( G ( m 1 + m 2 ) ) が出ます。これがケプラーの第三法則で、ケプラー自身の言明にはなかった ( m 1 + m 2 ) (m_1+m_2) ( m 1 + m 2 ) の補正が付きます。
逆二乗力には角運動量とエネルギーのほかにラプラス–ルンゲ–レンツベクトル という保存量があり、これが軌道が閉じる(近日点が動かない)理由です。
ヨハネス・ケプラーは、ティコ・ブラーエが残した火星の観測記録を十数年かけて解析し、三つの規則を取り出しました。
惑星は太陽を一つの焦点とする楕円を描く。
太陽と惑星を結ぶ線分が単位時間に掃く面積は一定である。
公転周期の二乗は軌道長半径の三乗に比例する。
これらは観測の要約 であって、説明ではありません。なぜ楕円なのか、なぜ面積なのか、なぜ二乗と三乗なのか。ケプラー自身はこの問いに答える手段を持っていませんでした。
ニュートンが与えた答えは、三つの規則が独立ではないというものです。第二法則は「力が太陽の方向を向いている」ことだけから出ます。力の強さがどう距離に依存するかには一切依りません。第一法則と第三法則は、そこにさらに「力が距離の二乗に反比例する」という条件を加えると出てきます。つまり、三つの経験則は逆二乗の中心力という一つの仮定にほぼ集約されるのです。
この記事ではその導出を最後まで実行します。使う道具は、ニュートン力学の基礎 で扱った運動方程式(第 2 法則(公理 3.3)[ニュートン力学の基礎] )と、多変数関数の微分と偏微分 程度の解析学だけです。後の章で扱うラグランジュ形式(同じ中心力運動を極座標のラグランジアンから扱う例が 例 5.3[ラグランジュ形式の力学] にあります)やネーターの定理は、ここで手を動かして得た保存則が、実は座標の取り方によらない対称性の帰結であることを教えてくれます。その視点を持って読み返せるように、どの保存則がどの仮定から来ているかを毎回明示します。
ノート
記号の約束です。二つの天体の質量を m 1 , m 2 m_1, m_2 m 1 , m 2 、位置ベクトルを r 1 , r 2 \boldsymbol{r}_1, \boldsymbol{r}_2 r 1 , r 2 とし、相対位置を r = r 1 − r 2 \boldsymbol{r} = \boldsymbol{r}_1 - \boldsymbol{r}_2 r = r 1 − r 2 、その大きさを r = ∣ r ∣ r = |\boldsymbol{r}| r = ∣ r ∣ 、単位ベクトルを r ^ = r / r \hat{\boldsymbol{r}} = \boldsymbol{r}/r r ^ = r / r と書きます。時間微分はドットで表します。
定義 2.1 (中心力 )
質点 1 が質点 2 から受ける力 F \boldsymbol{F} F が、ある実数値関数 f f f を用いて
F = f ( r ) r ^ , r = r 1 − r 2 , r = ∣ r ∣ , r ^ = r / r \boldsymbol{F} = f(r)\,\hat{\boldsymbol{r}},\qquad \boldsymbol{r} = \boldsymbol{r}_1 - \boldsymbol{r}_2,\quad r = |\boldsymbol{r}|,\quad \hat{\boldsymbol{r}} = \boldsymbol{r}/r F = f ( r ) r ^ , r = r 1 − r 2 , r = ∣ r ∣ , r ^ = r / r と書けるとき、この力を中心力 といいます。すなわち、力の向きが二質点を結ぶ直線に沿い、大きさが両者の距離だけで決まる場合です。f ( r ) < 0 f(r) < 0 f ( r ) < 0 のとき引力、f ( r ) > 0 f(r) > 0 f ( r ) > 0 のとき斥力です。
中心力は必ず保存力です。実際、U ( r ) = − ∫ r f ( s ) d s U(r) = -\int^r f(s)\,ds U ( r ) = − ∫ r f ( s ) d s とおけば f ( r ) = − U ′ ( r ) f(r) = -U'(r) f ( r ) = − U ′ ( r ) であり、球対称関数 U ( r ) U(r) U ( r ) の勾配は
∇ U ( r ) = U ′ ( r ) ∇ r = U ′ ( r ) r ^ \nabla U(r) = U'(r)\,\nabla r = U'(r)\,\hat{\boldsymbol{r}} ∇ U ( r ) = U ′ ( r ) ∇ r = U ′ ( r ) r ^
となるので F = − ∇ U \boldsymbol{F} = -\nabla U F = − ∇ U が成り立ちます(∇ r = r ^ \nabla r = \hat{\boldsymbol{r}} ∇ r = r ^ は r = x 2 + y 2 + z 2 r = \sqrt{x^2+y^2+z^2} r = x 2 + y 2 + z 2 を各成分で偏微分すればすぐに確かめられます。たとえば ∂ r / ∂ x = x / r \partial r/\partial x = x/r ∂ r / ∂ x = x / r です)。したがって力学的エネルギーが保存します(定理 7.5[ニュートン力学の基礎] )。
万有引力は f ( r ) = − G m 1 m 2 / r 2 f(r) = -Gm_1m_2/r^2 f ( r ) = − G m 1 m 2 / r 2 、点電荷間のクーロン力は f ( r ) = q 1 q 2 / ( 4 π ε 0 r 2 ) f(r) = q_1q_2/(4\pi\varepsilon_0 r^2) f ( r ) = q 1 q 2 / ( 4 π ε 0 r 2 ) 、等方調和振動子は f ( r ) = − μ ω 2 r f(r) = -\mu\omega^2 r f ( r ) = − μ ω 2 r で、いずれも中心力です。以下では逆二乗引力を
f ( r ) = − k r 2 , U ( r ) = − k r , k > 0 f(r) = -\frac{k}{r^2},\qquad U(r) = -\frac{k}{r},\qquad k > 0 f ( r ) = − r 2 k , U ( r ) = − r k , k > 0
と書きます。万有引力なら k = G m 1 m 2 k = Gm_1m_2 k = G m 1 m 2 です。
定理 2.2 (二体問題の分離 )
質量 m 1 , m 2 m_1, m_2 m 1 , m 2 の二質点が互いに 定義 2.1 の中心力を及ぼし合い、外力が働かないとする。全質量を M = m 1 + m 2 M = m_1+m_2 M = m 1 + m 2 、重心を R = ( m 1 r 1 + m 2 r 2 ) / M \boldsymbol{R} = (m_1\boldsymbol{r}_1+m_2\boldsymbol{r}_2)/M R = ( m 1 r 1 + m 2 r 2 ) / M 、換算質量を
μ = m 1 m 2 m 1 + m 2 \mu = \frac{m_1 m_2}{m_1 + m_2} μ = m 1 + m 2 m 1 m 2 とおくと、次が成り立つ。
R ¨ = 0 \ddot{\boldsymbol{R}} = \boldsymbol{0} R ¨ = 0 。すなわち重心は等速直線運動をする。
相対位置 r = r 1 − r 2 \boldsymbol{r} = \boldsymbol{r}_1 - \boldsymbol{r}_2 r = r 1 − r 2 は μ r ¨ = f ( r ) r ^ \mu\ddot{\boldsymbol{r}} = f(r)\hat{\boldsymbol{r}} μ r ¨ = f ( r ) r ^ を満たす。
全運動エネルギーは T = 1 2 M ∣ R ˙ ∣ 2 + 1 2 μ ∣ r ˙ ∣ 2 T = \tfrac12 M|\dot{\boldsymbol{R}}|^2 + \tfrac12\mu|\dot{\boldsymbol{r}}|^2 T = 2 1 M ∣ R ˙ ∣ 2 + 2 1 μ ∣ r ˙ ∣ 2 と分離する。
証明(定理 2.2) 作用反作用の法則(公理 3.4[ニュートン力学の基礎] )により、質点 1 が受ける力は F = f ( r ) r ^ \boldsymbol{F} = f(r)\hat{\boldsymbol{r}} F = f ( r ) r ^ 、質点 2 が受ける力は − F -\boldsymbol{F} − F です。運動方程式は
m 1 r ¨ 1 = F , m 2 r ¨ 2 = − F . m_1\ddot{\boldsymbol{r}}_1 = \boldsymbol{F},\qquad m_2\ddot{\boldsymbol{r}}_2 = -\boldsymbol{F}. m 1 r ¨ 1 = F , m 2 r ¨ 2 = − F . (1) 二式を辺々加えると m 1 r ¨ 1 + m 2 r ¨ 2 = 0 m_1\ddot{\boldsymbol{r}}_1 + m_2\ddot{\boldsymbol{r}}_2 = \boldsymbol{0} m 1 r ¨ 1 + m 2 r ¨ 2 = 0 です。左辺は M R ¨ M\ddot{\boldsymbol{R}} M R ¨ に等しいので R ¨ = 0 \ddot{\boldsymbol{R}} = \boldsymbol{0} R ¨ = 0 を得ます。
(2) 第一式を m 1 m_1 m 1 で、第二式を m 2 m_2 m 2 で割って差を取ると
r ¨ = r ¨ 1 − r ¨ 2 = F m 1 + F m 2 = ( 1 m 1 + 1 m 2 ) F = F μ \ddot{\boldsymbol{r}} = \ddot{\boldsymbol{r}}_1 - \ddot{\boldsymbol{r}}_2 = \frac{\boldsymbol{F}}{m_1} + \frac{\boldsymbol{F}}{m_2} = \left(\frac{1}{m_1}+\frac{1}{m_2}\right)\boldsymbol{F} = \frac{\boldsymbol{F}}{\mu} r ¨ = r ¨ 1 − r ¨ 2 = m 1 F + m 2 F = ( m 1 1 + m 2 1 ) F = μ F です。最後の等号は 1 / μ = 1 / m 1 + 1 / m 2 1/\mu = 1/m_1 + 1/m_2 1/ μ = 1/ m 1 + 1/ m 2 という換算質量の定義そのものです。両辺に μ \mu μ を掛ければ主張を得ます。
(3) 重心の定義から r 1 = R + ( m 2 / M ) r \boldsymbol{r}_1 = \boldsymbol{R} + (m_2/M)\boldsymbol{r} r 1 = R + ( m 2 / M ) r 、r 2 = R − ( m 1 / M ) r \boldsymbol{r}_2 = \boldsymbol{R} - (m_1/M)\boldsymbol{r} r 2 = R − ( m 1 / M ) r です。これを代入すると
T = 1 2 m 1 ∣ R ˙ + m 2 M r ˙ ∣ 2 + 1 2 m 2 ∣ R ˙ − m 1 M r ˙ ∣ 2 = 1 2 ( m 1 + m 2 ) ∣ R ˙ ∣ 2 + ( m 1 m 2 M − m 2 m 1 M ) ⟨ R ˙ , r ˙ ⟩ + 1 2 m 1 m 2 2 + m 2 m 1 2 M 2 ∣ r ˙ ∣ 2 . \begin{aligned}
T &= \tfrac12 m_1\left|\dot{\boldsymbol{R}} + \tfrac{m_2}{M}\dot{\boldsymbol{r}}\right|^2 + \tfrac12 m_2\left|\dot{\boldsymbol{R}} - \tfrac{m_1}{M}\dot{\boldsymbol{r}}\right|^2 \\
&= \tfrac12 (m_1+m_2)|\dot{\boldsymbol{R}}|^2 + \left(\tfrac{m_1m_2}{M} - \tfrac{m_2m_1}{M}\right)\langle \dot{\boldsymbol{R}}, \dot{\boldsymbol{r}}\rangle + \tfrac12\frac{m_1m_2^2 + m_2m_1^2}{M^2}|\dot{\boldsymbol{r}}|^2 .
\end{aligned} T = 2 1 m 1 R ˙ + M m 2 r ˙ 2 + 2 1 m 2 R ˙ − M m 1 r ˙ 2 = 2 1 ( m 1 + m 2 ) ∣ R ˙ ∣ 2 + ( M m 1 m 2 − M m 2 m 1 ) ⟨ R ˙ , r ˙ ⟩ + 2 1 M 2 m 1 m 2 2 + m 2 m 1 2 ∣ r ˙ ∣ 2 . 交差項の係数は 0 0 0 です。最後の項の係数は m 1 m 2 ( m 2 + m 1 ) / M 2 = m 1 m 2 / M = μ m_1m_2(m_2+m_1)/M^2 = m_1m_2/M = \mu m 1 m 2 ( m 2 + m 1 ) / M 2 = m 1 m 2 / M = μ なので、T = 1 2 M ∣ R ˙ ∣ 2 + 1 2 μ ∣ r ˙ ∣ 2 T = \tfrac12 M|\dot{\boldsymbol{R}}|^2 + \tfrac12\mu|\dot{\boldsymbol{r}}|^2 T = 2 1 M ∣ R ˙ ∣ 2 + 2 1 μ ∣ r ˙ ∣ 2 となります。
∎
この定理により、以後は μ r ¨ = f ( r ) r ^ \mu\ddot{\boldsymbol{r}} = f(r)\hat{\boldsymbol{r}} μ r ¨ = f ( r ) r ^ という一体問題 だけを考えれば十分です。重心系(R ˙ = 0 \dot{\boldsymbol{R}} = \boldsymbol{0} R ˙ = 0 となる慣性系)を取れば、二天体の実際の軌道は相対軌道 r ( t ) \boldsymbol{r}(t) r ( t ) を m 2 / M m_2/M m 2 / M 倍・− m 1 / M -m_1/M − m 1 / M 倍に縮小したものになります。太陽系では m 2 ≫ m 1 m_2 \gg m_1 m 2 ≫ m 1 なので μ ≈ m 1 \mu \approx m_1 μ ≈ m 1 となり、太陽はほぼ静止していると見なせます。しかし連星系のように質量が同程度の場合、両方の星が共通重心のまわりに相似な楕円を描きます。
ここからが本題です。まず、力が中心力であるという仮定「だけ」から何が出るかを見ます。
定理 3.1 (中心力の下での角運動量保存 )
μ r ¨ = f ( r ) r ^ \mu\ddot{\boldsymbol{r}} = f(r)\hat{\boldsymbol{r}} μ r ¨ = f ( r ) r ^ に従う運動に対し、角運動量
L = r × p = μ r × r ˙ \boldsymbol{L} = \boldsymbol{r}\times\boldsymbol{p} = \mu\,\boldsymbol{r}\times\dot{\boldsymbol{r}} L = r × p = μ r × r ˙ は時間によらず一定である。ここで f f f は任意の(連続な)関数でよい。
証明(定理 3.1) 積の微分法則をベクトル積に適用します。
d L d t = μ d d t ( r × r ˙ ) = μ ( r ˙ × r ˙ ) + μ ( r × r ¨ ) . \frac{d\boldsymbol{L}}{dt} = \mu\frac{d}{dt}(\boldsymbol{r}\times\dot{\boldsymbol{r}}) = \mu\,(\dot{\boldsymbol{r}}\times\dot{\boldsymbol{r}}) + \mu\,(\boldsymbol{r}\times\ddot{\boldsymbol{r}}). d t d L = μ d t d ( r × r ˙ ) = μ ( r ˙ × r ˙ ) + μ ( r × r ¨ ) . 第一項は同じベクトル同士のベクトル積なので 0 \boldsymbol{0} 0 です。第二項に 定理 2.2 (2) の運動方程式 μ r ¨ = f ( r ) r ^ \mu\ddot{\boldsymbol{r}} = f(r)\hat{\boldsymbol{r}} μ r ¨ = f ( r ) r ^ を代入すると
r × ( μ r ¨ ) = f ( r ) r × r ^ = f ( r ) r r × r = 0 \boldsymbol{r}\times(\mu\ddot{\boldsymbol{r}}) = f(r)\,\boldsymbol{r}\times\hat{\boldsymbol{r}} = \frac{f(r)}{r}\,\boldsymbol{r}\times\boldsymbol{r} = \boldsymbol{0} r × ( μ r ¨ ) = f ( r ) r × r ^ = r f ( r ) r × r = 0 となります。ここで使ったのは、力が r \boldsymbol{r} r に平行であるという中心力の定義(定義 2.1 )だけです。したがって d L / d t = 0 d\boldsymbol{L}/dt = \boldsymbol{0} d L / d t = 0 です。
∎
保存する理由を一言でいえば、中心力は原点まわりのトルク N = r × F \boldsymbol{N} = \boldsymbol{r}\times\boldsymbol{F} N = r × F を生まないからです。後の章で見るように、これは相対座標系の回転対称性 の帰結であり、対称性と保存則(ネーターの定理) の最も基本的な実例(例 4.4[対称性と保存則] )になっています。
系 3.2 (運動の平面性 )
L ≠ 0 \boldsymbol{L}\neq\boldsymbol{0} L = 0 ならば、運動はすべて原点を通り L \boldsymbol{L} L に垂直な一つの平面内で起こる。L = 0 \boldsymbol{L}=\boldsymbol{0} L = 0 ならば、運動は原点を通る一本の直線上に限られる。
証明(系 3.2) L = μ r × r ˙ \boldsymbol{L} = \mu\,\boldsymbol{r}\times\dot{\boldsymbol{r}} L = μ r × r ˙ はベクトル積なので、常に ⟨ L , r ⟩ = 0 \langle \boldsymbol{L}, \boldsymbol{r}\rangle = 0 ⟨ L , r ⟩ = 0 が成り立ちます。定理 3.1 より L \boldsymbol{L} L は定ベクトルですから、L ≠ 0 \boldsymbol{L}\neq\boldsymbol{0} L = 0 のとき ⟨ L , r ( t ) ⟩ = 0 \langle \boldsymbol{L}, \boldsymbol{r}(t)\rangle = 0 ⟨ L , r ( t )⟩ = 0 は「r ( t ) \boldsymbol{r}(t) r ( t ) が原点を通り法線 L \boldsymbol{L} L を持つ平面上にある」ことを意味します。これがすべての t t t で成り立つので運動は平面運動です。
L = 0 \boldsymbol{L}=\boldsymbol{0} L = 0 の場合は r × r ˙ = 0 \boldsymbol{r}\times\dot{\boldsymbol{r}} = \boldsymbol{0} r × r ˙ = 0 、すなわち r ˙ \dot{\boldsymbol{r}} r ˙ が常に r \boldsymbol{r} r に平行です。このとき r ≠ 0 \boldsymbol{r}\neq\boldsymbol{0} r = 0 である区間で r ^ \hat{\boldsymbol{r}} r ^ の微分を計算すると
d r ^ d t = r ˙ r − r ˙ r 2 r = 1 r ( r ˙ − r ˙ r ^ ) \frac{d\hat{\boldsymbol{r}}}{dt} = \frac{\dot{\boldsymbol{r}}}{r} - \frac{\dot r}{r^2}\boldsymbol{r} = \frac{1}{r}\left(\dot{\boldsymbol{r}} - \dot r\,\hat{\boldsymbol{r}}\right) d t d r ^ = r r ˙ − r 2 r ˙ r = r 1 ( r ˙ − r ˙ r ^ ) であり、r ˙ = λ r \dot{\boldsymbol{r}} = \lambda\boldsymbol{r} r ˙ = λ r と書けることから r ˙ = λ r \dot r = \lambda r r ˙ = λ r 、よって r ˙ − r ˙ r ^ = λ r − λ r r ^ = 0 \dot{\boldsymbol{r}} - \dot r\hat{\boldsymbol{r}} = \lambda\boldsymbol{r} - \lambda r\hat{\boldsymbol{r}} = \boldsymbol{0} r ˙ − r ˙ r ^ = λ r − λ r r ^ = 0 となります。つまり方向 r ^ \hat{\boldsymbol{r}} r ^ は一定で、運動は原点を通る直線上の一次元運動です。
∎
以後 L ≠ 0 \boldsymbol{L}\neq\boldsymbol{0} L = 0 とし、運動平面に極座標 ( r , θ ) (r,\theta) ( r , θ ) を取ります。この平面での位置は r = r e ^ r \boldsymbol{r} = r\,\hat{\boldsymbol{e}}_r r = r e ^ r と書け、単位ベクトルの微分は
d e ^ r d t = θ ˙ e ^ θ , d e ^ θ d t = − θ ˙ e ^ r \frac{d\hat{\boldsymbol{e}}_r}{dt} = \dot\theta\,\hat{\boldsymbol{e}}_\theta,\qquad \frac{d\hat{\boldsymbol{e}}_\theta}{dt} = -\dot\theta\,\hat{\boldsymbol{e}}_r d t d e ^ r = θ ˙ e ^ θ , d t d e ^ θ = − θ ˙ e ^ r
です(e ^ r = ( cos θ , sin θ ) \hat{\boldsymbol{e}}_r = (\cos\theta,\sin\theta) e ^ r = ( cos θ , sin θ ) 、e ^ θ = ( − sin θ , cos θ ) \hat{\boldsymbol{e}}_\theta = (-\sin\theta,\cos\theta) e ^ θ = ( − sin θ , cos θ ) を t t t で微分すればそのまま出ます)。これを用いて速度と加速度を書き下すと
r ˙ = r ˙ e ^ r + r θ ˙ e ^ θ , r ¨ = ( r ¨ − r θ ˙ 2 ) e ^ r + ( r θ ¨ + 2 r ˙ θ ˙ ) e ^ θ \dot{\boldsymbol{r}} = \dot r\,\hat{\boldsymbol{e}}_r + r\dot\theta\,\hat{\boldsymbol{e}}_\theta,\qquad
\ddot{\boldsymbol{r}} = (\ddot r - r\dot\theta^2)\,\hat{\boldsymbol{e}}_r + (r\ddot\theta + 2\dot r\dot\theta)\,\hat{\boldsymbol{e}}_\theta r ˙ = r ˙ e ^ r + r θ ˙ e ^ θ , r ¨ = ( r ¨ − r θ ˙ 2 ) e ^ r + ( r θ ¨ + 2 r ˙ θ ˙ ) e ^ θ
となります。角運動量の大きさは L = ∣ L ∣ = μ ∣ r e ^ r × ( r ˙ e ^ r + r θ ˙ e ^ θ ) ∣ = μ r 2 ∣ θ ˙ ∣ L = |\boldsymbol{L}| = \mu|r\hat{\boldsymbol{e}}_r \times (\dot r\hat{\boldsymbol{e}}_r + r\dot\theta\hat{\boldsymbol{e}}_\theta)| = \mu r^2|\dot\theta| L = ∣ L ∣ = μ ∣ r e ^ r × ( r ˙ e ^ r + r θ ˙ e ^ θ ) ∣ = μ r 2 ∣ θ ˙ ∣ です。θ \theta θ の向きを θ ˙ > 0 \dot\theta > 0 θ ˙ > 0 となるように選んで
L = μ r 2 θ ˙ L = \mu r^2\dot\theta L = μ r 2 θ ˙
と書きます。これが以下でくり返し使う関係式です。
定理 3.3 (ケプラーの第二法則(面積速度一定) )
中心力の下での平面運動において、原点と質点を結ぶ線分が時刻 t 0 t_0 t 0 から t t t までに掃く面積を A ( t ) A(t) A ( t ) とすると
d A d t = 1 2 r 2 θ ˙ = L 2 μ = 一定 \frac{dA}{dt} = \frac{1}{2}r^2\dot\theta = \frac{L}{2\mu} = \text{一定} d t d A = 2 1 r 2 θ ˙ = 2 μ L = 一定 である。とくに、等しい時間に掃く面積は等しい。ここでも f f f は任意の関数でよく、逆二乗性は使わない。
証明(定理 3.3) 微小時間 d t dt d t の間に動径が θ \theta θ から θ + d θ \theta + d\theta θ + d θ まで回り、長さが r r r から r + d r r + dr r + d r に変わるとします。この間に掃かれる領域は、二辺 r r r 、r + d r r+dr r + d r と中心角 d θ d\theta d θ を持つ細い扇形で、その面積は極座標の面積要素 d A = 1 2 r 2 d θ dA = \tfrac12 r^2\,d\theta d A = 2 1 r 2 d θ に O ( d r d θ ) O(dr\,d\theta) O ( d r d θ ) の高次項を加えたものです。厳密には、掃かれる領域の面積は重積分で
A = ∫ θ 0 θ 1 ∫ 0 r ( θ ) ρ d ρ d θ = ∫ θ 0 θ 1 r ( θ ) 2 2 d θ A = \int_{\theta_0}^{\theta_1}\!\!\int_0^{r(\theta)} \rho \,d\rho\,d\theta = \int_{\theta_0}^{\theta_1}\frac{r(\theta)^2}{2}\,d\theta A = ∫ θ 0 θ 1 ∫ 0 r ( θ ) ρ d ρ d θ = ∫ θ 0 θ 1 2 r ( θ ) 2 d θ と書けます(極座標のヤコビアンが ρ \rho ρ であることを使いました。重積分と累次積分 の 例 6.6[重積分と累次積分] を参照してください)。両辺を t t t で微分し、合成関数の微分法を用いると d A / d t = 1 2 r ( θ ) 2 θ ˙ dA/dt = \tfrac12 r(\theta)^2\,\dot\theta d A / d t = 2 1 r ( θ ) 2 θ ˙ です。
ここに 定理 3.1 から得た L = μ r 2 θ ˙ L = \mu r^2\dot\theta L = μ r 2 θ ˙ を代入すると r 2 θ ˙ = L / μ r^2\dot\theta = L/\mu r 2 θ ˙ = L / μ なので
d A d t = L 2 μ \frac{dA}{dt} = \frac{L}{2\mu} d t d A = 2 μ L となり、L L L と μ \mu μ が定数であることから面積速度は一定です。
∎
ケプラーが火星の観測から抽出した第二法則が、力の距離依存性を一切使わずに出てしまったことに注意してください。第二法則は「太陽が引く」という事実の幾何学的言い換えにすぎず、逆二乗則の証拠にはなりません。第二法則が破れて見えたら、それは力が中心力でない(太陽以外の天体の摂動がある、あるいは相対論的効果がある)という信号です。
flowchart TD
A["力が中心力:F = f(r) r̂"] --> B["トルク r × F = 0"]
B --> C["角運動量 L が保存"]
C --> D["運動は平面内(L ≠ 0)"]
C --> E["面積速度 L/2μ 一定=ケプラー第2法則"]
A --> F["F は保存力:U(r) が存在"]
F --> G["力学的エネルギー E が保存"]
C --> H["有効ポテンシャル U + L²/2μr² の一次元運動"]
G --> H
H --> I["f(r) = -k/r² を追加"]
I --> J["円錐曲線軌道=ケプラー第1・第3法則"] 中心力の運動で保存する量と、それぞれが依拠する仮定
保存則を二つ手に入れたので、自由度を落とします。
定義 4.1 (有効ポテンシャル )
角運動量の大きさ L L L を固定したとき、
U e f f ( r ) = U ( r ) + L 2 2 μ r 2 U_{\mathrm{eff}}(r) = U(r) + \frac{L^2}{2\mu r^2} U eff ( r ) = U ( r ) + 2 μ r 2 L 2 を有効ポテンシャル といいます。第二項 L 2 / ( 2 μ r 2 ) L^2/(2\mu r^2) L 2 / ( 2 μ r 2 ) を遠心障壁 と呼びます。
命題 4.2 (動径方向の一次元問題への帰着 )
中心力の下での平面運動において、力学的エネルギー
E = 1 2 μ ∣ r ˙ ∣ 2 + U ( r ) E = \frac{1}{2}\mu|\dot{\boldsymbol{r}}|^2 + U(r) E = 2 1 μ ∣ r ˙ ∣ 2 + U ( r ) は保存し、さらに
E = 1 2 μ r ˙ 2 + U e f f ( r ) E = \frac{1}{2}\mu\dot r^2 + U_{\mathrm{eff}}(r) E = 2 1 μ r ˙ 2 + U eff ( r ) と書ける。すなわち r ( t ) r(t) r ( t ) は、ポテンシャル U e f f U_{\mathrm{eff}} U eff の中を動く質量 μ \mu μ の質点の一次元運動と同じ方程式に従う。
証明(命題 4.2) まずエネルギー保存を確かめます。μ r ¨ = − ∇ U \mu\ddot{\boldsymbol{r}} = -\nabla U μ r ¨ = − ∇ U の両辺と r ˙ \dot{\boldsymbol{r}} r ˙ の内積を取ると
μ ⟨ r ¨ , r ˙ ⟩ = − ⟨ ∇ U , r ˙ ⟩ ⟺ d d t ( 1 2 μ ∣ r ˙ ∣ 2 ) = − d d t U ( r ( t ) ) \mu\langle\ddot{\boldsymbol{r}},\dot{\boldsymbol{r}}\rangle = -\langle\nabla U,\dot{\boldsymbol{r}}\rangle
\;\Longleftrightarrow\;
\frac{d}{dt}\left(\frac{1}{2}\mu|\dot{\boldsymbol{r}}|^2\right) = -\frac{d}{dt}U(r(t)) μ ⟨ r ¨ , r ˙ ⟩ = − ⟨ ∇ U , r ˙ ⟩ ⟺ d t d ( 2 1 μ ∣ r ˙ ∣ 2 ) = − d t d U ( r ( t )) となり(右辺は合成関数の微分法)、E E E の時間微分が 0 0 0 です。ここで使ったのは §2.1 で示した「中心力は保存力である」という事実です。
次に速度の分解 r ˙ = r ˙ e ^ r + r θ ˙ e ^ θ \dot{\boldsymbol{r}} = \dot r\,\hat{\boldsymbol{e}}_r + r\dot\theta\,\hat{\boldsymbol{e}}_\theta r ˙ = r ˙ e ^ r + r θ ˙ e ^ θ を用います。e ^ r \hat{\boldsymbol{e}}_r e ^ r と e ^ θ \hat{\boldsymbol{e}}_\theta e ^ θ は互いに直交する単位ベクトルなので
∣ r ˙ ∣ 2 = r ˙ 2 + r 2 θ ˙ 2 . |\dot{\boldsymbol{r}}|^2 = \dot r^2 + r^2\dot\theta^2 . ∣ r ˙ ∣ 2 = r ˙ 2 + r 2 θ ˙ 2 . 定理 3.1 の関係 θ ˙ = L / ( μ r 2 ) \dot\theta = L/(\mu r^2) θ ˙ = L / ( μ r 2 ) を代入すると r 2 θ ˙ 2 = r 2 ⋅ L 2 / ( μ 2 r 4 ) = L 2 / ( μ 2 r 2 ) r^2\dot\theta^2 = r^2\cdot L^2/(\mu^2r^4) = L^2/(\mu^2r^2) r 2 θ ˙ 2 = r 2 ⋅ L 2 / ( μ 2 r 4 ) = L 2 / ( μ 2 r 2 ) です。よって
E = 1 2 μ r ˙ 2 + L 2 2 μ r 2 + U ( r ) = 1 2 μ r ˙ 2 + U e f f ( r ) E = \frac{1}{2}\mu\dot r^2 + \frac{L^2}{2\mu r^2} + U(r) = \frac{1}{2}\mu\dot r^2 + U_{\mathrm{eff}}(r) E = 2 1 μ r ˙ 2 + 2 μ r 2 L 2 + U ( r ) = 2 1 μ r ˙ 2 + U eff ( r ) を得ます。定義 4.1 の定義を使いました。
∎
この帰着は強力です。r ˙ 2 ≥ 0 \dot r^2 \ge 0 r ˙ 2 ≥ 0 より運動は U e f f ( r ) ≤ E U_{\mathrm{eff}}(r)\le E U eff ( r ) ≤ E を満たす r r r の範囲に限られ、等号が成り立つ点(転回点 )で r ˙ = 0 \dot r = 0 r ˙ = 0 となります。逆二乗引力 U = − k / r U = -k/r U = − k / r の場合、
U e f f ( r ) = − k r + L 2 2 μ r 2 U_{\mathrm{eff}}(r) = -\frac{k}{r} + \frac{L^2}{2\mu r^2} U eff ( r ) = − r k + 2 μ r 2 L 2
は r → 0 + r\to 0^+ r → 0 + で + ∞ +\infty + ∞ (遠心障壁が勝つ)、r → ∞ r\to\infty r → ∞ で 0 − 0^- 0 − に近づき、途中でただ一つの極小を持ちます。極小の位置は U e f f ′ ( r ) = k / r 2 − L 2 / ( μ r 3 ) = 0 U_{\mathrm{eff}}'(r) = k/r^2 - L^2/(\mu r^3) = 0 U eff ′ ( r ) = k / r 2 − L 2 / ( μ r 3 ) = 0 より
r c = L 2 μ k , U e f f ( r c ) = − μ k 2 2 L 2 r_c = \frac{L^2}{\mu k},\qquad U_{\mathrm{eff}}(r_c) = -\frac{\mu k^2}{2L^2} r c = μ k L 2 , U eff ( r c ) = − 2 L 2 μ k 2
です。
エネルギー r U_eff E が負:楕円(束縛運動) E = 0:放物線 E が正:双曲線 円軌道 r_c = L²/μk 逆二乗引力の有効ポテンシャル。エネルギーの高さで軌道の型が決まる
図から読み取れることを整理します。E = U e f f ( r c ) = − μ k 2 / ( 2 L 2 ) E = U_{\mathrm{eff}}(r_c) = -\mu k^2/(2L^2) E = U eff ( r c ) = − μ k 2 / ( 2 L 2 ) のとき r r r は動けず円軌道 です。U e f f ( r c ) < E < 0 U_{\mathrm{eff}}(r_c) < E < 0 U eff ( r c ) < E < 0 のとき r r r は二つの転回点 r min , r max r_{\min}, r_{\max} r m i n , r m a x の間を往復し、軌道は原点から有限の範囲に留まります(束縛運動 )。E ≥ 0 E\ge 0 E ≥ 0 のときは転回点が一つだけで、質点は最接近後に無限遠へ去ります(非束縛運動 )。この分類がそのまま楕円・放物線・双曲線に対応することを §5 で確認します。
例 4.3 (べき乗則の中心力における円軌道の安定性 )
f ( r ) = − k / r n f(r) = -k/r^n f ( r ) = − k / r n (k > 0 k>0 k > 0 、n ≠ 1 n\neq 1 n = 1 )という引力を考えます。対応するポテンシャルは U ( r ) = − k / ( ( n − 1 ) r n − 1 ) U(r) = -k/\bigl((n-1)r^{n-1}\bigr) U ( r ) = − k / ( ( n − 1 ) r n − 1 ) で、有効ポテンシャルは
U e f f ( r ) = − k ( n − 1 ) r n − 1 + L 2 2 μ r 2 U_{\mathrm{eff}}(r) = -\frac{k}{(n-1)r^{n-1}} + \frac{L^2}{2\mu r^2} U eff ( r ) = − ( n − 1 ) r n − 1 k + 2 μ r 2 L 2 です。円軌道は U e f f U_{\mathrm{eff}} U eff の停留点に対応します。
U e f f ′ ( r ) = k r n − L 2 μ r 3 = 0 ⟺ r c 3 − n = L 2 μ k U_{\mathrm{eff}}'(r) = \frac{k}{r^n} - \frac{L^2}{\mu r^3} = 0 \iff r_c^{\,3-n} = \frac{L^2}{\mu k} U eff ′ ( r ) = r n k − μ r 3 L 2 = 0 ⟺ r c 3 − n = μ k L 2 より、n ≠ 3 n\neq 3 n = 3 なら各 L L L に対して円軌道半径 r c r_c r c がただ一つ定まります。安定性は二階微分の符号で判定します。
U e f f ′ ′ ( r ) = − n k r n + 1 + 3 L 2 μ r 4 . U_{\mathrm{eff}}''(r) = -\frac{nk}{r^{n+1}} + \frac{3L^2}{\mu r^4}. U eff ′′ ( r ) = − r n + 1 nk + μ r 4 3 L 2 . 停留点の条件 L 2 / μ = k r c 3 − n L^2/\mu = k\,r_c^{\,3-n} L 2 / μ = k r c 3 − n を第二項に代入すると 3 L 2 / ( μ r c 4 ) = 3 k r c − 1 − n 3L^2/(\mu r_c^4) = 3k\,r_c^{\,-1-n} 3 L 2 / ( μ r c 4 ) = 3 k r c − 1 − n なので
U e f f ′ ′ ( r c ) = ( 3 − n ) k r c n + 1 . U_{\mathrm{eff}}''(r_c) = \frac{(3-n)k}{r_c^{\,n+1}} . U eff ′′ ( r c ) = r c n + 1 ( 3 − n ) k . k > 0 k>0 k > 0 、r c > 0 r_c>0 r c > 0 なので、これが正になるのは n < 3 n < 3 n < 3 のときに限ります。つまり逆三乗より急な引力では円軌道が不安定 で、わずかな摂動で質点は中心へ落ち込むか無限遠へ飛び去ります。逆二乗力は n = 2 n=2 n = 2 なので U e f f ′ ′ ( r c ) = k / r c 3 > 0 U_{\mathrm{eff}}''(r_c) = k/r_c^3 > 0 U eff ′′ ( r c ) = k / r c 3 > 0 となり、円軌道は安定です。この安定性が惑星系が存続できる理由の一つです。
ここまでは r ( t ) r(t) r ( t ) の時間発展を追ってきました。しかし軌道の形 を知りたいなら、時間を消去して r r r を θ \theta θ の関数として求めるほうが早道です。鍵になるのが変数変換 u = 1 / r u = 1/r u = 1/ r です。
補題 5.1 (ビネの軌道方程式 )
L ≠ 0 \boldsymbol{L}\neq\boldsymbol{0} L = 0 とし、u ( θ ) = 1 / r ( θ ) u(\theta) = 1/r(\theta) u ( θ ) = 1/ r ( θ ) とおく。中心力 f ( r ) f(r) f ( r ) の下での軌道は
d 2 u d θ 2 + u = − μ L 2 u 2 f ( 1 u ) \frac{d^2u}{d\theta^2} + u = -\frac{\mu}{L^2u^2}\,f\!\left(\frac{1}{u}\right) d θ 2 d 2 u + u = − L 2 u 2 μ f ( u 1 ) を満たす。
証明(補題 5.1) L = μ r 2 θ ˙ L = \mu r^2\dot\theta L = μ r 2 θ ˙ より θ ˙ = L u 2 / μ \dot\theta = Lu^2/\mu θ ˙ = L u 2 / μ です。L ≠ 0 L\neq 0 L = 0 なので θ ˙ \dot\theta θ ˙ は符号を変えず、θ \theta θ を独立変数に取り直せます。
まず r ˙ \dot r r ˙ を θ \theta θ 微分で書き換えます。r = 1 / u r = 1/u r = 1/ u なので、合成関数の微分法により
r ˙ = d d t ( 1 u ) = − 1 u 2 d u d θ θ ˙ = − 1 u 2 d u d θ ⋅ L u 2 μ = − L μ d u d θ . \dot r = \frac{d}{dt}\left(\frac{1}{u}\right) = -\frac{1}{u^2}\frac{du}{d\theta}\dot\theta = -\frac{1}{u^2}\frac{du}{d\theta}\cdot\frac{Lu^2}{\mu} = -\frac{L}{\mu}\frac{du}{d\theta}. r ˙ = d t d ( u 1 ) = − u 2 1 d θ d u θ ˙ = − u 2 1 d θ d u ⋅ μ L u 2 = − μ L d θ d u . u 2 u^2 u 2 がきれいに約分される点がこの変換の要です。もう一度微分すると
r ¨ = − L μ d 2 u d θ 2 θ ˙ = − L μ d 2 u d θ 2 ⋅ L u 2 μ = − L 2 u 2 μ 2 d 2 u d θ 2 . \ddot r = -\frac{L}{\mu}\frac{d^2u}{d\theta^2}\dot\theta = -\frac{L}{\mu}\frac{d^2u}{d\theta^2}\cdot\frac{Lu^2}{\mu} = -\frac{L^2u^2}{\mu^2}\frac{d^2u}{d\theta^2}. r ¨ = − μ L d θ 2 d 2 u θ ˙ = − μ L d θ 2 d 2 u ⋅ μ L u 2 = − μ 2 L 2 u 2 d θ 2 d 2 u . 一方、§3 で求めた加速度の e ^ r \hat{\boldsymbol{e}}_r e ^ r 成分から、運動方程式の動径成分は
μ ( r ¨ − r θ ˙ 2 ) = f ( r ) \mu(\ddot r - r\dot\theta^2) = f(r) μ ( r ¨ − r θ ˙ 2 ) = f ( r ) です。ここで
r θ ˙ 2 = 1 u ( L u 2 μ ) 2 = L 2 u 3 μ 2 r\dot\theta^2 = \frac{1}{u}\left(\frac{Lu^2}{\mu}\right)^2 = \frac{L^2u^3}{\mu^2} r θ ˙ 2 = u 1 ( μ L u 2 ) 2 = μ 2 L 2 u 3 なので、代入して
μ ( − L 2 u 2 μ 2 d 2 u d θ 2 − L 2 u 3 μ 2 ) = f ( 1 / u ) ⟺ − L 2 u 2 μ ( d 2 u d θ 2 + u ) = f ( 1 / u ) \mu\left(-\frac{L^2u^2}{\mu^2}\frac{d^2u}{d\theta^2} - \frac{L^2u^3}{\mu^2}\right) = f(1/u)
\;\Longleftrightarrow\;
-\frac{L^2u^2}{\mu}\left(\frac{d^2u}{d\theta^2} + u\right) = f(1/u) μ ( − μ 2 L 2 u 2 d θ 2 d 2 u − μ 2 L 2 u 3 ) = f ( 1/ u ) ⟺ − μ L 2 u 2 ( d θ 2 d 2 u + u ) = f ( 1/ u ) を得ます。u ≠ 0 u\neq 0 u = 0 (r r r は有限)なので両辺を − L 2 u 2 / μ -L^2u^2/\mu − L 2 u 2 / μ で割れば主張の式になります。
∎
補題 5.1 の右辺は一般には u u u の非線形関数です。ところが f ( r ) = − k / r 2 f(r) = -k/r^2 f ( r ) = − k / r 2 、すなわち f ( 1 / u ) = − k u 2 f(1/u) = -ku^2 f ( 1/ u ) = − k u 2 のときに限り u 2 u^2 u 2 が約分されて右辺が定数 になります。逆二乗則が特別なのはこの一点です。
定理 5.2 (ケプラーの第一法則(軌道は円錐曲線) )
逆二乗引力 f ( r ) = − k / r 2 f(r) = -k/r^2 f ( r ) = − k / r 2 (k > 0 k>0 k > 0 )の下で、L ≠ 0 \boldsymbol{L}\neq\boldsymbol{0} L = 0 の運動の軌道は
r ( θ ) = ℓ 1 + e cos ( θ − θ 0 ) , ℓ = L 2 μ k , e ≥ 0 r(\theta) = \frac{\ell}{1 + e\cos(\theta-\theta_0)},\qquad \ell = \frac{L^2}{\mu k},\quad e \ge 0 r ( θ ) = 1 + e cos ( θ − θ 0 ) ℓ , ℓ = μ k L 2 , e ≥ 0 で与えられる。これは力の中心(原点)を一つの焦点とする円錐曲線であり、e < 1 e<1 e < 1 なら楕円、e = 1 e=1 e = 1 なら放物線、e > 1 e>1 e > 1 なら双曲線の一方の分枝である。とくに束縛運動(e < 1 e<1 e < 1 )では軌道は楕円で、力の中心はその焦点にある。
証明(定理 5.2) 補題 5.1 に f ( 1 / u ) = − k u 2 f(1/u) = -ku^2 f ( 1/ u ) = − k u 2 を代入します。
d 2 u d θ 2 + u = − μ L 2 u 2 ⋅ ( − k u 2 ) = μ k L 2 = 1 ℓ . \frac{d^2u}{d\theta^2} + u = -\frac{\mu}{L^2u^2}\cdot(-ku^2) = \frac{\mu k}{L^2} = \frac{1}{\ell}. d θ 2 d 2 u + u = − L 2 u 2 μ ⋅ ( − k u 2 ) = L 2 μ k = ℓ 1 . これは定数係数二階線形非同次常微分方程式です。特殊解は定数関数 u p = 1 / ℓ u_p = 1/\ell u p = 1/ ℓ で、同次方程式 u ′ ′ + u = 0 u'' + u = 0 u ′′ + u = 0 の一般解は C cos ( θ − θ 0 ) C\cos(\theta-\theta_0) C cos ( θ − θ 0 ) (C ≥ 0 C\ge 0 C ≥ 0 、θ 0 \theta_0 θ 0 は定数。定理 5.2[ニュートン力学の基礎] で角振動数を 1 1 1 、時間変数を θ \theta θ に読み替えたものです)ですから、一般解は
u ( θ ) = 1 ℓ + C cos ( θ − θ 0 ) u(\theta) = \frac{1}{\ell} + C\cos(\theta-\theta_0) u ( θ ) = ℓ 1 + C cos ( θ − θ 0 ) です。e = C ℓ ≥ 0 e = C\ell \ge 0 e = C ℓ ≥ 0 とおいて逆数を取ると
r ( θ ) = 1 u ( θ ) = ℓ 1 + e cos ( θ − θ 0 ) r(\theta) = \frac{1}{u(\theta)} = \frac{\ell}{1 + e\cos(\theta-\theta_0)} r ( θ ) = u ( θ ) 1 = 1 + e cos ( θ − θ 0 ) ℓ となります。以後 θ 0 = 0 \theta_0 = 0 θ 0 = 0 となるように θ \theta θ の基準を選びます。
これが円錐曲線であることを確かめます。r + e r cos θ = ℓ r + er\cos\theta = \ell r + er cos θ = ℓ を直交座標 x = r cos θ x = r\cos\theta x = r cos θ 、y = r sin θ y = r\sin\theta y = r sin θ で書くと x 2 + y 2 = ℓ − e x \sqrt{x^2+y^2} = \ell - ex x 2 + y 2 = ℓ − e x で、両辺を二乗して
x 2 + y 2 = ℓ 2 − 2 e ℓ x + e 2 x 2 ⟺ ( 1 − e 2 ) x 2 + 2 e ℓ x + y 2 = ℓ 2 . x^2 + y^2 = \ell^2 - 2e\ell x + e^2x^2 \iff (1-e^2)x^2 + 2e\ell x + y^2 = \ell^2 . x 2 + y 2 = ℓ 2 − 2 e ℓ x + e 2 x 2 ⟺ ( 1 − e 2 ) x 2 + 2 e ℓ x + y 2 = ℓ 2 . e < 1 e<1 e < 1 のとき x x x について平方完成すると
( 1 − e 2 ) ( x + e ℓ 1 − e 2 ) 2 + y 2 = ℓ 2 + e 2 ℓ 2 1 − e 2 = ℓ 2 1 − e 2 (1-e^2)\left(x + \frac{e\ell}{1-e^2}\right)^2 + y^2 = \ell^2 + \frac{e^2\ell^2}{1-e^2} = \frac{\ell^2}{1-e^2} ( 1 − e 2 ) ( x + 1 − e 2 e ℓ ) 2 + y 2 = ℓ 2 + 1 − e 2 e 2 ℓ 2 = 1 − e 2 ℓ 2 となり、両辺を右辺で割れば
( x + a e ) 2 a 2 + y 2 b 2 = 1 , a = ℓ 1 − e 2 , b = ℓ 1 − e 2 = a 1 − e 2 \frac{\bigl(x + a e\bigr)^2}{a^2} + \frac{y^2}{b^2} = 1,\qquad a = \frac{\ell}{1-e^2},\quad b = \frac{\ell}{\sqrt{1-e^2}} = a\sqrt{1-e^2} a 2 ( x + a e ) 2 + b 2 y 2 = 1 , a = 1 − e 2 ℓ , b = 1 − e 2 ℓ = a 1 − e 2 という楕円の標準形が得られます。中心は ( − a e , 0 ) (-ae, 0) ( − a e , 0 ) にあり、原点はそこから距離 a e ae a e だけ離れています。a 2 − b 2 = a 2 e 2 a^2 - b^2 = a^2e^2 a 2 − b 2 = a 2 e 2 なので、この距離はまさに焦点距離であり、原点は楕円の焦点の一つです。e = 1 e=1 e = 1 のときは x 2 x^2 x 2 の項が消えて y 2 = ℓ 2 − 2 ℓ x y^2 = \ell^2 - 2\ell x y 2 = ℓ 2 − 2 ℓ x という放物線、e > 1 e>1 e > 1 のときは 1 − e 2 < 0 1-e^2<0 1 − e 2 < 0 となり同様の変形で双曲線の標準形になります。
∎
定義 5.3 (軌道要素 )
定理 5.2 に現れた量を次のように呼びます。e e e を離心率 、ℓ = L 2 / ( μ k ) \ell = L^2/(\mu k) ℓ = L 2 / ( μ k ) を半直弦 (セミラタスレクタム)、楕円軌道における a = ℓ / ( 1 − e 2 ) a = \ell/(1-e^2) a = ℓ / ( 1 − e 2 ) を軌道長半径 、b = a 1 − e 2 b = a\sqrt{1-e^2} b = a 1 − e 2 を軌道短半径 といいます。θ = 0 \theta=0 θ = 0 で r r r が最小となる点を近点 (太陽まわりなら近日点)、θ = π \theta=\pi θ = π で r r r が最大となる点を遠点 といい、
r min = ℓ 1 + e = a ( 1 − e ) , r max = ℓ 1 − e = a ( 1 + e ) r_{\min} = \frac{\ell}{1+e} = a(1-e),\qquad r_{\max} = \frac{\ell}{1-e} = a(1+e) r m i n = 1 + e ℓ = a ( 1 − e ) , r m a x = 1 − e ℓ = a ( 1 + e ) です。とくに ℓ = a ( 1 − e 2 ) = b 2 / a \ell = a(1-e^2) = b^2/a ℓ = a ( 1 − e 2 ) = b 2 / a が成り立ちます。
楕円軌道の幾何。太陽は楕円の中心ではなく焦点にある
系 5.4 (離心率とエネルギーの関係 )
定理 5.2 の軌道について、力学的エネルギー E E E と角運動量の大きさ L L L の間に
e 2 = 1 + 2 E L 2 μ k 2 e^2 = 1 + \frac{2EL^2}{\mu k^2} e 2 = 1 + μ k 2 2 E L 2 が成り立つ。したがって E < 0 ⟺ e < 1 E<0 \iff e<1 E < 0 ⟺ e < 1 (楕円)、E = 0 ⟺ e = 1 E=0\iff e=1 E = 0 ⟺ e = 1 (放物線)、E > 0 ⟺ e > 1 E>0\iff e>1 E > 0 ⟺ e > 1 (双曲線)である。さらに楕円軌道では
a = − k 2 E a = -\frac{k}{2E} a = − 2 E k であり、軌道長半径はエネルギーだけで決まる。
証明(系 5.4) 定理 5.2 の証明中の一般解を u = 1 / ℓ + C cos θ u = 1/\ell + C\cos\theta u = 1/ ℓ + C cos θ (C = e / ℓ C = e/\ell C = e / ℓ )と書きます。補題 5.1 の証明で得た r ˙ = − ( L / μ ) d u / d θ \dot r = -(L/\mu)\,du/d\theta r ˙ = − ( L / μ ) d u / d θ より
r ˙ = − L μ ⋅ ( − C sin θ ) = L C μ sin θ . \dot r = -\frac{L}{\mu}\cdot(-C\sin\theta) = \frac{LC}{\mu}\sin\theta . r ˙ = − μ L ⋅ ( − C sin θ ) = μ L C sin θ . これを 命題 4.2 のエネルギー表式に入れます。U = − k / r = − k u U = -k/r = -ku U = − k / r = − k u に注意して
E = 1 2 μ r ˙ 2 + L 2 u 2 2 μ − k u = L 2 C 2 2 μ sin 2 θ + L 2 2 μ ( 1 ℓ + C cos θ ) 2 − k ( 1 ℓ + C cos θ ) = L 2 C 2 2 μ sin 2 θ + L 2 2 μ ( 1 ℓ 2 + 2 C cos θ ℓ + C 2 cos 2 θ ) − k ℓ − k C cos θ . \begin{aligned}
E &= \frac{1}{2}\mu\dot r^2 + \frac{L^2u^2}{2\mu} - ku \\
&= \frac{L^2C^2}{2\mu}\sin^2\theta + \frac{L^2}{2\mu}\left(\frac{1}{\ell}+C\cos\theta\right)^2 - k\left(\frac{1}{\ell}+C\cos\theta\right) \\
&= \frac{L^2C^2}{2\mu}\sin^2\theta + \frac{L^2}{2\mu}\left(\frac{1}{\ell^2} + \frac{2C\cos\theta}{\ell} + C^2\cos^2\theta\right) - \frac{k}{\ell} - kC\cos\theta .
\end{aligned} E = 2 1 μ r ˙ 2 + 2 μ L 2 u 2 − k u = 2 μ L 2 C 2 sin 2 θ + 2 μ L 2 ( ℓ 1 + C cos θ ) 2 − k ( ℓ 1 + C cos θ ) = 2 μ L 2 C 2 sin 2 θ + 2 μ L 2 ( ℓ 2 1 + ℓ 2 C cos θ + C 2 cos 2 θ ) − ℓ k − k C cos θ . ここで ℓ = L 2 / ( μ k ) \ell = L^2/(\mu k) ℓ = L 2 / ( μ k ) すなわち L 2 / ( μ ℓ ) = k L^2/(\mu\ell) = k L 2 / ( μ ℓ ) = k を使うと、cos θ \cos\theta cos θ に比例する項は L 2 μ C cos θ ℓ − k C cos θ = k C cos θ − k C cos θ = 0 \dfrac{L^2}{\mu}\dfrac{C\cos\theta}{\ell} - kC\cos\theta = kC\cos\theta - kC\cos\theta = 0 μ L 2 ℓ C cos θ − k C cos θ = k C cos θ − k C cos θ = 0 と打ち消し合います。また sin 2 θ + cos 2 θ = 1 \sin^2\theta + \cos^2\theta = 1 sin 2 θ + cos 2 θ = 1 より C 2 C^2 C 2 の項がまとまり、定数項は L 2 2 μ ℓ 2 − k ℓ = k 2 ℓ − k ℓ = − k 2 ℓ = − μ k 2 2 L 2 \dfrac{L^2}{2\mu\ell^2} - \dfrac{k}{\ell} = \dfrac{k}{2\ell} - \dfrac{k}{\ell} = -\dfrac{k}{2\ell} = -\dfrac{\mu k^2}{2L^2} 2 μ ℓ 2 L 2 − ℓ k = 2 ℓ k − ℓ k = − 2 ℓ k = − 2 L 2 μ k 2 です。結局
E = L 2 C 2 2 μ − μ k 2 2 L 2 . E = \frac{L^2C^2}{2\mu} - \frac{\mu k^2}{2L^2}. E = 2 μ L 2 C 2 − 2 L 2 μ k 2 . C = e / ℓ = e μ k / L 2 C = e/\ell = e\mu k/L^2 C = e / ℓ = e μ k / L 2 を代入すると L 2 C 2 / ( 2 μ ) = μ k 2 e 2 / ( 2 L 2 ) L^2C^2/(2\mu) = \mu k^2e^2/(2L^2) L 2 C 2 / ( 2 μ ) = μ k 2 e 2 / ( 2 L 2 ) なので
E = μ k 2 2 L 2 ( e 2 − 1 ) ⟺ e 2 = 1 + 2 E L 2 μ k 2 E = \frac{\mu k^2}{2L^2}\left(e^2 - 1\right)
\;\Longleftrightarrow\;
e^2 = 1 + \frac{2EL^2}{\mu k^2} E = 2 L 2 μ k 2 ( e 2 − 1 ) ⟺ e 2 = 1 + μ k 2 2 E L 2 を得ます。μ , k , L 2 \mu, k, L^2 μ , k , L 2 はすべて正なので E E E と e 2 − 1 e^2-1 e 2 − 1 の符号は一致し、型の判定が従います。
楕円の場合、定義 5.3 より a = ℓ / ( 1 − e 2 ) a = \ell/(1-e^2) a = ℓ / ( 1 − e 2 ) で、いま 1 − e 2 = − 2 E L 2 / ( μ k 2 ) 1-e^2 = -2EL^2/(\mu k^2) 1 − e 2 = − 2 E L 2 / ( μ k 2 ) ですから
a = L 2 μ k ⋅ μ k 2 − 2 E L 2 = − k 2 E a = \frac{L^2}{\mu k}\cdot\frac{\mu k^2}{-2EL^2} = -\frac{k}{2E} a = μ k L 2 ⋅ − 2 E L 2 μ k 2 = − 2 E k となります(E < 0 E<0 E < 0 なので a > 0 a>0 a > 0 です)。
∎
この系は実用上とても便利です。角運動量は軌道の「細さ」(離心率)だけを、エネルギーは軌道の「大きさ」(長半径)だけを決めます。円軌道は e = 0 e=0 e = 0 すなわち E = − μ k 2 / ( 2 L 2 ) E = -\mu k^2/(2L^2) E = − μ k 2 / ( 2 L 2 ) の場合で、これは §4 で求めた U e f f U_{\mathrm{eff}} U eff の最小値と一致します。二つの独立な導出が合致したことになります。
例 5.5 (恒星間天体の双曲線軌道 )
2017 年に発見された 1I/ʻOumuamua は、離心率 e ≈ 1.20 e \approx 1.20 e ≈ 1.20 、近日点距離 q ≈ 0.255 a u q \approx 0.255\ \mathrm{au} q ≈ 0.255 au の軌道で太陽系を通過しました。e > 1 e>1 e > 1 なので 系 5.4 により E > 0 E>0 E > 0 、すなわち太陽に束縛されていません。
r → ∞ r\to\infty r → ∞ となるのは分母が 0 0 0 になるとき、つまり cos θ ∞ = − 1 / e = − 0.833 \cos\theta_\infty = -1/e = -0.833 cos θ ∞ = − 1/ e = − 0.833 より θ ∞ = 146.4 ∘ \theta_\infty = 146.4^\circ θ ∞ = 146. 4 ∘ です。入射方向と射出方向のなす角(軌道の曲がり角)は 2 θ ∞ − 180 ∘ = 112.8 ∘ 2\theta_\infty - 180^\circ = 112.8^\circ 2 θ ∞ − 18 0 ∘ = 112. 8 ∘ で、太陽の重力によって進行方向が 113 113 113 度ほど曲げられた計算になります。
無限遠での速さ v ∞ v_\infty v ∞ を求めます。q = a ′ ( e − 1 ) q = a'(e-1) q = a ′ ( e − 1 ) (a ′ = ∣ a ∣ = k / ( 2 E ) a' = |a| = k/(2E) a ′ = ∣ a ∣ = k / ( 2 E ) は双曲線の実半軸)より a ′ = 0.255 / 0.20 = 1.275 a u = 1.907 × 10 11 m a' = 0.255/0.20 = 1.275\ \mathrm{au} = 1.907\times10^{11}\ \mathrm{m} a ′ = 0.255/0.20 = 1.275 au = 1.907 × 1 0 11 m です。E = 1 2 μ v ∞ 2 E = \tfrac12\mu v_\infty^2 E = 2 1 μ v ∞ 2 と E = k / ( 2 a ′ ) E = k/(2a') E = k / ( 2 a ′ ) (e > 1 e>1 e > 1 での 系 5.4 の符号違い版)から、μ ≈ m \mu\approx m μ ≈ m (天体の質量は太陽に比べて無視できる)として
v ∞ = G M ⊙ a ′ = 1.327 × 10 20 1.907 × 10 11 = 6.958 × 10 8 ≈ 2.64 × 10 4 m / s v_\infty = \sqrt{\frac{GM_\odot}{a'}} = \sqrt{\frac{1.327\times10^{20}}{1.907\times10^{11}}} = \sqrt{6.958\times10^{8}} \approx 2.64\times10^{4}\ \mathrm{m/s} v ∞ = a ′ G M ⊙ = 1.907 × 1 0 11 1.327 × 1 0 20 = 6.958 × 1 0 8 ≈ 2.64 × 1 0 4 m/s すなわち約 26 k m / s 26\ \mathrm{km/s} 26 km/s です。これは観測された値とよく一致し、この天体が太陽系外から来たことの根拠になりました。近日点での速さは v q = G M ⊙ ( 1 + e ) / q v_q = \sqrt{GM_\odot(1+e)/q} v q = G M ⊙ ( 1 + e ) / q より
v q = 1.327 × 10 20 × 2.20 3.815 × 10 10 = 7.65 × 10 9 ≈ 8.75 × 10 4 m / s v_q = \sqrt{\frac{1.327\times10^{20}\times 2.20}{3.815\times10^{10}}} = \sqrt{7.65\times10^{9}} \approx 8.75\times10^{4}\ \mathrm{m/s} v q = 3.815 × 1 0 10 1.327 × 1 0 20 × 2.20 = 7.65 × 1 0 9 ≈ 8.75 × 1 0 4 m/s で、約 87 k m / s 87\ \mathrm{km/s} 87 km/s に達します。
定理 6.1 (ケプラーの第三法則(調和の法則) )
逆二乗引力 f ( r ) = − k / r 2 f(r) = -k/r^2 f ( r ) = − k / r 2 の下での楕円軌道(e < 1 e<1 e < 1 )の公転周期 T T T と軌道長半径 a a a の間には
T 2 = 4 π 2 μ k a 3 T^2 = \frac{4\pi^2\mu}{k}\,a^3 T 2 = k 4 π 2 μ a 3 が成り立つ。とくに万有引力 k = G m 1 m 2 k = Gm_1m_2 k = G m 1 m 2 、μ = m 1 m 2 / ( m 1 + m 2 ) \mu = m_1m_2/(m_1+m_2) μ = m 1 m 2 / ( m 1 + m 2 ) の場合は
T 2 = 4 π 2 G ( m 1 + m 2 ) a 3 T^2 = \frac{4\pi^2}{G(m_1+m_2)}\,a^3 T 2 = G ( m 1 + m 2 ) 4 π 2 a 3 となる。比例係数は離心率に依らず、二天体の質量の和だけで決まる。
証明(定理 6.1) 定理 3.3 より面積速度は L / ( 2 μ ) L/(2\mu) L / ( 2 μ ) で一定です。一周する間に掃かれる面積は楕円全体の面積 π a b \pi ab π ab ですから、
T = π a b L / ( 2 μ ) = 2 π μ a b L . T = \frac{\pi ab}{L/(2\mu)} = \frac{2\pi\mu\,ab}{L}. T = L / ( 2 μ ) π ab = L 2 π μ ab . ここで L L L を軌道要素で書き換えます。定義 5.3 の ℓ = L 2 / ( μ k ) \ell = L^2/(\mu k) ℓ = L 2 / ( μ k ) と ℓ = b 2 / a \ell = b^2/a ℓ = b 2 / a を等置すると
L 2 = μ k ℓ = μ k b 2 a ⟹ L = b μ k a L^2 = \mu k \ell = \frac{\mu k b^2}{a} \;\Longrightarrow\; L = b\sqrt{\frac{\mu k}{a}} L 2 = μ k ℓ = a μ k b 2 ⟹ L = b a μ k です(L > 0 L>0 L > 0 と取りました)。これを上の T T T の式に代入すると b b b が約分されて
T = 2 π μ a b b a μ k = 2 π μ a a μ k = 2 π μ a 3 k . T = \frac{2\pi\mu\,ab}{b}\sqrt{\frac{a}{\mu k}} = 2\pi\mu a\sqrt{\frac{a}{\mu k}} = 2\pi\sqrt{\frac{\mu a^3}{k}} . T = b 2 π μ ab μ k a = 2 π μ a μ k a = 2 π k μ a 3 . 両辺を二乗すれば T 2 = 4 π 2 μ a 3 / k T^2 = 4\pi^2\mu a^3/k T 2 = 4 π 2 μ a 3 / k です。
万有引力の場合、μ / k = m 1 m 2 / ( m 1 + m 2 ) G m 1 m 2 = 1 G ( m 1 + m 2 ) \mu/k = \dfrac{m_1m_2/(m_1+m_2)}{Gm_1m_2} = \dfrac{1}{G(m_1+m_2)} μ / k = G m 1 m 2 m 1 m 2 / ( m 1 + m 2 ) = G ( m 1 + m 2 ) 1 なので、代入して第二式を得ます。
∎
例 6.3 (地球の公転周期を計算する )
太陽の重力定数積は G M ⊙ = 1.32712 × 10 20 m 3 / s 2 GM_\odot = 1.32712\times10^{20}\ \mathrm{m^3/s^2} G M ⊙ = 1.32712 × 1 0 20 m 3 / s 2 、地球の軌道長半径は a = 1.49598 × 10 11 m a = 1.49598\times10^{11}\ \mathrm{m} a = 1.49598 × 1 0 11 m です。地球質量は太陽の 3 × 10 − 6 3\times10^{-6} 3 × 1 0 − 6 倍なので 注意 6.2 の補正は無視して G ( m 1 + m 2 ) ≈ G M ⊙ G(m_1+m_2)\approx GM_\odot G ( m 1 + m 2 ) ≈ G M ⊙ とします。定理 6.1 より
T = 2 π a 3 G M ⊙ . T = 2\pi\sqrt{\frac{a^3}{GM_\odot}} . T = 2 π G M ⊙ a 3 . 順に計算します。
a 3 = ( 1.49598 × 10 11 ) 3 = 3.3479 × 10 33 m 3 , a^3 = (1.49598\times10^{11})^3 = 3.3479\times10^{33}\ \mathrm{m^3}, a 3 = ( 1.49598 × 1 0 11 ) 3 = 3.3479 × 1 0 33 m 3 , a 3 G M ⊙ = 3.3479 × 10 33 1.32712 × 10 20 = 2.5227 × 10 13 s 2 , \frac{a^3}{GM_\odot} = \frac{3.3479\times10^{33}}{1.32712\times10^{20}} = 2.5227\times10^{13}\ \mathrm{s^2}, G M ⊙ a 3 = 1.32712 × 1 0 20 3.3479 × 1 0 33 = 2.5227 × 1 0 13 s 2 , 2.5227 × 10 13 = 5.0226 × 10 6 s , T = 2 π × 5.0226 × 10 6 = 3.1559 × 10 7 s . \sqrt{2.5227\times10^{13}} = 5.0226\times10^{6}\ \mathrm{s},\qquad T = 2\pi\times5.0226\times10^{6} = 3.1559\times10^{7}\ \mathrm{s}. 2.5227 × 1 0 13 = 5.0226 × 1 0 6 s , T = 2 π × 5.0226 × 1 0 6 = 3.1559 × 1 0 7 s . これを日に直すと 3.1559 × 10 7 / 86400 = 365.3 3.1559\times10^{7}/86400 = 365.3 3.1559 × 1 0 7 /86400 = 365.3 日となり、実際の恒星年 365.256 365.256 365.256 日と有効数字 4 桁で一致します。逆二乗則という一つの仮定から、観測値がこの精度で再現されることを確認してください。
例 6.4 (等方調和振動子は「もう一つの」閉じた軌道 )
中心力 f ( r ) = − μ ω 2 r f(r) = -\mu\omega^2 r f ( r ) = − μ ω 2 r (ポテンシャル U = 1 2 μ ω 2 r 2 U = \tfrac12\mu\omega^2r^2 U = 2 1 μ ω 2 r 2 )を考えます。直交座標では運動方程式が x ¨ = − ω 2 x \ddot x = -\omega^2 x x ¨ = − ω 2 x 、y ¨ = − ω 2 y \ddot y = -\omega^2 y y ¨ = − ω 2 y と完全に分離するので、一般解は
x ( t ) = a cos ω t , y ( t ) = b sin ω t x(t) = a\cos\omega t,\qquad y(t) = b\sin\omega t x ( t ) = a cos ω t , y ( t ) = b sin ω t (初期条件で位相を調整)となり、軌道は ( x / a ) 2 + ( y / b ) 2 = 1 (x/a)^2 + (y/b)^2 = 1 ( x / a ) 2 + ( y / b ) 2 = 1 、すなわち中心が力の中心にある楕円 です。逆二乗力の場合(力の中心は焦点)と対比してください。
角運動量を確かめます。
L = μ ( x y ˙ − y x ˙ ) = μ ( a cos ω t ⋅ b ω cos ω t + b sin ω t ⋅ a ω sin ω t ) = μ a b ω L = \mu(x\dot y - y\dot x) = \mu\bigl(a\cos\omega t\cdot b\omega\cos\omega t + b\sin\omega t\cdot a\omega\sin\omega t\bigr) = \mu ab\omega L = μ ( x y ˙ − y x ˙ ) = μ ( a cos ω t ⋅ bω cos ω t + b sin ω t ⋅ aω sin ω t ) = μ abω で、確かに定数です(定理 3.1 と整合します)。周期は T = 2 π / ω T = 2\pi/\omega T = 2 π / ω で、振幅にまったく依りません。ケプラー問題の T ∝ a 3 / 2 T\propto a^{3/2} T ∝ a 3/2 とは異なる依存性です。
有界な軌道がすべて閉じる中心力は、逆二乗力 − k / r 2 -k/r^2 − k / r 2 と調和力 − μ ω 2 r -\mu\omega^2 r − μ ω 2 r の二つだけであることが知られています(ベルトランの定理)。この事実は、これら二つの力にだけ余分な保存量が存在することと表裏一体で、次節でその一方を見ます。
三次元の運動は自由度 3、つまり位相空間は 6 次元です。エネルギー E E E と角運動量 L \boldsymbol{L} L (3 成分)で保存量は 4 個ありますが、逆二乗力にはさらにもう一つ独立な保存量が存在します。
命題 7.1 (ラプラス–ルンゲ–レンツベクトルの保存 )
逆二乗引力 μ r ¨ = − k r 2 r ^ \mu\ddot{\boldsymbol{r}} = -\dfrac{k}{r^2}\hat{\boldsymbol{r}} μ r ¨ = − r 2 k r ^ の下で、p = μ r ˙ \boldsymbol{p} = \mu\dot{\boldsymbol{r}} p = μ r ˙ とおくとき
A = p × L − μ k r ^ \boldsymbol{A} = \boldsymbol{p}\times\boldsymbol{L} - \mu k\,\hat{\boldsymbol{r}} A = p × L − μ k r ^ は保存する。さらに A \boldsymbol{A} A は運動平面内にあって近点方向を向き、その大きさは ∣ A ∣ = μ k e |\boldsymbol{A}| = \mu k e ∣ A ∣ = μ k e である。
証明(命題 7.1) L \boldsymbol{L} L は 定理 3.1 により定ベクトルなので
d d t ( p × L ) = p ˙ × L = ( − k r 3 r ) × ( μ r × r ˙ ) \frac{d}{dt}(\boldsymbol{p}\times\boldsymbol{L}) = \dot{\boldsymbol{p}}\times\boldsymbol{L} = \left(-\frac{k}{r^3}\boldsymbol{r}\right)\times\left(\mu\,\boldsymbol{r}\times\dot{\boldsymbol{r}}\right) d t d ( p × L ) = p ˙ × L = ( − r 3 k r ) × ( μ r × r ˙ ) です(r ^ = r / r \hat{\boldsymbol{r}} = \boldsymbol{r}/r r ^ = r / r を使いました)。ベクトル三重積の公式 a × ( b × c ) = b ⟨ a , c ⟩ − c ⟨ a , b ⟩ \boldsymbol{a}\times(\boldsymbol{b}\times\boldsymbol{c}) = \boldsymbol{b}\langle\boldsymbol{a},\boldsymbol{c}\rangle - \boldsymbol{c}\langle\boldsymbol{a},\boldsymbol{b}\rangle a × ( b × c ) = b ⟨ a , c ⟩ − c ⟨ a , b ⟩ を適用すると
r × ( r × r ˙ ) = r ⟨ r , r ˙ ⟩ − r ˙ ∣ r ∣ 2 = r r ˙ r − r 2 r ˙ \boldsymbol{r}\times(\boldsymbol{r}\times\dot{\boldsymbol{r}}) = \boldsymbol{r}\langle\boldsymbol{r},\dot{\boldsymbol{r}}\rangle - \dot{\boldsymbol{r}}\,|\boldsymbol{r}|^2 = r\dot r\,\boldsymbol{r} - r^2\dot{\boldsymbol{r}} r × ( r × r ˙ ) = r ⟨ r , r ˙ ⟩ − r ˙ ∣ r ∣ 2 = r r ˙ r − r 2 r ˙ です。ここで ⟨ r , r ˙ ⟩ = r r ˙ \langle\boldsymbol{r},\dot{\boldsymbol{r}}\rangle = r\dot r ⟨ r , r ˙ ⟩ = r r ˙ は r 2 = ⟨ r , r ⟩ r^2 = \langle\boldsymbol{r},\boldsymbol{r}\rangle r 2 = ⟨ r , r ⟩ の両辺を微分して得られる関係です。よって
d d t ( p × L ) = − μ k r 3 ( r r ˙ r − r 2 r ˙ ) = μ k ( r ˙ r − r ˙ r 2 r ) = μ k d r ^ d t \frac{d}{dt}(\boldsymbol{p}\times\boldsymbol{L}) = -\frac{\mu k}{r^3}\left(r\dot r\,\boldsymbol{r} - r^2\dot{\boldsymbol{r}}\right) = \mu k\left(\frac{\dot{\boldsymbol{r}}}{r} - \frac{\dot r}{r^2}\boldsymbol{r}\right) = \mu k\,\frac{d\hat{\boldsymbol{r}}}{dt} d t d ( p × L ) = − r 3 μ k ( r r ˙ r − r 2 r ˙ ) = μ k ( r r ˙ − r 2 r ˙ r ) = μ k d t d r ^ となります(最後の等号は 系 3.2 の証明でも使った r ^ \hat{\boldsymbol{r}} r ^ の微分公式です)。したがって d A / d t = 0 d\boldsymbol{A}/dt = \boldsymbol{0} d A / d t = 0 です。
次に A \boldsymbol{A} A の向きと大きさを調べます。p × L \boldsymbol{p}\times\boldsymbol{L} p × L も r ^ \hat{\boldsymbol{r}} r ^ も L \boldsymbol{L} L と直交するので、A \boldsymbol{A} A は運動平面内のベクトルです。r \boldsymbol{r} r との内積を取ると、スカラー三重積の巡回性から
⟨ r , p × L ⟩ = ⟨ L , r × p ⟩ = ⟨ L , L ⟩ = L 2 \langle\boldsymbol{r}, \boldsymbol{p}\times\boldsymbol{L}\rangle = \langle\boldsymbol{L}, \boldsymbol{r}\times\boldsymbol{p}\rangle = \langle\boldsymbol{L},\boldsymbol{L}\rangle = L^2 ⟨ r , p × L ⟩ = ⟨ L , r × p ⟩ = ⟨ L , L ⟩ = L 2 なので
⟨ A , r ⟩ = L 2 − μ k r . \langle\boldsymbol{A},\boldsymbol{r}\rangle = L^2 - \mu k r . ⟨ A , r ⟩ = L 2 − μ k r . A \boldsymbol{A} A と r \boldsymbol{r} r のなす角を θ \theta θ とすれば左辺は ∣ A ∣ r cos θ |\boldsymbol{A}|\,r\cos\theta ∣ A ∣ r cos θ ですから、r r r について解いて
r = L 2 μ k + ∣ A ∣ cos θ = L 2 / ( μ k ) 1 + ( ∣ A ∣ / μ k ) cos θ r = \frac{L^2}{\mu k + |\boldsymbol{A}|\cos\theta} = \frac{L^2/(\mu k)}{1 + \bigl(|\boldsymbol{A}|/\mu k\bigr)\cos\theta} r = μ k + ∣ A ∣ cos θ L 2 = 1 + ( ∣ A ∣/ μ k ) cos θ L 2 / ( μ k ) を得ます。これを 定理 5.2 の軌道式と比べると ℓ = L 2 / ( μ k ) \ell = L^2/(\mu k) ℓ = L 2 / ( μ k ) 、e = ∣ A ∣ / ( μ k ) e = |\boldsymbol{A}|/(\mu k) e = ∣ A ∣/ ( μ k ) 、すなわち ∣ A ∣ = μ k e |\boldsymbol{A}| = \mu k e ∣ A ∣ = μ k e です。また θ = 0 \theta=0 θ = 0 、つまり A \boldsymbol{A} A の方向で r r r が最小になるので、A \boldsymbol{A} A は近点を指します。
∎
この計算は注目に値します。微分方程式を解かずに、保存量の内積を取っただけで軌道方程式が出てしまいました。物理的な意味も明快で、A \boldsymbol{A} A が定ベクトルであるということは近点の方向が動かない ということです。つまり軌道が閉じる理由が保存量として説明されます。
逆に、力が厳密な逆二乗からずれると A \boldsymbol{A} A は保存せず、近点はゆっくり回転します。水星の近日点移動(100 年あたり約 43 秒角の説明されない分)は、一般相対性理論が予言する − k / r 2 -k/r^2 − k / r 2 からのずれによるものでした。演習 8.3 では、ポテンシャルに 1 / r 2 1/r^2 1/ r 2 の補正項が加わると近点がどれだけ移動するかを実際に計算します。
演習 8.1 易
楕円軌道を描く惑星の近点での速さを v p v_p v p 、遠点での速さを v a v_a v a とする。近点距離 r min = a ( 1 − e ) r_{\min}=a(1-e) r m i n = a ( 1 − e ) 、遠点距離 r max = a ( 1 + e ) r_{\max}=a(1+e) r m a x = a ( 1 + e ) を用いて、比 v p / v a v_p/v_a v p / v a を離心率 e e e だけで表せ。また地球(e = 0.0167 e = 0.0167 e = 0.0167 )について比を数値で求めよ。
解答 近点と遠点では r r r が極値を取るので動径速度は r ˙ = 0 \dot r = 0 r ˙ = 0 です。したがって速度は e ^ θ \hat{\boldsymbol{e}}_\theta e ^ θ 方向だけを向き、速さは v = r ∣ θ ˙ ∣ v = r|\dot\theta| v = r ∣ θ ˙ ∣ になります。定理 3.1 より L = μ r 2 θ ˙ = μ r v L = \mu r^2\dot\theta = \mu r v L = μ r 2 θ ˙ = μ r v が両点で成り立ち、L L L は保存量なので
r min v p = r max v a ⟹ v p v a = r max r min = 1 + e 1 − e . r_{\min}v_p = r_{\max}v_a \;\Longrightarrow\; \frac{v_p}{v_a} = \frac{r_{\max}}{r_{\min}} = \frac{1+e}{1-e}. r m i n v p = r m a x v a ⟹ v a v p = r m i n r m a x = 1 − e 1 + e . 地球では ( 1 + 0.0167 ) / ( 1 − 0.0167 ) = 1.0167 / 0.9833 = 1.0340 (1+0.0167)/(1-0.0167) = 1.0167/0.9833 = 1.0340 ( 1 + 0.0167 ) / ( 1 − 0.0167 ) = 1.0167/0.9833 = 1.0340 です。近日点(1 月初旬)での公転速度は遠日点より約 3.4 % 3.4\% 3.4% 速く、実際の値では 30.29 k m / s 30.29\ \mathrm{km/s} 30.29 km/s と 29.29 k m / s 29.29\ \mathrm{km/s} 29.29 km/s にあたります。これがケプラーの第二法則の最も直接的な現れです。
演習 8.2 標準
地球の公転周期 T = 3.156 × 10 7 s T = 3.156\times10^{7}\ \mathrm{s} T = 3.156 × 1 0 7 s と軌道長半径 a = 1.496 × 10 11 m a = 1.496\times10^{11}\ \mathrm{m} a = 1.496 × 1 0 11 m 、および万有引力定数 G = 6.674 × 10 − 11 m 3 k g − 1 s − 2 G = 6.674\times10^{-11}\ \mathrm{m^3\,kg^{-1}\,s^{-2}} G = 6.674 × 1 0 − 11 m 3 k g − 1 s − 2 から太陽の質量を求めよ。地球質量は無視してよい。
解答 定理 6.1 の万有引力版 T 2 = 4 π 2 a 3 / ( G ( M ⊙ + m ⊕ ) ) T^2 = 4\pi^2a^3/\bigl(G(M_\odot+m_\oplus)\bigr) T 2 = 4 π 2 a 3 / ( G ( M ⊙ + m ⊕ ) ) で m ⊕ m_\oplus m ⊕ を無視すると
M ⊙ = 4 π 2 a 3 G T 2 . M_\odot = \frac{4\pi^2 a^3}{G T^2}. M ⊙ = G T 2 4 π 2 a 3 . 分子を計算します。a 3 = ( 1.496 × 10 11 ) 3 = 3.348 × 10 33 a^3 = (1.496\times10^{11})^3 = 3.348\times10^{33} a 3 = ( 1.496 × 1 0 11 ) 3 = 3.348 × 1 0 33 、4 π 2 = 39.478 4\pi^2 = 39.478 4 π 2 = 39.478 なので
4 π 2 a 3 = 39.478 × 3.348 × 10 33 = 1.3218 × 10 35 . 4\pi^2a^3 = 39.478\times3.348\times10^{33} = 1.3218\times10^{35}. 4 π 2 a 3 = 39.478 × 3.348 × 1 0 33 = 1.3218 × 1 0 35 . 分母は T 2 = ( 3.156 × 10 7 ) 2 = 9.960 × 10 14 T^2 = (3.156\times10^{7})^2 = 9.960\times10^{14} T 2 = ( 3.156 × 1 0 7 ) 2 = 9.960 × 1 0 14 より
G T 2 = 6.674 × 10 − 11 × 9.960 × 10 14 = 6.647 × 10 4 . GT^2 = 6.674\times10^{-11}\times9.960\times10^{14} = 6.647\times10^{4}. G T 2 = 6.674 × 1 0 − 11 × 9.960 × 1 0 14 = 6.647 × 1 0 4 . よって
M ⊙ = 1.3218 × 10 35 6.647 × 10 4 = 1.988 × 10 30 k g . M_\odot = \frac{1.3218\times10^{35}}{6.647\times10^{4}} = 1.988\times10^{30}\ \mathrm{kg}. M ⊙ = 6.647 × 1 0 4 1.3218 × 1 0 35 = 1.988 × 1 0 30 kg . 公表値 1.989 × 10 30 k g 1.989\times10^{30}\ \mathrm{kg} 1.989 × 1 0 30 kg と一致します。同じ手順で、惑星の衛星の周期と軌道半径からその惑星の質量が求まります。これが天体の質量を測る標準的な方法です。
演習 8.3 難
ポテンシャルが U ( r ) = − k r − β r 2 U(r) = -\dfrac{k}{r} - \dfrac{\beta}{r^2} U ( r ) = − r k − r 2 β (β \beta β は小さな定数)で与えられるとき、軌道が近点から次の近点まで進む間の θ \theta θ の増加量を求めよ。これが 2 π 2\pi 2 π からずれることを示し、β \beta β の一次までの近点移動角を計算せよ。
解答 力は f ( r ) = − U ′ ( r ) = − k r 2 − 2 β r 3 f(r) = -U'(r) = -\dfrac{k}{r^2} - \dfrac{2\beta}{r^3} f ( r ) = − U ′ ( r ) = − r 2 k − r 3 2 β 、すなわち f ( 1 / u ) = − k u 2 − 2 β u 3 f(1/u) = -ku^2 - 2\beta u^3 f ( 1/ u ) = − k u 2 − 2 β u 3 です。補題 5.1 に代入すると
d 2 u d θ 2 + u = − μ L 2 u 2 ( − k u 2 − 2 β u 3 ) = μ k L 2 + 2 μ β L 2 u . \frac{d^2u}{d\theta^2} + u = -\frac{\mu}{L^2u^2}\left(-ku^2 - 2\beta u^3\right) = \frac{\mu k}{L^2} + \frac{2\mu\beta}{L^2}u . d θ 2 d 2 u + u = − L 2 u 2 μ ( − k u 2 − 2 β u 3 ) = L 2 μ k + L 2 2 μ β u . u u u の項を左辺に移すと
d 2 u d θ 2 + ( 1 − 2 μ β L 2 ) u = μ k L 2 . \frac{d^2u}{d\theta^2} + \left(1 - \frac{2\mu\beta}{L^2}\right)u = \frac{\mu k}{L^2}. d θ 2 d 2 u + ( 1 − L 2 2 μ β ) u = L 2 μ k . γ 2 = 1 − 2 μ β / L 2 \gamma^2 = 1 - 2\mu\beta/L^2 γ 2 = 1 − 2 μ β / L 2 とおけば(β \beta β が小さいので γ 2 > 0 \gamma^2>0 γ 2 > 0 )、これは 定理 5.2 の証明とまったく同じ形の方程式で、一般解は
u ( θ ) = μ k γ 2 L 2 + C cos ( γ θ ) u(\theta) = \frac{\mu k}{\gamma^2L^2} + C\cos(\gamma\theta) u ( θ ) = γ 2 L 2 μ k + C cos ( γ θ ) です。u u u が最大(r r r が最小、すなわち近点)になるのは γ θ = 2 π n \gamma\theta = 2\pi n γ θ = 2 π n のときなので、連続する二つの近点の間の角度は
Δ θ = 2 π γ = 2 π 1 − 2 μ β / L 2 . \Delta\theta = \frac{2\pi}{\gamma} = \frac{2\pi}{\sqrt{1 - 2\mu\beta/L^2}} . Δ θ = γ 2 π = 1 − 2 μ β / L 2 2 π . β \beta β の一次まで展開すると ( 1 − x ) − 1 / 2 ≈ 1 + x / 2 (1-x)^{-1/2} \approx 1 + x/2 ( 1 − x ) − 1/2 ≈ 1 + x /2 (平均値の定理とテイラーの定理 の 定理 5.3[平均値の定理とテイラーの定理] を参照)より
Δ θ ≈ 2 π ( 1 + μ β L 2 ) = 2 π + 2 π μ β L 2 . \Delta\theta \approx 2\pi\left(1 + \frac{\mu\beta}{L^2}\right) = 2\pi + \frac{2\pi\mu\beta}{L^2}. Δ θ ≈ 2 π ( 1 + L 2 μ β ) = 2 π + L 2 2 π μ β . したがって一周ごとに近点が δ = 2 π μ β / L 2 \delta = 2\pi\mu\beta/L^2 δ = 2 π μ β / L 2 だけ進みます(β > 0 \beta>0 β > 0 なら軌道の回転と同じ向き)。β = 0 \beta=0 β = 0 なら γ = 1 \gamma=1 γ = 1 で Δ θ = 2 π \Delta\theta=2\pi Δ θ = 2 π 、つまり軌道は閉じ、命題 7.1 の A \boldsymbol{A} A の保存と整合します。一般相対論の水星近日点移動も、有効ポテンシャルに 1 / r 3 1/r^3 1/ r 3 項が加わることによる同種の効果です。
演習 8.4 難
命題 7.1 のベクトル A = p × L − μ k r ^ \boldsymbol{A} = \boldsymbol{p}\times\boldsymbol{L} - \mu k\hat{\boldsymbol{r}} A = p × L − μ k r ^ について
∣ A ∣ 2 = μ 2 k 2 + 2 μ E L 2 |\boldsymbol{A}|^2 = \mu^2k^2 + 2\mu E L^2 ∣ A ∣ 2 = μ 2 k 2 + 2 μ E L 2 を示せ。これと ∣ A ∣ = μ k e |\boldsymbol{A}| = \mu k e ∣ A ∣ = μ k e から 系 5.4 を再導出せよ。
解答 ∣ A ∣ 2 = ∣ p × L ∣ 2 − 2 μ k ⟨ p × L , r ^ ⟩ + μ 2 k 2 |\boldsymbol{A}|^2 = |\boldsymbol{p}\times\boldsymbol{L}|^2 - 2\mu k\langle\boldsymbol{p}\times\boldsymbol{L}, \hat{\boldsymbol{r}}\rangle + \mu^2k^2 ∣ A ∣ 2 = ∣ p × L ∣ 2 − 2 μ k ⟨ p × L , r ^ ⟩ + μ 2 k 2 を項ごとに計算します。
第一項。p \boldsymbol{p} p と L \boldsymbol{L} L は直交します(⟨ L , p ⟩ = ⟨ r × p , p ⟩ = 0 \langle\boldsymbol{L},\boldsymbol{p}\rangle = \langle\boldsymbol{r}\times\boldsymbol{p},\boldsymbol{p}\rangle = 0 ⟨ L , p ⟩ = ⟨ r × p , p ⟩ = 0 )から、∣ p × L ∣ = ∣ p ∣ ∣ L ∣ |\boldsymbol{p}\times\boldsymbol{L}| = |\boldsymbol{p}||\boldsymbol{L}| ∣ p × L ∣ = ∣ p ∣∣ L ∣ 、よって ∣ p × L ∣ 2 = p 2 L 2 |\boldsymbol{p}\times\boldsymbol{L}|^2 = p^2L^2 ∣ p × L ∣ 2 = p 2 L 2 です。
第二項。スカラー三重積の巡回性より
⟨ p × L , r ^ ⟩ = 1 r ⟨ p × L , r ⟩ = 1 r ⟨ L , r × p ⟩ = L 2 r . \langle\boldsymbol{p}\times\boldsymbol{L}, \hat{\boldsymbol{r}}\rangle = \frac{1}{r}\langle\boldsymbol{p}\times\boldsymbol{L},\boldsymbol{r}\rangle = \frac{1}{r}\langle\boldsymbol{L},\boldsymbol{r}\times\boldsymbol{p}\rangle = \frac{L^2}{r}. ⟨ p × L , r ^ ⟩ = r 1 ⟨ p × L , r ⟩ = r 1 ⟨ L , r × p ⟩ = r L 2 . 以上より
∣ A ∣ 2 = p 2 L 2 − 2 μ k L 2 r + μ 2 k 2 = 2 μ L 2 ( p 2 2 μ − k r ) + μ 2 k 2 . |\boldsymbol{A}|^2 = p^2L^2 - \frac{2\mu kL^2}{r} + \mu^2k^2 = 2\mu L^2\left(\frac{p^2}{2\mu} - \frac{k}{r}\right) + \mu^2k^2 . ∣ A ∣ 2 = p 2 L 2 − r 2 μ k L 2 + μ 2 k 2 = 2 μ L 2 ( 2 μ p 2 − r k ) + μ 2 k 2 . 括弧の中はまさに力学的エネルギー E = p 2 / ( 2 μ ) − k / r E = p^2/(2\mu) - k/r E = p 2 / ( 2 μ ) − k / r です(命題 4.2 )。よって ∣ A ∣ 2 = μ 2 k 2 + 2 μ E L 2 |\boldsymbol{A}|^2 = \mu^2k^2 + 2\mu EL^2 ∣ A ∣ 2 = μ 2 k 2 + 2 μ E L 2 です。
∣ A ∣ = μ k e |\boldsymbol{A}| = \mu k e ∣ A ∣ = μ k e を代入すると μ 2 k 2 e 2 = μ 2 k 2 + 2 μ E L 2 \mu^2k^2e^2 = \mu^2k^2 + 2\mu EL^2 μ 2 k 2 e 2 = μ 2 k 2 + 2 μ E L 2 、両辺を μ 2 k 2 \mu^2k^2 μ 2 k 2 で割って
e 2 = 1 + 2 E L 2 μ k 2 e^2 = 1 + \frac{2EL^2}{\mu k^2} e 2 = 1 + μ k 2 2 E L 2 となり、系 5.4 が再現されます。エネルギーと角運動量から離心率が決まるという事実が、保存量のノルムという形で自然に現れました。
H. Goldstein, C. Poole, J. Safko, Classical Mechanics , 3rd ed., Addison-Wesley, 2002 — 第 3 章「The Central Force Problem」。有効ポテンシャル、ビネの方程式、ラプラス–ルンゲ–レンツベクトル、ベルトランの定理を体系的に扱っています。
L. D. Landau, E. M. Lifshitz, Mechanics , 3rd ed., Butterworth-Heinemann, 1976 — 第 III 章「Integration of the equations of motion」。二体問題の帰着とケプラー問題を最短距離で扱う古典的な記述です。
原島鮮『力学 I』裳華房、1972 — 中心力と惑星運動の章。日本語で読める標準的な入門書です。
山本義隆『古典力学の形成 ニュートンからラグランジュへ』日本評論社、1997 — ケプラーの三法則から万有引力の逆二乗則が導かれた歴史的経緯を一次資料に即して追っています。
V. I. Arnold, Mathematical Methods of Classical Mechanics , 2nd ed., Springer, 1989 — 第 2 章。中心力場の運動を微分幾何的な視点から扱い、軌道が閉じる条件を論じています。
残された問題。 本文では軌道の「形」r ( θ ) r(\theta) r ( θ ) を完全に決めましたが、惑星が「いつ」どこにいるかは決めていません。定理 3.3 の面積速度一定則を積分すればよいのですが、d t = ( μ / L ) r ( θ ) 2 d θ dt = (\mu/L)\,r(\theta)^2 d\theta d t = ( μ / L ) r ( θ ) 2 d θ を θ \theta θ について直接積分しても閉じた形の逆関数が得られません。ここで役立つのが離心近点角 という補助変数です。
離心近点角の導入。 楕円軌道に対し、変数 ψ \psi ψ を
r = a ( 1 − e cos ψ ) r = a(1 - e\cos\psi) r = a ( 1 − e cos ψ )
で定義します。ψ = 0 \psi=0 ψ = 0 が近点(r = a ( 1 − e ) r = a(1-e) r = a ( 1 − e ) )、ψ = π \psi=\pi ψ = π が遠点(r = a ( 1 + e ) r = a(1+e) r = a ( 1 + e ) )に対応し、ψ \psi ψ は r r r の値域を過不足なく覆います。幾何学的には、楕円に外接する半径 a a a の円(補助円)に軌道上の点を垂直に射影したときの中心角にあたります。
時間の積分。 命題 4.2 のエネルギー式を r ˙ \dot r r ˙ について解き、系 5.4 の E = − k / ( 2 a ) E = -k/(2a) E = − k / ( 2 a ) と L 2 = μ k a ( 1 − e 2 ) L^2 = \mu k a(1-e^2) L 2 = μ k a ( 1 − e 2 ) (定義 5.3 の ℓ = a ( 1 − e 2 ) \ell = a(1-e^2) ℓ = a ( 1 − e 2 ) と ℓ = L 2 / μ k \ell = L^2/\mu k ℓ = L 2 / μ k から)を代入します。
r ˙ 2 = 2 μ ( E + k r ) − L 2 μ 2 r 2 = k μ a r 2 ( 2 a r − r 2 − a 2 ( 1 − e 2 ) ) = k μ a r 2 ( a 2 e 2 − ( r − a ) 2 ) . \dot r^2 = \frac{2}{\mu}\left(E + \frac{k}{r}\right) - \frac{L^2}{\mu^2r^2}
= \frac{k}{\mu a r^2}\left(2ar - r^2 - a^2(1-e^2)\right)
= \frac{k}{\mu a r^2}\left(a^2e^2 - (r-a)^2\right). r ˙ 2 = μ 2 ( E + r k ) − μ 2 r 2 L 2 = μ a r 2 k ( 2 a r − r 2 − a 2 ( 1 − e 2 ) ) = μ a r 2 k ( a 2 e 2 − ( r − a ) 2 ) .
最後の等号は 2 a r − r 2 − a 2 + a 2 e 2 = − ( r − a ) 2 + a 2 e 2 2ar - r^2 - a^2 + a^2e^2 = -(r-a)^2 + a^2e^2 2 a r − r 2 − a 2 + a 2 e 2 = − ( r − a ) 2 + a 2 e 2 という平方完成です。ここに r − a = − a e cos ψ r - a = -ae\cos\psi r − a = − a e cos ψ を入れると括弧の中は a 2 e 2 sin 2 ψ a^2e^2\sin^2\psi a 2 e 2 sin 2 ψ になります。一方 r ˙ = a e sin ψ ψ ˙ \dot r = ae\sin\psi\,\dot\psi r ˙ = a e sin ψ ψ ˙ なので、両辺を等置して
a 2 e 2 sin 2 ψ ψ ˙ 2 = k μ a r 2 a 2 e 2 sin 2 ψ ⟹ ψ ˙ = 1 r k μ a = 1 a ( 1 − e cos ψ ) k μ a a^2e^2\sin^2\psi\,\dot\psi^2 = \frac{k}{\mu a r^2}\,a^2e^2\sin^2\psi
\;\Longrightarrow\;
\dot\psi = \frac{1}{r}\sqrt{\frac{k}{\mu a}} = \frac{1}{a(1-e\cos\psi)}\sqrt{\frac{k}{\mu a}} a 2 e 2 sin 2 ψ ψ ˙ 2 = μ a r 2 k a 2 e 2 sin 2 ψ ⟹ ψ ˙ = r 1 μ a k = a ( 1 − e cos ψ ) 1 μ a k
を得ます。変数分離して近点通過時刻 t p t_p t p から積分すると
∫ 0 ψ ( 1 − e cos ψ ′ ) d ψ ′ = k μ a 3 ( t − t p ) , \int_0^{\psi}(1-e\cos\psi')\,d\psi' = \sqrt{\frac{k}{\mu a^3}}\,(t-t_p), ∫ 0 ψ ( 1 − e cos ψ ′ ) d ψ ′ = μ a 3 k ( t − t p ) ,
すなわち
ψ − e sin ψ = n ( t − t p ) , n = k μ a 3 = 2 π T \psi - e\sin\psi = n\,(t-t_p),\qquad n = \sqrt{\frac{k}{\mu a^3}} = \frac{2\pi}{T} ψ − e sin ψ = n ( t − t p ) , n = μ a 3 k = T 2 π
です(最後の等式は 定理 6.1 そのものです)。右辺 M = n ( t − t p ) M = n(t-t_p) M = n ( t − t p ) を平均近点角 と呼び、この関係式
M = ψ − e sin ψ M = \psi - e\sin\psi M = ψ − e sin ψ
をケプラー方程式 といいます。
数値解法。 ケプラー方程式は ψ \psi ψ について初等関数で解けません。実用上はニュートン法で解きます。g ( ψ ) = ψ − e sin ψ − M g(\psi) = \psi - e\sin\psi - M g ( ψ ) = ψ − e sin ψ − M とおくと g ′ ( ψ ) = 1 − e cos ψ ≥ 1 − e > 0 g'(\psi) = 1 - e\cos\psi \ge 1-e > 0 g ′ ( ψ ) = 1 − e cos ψ ≥ 1 − e > 0 なので、0 ≤ e < 1 0\le e<1 0 ≤ e < 1 では g g g は狭義単調増加で解が一意に存在し、ニュートン法は安定に収束します。
def solve_kepler ( M , e , tol= 1e-12 , max_iter= 60 ) :
""" ケプラー方程式 M = psi - e*sin(psi) をニュートン法で解く。 """
psi = M if e < 0.8 else np.pi
for _ in range ( max_iter ):
g = psi - e * np. sin ( psi ) - M
dpsi = - g / ( 1.0 - e * np. cos ( psi ))
def position ( t , a , e , n , t_p= 0.0 ) :
""" 時刻 t における近点からの角度 theta と動径 r を返す。 """
psi = solve_kepler ( n * (t - t_p) , e )
r = a * ( 1.0 - e * np. cos ( psi ))
theta = 2.0 * np. arctan2 ( np. sqrt ( 1 + e ) * np. sin ( psi / 2 ) ,
np. sqrt ( 1 - e ) * np. cos ( psi / 2 ))
最後の theta の式は、tan ( θ / 2 ) = ( 1 + e ) / ( 1 − e ) tan ( ψ / 2 ) \tan(\theta/2) = \sqrt{(1+e)/(1-e)}\,\tan(\psi/2) tan ( θ /2 ) = ( 1 + e ) / ( 1 − e ) tan ( ψ /2 ) という真近点角と離心近点角の関係を、象限を正しく扱えるように arctan2 で書いたものです。この関係は r = a ( 1 − e cos ψ ) r=a(1-e\cos\psi) r = a ( 1 − e cos ψ ) と r = a ( 1 − e 2 ) / ( 1 + e cos θ ) r = a(1-e^2)/(1+e\cos\theta) r = a ( 1 − e 2 ) / ( 1 + e cos θ ) を等置して cos θ \cos\theta cos θ を cos ψ \cos\psi cos ψ で表し、半角公式を使えば導けます。
天体暦の計算はこの手順の繰り返しです。観測から軌道要素 ( a , e , … ) (a, e, \ldots) ( a , e , … ) を決め、ケプラー方程式で任意時刻の位置を再現する。ケプラーが観測表から法則を読み取ったのとちょうど逆向きの操作を、私たちは毎日実行していることになります。