Two Body Kepler Propogation: Errors from using float / single precision

Question

I am looking into leveraging GPU's for two body orbit propagation. I am doing basic Kepler propagation, and only dealing with circular and elliptical orbits. I took a first shot at using single precision (8x speed for GPU's) and noticed some significant error in the orbit.

I wanted to know if it is known to be a problem dealing with single precision, or should it be ok? I am only concerned with error being in the region of 1km or so.

---

UPDATE

The precision error I am talking about is not to do with integration, like was mentioned I am using the closed form analytical method (very similar to the linked method) and using Newton-Rhapson integration technique for elliptical orbits.

I ended up finding out that the precision issue i was having was when calculating GMST and using that for the conversion from ECI to ECEF. I even tried using MJD, instead of JD, but it turns out if I just pre-compute the GMST for each time step using double precision, that solves my issue.

Answer 1

Single precision floating point uses 23 bit for the mantissa. Distances in meters may be up to 8,388,608 m with a resolution of a single meter. If you add small increments with a mean error of 0.5 m a thousand times, the error may accumulate to about 0.5 km.

If you need distances of up to 33,554,432 m you get a resolution of 4 m. The mean error will be 2 m and a thousand increments may produce an error of 2 km.

Distances of about 400,000 km like the Moon to Earth may be represented with a resolution of 64 m using single precision.

A small python program for demonstration. Using numpy arrays to force the use of single precision float.

import numpy as np

a=np.zeros((1), dtype='float32')
b=np.zeros((1), dtype='float32')
c=np.zeros((1), dtype='float32')
d=np.zeros((1), dtype='float32')

a[0]=b[0]=c[0]=d[0]=8E6

for i in range(0,1000):
    a[0] += 1.0
    b[0] += 0.75
    c[0] += 0.5
    d[0] += 0.25

print(a[0], b[0], c[0], d[0])

result:

8.001e+06 8.001e+06 8.0005e+06 8e+06

but without rounding errors the values should be:

8.001e+06 8.00075e+06 8.0005e+06 8.00025e+06

Adding 1.0 we get the expected result, 0.75 is rounded to 1.0, 0.5 rounded to 0.5 but 0.25 is rounded to 0.0

But it is a little bit more complicated. Single bit floating point numbers are normalized, that means the most significant bit of the fractional part is always 1. Since it is 1 for all nonzero numbers, it is not stored. But resolution is increased by using 24 instead of 23 bits. So distances up to 16,777,216 m may be represented with a resolution of 1 m.

Answer 2

I am looking into leveraging GPU's for two body orbit propagation. I am doing basic Kepler propagation, and only dealing with circular and elliptical orbits.

Just fyi 2-body Kepler orbits have closed form analytical solutions and require no numerical integration! You can get position (say $t(\theta)$ and $r(\theta)$ from simple equations, and get $\theta(t)$ (what we normally want) to numerical accuracy by solving using Newton's method with only a few iterations.

But maybe you are doing this as an exercise in numerical computing, in which case you can use those analytical solutions for comparison to see what your error is.

I took a first shot at using single precision (8x speed for GPU's) and noticed some significant error in the orbit.

There's a good chance that your error could be related to your numerical integration algorithm which you haven't mentioned, and you haven't mentioned what your "significant error" is with respect to. Even a simple RK45 will give good results if used skillfully and with variable step size. Canned integrators (ODE solvers) will generally perform well if you use them correctly. In extreme cases Symplectic integrators can improve results, for more on that see answers to What does “symplectic” mean in reference to numerical integrators, and does SciPy's odeint use them?

For more on this I recommend this answer to Calculating the planets and moons based on Newtons's gravitational force where I compare Euler Forward method to RK4 for example. Answers to How to calculate the planets and moons beyond Newtons's gravitational force? may be helpful especially if you want to try things beyond two-body Kepler orbits with Newtonian gravity only.

You may find some helpful information in the question need to understand better ODE solution accuracy vs numerical precision. I'll write an answer there someday if nobody else does.

I am only concerned with error being in the region of 1km or so.

I'm pretty sure if you are only propagating your Kepler orbits for a few orbits this might be possible depending on the situation and size of orbit (which you have not described) if you use the right integration step size. If your routine has an automatic step size then you need to keep an eye on the precision you've specified. If you are using a routine you've coded yourself then read up on step size.

Mathematically small step sizes are better but numerically the more steps the larger the roundoff error which is aggravated by your lower precision. So you will want to look hard at all decisions related to step size.