ケプラー問題の数値計算

0．これはなに

大学のプログラミング演習で提出したレポートです。去年のアドベントカレンダーでネタが本当になくなったときに使おうかと思っていましたが、Markdown + KaTeX をはてなマークダウンに変換するのが面倒すぎて没にしていたものです。ふと思い出して自動で変換してくれるスクリプトを探したところ存在したので、勢いのまま公開してしまいます。

次節からはレポートそのままです。

1．はじめに

今学期は航空宇宙工学科の「宇宙軌道力学」を履修しているので、その復習も兼ねて、2次元ケプラー問題の円軌道・楕円軌道を様々な方法で数値計算してみることにする。

一般の宇宙機の軌道設計においては、惑星の運動を所与とする近似を行ってもなお惑星間の相対運動によってポテンシャルに時間依存性が生じるため、あまりエレガントな計算法がない。対して、ケプラー問題にはハミルトニアン、角運動量、ルンゲ=レンツ・ベクトルの3つの保存量があり、しかも3つすべてを保存する構造保存解法の存在が知られている。更に、初等的な式変形だけで時刻から位置と速度を求める方程式（ケプラー方程式）を導くことができ、これを数値的に求根することによっても解ける。

このレポートでは、まずハミルトニアンを保存するシンプレクティック・オイラー法による計算を行う。シンプレクティック・オイラー法はルンゲ=レンツ・ベクトルを保存しないため、軌道が閉じないことがあることを見る。その後、より高精度に計算できる方法としてケプラー方程式による方法と、Minesaki-Nakamura [1] によるレヴィ=チヴィタ変換を使った構造保存解法を導入し、それらの挙動を見る。

2．ケプラー問題

定義

以下のハミルトニアン $H(q _ 1, q _ 2, p _ 1, p _ 2)$ に従って運動する質点の軌道を問うのがケプラー問題である。簡単のため、質点の質量を1、万有引力定数と中心天体質量の積を $\mu$ とした。

H = \frac{p_1^2+p_2^2}{2} - \frac{\mu}{\sqrt{q_1^2+q_2^2}}

正準方程式は、

\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix}q_k\\p_k\end{bmatrix}=\begin{bmatrix}p_k\\-\mu q_k\left(q_1^2+q_2^2\right)^{-3/2}\end{bmatrix}

となる。

また、角運動量ベクトル ${h}$ とルンゲ=レンツ・ベクトル $P$ は以下の表式で与えられ、それぞれ保存する。

h = \left(0,\, 0,\, q_1 p_2 - q_2 p_1\right) \\P = \left(p_2\left(q_1 p_2 - q_2 p_1\right)-\frac{\mu q_1}{\sqrt{q_1^2+q_2^2}},\,\, p_1\left(q_1 p_2 - q_2 p_1\right)-\frac{\mu q_2}{\sqrt{q_1^2+q_2^2}},\,\, 0\right)

解

軌道形状が円錐曲線になることはよく知られているが、ケプラー方程式のための準備も兼ねて、ここで証明過程を概観しておく。ハミルトン形式ではなくニュートン方程式で書いたほうが見通しがいい。

$r=\left[q _ 1, q _ 2\right]$ として、ニュートン方程式は以下。

\ddot{r} + \frac{\mu}{|r|^3} r=0

3つの保存量は次のようになる。それぞれニュートン方程式の両辺と $r$ の外積、$\dot{r}$ の内積、$h$ の外積を取れば導出できる。なおエネルギーはハミルトニアン関数ではないという意味を込めて $\varepsilon$ で書いた。

\begin{aligned}\varepsilon &= \frac{\dot{r}^2}{2} - \frac{\mu}{|r|}\\h &= r \times \dot{r} \\P &= \dot{r} \times h - \mu\frac{r}{|r|}\end{aligned}

ここで、$P$ と $r$ の内積を取り、両者の間の角度を $\nu$ と置く。

P \cdot r = |P||r|\cos{\nu}=r \cdot \left(\dot{r}\times h\right)-\mu\frac{r^2}{|r|}

最右辺の第1項はスカラー三重積で、適切に順番を変えれば $h ^ 2$ に等しいと分かる。よって、

|P||r|\cos{\nu}=h^2 -\mu |r|

これを $|r|$ について解き、$\rho = \frac{h ^ 2}{\mu},\, e=\frac{|P|}{\mu}$ とすれば、よく知られた円錐曲線の極座標表示を得る。

|r| = \frac{\rho}{1+e\cos{\nu}}

この表示からはルンゲ=レンツ・ベクトルについての直感的解釈が得られる。

$r$ と $P$ が同じ向きを向くとき $\cos{\nu}$ は最大になるので $|r|$ は最小となるが、このことから $P$ の方向は近心点を指していることがわかる。そのため、$\nu$ は真近点角と呼ばれる（後に真でない近点角も登場する）。さらに、$e=\frac{|P|}{\mu}$ であるから、 $P$ の大きさは離心率に比例している。

軌道形状が $h$ と $P$ のみで定まるということは、$h$ と $P$ から $\varepsilon$ を計算できるはずである^{[† 1]}。実際、$P ^ 2$ を計算することで、

\begin{aligned}P^2 &= \dot{r}^2h^2 - 2h^2\frac{\mu}{|r|}+\mu^2 \\P^2 - \mu^2 &= h^2\left(\dot{r}^2-2\frac{\mu}{|r|}\right) \\&= 2\varepsilon h^2 \\\therefore &\,\,\,\, e^2-1 = 2\frac{h^2}{\mu^2}\varepsilon\end{aligned}

という関係式を導ける。

3．シンプレクティック・オイラー法

$\tau$ を時間刻み幅として、前節の正準方程式を次のように離散化すると、各ステップ $\left(q _ k ^ {(i)}, p _ k ^ {(i)}\right)$ $→$ $\left(q _ k ^ {(i+1)}, p _ k ^ {(i+1)}\right)$ は正準変換となる。

\begin{bmatrix}q_k^{\left(i+1\right)} \\p_k^{\left(i+1\right)}\end{bmatrix}=\begin{bmatrix}q_k^{\left(i\right)} + \tau p_k^{\left(i\right)} \\p_k^{\left(i\right)} - \tau\mu q_k^{\left(i+1\right)}\left(\left(q_1^{\left(i+1\right)}\right)^2+\left(q_2^{\left(i+1\right)}\right)^2\right)^{-3/2}\end{bmatrix}

ハミルトニアンが $H(q, p) = K(p) + U(q)$ と書けるときこのハミルトニアンは可分であると言い、そのとき、このように $q _ k ^ {(i+1)}$ を決めたあとそれを使って $p _ k ^ {(i+1)}$ を決めるような離散化をすると一般に必ず正準変換となる。

正準変換になっているということは、この変換が表す運動があり、そこで保存するハミルトニアン $\tilde{H}$ がある。$\tilde{H}$ は $H$ とは一致しないが、両者のずれは $O(\tau)$ であることが知られている。

円軌道

$\mu=1$、初期条件 $\left(q _ 1 ^ {(0)}, q _ 2 ^ {(0)}, p _ 1 ^ {(0)}, p _ 2 ^ {(0)}\right)=(1, 0, 0, 1)$ とすると、厳密解は円軌道になる。シンプレクティック・オイラー法は円軌道のときにだけ軌道が閉じる。これは円軌道においてルンゲ=レンツ・ベクトルの大きさが0になるからである。

以下に $\tau=0.001$ で計算した1万ステップ分の軌道とハミルトニアンの時間変化を示す。

2枚目は縦の縮尺が非常に拡大されているため分かりづらくなっているが、ハミルトニアンの振幅は縦軸の上にも書いてあるように $10 ^ {-7}$ 程度である。

楕円軌道

円軌道のときと同様、$\mu=1$、$\tau=0.001$ とする。初期条件 $\left(q _ 1 ^ {(0)}, q _ 2 ^ {(0)}, p _ 1 ^ {(0)}, p _ 2 ^ {(0)}\right)=(1, 0, 0, 1)$ のうち $p _ 2 ^ {(0)}$ だけを少しずつ減少させていき近心点を中心天体に近づけると、それにつれて軌道形状の誤差は大きくなっていく。

$p _ 2 ^ {(0)}=0.5, 0.3, 0.18$ の3つの初期条件について3万ステップ分の軌道を以下に示した。ルンゲ=レンツ・ベクトルの図示も試みたが、そもそも軌道全体の誤差が小さい場合にも大きく振動したため、あまり意味がないと判断した。しかし、$p _ 2 ^ {(0)}$ が小さくなるにつれて楕円軌道の向きが少しずつずれていることから、軌道が安定しないという意味でもルンゲ=レンツ・ベクトルが保存していないことは見て取れる。

4．ケプラー方程式

楕円軌道の長軸半径を $a$ 、短軸半径を $b$ とする。また、楕円 $C$ 上の真近点角 $\nu$ の点 $A$ に対し、$C$ と長軸を共有する円 $C'$ に $A$ を長軸に垂直に投影したときの点 $A'$ を考えて、この点 $A'$ とルンゲ=レンツ・ベクトルとの角度を $E$ と置く。 $E$ を離心近点角と言う。

楕円の幾何学的な性質として、次が成立する。

\begin{aligned}|r|\cos{\nu} &= a\left(\cos{E}-e\right) \\|r|\sin{\nu} &= b\sin{E} = a\sqrt{1-e^2}\sin{E}\end{aligned}

この2式を微分し、適切な係数をかけてから引き算すると、$r ^ 2\dot{\nu}$ の表現が得られる。$r ^ 2\dot{\nu}$ は角運動量の極座標表示であるので、

r^2\dot{\nu} = a^2\sqrt{1-e^2}\left(1-e\cos{E}\right)\dot{E} = h

を得る。$E$ を左辺に集めて時間 $t$ で積分すると、

E-e\sin{E}=\frac{h}{a^2\sqrt{1-e^2}}\left(t-t_p\right)

これがケプラー方程式である。$t _ p$ は積分定数だが、左辺が 0 になるとき $t=t _ p$ なので、近心点通過時刻であるとわかる。$a$ は軌道形状の極座標表示から $a=\frac{\rho}{1-e ^ 2}$ と計算できて、$h, e, \rho$ はすべて初期条件だけから計算できる保存量だから、ケプラー方程式はこのまま求根できて、時刻 $t$ からその時の離心近点角 $E$ が求まる。

$E$ が求まれば、この節の最初の2式が質点の $q _ 1=|r|\cos{\nu}$ と $q _ 2=|r|\sin{\nu}$ を表しているのだから、その右辺に代入してやれば良い。

$\mu = 1$、$h=0.001$、初期条件 $\left(q _ 1 ^ {(0)}, q _ 2 ^ {(0)}, p _ 1 ^ {(0)}, p _ 2 ^ {(0)}\right)=(1, 0, 0, 0.18)$ で10万ステップ分のケプラー方程式を解いた。軌道形状を解析的に計算した上での数値計算なので軌道が閉じることは自明であり、普通に軌道を図示しても意義が薄いことから、プロットは $|r|$-$t$ グラフとし、シンプレクティック・オイラー法によるプロットに重ねて表示した。求根には Python の SciPy ライブラリの関数 scipy.optimize.newton を使用した。

反復数が多いため分かりづらいが、シンプレクティック・オイラー法の結果は遠心点が少しずつ遠ざかっている。また、位相もケプラー方程式に比べて $t=100$ 近くでは半分くらいずれてしまっている。

5．レヴィ=チヴィタ変換の方法

レヴィ=チヴィタ変換

この節の記述はほとんど Celletti [2] に基づいている。

まず次の正準変換 $\left(q _ 1, q _ 2\right) → \left(Q _ 1, Q _ 2\right)$ を行う。

\begin{aligned}q_1 &= {Q_1}^2 -{Q_2}^2 \\q_2 &= 2Q_1Q_2\end{aligned}

倍角公式の形になっているので、これは角度変数を $\nu → \nu/2$ としているのに等しい。こうすると軌道が平らに開かれるので、近心点高度が低ければ低いほど近心点付近での動き方は穏やかになる。また、長さは平方根になっているので、$\sqrt{{q _ 1} ^ 2+{q _ 2} ^ 2}={Q _ 1} ^ 2+{Q _ 2} ^ 2$ である。この値を後の計算のために $D/4$ と置いておく。新しいハミルトニアンは

H' = \frac{P_1^2+P_2^2}{2D} - \frac{\mu}{D/4}

となる。正準方程式は、

\frac{\mathrm{d}}{\mathrm{d}t}\begin{bmatrix}Q_k\\P_k\end{bmatrix}=\frac{1}{D}\begin{bmatrix}P_k\\8Q_i\left(\frac{P_1^2+P_2^2}{2D}-\frac{4\mu}{D}\right)\end{bmatrix}=\frac{1}{D}\begin{bmatrix}P_k\\8Q_i\varepsilon\end{bmatrix}

となる。ただし、最後に $H'$ と同じ形が出た部分をエネルギー $\varepsilon$ に置き換えた。

次に原点での特異性を除去するため、正準方程式を形式上で $D$ 倍したい。運動を不変に保つには、時間変数 $t$ を $\mathrm{d}s = (1/D)\mathrm{d}t$ なる $s$ に置き換えればいい。すると、最終的に解くべき正準方程式は次のようになる。

\frac{\mathrm{d}}{\mathrm{d}s}\begin{bmatrix}Q_k\\P_k\end{bmatrix}=\begin{bmatrix}P_k\\8\varepsilon Q_i\end{bmatrix}

これは角周波数 $\omega = \sqrt{-8\varepsilon}$ の調和振動子である^{[† 2]}。以上がレヴィ=チヴィタ変換である。

離散化

[1] によると、この調和振動子を陰的中点法で離散化して計算してからレヴィ=チヴィタ逆変換をかけてケプラー問題の解として読むと、ケプラー問題の3つの保存量がすべて保存される。それを確かめよう。

時間刻み幅を $\sigma$ とすると、陰的中点法の更新式は以下の連立一次方程式である。

Q_k^{(i+1)}- \frac{\sigma}{2}P_k^{(i+1)}=Q_k^{(i)}+\frac{\sigma}{2}P_k^{(i)} \\\,\\P_k^{(i+1)}- 4\sigma\varepsilon Q_k^{(i+1)}=P_k^{(i)}+4\sigma\varepsilon Q_k^{(i)}

解析的に解けるので解いてしまって、その答え

\begin{aligned}Q_k^{(i+1)} &= \frac{\left(1+2\sigma^2\varepsilon\right)Q_k^{(i)}+\sigma P_k^{(i)}}{1-2\sigma^2\varepsilon}\\P_k^{(i+1)} &= \frac{\left(1+2\sigma^2\varepsilon\right)P_k^{(i)}+8\sigma\varepsilon Q_k^{(i)}}{1-2\sigma^2\varepsilon}\end{aligned}

を数値計算することにする。$\mu=1$、$h=0.001$、初期条件 $\left(q _ 1 ^ {(0)}, q _ 2 ^ {(0)}, p _ 1 ^ {(0)}, p _ 2 ^ {(0)}\right)=(1, 0, 0, 0.2)$ としたときの10万ステップ分の軌跡を、シンプレクティック・オイラー法のものに重ねて表示したのが下図である。

レヴィ=チヴィタ変換による方法のプロットは散布図として設定したが、解が閉じていることや時間刻み幅が非常に小さいこともあり、ほぼつながった線に見えている。非常に高い精度で安定していることが分かる。一方、すでに見たように、シンプレクティック・オイラー法はこれほど近心点高度が低いとルンゲ=レンツ・ベクトルの向きが安定しない。

6．参考文献

[1] Minesaki, Y. and Nakamura, Y., A new discretization of the Kepler motion which conserves the Runge-Lenz vector, Phys. Lett. A, 306 (2002), 127-133.

[2] Celletti, A., Basics of regularization theory, in Chaotic Worlds: From Order to Disorder in Gravitational N-Body Dynamical Systems, Springer, Dordrecht (2006), 203-230.

7．脚注

† 1: ^ 3つの保存量が独立でなくなっているように見えるが、そうではない。ここでは軌道面が動かないことを前提して2次元平面で計算しているので、角運動量ベクトルの保存のうち方向の保存は問題設定中に織り込まれているのである。しかるにここで登場する $h$ は角運動量ベクトルの大きさであるから、$P$ と $\varepsilon$ が保存しても3次元空間で角運動量ベクトルの方向が保存するかどうかは以下の式のみからは分からない。
ただし、次のことは言える：2次元平面でケプラー問題を計算する数値解法について、それがハミルトニアンとルンゲ=レンツ・ベクトルを保存するなら、角運動量も必ず保存する。

† 2: ^ ケプラー問題で軌道が閉じているときエネルギーは負であるから、これでいい。