Iterating the hyperbolic version of Kepler's equation does not converge

93 views Asked by At

I'm creating an orbit simulator using Java. I have elliptical orbits working great, and now I'm trying to do hyperbolic orbits. The problem I'm having is that when the ship enters the orbit at periapsis, it moves very slowly at first, and quickly accelerates until it reaches infinite velocity at infinite altitude. Obviously, it should start at its fastest and go to zero, which is the opposite of what it's doing. The orbit path is rendering properly and the ship is following it, but the problem seems to be that Kepler's equation is growing exponentially instead of converging on a value of the hyperbolic anomaly, and I don't know why. I'm following my references precisely, and the equivalent code for e < 1 works perfectly. Here is the function that calculates the ship's position at time t, with n being the mean motion.

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 * sinh(H)) - M;}
    nu = 2 * atan((tanh(H / 2)) * sqrt((e + 1) / (e - 1))); //current true anomaly
    //position variables
    r = (a * (1.0 - (e * e))) / (1.0 + (e * cos(nu))); //current orbital radius
    phi = atan((e * sin(nu)) / (1.0 + (e * cos(nu)))); //current flight angle
    v = sqrt(GM * ((2.0 / r) - (1.0 / a))); //current orbital velocity
}

I've quadruple checked every equation, verified the values of the variables, and gone over each step with a fine toothed comb. I think the problem is that H = (e * sinh(H)) - M is growing exponentially with every iteration instead of converging, but I don't think that's what it's supposed to do. I wrote this simple test program to watch it work, and this is what it gave me.

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();
            }
        }
    }
}

Output:

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

I'm trying to teach this to myself, so I'm hoping it's something obvious like using the wrong equation or skipping a step.

1

There are 1 answers

1
MBobrik On BEST ANSWER

Trying to find a solution of X = f(X) by just iterating it as an assignment works only if d(f)/dX < 1 ... just use a different numerical method, like for example https://en.wikipedia.org/wiki/Newton%27s_method