下ごしらえ（ver.2）した万有引力の2体問題の運動方程式を Maxima で数値的に解く

無次元化された運動方程式

「万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ ver. 2」にまとめたように，軌道長半径 $a$ と周期 $P$ で無次元化した変数に対する運動方程式は以下のようになるのであった。

\begin{eqnarray}
\frac{d^2 X}{dT^2} = – 4\pi^2\frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{d^2 Y}{dT^2} = – 4\pi^2\frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}

無次元化された初期条件

無次元化された変数に対する初期条件は $T = 0$ で

\begin{eqnarray}
X(0) &=& 1-e \\
Y(0) &=& 0\\
V_x(0) &=& 0 \\
V_y(0) &=& 2\pi\sqrt{\frac{1+e}{1-e}}
\end{eqnarray}

この初期条件で数値計算すると，（規格化された）長半径 $1$，離心率 $e$ の楕円の軌道が得られるハズ。

じゃあ，そこまでわかっているのなら，なぜわざわざ数値計算するのかというと，それは楕円軌道の解が時間 $t$ の陽関数として与えられているわけではないから。時刻 $t$ のとき，どこにいるかが解析的にわかっていないので，それを知りたいから数値計算することになる。

Runge-Kutta 用に連立1階微分方程式系に

\begin{eqnarray}
\frac{dX}{dT} &=& F_1(V_x) = V_x \\
\frac{dY}{dT} &=& F_2(V_y) = V_y \\
\frac{dV_x}{dT} &=& F_3(X, Y)
= -4\pi^2 \frac{X}{\left(X^2 + Y^2\right)^{\frac{3}{2}}} \\
\frac{dV_y}{dT} &=& F_4(X, Y)
= -4\pi^2 \frac{Y}{\left(X^2 + Y^2\right)^{\frac{3}{2}}}
\end{eqnarray}

rk() で連立常微分方程式を数値的に解く

ちょうど1周期 $T=t/P = 1$ までいくと，$T=0$ の時の初期値になっているはずだから，数値誤差を確認する。

rksol[length(rksol)] はリスト rksol の最後の項，rksol[1] は最初の項。

rksol は順番に [T, x, y, v, w] の5つの項からなるから，rksol[1] の x の値は rksol[1][2] などとして参照する。

一定時間間隔ごとの位置をグラフに

plot2d() でグラフ作成

draw2d() でグラフ作成

ケプラー方程式を数値的に解いた場合の解との比較

「ケプラー方程式を数値的に解いてケプラーの第2法則を視覚的に確認する」にまとめているように，

\begin{eqnarray}
x &=& r \cos\varphi = a(\cos u – e)\\
y &=& r \sin\varphi = a \sqrt{1-e^2} \sin u\\
\frac{2\pi}{P} t &=& u – e \sin u
\end{eqnarray}

だから，無次元化した変数では

\begin{eqnarray}
X_u(u) &=& \frac{x}{a} = \cos u – e\\
Y_u(u) &=& \frac{y}{a} = \sqrt{1-e^2} \sin u\\
2\pi T &=& u – e \sin u, \quad T = \frac{t}{P}
\end{eqnarray}

ケプラー方程式の数値解

$$g(u) \equiv u – e \sin u$$

として $T_i = \frac{i}{N_{\rm div}}$ のときにケプラー方程式 $g(u_i) = 2\pi T_i $ を数値的に解いて $u_i$ を求める。

一定時間間隔ごとの位置

2つの解の比較

運動方程式を数値的に解いた解である xy_rk とケプラー方程式を数値的に解いた解である xy_kep は，当然ながら数値誤差の範囲内で一致しているはずである。

2つのリストの引き算をして差を調べる。37組全て表示すると冗長なので，代表して3組ほど。

念のために2つの数値解をグラフにすると，当然ながら重なっていることがわかる。