万有引力の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\) 軸方向にとると,運動は \(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}}}$$

となることもわかってしまっているのであった。