水星の近日点移動(8)

前回
$\Large \left(\frac{{\rm d}x}{{\rm d}\psi}\right)^2 = 1 - x^2 + \beta x^3 \hspace{560pt}(1)$
の近似解を求めて
$\Large x= \cos(\Omega\psi)+\frac{1}{4}\beta\{3-\cos(2\Omega\psi)\} +\frac{1}{64}\beta^2\{37\cos(\Omega\psi)+3\cos(3\Omega\psi)\}+O(\beta^3)\hspace{52pt}(2)\\ \Omega = 1-\frac{15}{16}\beta^2\hspace{600pt}\hspace{43pt}(3)$
となるところまで計算した。

ところで、前回の (3) のところで変数 ξ を導入したとき ξ≧0 とした。x で考えると、これは x が最大値から最初に 0 になるところまで、およそ4分の1周期しか計算できていないということになる。
(1) の振り子の運動は周期的で時間反転に対する対称性をもっているはずだが、それに対応して (2) は $\frac{2\pi}{\Omega}$ の周期を持ち、ψ について偶関数なので、まず x≧0 のところは(今の近似の精度で)正しいことが分かる。
x の方向を反転し(つまり x → -x とし)、同時に β の符号を反転させると (1) はもとと同じ形になる。ψ を半周期ずらせば(つまり $\psi \rightarrow \psi+\frac{\pi}{\Omega}$ とすれば)、x の最大値だったところが最小値になるので、結局、 $x \rightarrow -x,\;\beta \rightarrow -\beta,\;\psi \rightarrow \psi+\frac{\pi}{\Omega}$ の変換を同時に行って (2) が形を変えなければ、(2) の解は全域で正しいと言える。(2) を見れば実際にそうなっていることが分かる。

(2) から水星の近日点移動(5) の (3) 式
$\Large \left(\frac{{\rm d}y}{{\rm d}\phi}\right)^2 = e^2 - (y-1)^2 + \alpha y^3\hspace{514pt}(4)$
の解を求める。

いろいろ思いつくままに書いてきたので、無計画に建て増しを重ねた家みたいなことになってるが、水星の近日点移動(5) の (5) 式、水星の近日点移動(6) の (2)(3)(5) 式から
$\Large y = y_0+\frac{\gamma}{\omega}x \hspace{600pt}\hspace{50pt}(5)\\ \psi = \omega\phi \hspace{600pt}\hspace{85pt}(6)\\ \beta=\frac{\gamma}{\omega^3}\alpha \hspace{600pt}\hspace{80pt}(7)$
で、ここに出てくる $y_0\;\omega,\;\gamma$ は水星の近日点移動(6) の (13)(15)、水星の近日点移動(5) の (7) を見ると
$\Large y_0=1 + \frac{3}{2}\alpha+\frac{9}{2}\alpha^2 + O(\alpha^3)\hspace{520pt}(8) \\ \frac{1}{\omega} = 1+\frac{3}{2}\alpha+\frac{45}{8}\alpha^2+O(\alpha^3)\hspace{521pt}(9) \\ \gamma^2 = e^2-(y_0-1)^2 + \alpha{y_0}^3\hspace{527pt}(10)$
で、γ＞0 だった。

(10) に (8) を代入して
$\Large \gamma^2=e^2+\alpha+\frac{9}{4}\alpha^2+O(\alpha^3)$
これと $\sqrt{1+\epsilon}=1+\frac{1}{2}\epsilon-\frac{1}{8}\epsilon^2+O(\epsilon^3)$ から
$\Large \gamma = e\sqrt{1+\frac{1}{e^2}\alpha + \frac{9}{4e^2}\alpha^2} +O(\alpha^3)\\ =e\left\{1 + \frac{1}{2}\left(\frac{1}{e^2}\alpha + \frac{9}{4e^2}\alpha^2\right)-\frac{1}{8}\left(\frac{1}{e^2}\alpha\right)^2\right\}+O(\alpha^3)\\=e + \frac{1}{2e}\alpha+ \frac{9e^2-1}{8e^3}\alpha^2+O(\alpha^3)$
これと (9) から
$\Large \frac{\gamma}{\omega}=e+\frac{1+3e^2}{2e}\alpha+\frac{-1+15e^2+45e^4}{8e^3}\alpha^2+O(\alpha^3)\hspace{377pt}(11)$
これと (5)(8) から
$\Large y=1+\frac{3}{2}\alpha+{9}{2}\alpha^2 +\left(e+\frac{1+3e^2}{2e}\alpha+\frac{-1+15e^2+45e^4}{8e^3}\alpha^2\right)x+O(\alpha^3)\hspace{201pt}(12)$
(7)(9)(11) から
$\Large \beta=e\alpha + \frac{1+9e^2}{2e}\alpha^2+O(\alpha^3)\hspace{512pt}(13)$
(9) と $\frac{1}{1+\epsilon}=1-\epsilon+\epsilon^2+O(\epsilon^3)$ から
$\Large \omega=1 - \frac{3}{2}\alpha- \frac{27}{8}\alpha^2+O(\alpha^3)$
これと (3)(13) から
$\Large \Omega\psi = \omega'\phi\hspace{600pt}\hspace{57pt}(14)$
ただし
$\Large \omega' = 1 - \frac{3}{2}\alpha - \frac{3(18+5e^2)}{16}\alpha^2+O(\alpha^3)\hspace{450pt}(15)$
とした。
あとは (2)(12)(13)(14) から (4) の解は α の 2次までで
$\Large y = 1+e\cos(\omega'\phi) + \alpha\left\{\frac{3(2+e^2)}{4} + \frac{1+3e^2}{2e}\cos(\omega'\phi)-\frac{e^2}{4}\cos(2\omega'\phi)\right\}\\\;\;\; + \alpha^2\left\{\frac{3(7+6e^2)}{4} + \frac{-8+120e^2+360e^4+37e^6}{64e^3}\cos(\omega'\phi) - \frac{1+6e^2}{4}\cos(2\omega'\phi)+\frac{3e^3}{64}\cos(3\omega'\phi)\right\}\\\;\;\;+O(\alpha^3)\hspace{600pt}\hspace{50pt}(16)$

(4) の厳密解も求めておこう。
水星の近日点移動(6) の後半と同じようなことなので多少省略して書く。
(4) の右辺の零点を ξ, η, ζ (ξ＜η＜ζ) とすると、y は ξ≦y≦η の範囲で動く。(16) に合わせて φ=0 で y が最大値(=η)を取るように初期条件を選べば
$\Large \phi = \int_y^\eta \frac{{\rm d}y}{\sqrt{e^2-(y-1)^2+\alpha y^3}} = \frac{1}{\sqrt{\alpha}}\int_y^\eta \frac{{\rm d}y}{\sqrt{(y-\xi)(y-\eta)(y-\zeta)}}$
$\int_y^\eta f{\rm d}y = \int_\xi^\eta f{\rm d}y -\int_\xi^y f{\rm d}y$
だから ∫dx/√{(x-ξ)(x-η)(x-ζ)} の (22)(23) から
$\Large \phi = \sqrt{\frac{2k}{\alpha\Delta}}\,\left\{K(k)-F\left(\frac{(y-\mu)/\Delta-k}{1-k(y-\mu)/\Delta},\,k\right)\right\}$
これはヤコビの楕円関数 sn を使えば y について解ける。sn は第一種完全楕円積分 F の逆関数で、 $u = F(z,k)$ のとき $\Large z = {\rm sn}(u,k)$ 。途中省略して結果だけ書けば、s を
$\Large s = {\rm sn}\left(K(k)-\sqrt{\frac{\alpha\Delta}{2k}}\phi,\, k\right)$
として
$\Large y = \mu +\Delta\frac{k+s}{1+ks}\hspace{600pt}\hspace{25pt}(17)$

下は高精度計算サイトのフリー計算で、厳密解 (17) と α の2次までの近似解 (16) を出力するプログラム。
$\alpha$ の値を適当に動かしてやると誤差がちゃんと $\alpha^3$ 程度になっていることが分かる。

phi = 1; 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));
s = jacobiSN(ellipticK(k) - sqrt(alpha*Delta/2/k)*phi, k);
y = mu + Delta*(k+s)/(1+k*s);

omega1 = 1 - 3/2*alpha - 3*(18+5*e0^2)/16*alpha^2;
c = cos(omega1*phi); c2 = cos(2*omega1*phi); c3 = cos(3*omega1*phi);
y2 = 1 + e0*c + alpha*(3*(2+e0^2)/4 + (1+3*e0^2)/2/e0*c - e0^2/4*c2)
     + alpha^2*(3*(7+6*e0^2)/4 + (-8+120*e0^2+360*e0^4+37*e0^6)/64/e0^3*c
     - (1+6*e0^2)/4*c2 + 3*e0^3/64*c3);
println(y, y2);

他との比較のため、最後に水星の近日点移動(5) の (1) 式
$\Large \left(\frac{{\rm d}u}{{\rm d}\phi}\right)^2 = \frac{e^2}{l^2} - \left(u-\frac{1}{l}\right)^2+r_{\rm g} u^3$
の摂動論による近似解を書いておく。 $u=\frac{y}{l},\;\alpha=\frac{r_{\rm g}}{l}$ なので (16) より
$\Large u = \frac{1+e\cos(\omega'\phi)}{l} + \frac{r_{\rm g}}{l^2}\left\{\frac{3(2+e^2)}{4} + \frac{1+3e^2}{2e}\cos(\omega'\phi)-\frac{e^2}{4}\cos(2\omega'\phi)\right\}\\\;\;\; + \frac{r_{\rm g}^2}{l^3}\left\{\frac{3(7+6e^2)}{4} + \frac{-8+120e^2+360e^4+37e^6}{64e^3}\cos(\omega'\phi) - \frac{1+6e^2}{4}\cos(2\omega'\phi)+\frac{3e^3}{64}\cos(3\omega'\phi)\right\}+O\left(\frac{{r_{\rm g}}^3}{l^3}\right)\vspace{40pt}$
ただし ω' は (15) で与えられる。
軌道の動径 r の式にしたければ $r=\frac{1}{u}$ だったことを思い出せばよい。