「万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ」の改良版。変数の規格化をもう少しシンプルに。
万有引力の2体問題（は結局1体問題に帰着するんだけど）の運動方程式を数値的に解く前の下ごしらえ。闇雲な初期条件からはじめるのではなく，そもそも数値計算しなくても楕円になることはわかっているのだから，ルンゲ・クッタ法なりを使って数値的に解く前に，それなりの下ごしらえをしておこう。
運動方程式
万有引力の2体問題は相対位置ベクトル \(\boldsymbol{r}\) と全質量 \(M\) を使って書くと結局1体問題に帰着して
$$\frac{d^2 \boldsymbol{r}}{dt^2} = – \frac{GM}{r^3} \boldsymbol{r}$$
保存量（運動の定数）
運動方程式から得られる保存量は2つ。（単位質量当たりの）角運動量 \(\boldsymbol{\ell}\)
$$\boldsymbol{\ell} \equiv \boldsymbol{r} \times \dot{\boldsymbol{r}} = \mbox{const.}$$
と（単位質量当たりの）全エネルギー \(\varepsilon\)
$$\varepsilon \equiv \frac{1}{2} \dot{\boldsymbol{r}} \cdot\dot{\boldsymbol{r}} – \frac{GM}{r}$$
一定のベクトル \(\boldsymbol{\ell}\) の向きを \(z\) 軸方向にとると，
$$ \boldsymbol{\ell} = (0, 0, \ell)$$
運動は \(xy\) 平面上に制限されて
$$\boldsymbol{r} = (x, y, 0)$$とおける。
保存量（運動の定数）と楕円軌道要素
数値計算をするまでもなく，この解は楕円軌道になることがわかっている。$xy$ 平面上の極座標 $r, \phi$ であらわすと，
$$x = r \cos\phi, \quad y = r \sin \phi$$
$$r = \frac{a (1-e^2)}{1 + e \cos\phi}$$
ここで，$a$ は軌道長半径，$e$ は離心率。$\phi = 0 $ で $r = r_{\rm min} = a (1-e)$ なる近点距離になるように積分定数を選んでいる。
運動の定数である \(\varepsilon\) や \(\ell\) は，以下のように楕円軌道要素を使って表すことができる。
\begin{eqnarray}
\varepsilon &=& – \frac{GM}{2 a} \\
\ell &=& \sqrt{GM a (1-e^2)}
\end{eqnarray}
軌道長半径 $a$ と周期 $P$ の間にはケプラーの第3法則が成り立ち，
$$ \frac{a^3}{P^2} = \frac{GM}{4 \pi^2}, \quad \therefore\ \ P = 2 \pi \frac{a \sqrt{a}}{GM}$$
無次元化
この系に特徴的な長さである軌道長半径 \(a\) および時間スケールである周期 \(P\) でこの系を以下のように規格化する。
\begin{eqnarray}
X &\equiv& \frac{x}{a} \\
Y &\equiv& \frac{y}{a} \\
T &\equiv& \frac{t}{P}
\end{eqnarray}
無次元化された運動方程式は
\begin{eqnarray}
\frac{d^2 X}{dT^2} = – 4\pi^2 \frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \tag{1}\\
\frac{d^2 Y}{dT^2} = – 4\pi^2 \frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \tag{2}
\end{eqnarray}
初期条件
初期条件を \( t = 0\) で
\begin{eqnarray}
x(0) &=& r_{\rm min} = a(1-e) \\
y(0) &=& 0 \\
v_x(0) &=& 0 \\
v_y(0) &=& v_{y0}
\end{eqnarray}
とする。ここで $v_{y0}$ は $t = 0$ での全エネルギーの式から以下のようになる。
\begin{eqnarray}
\varepsilon = – \frac{GM}{2a} &=& \frac{1}{2} \left(v_{y0}\right)^2 – \frac{GM}{a(1-e)} \\
\therefore\ \ v_{y0} &=& \sqrt{\frac{GM (1+e)}{a(1-e)}}
\end{eqnarray}
無次元化された初期条件
$T = t/P = 0$ のとき
\begin{eqnarray}
X(0) &=& \frac{x(0)}{a} = 1-e \\
Y(0) &=& \frac{y(0)}{a} = 0 \\
V_x(0) &=& \frac{P}{a} v_x(0) = 0 \\
V_y(0) &=& \frac{P}{a} v_{y0} = 2 \pi \sqrt{\frac{1+e}{1-e}}
\end{eqnarray}
この初期条件で，式 $(1)$, $(2)$ を数値的に解くと，（規格化された）長半径 $1$，離心率 $e$ の楕円の軌道が得られるハズ。
じゃあ，そこまでわかっているのなら，なぜわざわざ数値計算するのかというと，それは楕円軌道の解が時間 $t$ の陽関数として与えられているわけではないから。時刻 $t$ のとき，どこにいるかが解析的にわかっていないので，それを知りたいから数値計算することになる。