一种实用的开普勒方程求解方法及其 C 语言实现

本文介绍一种实用的开普勒方程求解方法，由于其算法简单、高效且利于编码实现，因而具有较好的实用性。

开普勒方程是航天动力学基础方程，是开普勒定律的数学描述。由于开普勒方程属于超越方程，因而无法通过解析法精确求解，这一问题在历史上困扰全世界数学家们达 200 年之久，直至牛顿迭代方法的出现。

本文将介绍一种实用的开普勒方程求解方法，并采用 C 语言实现其算法。该方法出自美国海军天文台 Marc A. Murison 名为 A Practical Method for Solving the Kepler Equation 的文章，由于其算法简单、高效且利于编码实现，因而具有较好的实用性。

1. 开普勒方程描述

开普勒方程可表达为：

$f(E(t))=E(t)-esinE(t)-M(t)=0$

它描述了线性时间 $t$ 下偏近角 $E(t)$ 和平近角 $M(t)$ 之间的非线性转换关系。式中， $e$ 为轨道偏心率， $M(t) = n(t-\tau)$ 为平近角。其中， $n = \sqrt{\mu / a^3}$ 表示二体模型中的平运动， $\tau$ 为天体（或卫星）飞越近地点后累积的时长， $a$ 为轨道半长轴。

实际计算中，有两种情况需要求解开普勒方程，一是已知真近角 $\theta$ ，求平近角 $M$ ；另外则是已知平近角 $M$ ，求真近角 $\theta$ 。可见偏近角 $E$ 只是过程量，真近角 $\theta$ 和 $M$ 才是最终想要的结果量。在航天动力学中偏近角 $E$ 本就是不真实的虚拟量，但它与真近角 $\theta$ 存在直观的几何关系，如图 1 所示。当然平近角 $M$ 也是不真实的虚拟量，它必须通过偏近角和开普勒方程与真近角实现数值上的转换。

真近角和偏近角之间的几何关系

图 1: 真近角和偏近角之间的几何关系

由图 1 内部椭圆（真实运动轨迹）几何关系可得：

$r = a(1 - EcosE)$

由图 1 外部圆（假想的平运动）几何关系可得：

$acosE = rcos \theta + ae$

于是可得真近角和偏近角之间的相互转换关系为：

$tan \theta = \frac{sin \theta}{cos \theta} = \frac{\sqrt{1-e^2} sinE}{cosE-e} \\ tanE = \frac{sinE}{cosE} = \frac{\sqrt{1-e^2} sin \theta}{e+cos \theta}$

开普勒方程属于超越方程，故需采用数值法才能求得精确值。理想情况下，一个数值方法是否实用，应该看它是否同时具备以下两点：

计算 $E$ 的 CPU 耗时是否最小化;
程序的复杂度是否最小化。

这两个需求往往是互相制约的，但幸运的是通过分析研究可以找到一个折中的方法来解决这个矛盾。开普勒方程只能通过迭代法求解，任何一个迭代过程的设计都有两个任务：

第一个任务是设计迭代循环中的逼近算法。这个逼近算法需要重复执行直到结果达到令人满意的精度，一般说来迭代公式所用的阶数越高，迭代所需的次数就越少。然而，更高阶的迭代公式将使得公式的表达式变得十分复杂，这将极大地耗费 CPU 的运行时间。因此，无论我们选择何种算法，一个恰当（通常是比较低）的阶数会使得 CPU 的耗时最短；
第二个任务是选择一个迭代循环的初始值。初始值选的越精确，循环收敛的越快。选择初始值的方法不需要和迭代方法一样，哪怕两者极其相近。类似于迭代逼近算法，对于一个初值确定的算法，将有一个理想的阶数能够最大限度地减少 CPU 的耗时。

接下来将给出一个非常简单的初值确定方法，紧接着是一个快速迭代算法。最后，直接给出算法性能分析结果及其伪代码，并用 C 语言加以实现之。

2. 初值确定方法

由于必须用迭代法求解开普勒方程，因而循环迭代的初始值越精确，迭代效果就越好，至少在迭代表达式复杂得令人反感之前应该是如此。将开普勒方程写成：

$E = M + esinE$

在 $e=0$ 的极限情况下，有 $E=M$ ，这是最简单的初始值近似。因此上面的方程可改写成如下简单的迭代近似公式：

$E_k = M + esinE_{k-1}$

其中初始条件 $E_0 = M$ 。我们可以对上述递归表达式进行反复迭代，直至得到足够高的偏心率阶数。例如，三阶近似公式为：

$E = M + esinM + e^2sinMcosM + \frac{1}{2}\ e^3sinM(3cos^2M-1)$

程序化后的三阶伪代码如下：

KeplerStart3 := proc(e,M)
    local t33, t35, t34;
    t34 := eˆ2;
    t35 := e*t34;
    t33 := cos(M);
    return M+(-1/2*t35+e+(t34+3/2*t33*t35)*t33)*sin(M);
end proc;

3. 循环迭代方法

循环迭代的最终目的是要得到收敛的结果。具体到开普勒方程，当给定一个带有误差的 $E$ 时，必须找到一个迭代公式，使其每次返回一个误差更小的近似值，同时它也必须收敛。基于此，按如下方式改写开普勒方程：

$f(x) = x - esinx - M$

式中， $f(x)=0$ 的解是 $x=E$ 。令 $\epsilon \equiv x-E$ 为 $x$ 逼近 $E$ 时存在的误差，将 $f(x)$ 在 $x=E$ 处进行泰勒展开，于是得到：

$f(E) = f(x-\epsilon) = x-esinx-M-(1-ecosx) \epsilon + \frac{1}{2}\ \epsilon ^2 esinx - \frac{1}{6}\ \epsilon ^3 ecosx + ...$

式中，假设 $\epsilon$ 充分小。若只考虑一阶展开，可得：

$\epsilon = \frac{x-esinx-M}{1-ecosx}$

上式可作为一阶迭代方案的的核心算法。假设一开始猜测 $x = x_0$ ，那么 $x_1 = x_0 + \epsilon$ 比 $x_0$ 更接近于 $E$ 。于是得到一阶迭代方程：

$\epsilon_{n+1} = \frac{x_n-esinx_n-M}{1-ecosx_n}$

其中初始值 $x_0$ 的取值将在后面讨论。由上式可以得到单步一阶迭代法来估计 $E_{n+1} = E_n - \epsilon_n$ ，对应的伪代码如下：

1
2
3

eps1 := proc(e,M,x)
    return (x-e*sin(x)-M)/(1-e*cos(x));
end proc;

现在把二阶项考虑进去，方程展开后写成 Horner 形式：

$f(x-\epsilon) = x-esinx-M- \left( 1 - ecosx - \frac{1}{2}\ esinx \cdot \epsilon \right) \epsilon$

同理，令 $f(x-\epsilon)=0$ ，整理得到如下形式：

$\epsilon = \frac{x-esinx-M}{1-ecosx-\frac{1}{2}\ \epsilon esinx}$

类比单步一阶形式，创建单步二阶迭代方程：

$\epsilon_{n+1} = \frac{x_n-esinx_n-M}{1-ecosx_n-\frac{1}{2}\ \epsilon_n esinx_n}$

重点来了，此处可以将 $\epsilon_n$ 用一阶迭代公式替换，创建一个两步迭代过程，于是得到新的迭代方程：

$\epsilon_{n+1} = \frac{x_n-esinx_n-M}{1-ecosx_n-\frac{1}{2}\ esinx_n \frac{x_n-esinx_n-M}{1-ecosx_n}}$

上述方程经优化后的伪代码为：

eps2 := proc(e,M,x)
    localt1,t2,t3;
    t1 := -1+e*cos(x);
    t2 := e*sin(x);
    t3 := -x+t2+M;
    return t3/(1/2*t3*t2/t1+t1);
end proc;

更进一步， $f(E)=f(x-\epsilon)$ 三阶展开方程的 Horner 形式为：

$f(x-\epsilon) = x-esinx-M- \left( 1 - ecosx - \left( \frac{1}{2}\ esinx - \frac{1}{6}\ ecosx \cdot \epsilon \right) \epsilon \right) \epsilon$

同理，令 $f(x-\epsilon)=0$ ，整理得到如下形式：

$\epsilon_{n+1} = \frac{x_n-esinx_n-M}{1-ecosx_n-\frac{1}{2}\ \left( esinx_n - \frac{1}{3}\ ecosx \cdot \epsilon_n \right) \epsilon_n}$

重点又来了，此处可以先将 $\epsilon_n$ 用二阶迭代公式替换，再将 $\epsilon_n$ 用一阶迭代公式替换，创建一个三步迭代过程，于是得到新的迭代方程，由于方程过于复杂，直接给出优化后的伪代码：

eps3 := proc(e,M,x)
    local t1, t2, t3, t4, t5, t6;
    t1 := cos(x);
    t2 := -1+e*t1;
    t3 := sin(x);
    t4 := e*t3;
    t5 := -x+t4+M;
    t6 := t5/(1/2*t5*t4/t2+t2);
    return t5/((1/2*t3 - 1/6*t1*t6)*e*t6+t2);
end proc;

可以继续利用这种方式到更高阶的形式，本文的推导止于此，后文我们会看到三阶迭代已经处于 CPU 时耗最优区间。

4. 实用的开普勒方程求解方法

数值测试，令算法执行的 CPU 耗时为迭代阶数 (Niter) 和初始值逼近阶数 (Nstart) 的二元函数，绘制求解 16000 个开普勒方程的 CPU 耗时等高线，其中 $e$ 和 $M$ 在 $\{ \mathbb{R} \times \mathbb{R}: e \in [0,1), M \in [0, \pi] \}$ 等间隔 $400 \times 400$ 网格域中选取。结果表明，不论是对于初值算法，还是迭代算法，三阶形式都接近最优，如图 2 所示。

真近角和偏近角之间的几何关系

图 2: 真近角和偏近角之间的几何关系

由此，给出优化后的三阶迭代和初始值方法求解开普勒方程的伪代码：

KeplerSolve := proc( e, M, tol:=1.0e-14 )
    local dE, E, E0, Mnorm, count;
    global Estart3, eps3;
    Mnorm := fmod(M,2*Pi);
    E0 := KeplerStart3(e,Mnorm);
    dE := tol + 1;
    count := 0;
    while dE > tol do
        E := E0 - eps3(e,Mnorm,E0);
        dE := abs(E-E0);
        E0 := E;
        count := count + 1;
        if count=100 then
            print “太令人惊讶了，KeplerSolve 竟然不收敛！”;
            break;
        end if;
    end do;
    return E;
end proc;

5. C 语言实现

根据优化后的实用开普勒方程求解伪代码，可编写其计算机语言实现代码。下面给出三个主要函数的 C 语言实现，这些代码已经过充分测试，可放心使用。函数中的 fmod 函数为取模函数，用于处理平近角周期问题，请读者自行补脑。至此，码完收工，拜了个拜！

/* ---------------------------------------------------------------------
**
** 函数名称: KeplerStart3
** 函数功能: 开普勒方程三阶初值估计
** 输入参数: ecc -> 偏心率 [N/A]
**           M   -> 平近角 [rad]
** 输出参数: 无
** 返回值  : E   -> 偏近角 [rad]
**
** ---------------------------------------------------------------------*/
double KeplerStart3(double ecc, double M){

    double t34, t35, t33, E;

    if(ecc < 0 || ecc >= 1){
        printf("!!!错误的偏心率，该程序只适合椭圆轨道(0 =< ecc < 1)!!!\n");
        return -1;
    }

    if(M >= D2PI || M < 0){
        M = fmod(M, 0, D2PI);
    }

    t34 = ecc * ecc;
    t35 = ecc * t34;
    t33 = cos(M);

    E = M + (-0.5*t35 + ecc + (t34 + 3.0/2.0*t33*t35) * t33) * sin(M);

    return E;
}

/* ---------------------------------------------------------------------
**
** 函数名称: eps3
** 函数功能: 开普勒方程三阶误差估计
** 输入参数: ecc -> 偏心率             [N/A]
**           M   -> 平近角             [rad]
**           x   -> 偏近角过程值       [rad]
** 输出参数: 无
** 返回值  : E_error -> 偏近角估计误差 [rad]
**
** ---------------------------------------------------------------------*/
double eps3(double ecc, double M, double x){

    double t1, t2, t3, t4, t5, t6, E_error;

    if(ecc < 0 || ecc >= 1){
        printf("!!!错误的偏心率，该程序只适合椭圆轨道(0 =< ecc < 1)!!!\n");
        return -1;
    }

    if(M >= D2PI || M < 0){
        M = fmod(M, 0, D2PI);
    }

    t1 = cos(x);
    t2 = -1 + ecc * t1;
    t3 = sin(x);
    t4 = ecc * t3;
    t5 = -x + t4 + M;
    t6 = t5 / (0.5 * t5 * t4 / t2 + t2);

    E_error = t5 / ((0.5*t3 - 1.0/6.0*t1*t6) * ecc * t6 + t2);

    return E_error;
}

/* ---------------------------------------------------------------------
**
** 函数名称: KeplerSolve
** 函数功能: 求解开普勒方程
** 输入参数: ecc -> 偏心率 [N/A]
**           M   -> 平近角 [rad]
** 输出参数: 无
** 返回值  : E   -> 偏近角 [rad]
**
** ---------------------------------------------------------------------*/
double KeplerSolve(double ecc, double M){

    double E0, dE, E, Mnorm;
    double tol = 1.0e-20;
    int    iter_count = 0;

    if(ecc < 0 || ecc >= 1){
        printf("!!!错误的偏心率，该程序只适合椭圆轨道(0 =< ecc < 1)!!!\n");
        return -1;
    }

    if(M >= D2PI || M < 0){
        M = fmod(M, 0, D2PI);
    }

    Mnorm = fmod(M, 0, D2PI);

    E0 = KeplerStart3(ecc, Mnorm);
    dE = tol + 1;

    while(dE > tol){
        E  = E0 - eps3(ecc, Mnorm, E0);
        dE = fabs(E - E0);
        E0 = E;

        iter_count = iter_count + 1;

        if(iter_count > 10){
            printf("太令人惊讶了，KeplerSolve 竟然不收敛！\n");
            return -1;
        }

        // printf("迭代次数：iter_count = %d\n", iter_count);
    }

    return E;
}