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

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

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

ケプラー方程式

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

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

In [1]:
# 最近の gnuplot にはベッセル関数 Jn(x) もある
help besjn
Out[1]:
 関数 `besjn(n,x)` は引数の Jn ベッセル関数 (n 次の第 1 種円柱関数 Jn、
 n 次ベッセル関数) の値を返します。引数 x はラジアンで与えます。

以下のようにして，三項演算子を使って再帰的に定義します。
\begin{eqnarray}
u(0, e, \omega t) &=& \omega t \\
u(N, e, \omega t) &=& \omega t+ \sum_{n=1}^{N} \frac{2}{n} J_n(n e) \sin(n\, \omega t) \\
&=& u(N-1, e, \omega t) + \frac{2}{N} J_N(N e) \sin(N\, \omega t)
\end{eqnarray}

In [2]:
# 三項演算子を使った再帰定義による sum2 の定義例
# i**2 を i=1 から N まで足す例
# sum2(N) = (N==1 ? 1**2 : sum2(N-1) + N**2)

Ubes(N, e, omegat) = (N == 0 ? \
   omegat \
   : Ubes(N-1, e, omegat) + 2./N * besjn(N, N*e)*sin(N*omegat))
In [3]:
# たとえば N=10, e=0.5, omegat = 1.0 の値 
print Ubes(10, 0.5, 1)
Out[3]:
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)$ を三項演算子を使って再帰的に定義します。

\begin{eqnarray}
N = 0 \ \ &\Rightarrow \ \ u_0 =& \omega t\\
N > 0 \ \ &\Rightarrow \ \ u_{N} =& \omega t + e \sin u_{N-1} \\
\end{eqnarray}

In [4]:
# gnuplot の三項演算子を使った関数の定義例 
# N==0 なら f(N)=A, そうでないなら f(N)=B 
# f(N) = (N==0 ? A : B) 

Uite(N, e, omegat) = (N == 0 ? \
  omegat \
  : omegat + e*sin(Uite(N-1, e, omegat)))
In [5]:
# たとえば N=10, e=0.5, omegat = 1.0 の値 
print Uite(10, 0.5, 1)
Out[5]:
1.49870113351784

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

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

In [6]:
f(e, omt, u) = u - e*sin(u) - omt

# 再帰定義による解
# x には探索の初期値を入れる
# h は解の収束条件である微小量と数値微分をするときに使う微小量を兼ねる。
newton_f(e, omt, x, h) = \
  (abs(f(e, omt, x)) > h ? \
   newton_f(e,omt,x-f(e,omt,x)*2.*h/(f(e,omt,x+h)-f(e,omt,x-h)),h) \
   : x)

Unum(e, omegat) = newton_f(e, omegat, omegat, 1E-15)
In [7]:
# たとえば e=0.5, omegat = 1.0 の値
print Unum(0.5, 1)
Out[7]:
1.49870113351785

数値解と近似解との比較

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

do for [i=1: 9] {
  omti = pi/10*i
  print sprintf("%7.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))
}
Out[8]:
omega t    Unum               Uite-Unum          Ubes-Unum
0.31416  0.593999023813606 -0.000048370871787  0.000205334548065
0.62832  1.065940683889791 -0.000000479808531 -0.000233511559215
0.94248  1.438080909968085 -0.000000000003061  0.000199446542732
1.25664  1.748741781633489  0.000000000005287 -0.000159160439054
1.57080  2.020979938089771 -0.000000056633643  0.000123562038968
1.88496  2.268208852924499 -0.000003478720316 -0.000093184551056
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 [9]:
# RMSite
tmpsum = 0
do for [i=1: 9] {
  omti = pi/10*i
  tmpsum = tmpsum + (Uite(N, e, omti)-Unum(e, omti))**2
}
print sprintf("Uite: %18.15e", sqrt(1./9 * tmpsum))

# RMSbes
tmpsum = 0
do for [i=1: 9] {
  omti = pi/10*i
  tmpsum = tmpsum + (Ubes(N, e, omti)-Unum(e, omti))**2
}
print sprintf("Ubes: %18.15e", sqrt(1./9 * tmpsum))
Out[9]:
Uite: 4.148349447013503e-05
Ubes: 1.462591053870202e-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 [10]:
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))
In [11]:
%gnuplot inline svg size 600,480 fixed enhanced font 'Noto Sans CJK JP,14'
In [12]:
set parametric
set xrange [-9:3]
set yrange [-6:6]
set size square
set zeroaxis

plot [t=0:2*pi] x(5, 0.6, t), y(5, 0.6, t) lc "blue" notitle

Out[12]:
	dummy variable is t for curves, u/v for surfaces
In [13]:
set term svg font 'Noto Sans CJK JP,14'
set output "gn-fig1.svg"
replot
Out[13]:
Terminal type is now 'svg'
Options are 'size 600,480 fixed enhanced font 'Noto Sans CJK JP,14' butt dashlength 1.0 '
In [14]:
%gnuplot inline svg size 600,480 fixed enhanced font 'Noto Sans CJK JP,14'
In [15]:
array X[36]
array Y[36]

do for [i=1: 36] {
  omegati = pi/18 * i
  X[i] = x(5, 0.6, omegati)
  Y[i] = y(5, 0.6, omegati)
}
In [16]:
plot Y using (X[$1]):2 pt 7 ps 0.8 lc "blue" notitle

In [17]:
set term svg font 'Noto Sans CJK JP,14'
set output "gn-fig2.svg"
replot
Out[17]:
Terminal type is now 'svg'
Options are 'size 600,480 fixed enhanced font 'Noto Sans CJK JP,14' butt dashlength 1.0 '