(62) 水星の近日点移動 - 読み書きプログラミング

一般相対論的な重力場での運動の場合、軌道は楕円にならず、近日点は移動します。これを導くため、Schwarzschild 解から出発します。Schwarzschild 解は以下のような重力場の方程式の球対称解です。

$ds^2=(c^2-\frac{2GM}{r})dt^2-\frac{1}{c^2-\frac{2GM}{r}}dr^2-r^2d\theta^2$

この計量で測地線が満たさなければいけない方程式は以下になります。

$\frac{d^2u}{d\theta^2}+u=a+\frac{3GM}{c^2}u^2$

ここで、

$a=\frac{GM}{h_1^2}, h_1&=r^2\frac{d\theta}{ds}$

歳差運動を抽出するため、以下の変数変換を行います。

$x=(1+bg)\theta, b=\frac{3GM}{ac^2}$

変数変数後の微分方程式は、以下のようになります。

$(1+bg)^2\frac{d^2u}{dx^2}+u=a(1+bu^2)$

uは0次ではNewton項を示すはずですから、以下のように展開します。

$u(x)=a-A\sin x+bu_1(x)+...$

1次の摂動項 $u_1(x)$ について微分方程式を解きましょう。

/* 9.4.20m */
u : a - A*sin(x) + b*u1(x)$
rel : taylor((1 + b*g)^2*diff(u, x, 2) + u, b, 0, 1) = taylor(a*(1 + b*u^2), b, 0, 1);

$a+\left( 2\,g\,\mathrm{sin}\left( x\right) \,A+\frac{{d}^{2}}{d\,{x}^{2}}\,\mathrm{u1}\left( x\right) +\mathrm{u1}\left( x\right) \right) \,b+...=a+\left( a\,{\mathrm{sin}\left( x\right) }^{2}\,{A}^{2}-2\,{a}^{2}\,\mathrm{sin}\left( x\right) \,A+{a}^{3}\right) \,b+...$

0次(b=0)が合っていることは目で確認できるので、1次の式を導きましょう。

/* 9.4.22m */
rel2 : subst(0, b, diff(rel, b));

$2\,g\,\mathrm{sin}\left( x\right) \,A+\frac{{d}^{2}}{d\,{x}^{2}}\,\mathrm{u1}\left( x\right) +\mathrm{u1}\left( x\right) =a\,{\mathrm{sin}\left( x\right) }^{2}\,{A}^{2}-2\,{a}^{2}\,\mathrm{sin}\left( x\right) \,A+{a}^{3}$

これは角速度1の調和振動子にsin(x)の関数の強制力が加わったような運動方程式の形をしています。sin(x)の一次の項は「共振力」なので摂動項の安定性としてはのぞましくありません。従って安定な摂動解のためには、パラメータに
$g=-{a}^{2}$ という条件が付くことになります。この時、微分方程式は以下になります。

$a\,{\mathrm{sin}\left( x\right) }^{2}\,{A}^{2}-\frac{{d}^{2}}{d\,{x}^{2}}\,\mathrm{u1}\left( x\right) -\mathrm{u1}\left( x\right) +{a}^{3}=0$

これをMaximaに解かせましょう。

desolve(a*sin(x)^2*A^2-'diff(u1(x),x,2)-u1(x)+a^3, u1(x));

$\mathrm{u1}\left( x\right) =-\frac{\mathrm{cos}\left( x\right) \,\left( 2\,a\,{A}^{2}+3\,{a}^{3}-3\,\mathrm{u1}\left( 0\right) \right) }{3}+\frac{a\,{A}^{2}+2\,{a}^{3}}{2}+\frac{a\,\mathrm{cos}\left( 2\,x\right) \,{A}^{2}}{6}+\mathrm{sin}\left( x\right) \,\left( \left. \frac{d}{d\,x}\,\mathrm{u1}\left( x\right) \right|_{x=0}\right)$

適当な初期条件の下では以下になることが見て取れます。

$\mathrm{u1}\left( x\right) =\frac{a\,{A}^{2}+2\,{a}^{3}}{2}+\frac{a\,\mathrm{cos}\left( 2\,x\right) \,{A}^{2}}{6}$

以上で、摂動解全体は

$r=\frac{1}{a-A\sin((1-ba^2)\theta)+bu_1}$

であることがわかりました。u1が安定な解を持つための条件が、sinの中のθの係数に反映されており、これが近日点移動をもたらします。

では、近日点移動を示す軌道をMaximaに秒がさせてみましょう。デフォルメのため、太陽の質量は10の6乗倍になっています。

/* 9.4.27m */
G : 6.67e-11$
M : 1.99e30*10^6$
c : 3e8$
a : 1.8e-11$
e : 0.205$
b : 3*G*M/(c^2*a)$
u1(x) := a^3*(1 + e^2/2 + e^2/6*cos(2*x))$
rad(the) := block([x],
  x : (1 - a^2*b)*the,
  return(1/(1 - e*sin(x) + b/a*u1(x)))
  )$
plot2d([parametric, rad(th)*cos(th), rad(th)*sin(th), [th, 0, 10*%pi], [nticks, 500]], [x, -2, 2], [y, -1.5, 1.5])$

参考文献

Crandall, Mathematica―理工系ツールとしての (アジソンウェスレイ・トッパン情報科学シリーズ) p.298-p.302

あとがき

Crandallシリーズは一応これで終わりになります。
東北大震災の時に、何か継続して努力する日常を求めて始めたことでしたが、意外と早くフィナーレを迎えました。

Crandallさんのこの書籍は、正直そんなに読んでいませんでした。Mathematicaも身近にありませんでしたから。
高価なMathematicaを使った内容が、今ではそのほとんどすべてが、オープンソースのソフトウエアの一部の機能でできてしまうことがわかりました。

これだけ教育が行き届いても、数学の道具としての一面を活用する人たちは極限られた人たちだけというのが現状です。書籍に数式を載せるのは専門書だけというコンセンサスができてしまっています。
一方で、表計算を使う人はかなり大勢と想像します。数字を扱うことが重要になっている。適したツールがあれば利用者は増える。

iPadがコンピュータのイメージを変えてくれました。変えたというよりは、KayさんのDynabookのイメージに少し戻してくれたというべきでしょうか。
こういうハードウエアやプラットフォームの変化の中で、数式を日常的に活用できるようなブレイクスルーが起こりそうな予感があります。それに貢献できる機会がもしあれば幸せです。

62回のおつきあい、ありがとうございました。Crandallシリーズは、記事の加筆作業に移行します。