Maxima でケプラー方程式の近似解と数値解

ベッセル関数によるケプラー方程式のフーリエ級数解

ケプラー方程式

$$u – e \sin u = \omega t$$

のフーリエ級数解は，ベッセル関数 $J_n(x)$ を使って以下のように書ける。

$$u = \omega t+ \sum_{n=1}^{\infty} \frac{2}{n} J_n(n e) \sin(n\, \omega t)$$

現実世界では無限大まで足し算はできないので，以下のように $N$ までの和をとることにする。

$$u(N, e, \omega t) = \omega t+ \sum_{n=1}^{N} \frac{2}{n} J_n(n e) \sin(n\, \omega t)$$

これを Maxima の関数として定義する。

In [1]:
Ubes(N, e, omegat):= 
  omegat + sum(2/n * bessel_j(n, n*e) * sin(n * omegat), n, 1, N);
Out[1]:
\[\tag{${\it \%o}_{1}$}{\it Ubes}\left(N , e , {\it omegat}\right):={\it omegat}+{\it sum}\left(\frac{2}{n}\,{\it bessel\_j}\left(n , n\,e\right)\,\sin \left(n\,{\it omegat}\right) , n , 1 , N\right)\]
In [2]:
/* たとえば N=10, e=0.5, omegat = 1.0 の値 */
Ubes(10, 0.5, 1.0);
Out[2]:
\[\tag{${\it \%o}_{2}$}1.49885975062147\]

逐次近似法によるケプラー方程式の近似解

離心率 $e$ は $0 \leq e < 1$ であることから，ケプラー方程式
$$u – e \sin u = \omega t$$
を以下のように逐次近似的に解いてみます。

\begin{eqnarray}
u &=& \omega t + e \sin u \\
u_0 &=& \omega t\\
u_1 &=& \omega t + e \sin u_0 = \omega t + e \sin \omega t\\
u_2 &=& \omega t + e \sin u_1 =\omega t + e \sin\left(\omega t + e \sin \omega t\right) \\
u_3 &=& \omega t + e \sin u_2 =\omega t + e \sin\left\{\omega t + e \sin\left(\omega t + e \sin \omega t\right)\right\} \\
&\vdots&\\
u_{n} &=& \omega t + e \sin u_{n-1} = \dots\\
\end{eqnarray}

$n$ が大きくなると，入れ子になっている項がどんどん増殖していきますが，$u_3$ のあたりまでは，近似的に $u$ は $t$ の陽関数としてあらわされているなぁ… という見た目がします。

上の式にそって，以下のように関数 $u(n, e, \omega t)$ を再帰的に定義します。

In [3]:
Uite(N, e, omegat):=
block(  
  if N = 0 then
    omegat
  else
    omegat + e*sin(Uite(N-1, e, omegat))
)$
In [4]:
/* たとえば N=10, e=0.5, omegat = 1.0 の値 */
Uite(10, 0.5, 1.0);
Out[4]:
\[\tag{${\it \%o}_{4}$}1.498701133517836\]

方程式の数値解法によるケプラー方程式の数値解

$\omega t$ を与えたときに，ケプラー方程式を数値的に解いて $u(\omega t)$ を求める。ケプラー方程式は超越方程式であるので，find_root() 関数を使って数値に解く。

In [5]:
Unum(e, omegat):= 
  block(find_root(u - e*sin(u) = omegat, u, 0, 2*%pi));
Out[5]:
\[\tag{${\it \%o}_{5}$}{\it Unum}\left(e , {\it omegat}\right):=\mathbf{block}\;\left({\it find\_root}\left(u-e\,\sin u={\it omegat} , u , 0 , 2\,\pi\right)\right)\]
In [6]:
/* たとえば e=0.5, omegat = 1.0 の値 */
Unum(0.5, 1.0);
Out[6]:
\[\tag{${\it \%o}_{6}$}1.498701133517848\]

数値解と近似解との比較

In [7]:
omt: makelist(float(%pi/10*i), i, 1, 9)$
In [8]:
/* 数値解 Unum と逐次近似解 Uite, フーリエ級数解 Ubes の差 */
e: 0.5$
N: 10$
printf(true," omega t    Unum               Uite-Unum          Ubes-Unum ~%")$

for i thru length(omt) do(
  printf(true, "~8,5f ~18,15f ~18,15f ~18,15f~%", 
                omt[i],
                Unum(e, omt[i]), 
                Uite(N, e, omt[i])-Unum(e, omt[i]), 
                Ubes(N, e, omt[i])-Unum(e, omt[i]))
)$

 omega t    Unum               Uite-Unum          Ubes-Unum 
 0.31416  0.593999023813608 -0.000048370871788  0.000205334548063
 0.62832  1.065940683889791 -0.000000479808530 -0.000233511559215
 0.94248  1.438080909968085 -0.000000000003061  0.000199446542732
 1.25664  1.748741781633489  0.000000000005287 -0.000159160439054
 1.57080  2.020979938089770 -0.000000056633642  0.000123562038969
 1.88496  2.268208852924498 -0.000003478720316 -0.000093184551055
 2.19911  2.498822425235399 -0.000028291765788  0.000066804145119
 2.51327  2.718544855625697 -0.000076315203530 -0.000043154412396
 2.82743  2.931640124182721 -0.000080693322920  0.000021180282106

数値解との差の二乗平均平方根

いわゆる RMS (root mean square)

$$RMS \equiv \sqrt{\frac{1}{n} \sum_{i=1}^n \left( U(\omega t_i)- U_{\rm num}(\omega t_i)\right)^2}$$

In [9]:
RMSite(N, e):= 
  sqrt(1/length(omt) * 
       sum((Uite(N, e, omt[i])-Unum(e, omt[i]))**2, i, 1, length(omt)))$

RMSbes(N, e):=
  sqrt(1/length(omt) * 
       sum((Ubes(N, e, omt[i])-Unum(e, omt[i]))**2, i, 1, length(omt)))$
In [10]:
/* 浮動小数点表示の書式指定 */
printf(true, "~,15e~%", RMSite(10, 0.5))$
printf(true, "~,15e~%", RMSbes(10, 0.5))$

4.148349447033673e-5
1.462591053867549e-4

グラフ

$a = 5, \ e = 0.6$ の場合の楕円のグラフを描く。

\begin{eqnarray}
x(a, e, \omega t) &=& a (\cos u(\omega t) -e ) \\
y(a, e, \omega t) &=& a \sqrt{1 – e^2} \sin u(\omega t) \\
\end{eqnarray}

In [11]:
x(a, e, omegat):= 
  a*(cos(Uite(10, e, omegat)) - e);
y(a, e, omegat):= 
  a*sqrt(1-e**2)*sin(Uite(10, e, omegat));
Out[11]:
\[\tag{${\it \%o}_{16}$}x\left(a , e , {\it omegat}\right):=a\,\left(\cos {\it Uite}\left(10 , e , {\it omegat}\right)-e\right)\]
Out[11]:
\[\tag{${\it \%o}_{17}$}y\left(a , e , {\it omegat}\right):=a\,\sqrt{1-e^2}\,\sin {\it Uite}\left(10 , e , {\it omegat}\right)\]
In [12]:
/* plot2d 版 */
plot2d([parametric, x(5, 0.6, t), y(5, 0.6, t), [t, 0, 2*%pi]], 
       [x,-9, 3], [y, -6, 6], [same_xy])$

In [13]:
/* draw2d 版 */
draw2d(
  /* 滑らかに曲線を描くように */
  nticks = 300, 
  /* 表示範囲 */
  xrange = [-9, 3], yrange = [-6, 6],
  /* 縦横比。*/
  proportional_axes = xy,
  /* x 軸 y 軸の表示 */
  xaxis = true, yaxis = true, 

  parametric(x(5, 0.6, omegat), y(5, 0.6, omegat), omegat, 0, 2*%pi)
)$

In [14]:
pxy: makelist([x(5, 0.6, float(%pi/18*i)), 
               y(5, 0.6, float(%pi/18*i))], i, 1, 36)$
In [15]:
/* plot2d 版 */
plot2d([discrete, pxy], 
       [style, [points, 1.5]],
       [x,-9, 3], [y, -6, 6], [same_xy])$

In [16]:
/* draw2d 版 */
draw2d(
  /* 表示範囲 */
  xrange = [-9, 3], yrange = [-6, 6],
  /* 縦横比。*/
  proportional_axes = xy,
  /* x 軸 y 軸の表示 */
  xaxis = true, yaxis = true, 

  point_type = 7,
  point_size = 1,
  points(pxy)
)$