Java轨道模拟器：双曲线开普勒方程迭代不收敛问题

我用Java开发轨道模拟器，椭圆轨道（e<1）运行正常，但双曲线轨道（e>1）出现两个核心问题：

• 飞船从近心点进入轨道后初始速度极慢，随后迅速加速至无穷，和预期的“从最快速度逐渐降低”完全相反
• 求解双曲线异常H的开普勒方程迭代过程发散，H值指数增长至无穷，导致后续速度计算错误

以下是当前计算飞船t时刻位置的核心代码：

if(e > 1) {
    M = n * t; //current mean anomaly
    //iterate Kepler's equation to converge on a value of the current hyperbolic anomaly
    for(int i = 0; i < KEPLER; i++) {H = (e * Math.sinh(H)) - M;}
    nu = 2 * Math.atan((Math.tanh(H / 2)) * Math.sqrt((e + 1) / (e - 1))); //current true anomaly
    //position variables
    r = (a * (1.0 - (e * e))) / (1.0 + (e * Math.cos(nu))); //current orbital radius
    phi = Math.atan((e * Math.sin(nu)) / (1.0 + (e * Math.cos(nu)))); //current flight angle
    v = Math.sqrt(GM * ((2.0 / r) - (1.0 / a))); //current orbital velocity
}

我排查后确定问题出在H = (e * Math.sinh(H)) - M这个迭代公式上。测试程序的输出显示H值迅速发散：

测试代码：

public class Main {
    public static void main(String[] args) throws java.io.IOException {
        double E = 0, H = 0, E0 = 0, H0 = 0;
        int k = 0, k1 = 11;
        double M = 0.04;
        double e1 = 0.1;
        double e2 = 2.0;
        while(k < k1) {
            E = M + (e1 * Math.sin(E));
            H = (e2 * Math.sinh(H)) - M;
            System.out.printf("k = %d,  E = %.6f,  H = %.6f\n", k, E, H);
            System.out.printf("       dE = %.6f, dH = %.6f\n", E - E0, H - H0);
            E0 = E;
            H0 = H;
            k++;
            try {
                Thread.sleep(100);
            } catch (InterruptedException e) {
                e.printStackTrace();
            }
        }
    }
}

输出结果：

k = 0,  E = 0.040000,  H = -0.040000
       dE = 0.040000, dH = -0.040000
k = 1,  E = 0.043999,  H = -0.120021
       dE = 0.003999, dH = -0.080021
k = 2,  E = 0.044398,  H = -0.280619
       dE = 0.000400, dH = -0.160598
k = 3,  E = 0.044438,  H = -0.608634
       dE = 0.000040, dH = -0.328014
k = 4,  E = 0.044442,  H = -1.333825
       dE = 0.000004, dH = -0.725191
k = 5,  E = 0.044443,  H = -3.572066
       dE = 0.000000, dH = -2.238241
k = 6,  E = 0.044443,  H = -35.601966
       dE = 0.000000, dH = -32.029899
k = 7,  E = 0.044443,  H = -2895592024320601.000000
       dE = 0.000000, dH = -2895592024320565.500000
k = 8,  E = 0.044443,  H = -Infinity
       dE = 0.000000, dH = -Infinity
k = 9,  E = 0.044443,  H = -Infinity
       dE = 0.000000, dH = NaN
k = 10,  E = 0.044443,  H = -Infinity
       dE = 0.000000, dH = NaN

问题根源与修正方案 #

1. 开普勒方程的形式错误 #

双曲线轨道的开普勒方程是：
$M = e \sinh H - H$
你当前的迭代式H = e*Math.sinh(H) - M是把方程错误变形为H = e sinhH - M，属于不动点迭代，但这种形式对e>1的情况完全不收敛：

• 不动点迭代收敛的条件是迭代函数的导数绝对值小于1，这里迭代函数$g(H)=e\sinh H - M$的导数是$g’(H)=e\cosh H$
• 因为$\cosh H \geq 1$且e>1，所以$|g’(H)| >1$，迭代必然发散

2. 正确的迭代方法：牛顿-拉夫逊法 #

必须使用牛顿-拉夫逊法求解双曲线开普勒方程，迭代式推导如下：
定义函数：$f(H) = e\sinh H - H - M$
求导得：$f’(H) = e\cosh H - 1$
牛顿迭代式：$H_{n+1} = H_n - \frac{e\sinh H_n - H_n - M}{e\cosh H_n - 1}$

3. 修正后的迭代代码 #

替换原来的for循环，使用牛顿法迭代：

if(e > 1) {
    M = n * t; //current mean anomaly
    // 牛顿法迭代求解双曲线异常H
    double H = 0.0; // 初始值，近心点时刻H=0
    for(int i = 0; i < KEPLER; i++) {
        double sinhH = Math.sinh(H);
        double coshH = Math.cosh(H);
        double f = e * sinhH - H - M;
        double fPrime = e * coshH - 1;
        H -= f / fPrime;
        // 可选：加入收敛判断，当|f|小于阈值时提前退出
        if(Math.abs(f) < 1e-10) break;
    }
    nu = 2 * Math.atan(Math.tanh(H / 2) * Math.sqrt((e + 1) / (e - 1))); //current true anomaly
    //position variables
    r = (a * (1.0 - (e * e))) / (1.0 + (e * Math.cos(nu))); //current orbital radius
    phi = Math.atan((e * Math.sin(nu)) / (1.0 + (e * Math.cos(nu)))); //current flight angle
    v = Math.sqrt(GM * ((2.0 / r) - (1.0 / a))); //current orbital velocity
}

4. 测试验证 #

用你的测试参数（M=0.04，e=2.0）运行牛顿法，H会快速收敛到正确值（约0.040011），不会发散，后续的速度计算也会符合预期：从近心点的最快速度逐渐降低（趋近于双曲线的逃逸速度）。

内容的提问来源于stack exchange，提问作者Joshua Harvey