N-body simulation still losing precision despite using arbitrary-precision arithmetic and symplectic integrator

Question

A while ago I asked 2-body orbits expanding in nbody simulation?. One thing that I noticed was that JavaScript does not have arbitrary precision, and that I used RK4 as my integrator, while it is accurate short-term, there is an energy drift, causing orbits to decay/expand.

I've scrapped the old concept and have developed a new n-body simulation in Python, with the arbitrary-precision module mpmath with the precision set to 512 decimal digits, as well as incorporating the Verlet algorithm from the previous JavaScript engine.

I tested my code with a testmass (of mass 0) in a circular, zero-inclination orbit, orbiting a Jupiter-mass object at a distance of $10\ 000 \ 000$ meters (at coordinates (10000000,0,0)), with a timestep of $h=0.1$. After running the simulation for 2296 steps (a quarter of an orbit or about 45 in-simulation seconds), I checked the distance of the test mass from the planet. The initial distance was obviously $10\ 000\ 000$ meters, while the new distance was $10\ 000\ 029.529$ meters, which is quite a significant deviation. Am I doing something wrong?

Note: Source code is available if needed.

Answer 1

There are various possible causes for an orbit simulator to behave poorly, even when using an excellent arbitrary precision package like mpmath. When supplying numerical values, you need to pass them to mpmath as strings or integers, not floats, so it can convert them properly. Please see the mpmath docs for further details.

Another possible cause of problems is the catastrophic cancellation which can occur when combining numbers of wildly different sizes. You can often minimize that problem by carefully organising the order of computations, and by working in a scale where the range of magnitudes of your numbers isn't too large. Of course, when using a library like mpmath you can simply increase the precision, although that will slow things down.

The standard thing to do when a correctly coded orbit sim is giving unacceptably high errors is to make the time step smaller. However, using a huge number of time steps increases the cumulative rounding errors, so if you make the time step too small the result will actually be worse, as David Hammen mentioned earlier. Of course, when using mpmath you can get around that problem by increasing the precision.

---

For reference, here's some Python code that does a simple 1-body orbit simulation using the synchronised form of Leapfrog integration, which is actually the same as the velocity Verlet algorithm. It uses plain Python floating-point arithmetic, not mpmath. The Python float type uses the IEEE-754 binary64 format, which has 53 bits of precision (almost 16 decimal digits).

This program uses natural units, i.e., the gravitational parameter of the central body is $GM=\mu=4\pi^2$, so a circular orbit of radius 1 space unit has a period of 1 time unit.

For a simple circular orbit, using 10,000 time steps per orbit, the maximum error in the radius is under 2 parts in 10,000,000. A smaller time step is required for eccentric orbits.

''' Simple orbit sim using the
    synchronised Leapfrog algorithm
    Written by PM 2Ring 2017.05.20
'''

from math import pi

# Gravitational parameter of central body in "natural" units,
# i.e. where R**3 == T**2
mu = 4.0 * pi * pi

# Speed of a body in a circular orbit of radius r
def speed(r):
    return (mu / r) ** 0.5

# Acceleration vector due to central gravity at (x, y)
def acc(x, y):
    a = -mu / (x * x + y * y) ** 1.5
    return a * x, a * y

class Body(object):
    ''' An orbiting body '''
    def __init__(self, x, y, vx, vy, delta_time):
        # Current position
        self.x, self.y = x, y

        # Current velocity
        self.vx, self.vy = vx, vy

        # Time step
        self.delta_time = delta_time

        # Current acceleration due to central gravity
        self.ax, self.ay = acc(x, y)

        self.points = [(x, y)]

    def update(self):
        ''' Update the body's orbit parameters using
            the synchronised leapfrog algorithm
        '''
        dt = self.delta_time
        dt2 = dt * dt

        # Update position
        x, y = self.x, self.y
        x += self.vx * dt + 0.5 * self.ax * dt2
        y += self.vy * dt + 0.5 * self.ay * dt2

        # Update velocity using mean acceleration
        ax, ay = acc(x, y)
        self.vx += 0.5 * (self.ax + ax) * dt
        self.vy += 0.5 * (self.ay + ay) * dt

        self.ax, self.ay = ax, ay
        self.x, self.y = x, y
        lx, ly = self.points[-1]
        if abs(x - lx) > 1 or abs(y - ly) > 1:
            self.points.append((x, y))
        return x, y

# Radius and velocity of a circular orbit
radius = 100.0
print('Semi-major axis:', radius)

v = speed(radius)
period = radius ** 1.5
print('Period:', period)

steps = 10000
delta_time = period / steps
istep = steps // 100

body = Body(radius,0, 0,v, delta_time)

print("step radius x y")
for i in range(1, steps + 2):
   x, y = body.update()
   if i % istep == 0: print(i, (x*x + y*y)**0.5, x, y)

# Plot
P = points(body.points, size=1, color="blue")
P.show(aspect_ratio=1)

That's mostly plain Python code, apart from the last few lines which use Sage to do plotting (via matplotlib).

Here's the output. As you can see, the maximum error occurs at the halfway point.

Semi-major axis: 100.0
Period: 1000.0
step radius x y
100 100.00000001947542 99.8026728363567 6.279052365941193
200 100.00000007782481 99.2114701058467 12.533324180026883
300 100.00000017481786 98.22872501631164 18.738132688008204
400 100.00000031007181 96.85831601489033 24.86899034488533
500 100.0000004830529 95.10565148152806 30.901701456144753
600 100.00000069307836 92.97764838452808 36.812457667191076
700 100.0000009393194 90.48270498238338 42.577931924118325
800 100.00000122080417 87.63066767962191 48.17537053499955
900 100.00000153642179 84.43279216747251 53.582682968369106
1000 100.00000188492673 80.90169900271201 58.7785290345032
1100 100.0000022649435 77.05132380000519 63.74240310543257
1200 100.00000267497244 72.89686223430738 68.45471504130965
1300 100.00000311339528 68.4547100703815 72.89686750374857
1400 100.0000035784818 63.74239845610735 77.05132935101909
1500 100.00000406839649 58.77852473495271 80.90170482543533
1600 100.00000458120591 53.582679050660076 84.43279825988925
1700 100.00000511488614 48.17536703380712 87.63067404816206
1800 100.00000566733111 42.577928875365075 90.48271164233866
1900 100.00000623636048 36.812455106635035 92.97765536027522
2000 100.0000068197286 30.90169941794281 95.10565880655513
2100 100.00000741513318 24.86898886015765 96.85832373162384
2200 100.00000802022444 18.73813178343009 98.22873317574528
2300 100.00000863261437 12.53332387647376 99.21147876697717
2400 100.00000924988619 6.279052677211943 99.8026820654341
2500 100.00000986960366 9.316423505156823e-07 100.00000986960366
2600 100.00001048932121 -6.279050817680861 99.80268342431145
2700 100.00001110659296 -12.533322029108925 99.21148149399528
2800 100.00001171898288 -18.738129959342203 98.22873728916358
2900 100.00001232407413 -24.868987074894257 96.8583292581682
3000 100.00001291947869 -30.901697693138058 95.10566578063481
3100 100.00001350284674 -36.81245347153712 92.97766382296042
3200 100.00001407187602 -42.577927368209636 90.48272164011684
3300 100.0000146243209 -48.17536570300814 87.63068563149452
3400 100.00001515800102 -53.582677955791624 84.43281148159379
3500 100.00001567081028 -58.7785239475003 80.9017197389242
3600 100.00001616072485 -63.74239805996919 77.05134600842264
3700 100.00001662581123 -68.45471016210631 72.89688595395924
3800 100.0000170642339 -72.8968629230579 68.45473532798135
3900 100.00001747426262 -77.0513252072431 63.742425264971644
4000 100.00001785427928 -80.90170126160325 58.77855309408762
4100 100.00001820278406 -84.43279542201651 53.58270894403392
4200 100.00001851840162 -87.6306720835148 48.17539842982273
4300 100.00001879988628 -90.48271069762764 42.577961726543805
4400 100.00001904612719 -92.97765557980723 36.81248934952026
4500 100.00001925615258 -95.10566033036771 30.9017349732278
4600 100.00001942913354 -96.85832669363184 24.869025633039456
4700 100.00001956438729 -98.22873770192699 18.738169664187293
4800 100.00001966138026 -99.21148497362778 12.533362741262588
4900 100.00001971972958 -99.80269005751201 6.279092389071799
5000 100.00001973920497 -100.00001973919642 4.13416988284121e-05
5100 100.00001971972962 -99.80269524924384 -6.279009868830847
5200 100.00001966138022 -99.21149533660193 -12.533280709847839
5300 100.0000195643871 -98.2287531952455 -18.738088445339077
5400 100.00001942913315 -96.85834725614961 -24.868945547291286
5500 100.00001925615217 -95.10568588093415 -30.901656336641377
5600 100.00001904612671 -92.97768601758592 -36.812412472438055
5700 100.00001879988581 -90.48274590249471 -42.577886912364384
5800 100.00001851840106 -87.6307119165327 -48.17532597380332
5900 100.00001820278348 -84.43283972598266 -53.582639132125166
6000 100.00001785427867 -80.90174986167018 -58.778486201805165
6100 100.00001747426198 -77.05137791160872 -63.74236155630852
6200 100.00001706423315 -72.89691952372205 -68.45467505436653
6300 100.00001662581037 -68.4547704356921 -72.89682935326505
6400 100.00001616072397 -63.7424617686046 -77.05129330402596
6500 100.00001567080935 -58.77859083975651 -80.9016711388251
6600 100.00001515800002 -53.582747767675556 -84.4327671775946
6700 100.00001462431985 -48.175438159004266 -87.63064579844279
6800 100.0000140718749 -42.578002182367385 -90.4826864352153
6900 100.00001350284549 -36.81253034859931 -92.97763338514665
7000 100.00001291947744 -30.90177632970622 -95.10564023003296
7100 100.00001232407287 -24.86906716062583 -96.8583086956148
7200 100.00001171898157 -18.738211178175494 -98.22872179580955
7300 100.0000111065917 -12.533404060510405 -99.21147113098588
7400 100.00001048931996 -6.279133337910164 -99.80267823254468
7500 100.00000986960247 -8.175174521653228e-05 -100.00000986956906
7600 100.00000924988495 6.278970156979581 -99.8026872571317
7700 100.00000863261317 12.53324184506617 -99.21148912991777
7800 100.0000080202232 18.73805056458767 -98.228748669031
7900 100.00000741513185 24.868908774414006 -96.85834429410973
8000 100.00000681972718 30.90162078135972 -95.10568435709045
8100 100.00000623635901 36.81237822955507 -92.97768579802373
8200 100.00000566732959 42.577854061186905 -90.48274684717644
8300 100.0000051148845 48.17529457778811 -87.63071388115166
8400 100.00000458120421 53.582609238750955 -84.43284256382812
8500 100.00000406839483 58.7784578426692 -80.90175342547614
8600 100.0000035784801 63.74233474744264 -77.05138205535967
8700 100.00000311339362 68.45464979676468 -72.89692410438892
8800 100.00000267497069 72.89680563361073 -68.45477531487283
8900 100.00000226494173 77.05127109560578 -63.74246681404648
9000 100.00000188492496 80.90165040261006 -58.77859592673901
9100 100.00000153642003 84.4327478634704 -53.582752780233704
9200 100.00000122080235 87.63062784656721 -48.17544299097739
9300 100.00000093931756 90.48266977747897 -42.578006738258736
9400 100.00000069307654 92.97761794671167 -36.812534544236804
9500 100.000000483051 95.10562593092368 -30.901780092697216
9600 100.00000031006994 96.85829545233462 -24.869070430601987
9700 100.0000001748159 98.22870952295547 -18.73821390682732
9800 100.00000007782285 99.21145974283534 -12.533406211414873
9900 100.00000001947349 99.80266764458823 -6.279134886157624
10000 99.99999999999805 99.99999999996388 -8.268337527927994e-05

Here's a live version running on the SageMathCell server.

---

As Wikipedia mentions, there are higher order versions of Leapfrog, discovered independently in the 1980s & 1990s by several people (Forest & Ruth, Yoshida, and Candy & Rozmus). The 4th order version is certainly superior to plain Leapfrog or Verlet, but the improvements vs the increased number of calculations diminish with the higher orders, so there's not much point going higher than 4th or 8th order.

Here's a version of the above code using 4th order Yoshida coefficients.

''' Simple orbit sim
    4th order Leapfrog with Yoshida coefficients
    Written by PM 2Ring 2017.05.20
    Added Yoshida 2022.02.12
'''

from math import pi, radians, cos, sin

# Gravitational parameter of central body in "natural" units,
# i.e. where R**3 == T**2
mu = 4.0 * pi * pi

# Speed of a body in a circular orbit of radius r
def speed(r):
    return (mu / r) ** 0.5

# Acceleration vector due to central gravity at (x, y)
def acc(x, y):
    a = -mu / (x * x + y * y) ** 1.5
    return a * x, a * y

# Yoshida 4th order coefficients
def yosh4():
    q = 2 ** (1 / 3)
    w1 = 1 / (2 - q)
    w0 = -q * w1
    c1 = w1 / 2
    c2 = (w0 + w1) / 2
    return (c1, c2, c2, c1), (w1, w0, w1)

class Body(object):
    ''' An orbiting body '''
    def __init__(self, x, y, vx, vy, delta_time):
        # Current position
        self.x, self.y = x, y

        # Current velocity
        self.vx, self.vy = vx, vy

        # The Yoshida coefficients, multiplied by the time step
        yc, yd = yosh4()
        self.yc = [k * delta_time for k in yc]
        self.yd = [k * delta_time for k in yd]

        self.points = [(x, y)]

    def update(self):
        ''' Update the body's orbit parameters using
            the 4th order synchronised leapfrog algorithm
        '''
        x, y = self.x, self.y
        vx, vy = self.vx, self.vy

        # 4th order Leapfrog integration
        for c, d in zip(self.yc, self.yd):
            x += vx * c
            y += vy * c
            ax, ay = acc(x, y)
            vx += ax * d
            vy += ay * d
        c = self.yc[-1]
        x += vx * c
        y += vy * c

        self.vx, self.vy = vx, vy
        self.x, self.y = x, y

        lx, ly = self.points[-1]
        if abs(x - lx) > 1 or abs(y - ly) > 1:
            self.points.append((x, y))
        return x, y

# Radius and velocity of a circular orbit
radius = 100.0
print('Semi-major axis:', radius)

v = speed(radius)
period = radius ** 1.5
print('Period:', period)

steps = 1000
delta_time = period / steps
print('Step:', delta_time)
istep = steps // 100

# Initial direction. Use 0 for a circular orbit
th = radians(0)
vx, vy = v * sin(th), v * cos(th)

body = Body(radius,0, vx,vy, delta_time)
for i in range(1, steps + 2):
   x, y = body.update()
   if i % istep == 0: print(i, (x*x + y*y)**0.5, x, y)

P = points(body.points, size=1, color="blue")
P.show(aspect_ratio=1)

At only 1000 steps per orbit, the maximum error in the radius has dropped to 1.1 parts per billion. Here's the SageMathCell link.

Here's a plot of 10 periods of an eccentric orbit created using that code. It uses the same time step of 1.0, with an initial angle of 60° to the vertical.

The radius error has jumped up to 5 parts in ten million, per orbit, although that isn't visible at this scale.

If we increase the time step to 2.0, the error goes to 6 parts per million, and it starts to become obvious in the plot.

Answer 2

So there are a number of things that can happen here under the proverbial hood that can cause this kind of error to come up.

First off, round off error is going to be at the heart of all your problems. Let's say your machine precision is about $10^{-16}$, which it is if you're using 64-bit floating-point representation (this is pretty standard across the board, you have to go out of your way to change this). Let's say you're solving the ODE $\frac{d^2 r}{dt^2}=\frac{GM}{r^2}$ and you break it up into individual components, etc. Your distance is 10 000 000 or $10^7$, which is then being squared so you're working with $10^{14}$ on the bottom of that equation, $1.9 \times 10^{27}$ kg on the top, and all of that multiplied by $6.67 \times 10^{-11}$. As you can guess, these massively different numbers can start to create problems with round-off errors (if this isn't clear, let me know and I can add a more concrete example). You will start to lose effective sig figs, and that error can lead to things like you're seeing.

Furthermore, not all solvers are created equal. Symplectic solvers conserve energy, which is great, but there is no free lunch in computational physics; they come at a cost. The exact nature of that cost is not clear to me, but I have always been advised to avoid them unless energy conservation on a granular level is important (my current research is with N-body relativistic objects and my advisor has advised us against symplectic integrators, he gave us a reason but I can’t quite remember what it was, I’ll check into this)

Beyond that, not all solvers are created equally. If you can find an implicit solver, it'll be slower but more accurate, and typically the higher-order the better. Some solvers have capabilities to deal with stiffness and other things while others don't; it starts becoming more of a dark art than a science at a certain level because the reasons some work more than others become murkier and murkier.

tl;dr, try a bunch of different solvers and scale your units closer to unity (e.g., use mass as 1 to reprint 1 Jupiter mass, use something near unity for distance, and rescale G accordingly. Mess with it till G gets within a couple of orders of magnitude to unity as well)

Answer 3

What follows is too long to be a comment. Moreover, since comments can be deleted or moved to chat, I'm making this an answer.

Every numerical integrator based on advancing Cartesian state from one time step to the next has an optimal step size that depends on

• The nature of the problem at hand,
• The units used to represent the problem at hand,
• The mechanisms used to represent the real numbers, and
• The mechanisms used to advance state from one time step to the next.

I'll look at the last two items, with a focus on integration techniques that advance Cartesian position and velocity based on acceleration.

Numerically integrating the position of a vehicle in space inherently is a second order ordinary differential equation (ODE) problem. Some numerical integration techniques (e.g., canonical RK4) can only solve first order ODEs. A second order ODE can be solved using a first order ODE solver by adding auxiliary equations, but doing so tosses geometry out the window. (That this is a second order ODE is geometry.) In addition tossing geometry, techniques such as RK4 also toss the conservation laws out the window.

Every numerical integration technique that advances Cartesian velocity and position based on acceleration faces two basic challenges. If the time step is extremely small, velocity and/or position will be frozen. The problem with extremely small time steps is that $1+10^{-16}$ is exactly equal to one in double precision arithmetic. On the other side of time step size, errors inherent in the technique itself arise if the time step is made to be extremely large. No numerical integrator will perform well with a ten orbit step size.

For a given problem, a given set of units, a given representation technique of the reals, and a given integration technique, there is an optimal step size that lies somewhere between too small and too large. This optimal step size might well not be constant. For example, for an elliptical orbit, it's best to take smaller steps near periapsis and larger steps near apoapsis.

This optimal step size is somewhat easy to find for circular orbits as this makes the optimal step size constant. What one will see is a bathtub error curve. Ideally, integrators should swim close to the deep end of this bathtub.