楕円軌道上の時刻ごとの位置を Maxima で数値的に求める

「楕円軌道上の時刻ごとの位置を求めるための下ごしらえ」で下ごしらえした式を Maxima で数値的に解く。

無次元化された式

「楕円軌道上の時刻ごとの位置を求めるための下ごしらえ」で下ごしらえした式（無次元化済み）。

\begin{eqnarray}
\frac{d\phi}{dT} &=& f(\phi, e) = 2 \pi \frac{(1+e \cos\phi)^2}{(1-e^2)^{3/2}} \tag{1}\\
R(\phi) &\equiv& \frac{r}{a} = \frac{1-e^2}{1+e\cos\phi} \tag{2}\\
X(\phi) &\equiv& \frac{x}{a} =R(\phi) \cos\phi \tag{3}\\
Y(\phi) &\equiv& \frac{y}{a} =R(\phi) \sin\phi \tag{4}
\end{eqnarray}

In [1]:
/* 微分方程式の右辺 */
f(phi, e):= 2*%pi*(1+e*cos(phi))**2/(1 - e**2)**(3/2);

R(phi, e):= (1 - e**2)/(1 + e*cos(phi));
X(phi, e):= R(phi, e) * cos(phi);
Y(phi, e):= R(phi, e) * sin(phi);
Out[1]:
\[\tag{${\it \%o}_{1}$}f\left(\varphi , e\right):=\frac{2\,\pi\,\left(1+e\,\cos \varphi\right)^2}{\left(1-e^2\right)^{\frac{3}{2}}}\]
Out[1]:
\[\tag{${\it \%o}_{2}$}R\left(\varphi , e\right):=\frac{1-e^2}{1+e\,\cos \varphi}\]
Out[1]:
\[\tag{${\it \%o}_{3}$}X\left(\varphi , e\right):=R\left(\varphi , e\right)\,\cos \varphi\]
Out[1]:
\[\tag{${\it \%o}_{4}$}Y\left(\varphi , e\right):=R\left(\varphi , e\right)\,\sin \varphi\]

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

In [2]:
/* T の初期値 */
T0: 0$
/* phi の初期値 */
phi0: 0$
/* T の終了値 */
Tend: 1$
/* 分割数 */
ndiv: 100$
Ndiv: 36$
N: Ndiv * ndiv$
/* 刻み幅 h は以下のようにして計算される。*/
h: float((Tend - T0)/N)$

/* 1階上微分方程式を Runge-Kutta 法で解く */
/* e = 0.6 */
e: 0.6$
rkphi: rk(f(phi, e), phi, phi0, [T, T0, Tend, h])$

ちょうど1周期 $T = T_{\rm end}$ たてば $\phi = 2 \pi$ になっているはずだから，数値誤差を確認する。

$$\phi(T_{\rm end}) – 2 \pi$$

In [3]:
rkphi[length(rkphi)][2]- float(2*%pi);
Out[3]:
\[\tag{${\it \%o}_{14}$}-2.303579549334245 \times 10^{-11}\]

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

一定時間間隔 $\displaystyle \Delta T = \frac{T_{\rm end}}{N_{\rm div}}$ ごとの位置：

In [4]:
phii: makelist(rkphi[i][2], i, 1, length(rkphi), ndiv)$

pos: makelist([X(phii[i], e), Y(phii[i], e)], i, length(phii))$

plot2d() でグラフ作成

In [5]:
plot2d([discrete, pos], [style, [points,1]], 
       [same_xy], [x, -2, 1], [y, -1.5, 1.5])$

draw2d() でグラフ作成

In [6]:
draw2d(
  /* 点を塗りつぶした丸に */
  point_type = 7, 
  /* 表示範囲 */
  xrange = [-2, 1], yrange = [-1.5, 1.5],
  /* 縦横比 */
  proportional_axes=xy,
  /* x 軸 y 軸 */
  xaxis = true, yaxis = true,   
  
  /* 点の大きさ */
  point_size = 0.5, 
  points(pos)
)$

参考：ケプラー方程式の数値解との比較

「ケプラー方程式を数値的に解いてケプラーの第2法則を視覚的に確認する」にまとめているように，楕円軌道の $x, y$ を軌道長半径 $a$ で規格化した $X_u, Y_u$ を，離心近点離角 $u$ を使って以下のようにあらわすことができる。

\begin{eqnarray}
X_u &\equiv& \frac{x}{a} = \cos u – e\\
Y_u &\equiv& \frac{y}{a} = \sqrt{1-e^2} \sin u
\end{eqnarray}

離心近点離角 $u$ と周期 $P$ で規格化された時間 $T$ との関係は，以下のケプラー方程式を数値的に解いて得られる。

\begin{eqnarray}
2\pi T &=& u – e \sin u, \quad T \equiv \frac{t}{P}
\end{eqnarray}

ケプラー方程式の数値解

$$g(u_i) \equiv u_i – e \sin u_i = \frac{2\pi}{T} t_i = \frac{2\pi}{N_{\rm div}} \times i$$

を $u_i$ について find_root() 関数で数値的に解きます。

In [7]:
g(u):= u - e*sin(u);
Out[7]:
\[\tag{${\it \%o}_{21}$}g\left(u\right):=u-e\,\sin u\]
In [8]:
for i:0 thru Ndiv do 
    ui[i]: find_root(g(u) = 2*%pi/Ndiv * i, u, 0, 2*%pi)$

一定時間間隔ごとの位置

In [9]:
Xu(u, e):= (cos(u) - e);
Yu(u, e):= sqrt(1-e**2)*sin(u);
Out[9]:
\[\tag{${\it \%o}_{23}$}{\it Xu}\left(u , e\right):=\cos u-e\]
Out[9]:
\[\tag{${\it \%o}_{24}$}{\it Yu}\left(u , e\right):=\sqrt{1-e^2}\,\sin u\]
In [10]:
posu: makelist([Xu(ui[i], e), Yu(ui[i], e)], i, 0, Ndiv)$

2つの解の比較

$\phi$ に関する1階常微分方程式を数値的に解いた解である pos と，ケプラー方程式を数値的に解いた解である posu は，当然ながら数値誤差の範囲で一致してるはずである。

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

In [11]:
/* 2つの数値解 X, Y の差 */
(pos-posu)[10];
(pos-posu)[20];
(pos-posu)[30];
Out[11]:
\[\tag{${\it \%o}_{26}$}\left[ 1.530775506353166 \times 10^{-12} , 7.017719738655614 \times 10^{-13} \right] \]
Out[11]:
\[\tag{${\it \%o}_{27}$}\left[ -1.576516694967722 \times 10^{-13} , 1.149330630667578 \times 10^{-12} \right] \]
Out[11]:
\[\tag{${\it \%o}_{28}$}\left[ -1.973976537783528 \times 10^{-12} , 3.774758283725532 \times 10^{-13} \right] \]

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

In [12]:
draw2d(
  /* 点を塗りつぶした丸に */
  point_type = 7, 
  /* 表示範囲 */
  xrange = [-2, 1], yrange = [-1.5, 1.5],
  /* 縦横比 */
  proportional_axes=xy,
  /* x 軸 y 軸 */
  xaxis = true, yaxis = true,   

  point_size = 0.5, color = red, key = "常微分方程式の数値解",
  points(pos), 
  
  point_size = 0.2, color = blue, key = "ケプラー方程式の数値解",
  points(posu)
)$