「万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ」で下ごしらえした運動方程式を Maxima で Runge-Kutta 法を使って数値的に解くという話。
無次元化された運動方程式
「万有引力の2体問題の運動方程式を数値的に解く前の下ごしらえ」にまとめたように，系に特徴的な長さと時間スケールで無次元化した変数に対する運動方程式は以下のようになるのであった。
\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$ で
\begin{eqnarray}
X &=& 1 \\
Y &=& 0\\
V_x &=& \frac{dX}{dT} = 0 \\
V_y &=& \frac{dY}{dT} = \frac{v_{y0}}{v_0} = k \ \ ( 0 < k < \sqrt{2})
\end{eqnarray}
楕円の軌道要素（長半径，離心率，周期）
このとき，数値計算する前にもう，軌道は規格化された長半径が
$$ A \equiv \frac{a}{x_0} = \frac{1}{2 – k^2} $$
離心率が
$$e = |k^2 – 1|$$
の楕円になることがわかってしまう。また，規格化された周期が
$$P \equiv \frac{p}{\tau} = 2 \pi A^{\frac{3}{2}} = \frac{2 \pi}{\left(2 – k^2\right)^{\frac{3}{2}}}$$
となることもわかってしまっているのであった。
じゃあ，そこまでわかっているのなら，なぜわざわざ数値計算するのかというと，それは楕円軌道の解が時間 $t$ の陽関数として与えられているわけではないから。時刻 $t$ のとき，どこにいるかが解析的にわかっていないので，それを知りたいから数値計算することになる。
Runge-Kutta 用に連立1階微分方程式系に
全ての変数は無次元化されいること忘れないようにして，Maxima のコードにする際には，（見やすさの観点から）小文字で書いてみます。
\begin{eqnarray}
\frac{dx}{dt} &=& F_1(x, y, v, w, t) \Rightarrow F_1(v) = v \\
\frac{dy}{dt} &=& F_2(x, y, v, w, t) \Rightarrow F_2(w) = w \\
\frac{dv}{dt} &=& F_3(x, y, v, w, t) \Rightarrow F_3(x, y)
= -\frac{x}{\left(x^2 + y^2\right)^{\frac{3}{2}}} \\
\frac{dw}{dt} &=& F_4(x, y, v, w, t) \Rightarrow F_4(x, y)
= -\frac{y}{\left(x^2 + y^2\right)^{\frac{3}{2}}}
\end{eqnarray}
F1(v):= v;
F2(w):= w;
F3(x, y):= -x/(x**2 + y**2)**(3/2);
F4(x, y):= -y/(x**2 + y**2)**(3/2);
P(k):= float(2*%pi/(2 - k**2)**(3/2));
Runge-Kutta 法で解く
初期条件 $k = 1$ で円軌道となることを確認
まずは，$k = 1$ という初期条件で円軌道になることを確認。1周期を $N$ 等分し，$t = P(1) = 2 \pi$ まで計算したときの答えが $x_0 = 1, y_0 = 0$ になっているか。
fpprintprec: 5$ /* 浮動小数点表示は有効数字 5 桁に */
N: 36*20$
h: P(1)/N$
solc:
rk([F1(v), F2(w), F3(x, y), F4(x, y)],
[x, y, v, w], [1, 0, 0, 1],
[t, 0, P(1), h])$
print("N = ", N, " h = ", h)$
print("x の誤差 ", 1 - solc[length(solc)][2])$
print("y の誤差 ", 0 - solc[length(solc)][3])$
ちなみに，Maxima の rk() は4次の Runge-Kutta 法で解くので，誤差は $h^4$ の程度であることが予想される。
/* 4次の Runge-Kutta 法の誤差は h**4 の程度 */
h**4;
$e = 0.6$ となる初期条件で数値計算
初期条件 $k$ と楕円の離心率 $e$ との関係は $e = |k^2 – 1|$ であったから…
k: sqrt(1 + 0.6)$
N: 36*80$
h: P(k)/N$
sole:
rk([F1(v), F2(w), F3(x, y), F4(x, y)],
[x, y, v, w], [1, 0, 0, k],
[t, 0, P(k), h])$
print("N = ", N, " h = ", h)$
print("x の誤差 ", 1 - sole[length(sole)][2])$
print("y の誤差 ", 0 - sole[length(sole)][3])$
/* 4次の Runge-Kutta 法の誤差は h**4 の程度 */
h**4;
ケプラー方程式を数値的に解いた場合の解との比較
「ケプラー方程式を数値的に解いてケプラーの第2法則を視覚的に確認する」にまとめているように
/* 離心率 */
e: 0.6$
/* 分割数 */
N: 36$
/* ケプラー方程式を find_root で数値的に解く */
g(u):= u - e*sin(u)$
for i:0 thru N do
u[i]: find_root(g(u) = 2*%pi/N * i, u, 0, 2*%pi)$
/* t = 0 で xu = 1 ==> a = 1/(1-e) */
xu(u):= (cos(u) - e)/(1-e)$
yu(u):= sqrt(1-e**2)*sin(u)/(1-e)$
xy_kep: makelist([xu(u[i]), yu(u[i])], i, 0, N)$
/* Runge-Kutta 法で求めた解を間引く */
xy_rk: makelist([sole[i][2], sole[i][3]], i, 1, length(sole), 80)$
/* Runge-Kutta 法の解と，ケプラー方程式の数値解を使った解との差を */
/* 行列的表示で */
apply('matrix, xy_rk - xy_kep);
一定時間間隔ごとの位置をグラフに
draw2d(
title = "一定時間間隔ごとの位置",
/* 滑らかに曲線を描くように */
nticks = 100,
/* 表示範囲 */
xrange = [-4.5, 1.5], yrange = [-3, 3],
/* 縦横比。*/
proportional_axes = xy,
/* x 軸 y 軸の表示 */
xaxis = true, yaxis = true,
/* 軸のきざみ目盛の非表示 */
xtics = false, ytics = false,
/* 楕円*/
line_width = 1,
color = red,
parametric(xu(u), yu(u), u, 0, 2*%pi),
/* 周期の 1/N ごとの位置 */
point_type = 7,
point_size = 0.7,
color = blue,
points(xy_rk),
/* 原点 */
point_type = 7,
point_size = 1,
color = black,
points([[0, 0]])
)$
