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

「Maxima でケプラー方程式の近似解と数値解」の SymPy 版。

In [1]:
from sympy import *
from sympy.abc import *
from sympy import pi

# グラフの大きさの設定，inline svg 表示
import matplotlib.pyplot as plt
plt.rcParams['figure.figsize'] =8,6
plt.rcParams['font.size'] = 12
from IPython.display import set_matplotlib_formats
%matplotlib inline
set_matplotlib_formats('svg')

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

ケプラー方程式

$$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)$$

これを関数 Ubes(N, e, omegat) として定義する。

In [2]:
omegat = Symbol('omegat')

def Ubes(N, e, omegat):
    return(omegat + 
           summation(2/n * besselj(n, n*e) * sin(n * omegat), (n, 1, N))
          )
In [3]:
Ubes(2, e, omegat)
Out[3]:
$\displaystyle omegat + 2 \sin{\left(omegat \right)} J_{1}\left(e\right) + \sin{\left(2 omegat \right)} J_{2}\left(2 e\right)$
In [4]:
# たとえば N=10, e=0.5, omegat = 1.0 の値 
print(float(Ubes(10, 0.5, 1.0)))

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$ の陽関数としてあらわされているなぁ… という見た目がします。

上の式にそって，以下のように関数 Uite(N, e, omegat) を再帰的に定義します。

In [5]:
def Uite(N, e, omegat):
    if N ==0:
        return(omegat)
    else:
        return(omegat + e*sin(Uite(N-1, e, omegat)))
In [6]:
Uite(3, e, omegat)
Out[6]:
$\displaystyle e \sin{\left(e \sin{\left(e \sin{\left(omegat \right)} + omegat \right)} + omegat \right)} + omegat$
In [7]:
# たとえば N=10, e=0.5, omegat = 1.0 の値
print(float(Uite(10, 0.5, 1.0)))

1.4987011335178357

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

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

In [8]:
def Unum(e, omegat):
    return(nsolve(u-e*sin(u)-omegat, u, omegat, prec=16))
In [9]:
# たとえば e=0.5, omegat = 1.0 の値
Unum(0.5, 1.0)
Out[9]:
$\displaystyle 1.498701133517848$

数値解と近似解との比較

In [10]:
# 数値解 Unum と逐次近似解 Uite, フーリエ級数解 Ubes の差
e = 0.5
N = 10

print(" omega t    Unum               Uite-Unum          Ubes-Unum ")
for i in range(1, 10):
    omti = float(pi/10*i)
    print('%8.5f %18.15f %18.15f %18.15f' % 
          (omti, 
           Unum(e, omti), 
           Uite(N, e, omti) - Unum(e, omti), 
           Ubes(N, e, omti) - Unum(e, omti)
          )
         )

 omega t    Unum               Uite-Unum          Ubes-Unum 
 0.31416  0.593999023813608 -0.000048370871788  0.000205334548064
 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.000043154412395
 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 [11]:
def RMSite(N, e):
    tmpsum = 0
    for i in range(1, 10):
        omti = float(pi/10*i)
        tmpsum = tmpsum + (Uite(N, e, omti)-Unum(e, omti))**2
    return(sqrt(1/9 * tmpsum))

def RMSbes(N, e):
    tmpsum = 0
    for i in range(1, 10):
        omti = float(pi/10*i)
        tmpsum = tmpsum + (float(Ubes(N, e, omti))-Unum(e, omti))**2
    return(sqrt(1/9 * tmpsum))
In [12]:
print('%18.15e' % RMSite(10,0.5))
print('%18.15e' % RMSbes(10,0.5))

4.148349447036468e-05
1.462591053868229e-04

グラフ

$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 [13]:
def x(a, e, omegat):
    return(a*(cos(Uite(10, e, omegat)) - e))

def y(a, e, omegat):
    return(a*sqrt(1-e**2)*sin(Uite(10, e, omegat)))
In [14]:
# 縦横比。
import matplotlib.pyplot as plt
# デフォルトでは plt.rcParams['figure.figsize'] = (6.0, 4.0) だった
plt.rcParams['figure.figsize'] = (5, 5)

plot_parametric(
    x(5, 0.6, t), y(5, 0.6, t), (t, 0, 2*pi), 
    xlim=(-9, 3), ylim=(-6, 6));

In [15]:
omt = []
for i in range(1,37):
    omt.append(float(pi/18*i))

X = []
for omti in omt:
    X.append(x(5, 0.6, omti))

Y = []
for omti in omt:
    Y.append(y(5, 0.6, omti))
In [16]:
plt.xlim(-9, 3)
plt.ylim(-6, 6)
plt.plot(X, Y, 'bo');