名著「天体と軌道の力学」（品切・重版未定で Amazon ではかなりの高値）の「2.7 ケプラー方程式の解法」には「ケプラー方程式の解法には実にさまざまな方法がある」と書いている。ここでは，Maxima-Jupyter で数値的に解く方法（別件で既に紹介済み）と，逐次近似法による極めて素朴な近似解法について，Maxima における関数の再帰的定義の練習問題もかねてメモしておく。
ケプラー方程式
ケプラー方程式
$$ u – e \sin u = \frac{2\pi}{T} t \equiv \omega t$$
について，$u$ を $t$ の陽関数として近似的に表すという話。
Maxima における関数の再帰的定義の練習問題として。
find_root() 関数による数値的解法
ケプラー方程式は超越方程式であるため，$u$ について解いて $t$ の陽関数として表すことはできない。そこで，Maxima の find_root() 関数を使って数値的に解く。
周期 $T$ を $N$ 等分し，
$$t_i = \frac{T}{N} \times i, \quad (i = 0, 1, \dots, N)$$
に対して，$u_i = u(t_i)$ を数値的に求める。
/* 離心率 */
e: 6/10$
/* 分割数 */
N: 36$
$$g(u_i) \equiv u_i – e \sin u_i = \frac{2\pi}{T} t_i = \frac{2\pi}{N} \times i$$を $u_i$ について 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)$
find_root() 関数で求めた u[i] をあらためて g(u) に代入して，誤差
$$\left|g(u_i) – \frac{2\pi}{N} \times i\right|$$
を表示すると，以下のように誤差は $10^{-15}$ より小さいことがわかる。
err: makelist(abs(g(u[i]) - float(2*%pi/N * i)), i, 0, N)$
lmax(err);
逐次近似法による素朴な近似的解法
離心率 $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} \\
\end{eqnarray}
$n$ が大きくなると，入れ子になっている項がどんどん増殖していきますが，$u_3$ のあたりまでは，近似的に $u$ は $t$ の陽関数としてあらわされているなぁ… という見た目がします。
上の式にそって，以下のように関数 $U(n, e, \omega t)$ を再帰的に定義します。
U(n, e, omegat):=
block(
if n = 0 then
omegat
else
omegat + e*sin(U(n-1, e, omegat))
)$
この定義から推測されるように，$U(n, e, \omega t)$ は $e^n$ 程度の精度と考えられるので，今回の例のように $e = 0.6$ だと，例えば誤差を $10^{-15}$ より小さくしようなどと考えると，$n$ の値をかなり大きくしなければなりません。
e;
float(e**68);
以下では例として $n = 70$ として逐次近似法による近似解 $U(70, e, \omega t)$ と数値解 u[i] との誤差の最大値を表示しています。
err: makelist(abs(u[i] - U(70, e, float(2*%pi/N * i))), i, 1, N)$
lmax(err);