Adjusting Keplerian orbits for thrust with numerical stability

Question

I'm writing a mod for a game that models orbital physics (Kerbal Space Program, or KSP). I'm attempting to model the effects of thrust on spacecraft in certain states where the game only models them as a Keplerian orbit (The ship is 'on rails'), and I'm having some difficulty with numerical stability - even if I don't have any thrust applied I end up constantly increasing eccentricity for the orbit.

In the 'on rails' state, a spacecraft in the game follows a Keplerian orbit exactly. The orbit is calculated when the spacecraft is put on rails, and then when it goes off rails the game calculates the expected position and velocity at the current time and sets up the spacecraft to do normal Newtonian things from there.

At the moment, all my code does is take the calculated position and velocity of the spacecraft at the current time (using the game's own code) and then call a function already present in the game to calculate a Keplerian orbit for an object with that position and velocity around the currently-dominant body - that is, it should be a no-op. I'm actually seeing the orbit's eccentricity drift rapidly, certainly quickly enough that this won't work for gameplay purposes. This code is being called every physics timestep in the game, which is rather often - something like 50 times a second - so I suspect what I'm seeing is that the code in the game for calculating position and velocity from orbit and vice-versa isn't stable. I've verified that in actual gameplay taking a vessel off rails and then putting it back on rails can cause observable orbital differences, and I've confirmed that the calculated position/velocity of the craft change after I recalculate the orbit. Relative error in position is about $10^{-8}$, relative error in velocity is about $10^{-2}$.

I'm not sure what the best solution is here, and I only have very limited experience with numerical computation. I'd like to know what a good approach to this problem might look like. As far as I can see there are essentially three classes of solution:

• Calculate a new orbit directly from the old orbital parameters and the thrust. I don't know how plausible calculating a delta on a Keplerian orbit is, though, and I can't find much online about doing that.
• Write code with better precision or better numerical stability for doing orbit -> pos/vee -> orbit conversion. I don't know how plausible any stability is when you're doing that fifty times a second, though. I've found some code for calculating that, but I don't know where it sits on the fast/stable/good spectrum.
• Keep my own position/velocity values for the ship, update them continually, recalculate the orbit at intervals just so the player gets some feedback. That would require implementing a general two-body problem solution, which sounds like a lot of work.

What's a good option here?

Answer 1

Relative error of $10^{-8}$ for what should be an identity operation seems okay (but not great), but $10^{-2}$ is absolutely terrible. With a relative error per time step of $\epsilon$, the relative error after time $t$ at $50$ time steps per second will be $$ (1+\epsilon)^{50t} \approx e^{50\epsilon t}, $$ so the time-scale for accumulation of errors is $0.02/\epsilon$. For $\epsilon=10^{-2}$, this is $2$ seconds. From your drift in eccentricity it sound like that relative error is not only large, but also biased in one direction, and that bias is what would cause the drift. It's not completely clear whether this is specifically numerical stability, or if this is just an inaccurate computation method being used.

0. Implementing your own functions to compute position-velocty $\Leftrightarrow$ orbit parameters transformations is certainly feasible. A good textbook on orbital dynamics would explain the necessary formulas; there is also this helpful site. A book such as Accuracy and Stability of Numerical Algorithms by Higham will explain what numerical stability is and how to design a stable algorithm for these computations.

The main difficulty is covering all the corner cases correctly, since the equations are a little more complicated for inclined orbits. If I remember correctly it's all just coordinate transformations, with just one non-linear Kepler equation that is quite easy to solve.

Also note that adding a small amount of acceleration at the beginning of every time step may or may not be accurate enough, this would need to be experimentally checked. This is only first-order accurate (or second-order if the thrust vanishes at end-points of integration), so if it's not enough, then writing your own integrator is the way to go.

1. Implementing your own ODE integrator is quite feasible, and easier than it sounds. The integrator needs to be symplectic, so something like the Stormer-Verlet method is a good start: $$ \dot x = v, \qquad \dot v = f(x), $$ $$ v_{\frac12} = v_0 + \tfrac{h}{2} f(x_0), \qquad x_1 = x_0 + h v_{\frac12}, \qquad v_1 = v_{\frac12} + \tfrac{h}{2} f(x_1). $$ However, you probably need much more accuracy that just one Stormet-Verlet step per in-game time-step, so it's probably worth looking at a higher-order method.

The relevant references are Geometric Numerical Integration by Hairer, Lubich, Wanner, and Solving Ordinary Differential Equations by Hairer, Norsett, Wanner. According to that, implicit Gauss methods (Butcher-Kuntzmann methods) are symplectic, so for example if you don't mind solving a system of 6 nonlinear equations, then the 6th-order Butcher-Kuntzmann method might be worth trying; it's $\nu=3$ on p.56 of Butcher Implicit Runge-Kutta Processes (1964).