0

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.

  • I'm not sure I understand what your difficulty is, but I don't get any issues. From your data I computed `e = (M0+F0)/sinh(F0) = 1.0128321944527299`. Then with `M = M0; F = M; ` and your `for` loop for Newton-Raphson (although it should be something on the lines of `while(abs(Fnext-F) > eps)... `), I get very fast (4 iterations) the solution `F = 0.04402263120230289` – kikon Mar 26 '23 at 23:38
  • That's odd, I have e = 4.719... maybe something else is off somewhere – Joshua Harvey Mar 27 '23 at 00:07
  • The equation I'm using for e is 4.30 from [here](http://braeunig.us/space/orbmech.htm#launch). Is this maybe only valid for elliptical cases? – Joshua Harvey Mar 27 '23 at 00:14
  • That's okay but then for `e=4.719` and `F0 = 0.04402263120230271` you get `M0 = 0.1637872725920228` which also works out great to get `F` by solving Keplers equation by your Newton-Raphson iteration – kikon Mar 27 '23 at 00:17
  • I don't think it only works for elliptic orbits, it's just [eq. 112](https://orbital-mechanics.space/the-orbit-equation/the-orbit-equation.html#equation-eq-vector-orbit-equation) for the launch geometry. But that doesn't matter for the problem of solving Kepler's equation as long as `e>1` you get an `M` from `F` directly and the same `M` from that `F` by applying Newton-Raphson – kikon Mar 27 '23 at 00:37
  • Hmm... it does seem that there are some errors at large values of nu, but you're right, it is working most of the time. But it helps narrow down the problem, it seems that there's something screwy with the way the orbits are being initialized. – Joshua Harvey Mar 27 '23 at 01:44
  • If Newton doesn't converge, you may increase the number of iterations; it would be a great advantage if you changed the fixed number of iterations to a while loop logic so you could detect if the method appeared to have converged or not. – kikon Mar 27 '23 at 01:51
  • That makes sense, but what would I test against for convergence? I get what you mean with `while(abs(Fnext - F) > ?)` but I'm not sure what the ? should be. – Joshua Harvey Mar 28 '23 at 00:09
  • The value of `eps` is the precision you want for the solution. I added a more detailed description of this as an answer since it doesn't fit here. Please let me know if some parts of it don't make sense to you. – kikon Mar 28 '23 at 09:46

1 Answers1

0

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.

kikon
  • 3,670
  • 3
  • 5
  • 20