Newton-Raphson method on the hyperbolic Kepler's equation

Question

I'm here once again because I can't figure this out. I'm building an orbit simulator and currently working on placing the ship on a hyperbolic trajectory upon entering the SoI of a body. (I'm using patched conics in 2D.) I'm having an issue with the math though, where the orbit is calculating all the correct parameters, but the ship is ending up in the wrong spot.

I tracked the problem down to the point where the current hyperbolic anomaly (F) is calculated from the current mean anomaly (M). Based on a previous question I asked, I'm using the Newton-Raphson method (based on this) to find F:

for(int i = 0; i < 10; i++) {F = F - (((e * sinh(F)) - F - M) / ((e * cosh(F)) -1));}

The problem is that I'm not getting a symmetrical result from M -> F as from F -> M. A ship that has

nu0 = -0.5346949277282228
F0 = 0.04402263120230271
M0 = 5.793100753021599E-4
gets
M = 5.793100753021599E-4 (good)
F = 0.01027520200339216 (wrong)
nu = 0.1276522417546593 (also wrong)

It ends up at the wrong point on the orbit and nothing else is right.

To try to narrow down the problem, I graphed the equations to visualize what they were doing. In this graph, I have both the hyperbolic equation and my method of solving it. The first graph, M from E, does what I would expect: smoothly curves from (-pi, -infinity) to (+pi, +infinity). That's what the proper shape of a hyperbolic orbit is. I would expect the Newton method to give a perfect inverse, going from (-infinity, -pi) to (+pi, +infinity). But that's not what it does, it has a couple weird humps near 0,0, but otherwise goes from (-infinity, -infinity) to (+infinity, +infinity). The asymptote is sloped, not horizontal.

I also did the same thing with the elliptical case, and it produced a perfect inverse equation, exactly as I expected. But the hyperbolic equivalent does not, and I am completely lost as to why. I've tried numerous different forms of the equation, but they all give this same shape. I've tried different starting guesses and parameters, but none give me the mirrored function I want.

Am I doing something wrong? Is this actually right and I'm just misleading myself? I've taken calculus but not enough to diagnose this myself. Hopefully it's something simple that I'm doing wrong.

Answer 1

The Newton-Raphson solver you wrote for Kepler's equation is correct. Still, I recommend using a fixed-point verification procedure as a test of convergence. This is a very common pattern in the implementation of approximate methods.

The idea is that we run the loop until the solution doesn't change anymore. That means we found a fixed point. The iterated function applied to the current value produces the same value.

Note that a fixed point might not necessarily be the solution - it requires a mathematical proof of about the relation between the two, but: a) it is so in a very large set of cases; b) it's the best we can do, and c) we can verify the equation and make sure it is indeed a solution what we found.

Also, it is impractical and dangerous to check if the value is exactly equal to the next value, as in general comparison between floating point numbers should not be done by x == y but abs(x-y) < eps where eps is a very small number. The value of eps is related to the number of decimals we want the solution to be computed precisely: if eps = 5e-7 it means that x and y are equal up to 6 decimals, if eps = 5e-11, they agree up to 10 decimals, etc.

We also need to set an NMAX - the maximum number of iterations after which we give up. It can be a large number, as in practice if the equation is well-behaved (as this is), the number of iterations never gets close to that value.

With these, your code for the Newton-Raphson method applied to Kepler's equation in the hyperbolic case can be written in java:

public static double solveKepler(double e, double M, 
        double epsFixedPoint, int NMAX, double epsEquation) 
        throws SolverFailedExeption{
    double FNext = M,
        F = FNext + 1000 * epsFixedPoint, // make sure the initially F-FNext > eps
        i = 0;
    while(abs(FNext - F) > epsFixedPoint && i < NMAX){
        F = FNext;
        FNext = F - (e * sinh(F) - F - M) / (e * cosh(F) -1);
        i++;
    }
    if(abs(M + F - e *sinh(F)) > epsEquation){ // verify the equation, not the fixed point
        throw new SolverFailedExeption("Newton-Raphson method for Kepler equation failed");
    }
    return FNext;
}

public static double solveKepler(double e, double M) throws SolverFailedExeption{
    return Main.solveKepler(e, M, 5e-7, 10000, 5e-7);
}

Some might prefer the equivalent implementation of the loop by a for and a break.

Here's the link for the full code of this example.