万有引力の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\) 軸方向にとると，運動は \(xy\) 平面上に制限されて
$$\boldsymbol{r} = (x, y, 0)$$とおける。
初期条件で表される保存量（運動の定数）と楕円軌道要素
数値計算をする際，\(t = 0\) で \(x = x_0, \ y = 0\) という \(x\) 軸上から，\(v_x = 0\) , \( v_y = v_{y0}\) という\(y\) 軸方向への初速度ではじめることにする。
初速度 \(v_{y0}\) が以下の式で与えられる \(v_0\) のとき，軌道は円になる。
$$\frac{m v_0^2}{x_0} = \frac{GM m}{x_0^2}, \quad \therefore\ \ v_0 = \sqrt{\frac{GM}{x_0}}$$
なので，初速度 \(v_{y0}\) を \(v_0\) の定数倍というふうにおく。
$$v_{y0} = k v_0 = k \sqrt{\frac{GM}{x_0}}$$
運動の定数である \(\varepsilon\) や \(\ell = |\boldsymbol{\ell}|\) は，以下のように初期条件を使って表すことができる。
\begin{eqnarray}
\varepsilon &=& \frac{1}{2} \left(v_{y0}\right)^2 – \frac{GM}{x_0} \\
&=& – \frac{GM (2 – k^2)}{2 x_0} \\
\ell &=& x_0 v_{y0} = k x_0 v_0 \\
&=& k x_0 \sqrt{\frac{GM}{x_0}}
\end{eqnarray}
楕円軌道の長半径 \(a\) および離心率 \(e\) も，以下のように初期条件を使って表すことができる。
\begin{eqnarray}
a &=& \frac{GM}{2 |\varepsilon|} = \frac{x_0}{2 – k^2} \\
e &=& \sqrt{1 – \frac{2 |\varepsilon| \ell^2}{(GM)^2} } \\
&=& \sqrt{1 – k^2 (2-k^2)} = |k^2 – 1|
\end{eqnarray}
束縛軌道 \(\varepsilon < 0\) とするので，
$$ 0 < k < \sqrt{2}$$
無次元化
この系に特徴的な長さ \(x_0\) および時間スケール \(\displaystyle \tau \equiv \frac{x_0}{v_0} = \sqrt{\frac{x_0^3}{GM}}\) でこの系を以下のように規格化する。
\begin{eqnarray}
X &\equiv& \frac{x}{x_0} \\
Y &\equiv& \frac{y}{x_0} \\
T &\equiv& \frac{t}{\tau} = \frac{v_0 t}{x_0}
\end{eqnarray}
無次元化された運動方程式は
\begin{eqnarray}
\frac{d^2 X}{dT^2} = – \frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{d^2 Y}{dT^2} = – \frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}
無次元化された初期条件と楕円軌道要素
無次元化された変数に対する初期条件は \( T = 0\) で
$$X = 1, \quad Y = 0, \quad V_x = \frac{dX}{dT} = 0, \quad V_y = \frac{dY}{dT} = \frac{v_{y0}}{v_0} = k \ \ ( 0 < k < \sqrt{2})$$
このとき，数値計算する前にもう，軌道は規格化された長半径が
$$ A \equiv \frac{a}{x_0} = \frac{1}{2 – k^2} $$
離心率が
$$e = |k^2 – 1|$$
の楕円になることがわかってしまう。また，ケプラーの第3法則から，周期を \(p\) とすると
$$ \frac{a^3}{p^2} = \frac{GM}{4\pi^2}$$
より，規格化された周期が
$$P \equiv \frac{p}{\tau} = 2 \pi A^{\frac{3}{2}} = \frac{2 \pi}{\left(2 – k^2\right)^{\frac{3}{2}}}$$
となることもわかってしまっているのであった。
$k$ の範囲
なお，初速度の大きさを表す $k$ は束縛軌道であるという条件から $ 0 < k < \sqrt{2}$ であるが，特に断らない限り
$$ 1 \le k < \sqrt{2}$$
とする。この範囲にすることで初期位置が近点になる。なぜ近点になるかというと，$k=1$ で円軌道であり，$k > 1$ とすれば円軌道となる初速度よりも大きい初速度で $x$ 軸から打ち出されるので，反対側でより膨れる位置で $x$ 軸を横切ることが予想されるから。
長半径 $a$ を一定に保ちながら離心率 $e$ を変化させる初期条件
上記のように，無次元化された変数に対する初期条件を \( T = 0\) で
$$X = 1, \quad Y = 0, \quad V_x = \frac{dX}{dT} = 0, \quad V_y = \frac{dY}{dT} = k$$
とすると，初速度 $k$ の取り方しだいで楕円軌道になるのであるが，その際，離心率だけでなく軌道長半径 $a$ も $k$ に依存して変化し，したがって周期も $k$ に依存するようになる。
それはそれでおもしろいのであるが，もう少し工夫すると，軌道長半径 $a$ を一定に保ちながら離心率 $e$ が変化するような初期条件を設定することができる。$a$ が一定のままであるから周期も変化しない。（規格化された周期 $P$ は $P = 2 \pi$ のまま。）
初期条件を
$$ x_i = f x_0, \quad y_i = 0, \quad v_{x 0} = 0, \quad v_{y0} = k v_0$$
とおく。これは無次元化された変数に対しては
$$X = f, \quad Y = 0, \quad V_x = \frac{dX}{dT} = 0, \quad V_y = \frac{dY}{dT} = k$$
としたことに対応する。運動の定数 $\varepsilon$ が長半径 $a$ で書けることから，
\begin{eqnarray}
\varepsilon \equiv – \frac{GM}{2 a} &=& \frac{1}{2} v_{y0}^2 – \frac{GM}{x_i} \\
&=& \frac{1}{2} (k v_0)^2 – \frac{GM}{f x_0} \\
&=& \frac{1}{2} k^2 \frac{GM}{x_0} – \frac{GM}{f x_0} \\
&=& \frac{GM}{x_0} \left( \frac{k^2}{2} – \frac{1}{f}\right)
\end{eqnarray}
円軌道の初期条件 $k = 1, \ f = 1$ のときは $\displaystyle \frac{a}{x_0} = 1$ となる。 $k > 1$ のときも $\displaystyle \frac{a}{x_0} = 1$ となるようにするには，
\begin{eqnarray}
– \frac{1}{2} &=& \frac{k^2}{2} – \frac{1}{f} \\
\therefore\ \ f &=& \frac{2}{k^2 + 1}
\end{eqnarray}
とすればよい。このときの離心率を求めるにはまず，運動の定数 $\ell$ を初期条件から求めて，
\begin{eqnarray}
\ell &\equiv& x \frac{dy}{dt} – y \frac{dx}{dt} \\
&=&x_i \times v_{y0} \\
&=& f x_0 \times k v_0 \\
&=& \frac{2 k}{k^2 + 1} x_0 \sqrt{\frac{GM}{x_0}}
\end{eqnarray}
したがって
\begin{eqnarray}
e &=& \sqrt{1 – \frac{2 |\varepsilon| \ell^2}{(GM)^2}} \\
&=& \sqrt{1 – \frac{2}{(GM)^2} \frac{GM}{2a} \left(\frac{2k}{k^2 + 1}\right)^2 x_0^2 \frac{GM}{x_0} } \\
&=& \sqrt{1 – \left(\frac{2k}{k^2 + 1}\right)^2} \\
&=& \frac{|k^2 – 1|}{k^2 + 1}
\end{eqnarray}
$k$ の範囲
この場合も，初期位置が近点となるように $ k \ge 1$ とすればよい。
$$\lim_{k \rightarrow \infty} e = 1$$
なので，$k$ の上限は実質上ないことになる。
$e$ の式を $k$ について解くと，
\begin{eqnarray}
e &=&\frac{k^2 – 1}{k^2 + 1} \\
\therefore\ \ k &=& \sqrt{\frac{1+e}{1-e}}
\end{eqnarray}