水星の近日点移動(6)

前回までで一応近日点移動の値は求まったのだが、前回の (13) 式
$\Large \left(\frac{{\rm d}y'}{{\rm d}\phi}\right)^2 = \gamma^2 - \omega^2 y'^2 + \alpha y'^3\hspace{287pt}(1)$
で、 $\alpha y'^3$ の項を落とすという多少乱暴というか謎な方法を採用していた。
今回はこの項を落とさずに摂動論で近日点移動を求めてみる。

式を見やすくするために変数 x, ψ を導入して
$\Large \phi=\frac{1}{\omega}\psi\hspace{464pt}(2) \\ y' = \frac{\gamma}{\omega}x \hspace{463pt}(3)$
と変換すると
$\Large \left(\frac{{\rm d}x}{{\rm d}\psi}\right)^2 = 1 - x^2 + \beta x^3 \hspace{340pt}(4)$
ここで
$\Large \beta = \frac{\gamma}{\omega^3}\alpha \hspace{455pt}(5)$
とした。
(4) は昔書いた非線形振り子の記事の n=3 の場合で、ψ を時間と見れば x=0 から x が最大値になるまでにかかる時間は、既に計算してあるように
$\Large \frac{\pi}{2} + \beta + \frac{15\pi}{32}\beta^2+O(\beta^3) \hspace{330pt}(6)$
ただし、あっちの記事では省略していた $O(\beta^3)$ の項を復活させた。
x=0 から x が最小値になるまでにかかる時間は、(4) の形を見れば β → -β と置き換えればよいことが分かる。
$\Large \frac{\pi}{2} - \beta + \frac{15\pi}{32}\beta^2+O(\beta^3) \hspace{330pt}(7)$
(6)(7) を足せば半周期ぶんの時間だから、1周期はその2倍で、結局 (4) の振動の周期 Ψ は
$\Large \Psi = 2\pi + \frac{15\pi}{8}\beta^2+O(\beta^3) = 2\pi\left(1+\frac{15}{16}\beta^2\right)+O(\beta^3)\hspace{84pt}(8)$
となって β の1次の項が相殺する。

(8) を使って (1) の周期を α の2次の精度で計算してみよう。
前回の (8)(9)(10)(12) を少し変形して
$\Large \omega^2 = 1-3\alpha y_0\hspace{405pt}(9) \\ y_0 = 1 + \frac{3}{2}\alpha{y_0}^2\hspace{399pt}(10) \\ y_0 = 1 + O(\alpha)\hspace{410pt}(11) \\ \gamma = e+O(\alpha)\hspace{415pt}(12)$
(10) の右辺に (11) を代入。
$\Large y_0=1 + \frac{3}{2}\alpha + O(\alpha^2)$
これをもう一度 (10) の右辺に代入する*1。
$\Large y_0=1 + \frac{3}{2}\alpha\left(1+\frac{3}{2}\alpha\right)^2 + O(\alpha^3) = 1 + \frac{3}{2}\alpha+\frac{9}{2}\alpha^2 + O(\alpha^3)\hspace{50pt}(13)$
これを (9) に代入。
$\Large \omega^2 = 1-3\alpha - \frac{9}{2}\alpha^2+O(\alpha^3)\hspace{292pt}(14)$
$\Large \frac{1}{\omega} = (\omega^2)^{-\frac{1}{2}} = 1-\frac{1}{2}\left(-3\alpha-\frac{9}{2}\alpha^2\right)+\frac{3}{8}(-3\alpha)^2 + O(\alpha^3) \\= 1+\frac{3}{2}\alpha+\frac{45}{8}\alpha^2+O(\alpha^3)\hspace{310pt}(15)$
上では $|\epsilon|$ が小さいときに成り立つ式 $(1+\epsilon)^{-\frac{1}{2}}=1-\frac{1}{2}\epsilon+\frac{3}{8}\epsilon^2+O(\epsilon^3)$ を使った。
(5)(12)(15) より
$\Large \beta=e\alpha+O(\alpha)^2$
これを (8) に代入。
$\Large \Psi = 2\pi\left(1 + \frac{15e^2}{16}\alpha^2\right)+O(\alpha^3)\hspace{271pt}(16)$

Ψ が (4) の周期なので (1) の周期*2 Φ は、(3) より
$\Large \Phi = \frac{1}{\omega}\Psi$
これに (15)(16) を代入して
$\Large \Phi = 2\pi\left(1 + \frac{3}{2}\alpha + \frac{15(6+e^2)}{16}\alpha^2\right)+O(\alpha^3)\hspace{192pt}(17)$
1公転当たりの近日点移動 δφ は、前回やったようにこれから 2π を引けばいい。
$\Large \delta\phi = \Phi-2\pi = 3\pi\alpha + \frac{15\pi(6+e^2)}{8}\alpha^2+O(\alpha^3)\hspace{150pt}(18)$

前回求めた値は
$\Large \delta\phi = 3\pi\alpha+O(\alpha^2)\hspace{360pt}(19)$
で、(18) と比べてみればこれが α の1次までの精度で正しいことが分かる。
(1) から $\alpha y'^3$ の項を落として α の 1次まで正しい結果が得られた理由はほぼ明らかで、(4) の $\beta x^3$ の項が、(8) を見ればわかるとおり、β の1次(つまり α の1次) では周期に影響しないからだ。

ついでに δφ の厳密な値を求めておこう。
前回の (3)
$\Large \left(\frac{{\rm d}y}{{\rm d}\phi}\right)^2 = e^2 - (y-1)^2 + \alpha y^3\hspace{288pt}(20)$
(これを $y=y_0+y'$ で変数変換したのが今回の (1)) は変数分離形の微分方程式で、解は
$\Large \phi = \int \frac{{\rm d}y}{\sqrt{e^2-(y-1)^2+\alpha y^3}}$
(20) の右辺の3個の零点を ξ, η, ζ (ξ＜η＜ζ) とすると、ξ から η まで積分したものが半周期だから
$\Large \frac{1}{2}\Phi = \int_\xi^\eta \frac{{\rm d}y}{\sqrt{e^2-(y-1)^2+\alpha y^3}}$
よって
$\Large \delta\phi = \Phi-2\pi = 2\int_\xi^\eta \frac{{\rm d}y}{\sqrt{e^2-(y-1)^2+\alpha y^3}} - 2\pi = \frac{2}{\sqrt{\alpha}}\int_\xi^\eta \frac{{\rm d}y}{\sqrt{(y-\xi)(y-\eta)(y-\zeta)}}-2\pi$
この積分は先日の ∫dx/√{(x-ξ)(x-η)(x-ζ)} の記事で書いたように楕円積分を使って書けて、結果は
$\Large \delta\phi = 4\sqrt{\frac{2k}{\alpha\Delta}}K(k) - 2\pi\hspace{330pt}(21)$
となる。ただし、K(k) は第一種完全楕円積分、Δ と k は ξ, η, ζ から求められる値。

e と α を水星についての値 $e=0.2056,\;\alpha = 5.325\times 10^{-8}$ として(近似値ではなくて正確にこの値とする)、(21)(18)(19) の値を計算すると
$\Large 4\sqrt{\frac{2k}{\alpha\Delta}}K(k) - 2\pi\hspace{13pt} = 5.01869527334098\times 10^{-7} \\ 3\pi\alpha + \frac{15\pi(6+e^2)}{8}\alpha^2 = 5.01869527334072\times 10^{-7} \\ 3\pi\alpha\hspace{123pt} = 5.0186943\times 10^{-7}$
となる。

数値計算には CASIO の高精度計算サイトのフリー計算を使わせていただいた。
「計算式」のところに下のプログラムを貼り付けて計算させると上の値が出る。
これを見れば (21) で端折った計算がだいたい分かるだろう。
(20) の右辺を $y=Y+\frac{1}{3\alpha}$ と変数変換して $Y^3-3pY-2q$ の形にして、3次方程式 $Y^3-3pY-2q=0$ が3実根を持つ場合の解の公式を前半で使っている。

e0 = 0.2056; alpha = 5.325E-8;

p = (1-6*alpha)/9/alpha^2;
q = (1-9*alpha)/27/alpha^3 + (1-e0^2)/2/alpha;
r = sqrt(p); theta = acos(q/r^3)/3;
xi = 1/3/alpha + 2*r*cos(theta+2*pi/3);
eta = 1/3/alpha + 2*r*cos(theta-2*pi/3);
Zeta = 1/3/alpha + 2*r*cos(theta);

Delta = (eta-xi)/2; mu = (xi+eta)/2;
zeta1 = (Zeta-mu)/Delta; k = 1/(zeta1 + sqrt(zeta1^2-1));
deltaphi = 4*sqrt(2*k/alpha/Delta)*ellipticK(k) - 2*pi;
deltaphi2 = 3*pi*alpha + 15*pi*(6+e0^2)/8*alpha^2;
deltaphi1 = 3*pi*alpha;
println(deltaphi, deltaphi2, deltaphi1);

*1:今の計算には $y_0$ は α の1次まで求めれば十分だが、後で使うと思うので2次まで計算している。

*2: (1) 式も振動子と見て「周期」という言葉を使う。惑星の公転周期とはもちろん別物。