H..1.5 Kepler の方程式を数値計算で解く

Kepler の方程式は文献に現れた初めての超越方程式であり、 これまでに何百という解法が考えられてきたそうである (平凡社百科事典の記述 (by 堀 源一郎) による)。

私の記憶が確かならば、 中野主一氏の天文計算プログラムでも、 反復法で解いていた。 今度、チェックしよう。

$M$ (平均近点角), $e$ (離心率) が与えられているとき、

$$E=M+e \sin E$$

を満す $E$ (離心近点角) を求めよ。

反復法で解く場合、初期値が必要になるが、 よほど離心率が大きく ($1$ に近く) ない限り、 $M$ を初期値に選ぶので良いと思われる。

$M=0.8$, $e=0.2$ とすると、 $E=0.9643338876952228\cdots$. $E=0.964333887695222658$b

Mathematicaで

EE=2/10
m=8/10
FindRoot[x==m+EE*Sin[x],{x,1},WorkingPrecision->100]
0.96433388769522264457506271899845052028615328006880 (改行)
  70848686011389985763852645291991226699739659693633

[

l]Cで

/*
 * solve-kepler.c
 */

#include <stdio.h>
#include <math.h>

int verbose = 1;

int maxcount = 0;
double M, e;

double f(double E)
{
  return E - e * sin(E) - M;
}

double df(double E)
{
  return 1.0 - e * cos(E);
}

typedef double func(double);

#define MAX_ITER 100

// Newton法で f(x)=0 の x0 の近くの解を求める
double newton(func f, double x0, double eps)
{
  int i;
  double x, dx;
  x = x0;
  for (i = 0; i <= MAX_ITER; i++) {
    dx = f(x) / df(x);
    x -= dx;
    if (fabs(dx) < eps)
      break;
  }
  if (i > maxcount)
    maxcount = i;
  /* 簡単に報告 */
  if (fabs(dx) > eps)
    printf("%d回で収束しなかった。dx=%g\n", MAX_ITER, dx);
  else if (verbose)
    printf("x0=%g, count=%d, dx=%g\n", x0, i, dx);
  /* 最後にもう1回反復 */
  x -= f(x) / df(x);
  return x;
}

int main(void)
{
  double E;
  int i, n;
  double pi, dM;
  pi = 4.0 * atan(1.0);
  /* 区間 [0,2π] を 100 等分 */
  n = 100;
  dM =2 * pi / n;
  /* 離心率 */
  printf("e="); scanf("%lf", &e);
  for (i = 0; i <= n; i++) {
    M = i * dM;
    /* M = E - e * sin(E) を解く */
    /* f(E) = E - e * sin(E) - M */
    E = newton(f, M, 1e-14);
    printf("%f %f\n", M, E);
  }
  printf("max count=%d\n", maxcount);

  M = 0.8; e = 0.2;
  E = newton(f, M, 1e-12);
  printf("e=%g, M=%g, E=%17.15f\n", e, M, E);
  return 0;
}