Is Leapfrog the best choice for N-body simulations?

[iLoveBigData]

Hey guys!

Working on a side project using machine learning and n-body simulations. I could just use the library Rebound (http://rebound.readthedocs.io/en/latest/) but I decided to first try building my own simulator before taking shortcuts.

Simulation uses the masses of the bodies (from Sun to Pluto) in kgs, gets the positions and velocities for a certain date from NASA's Horizons (positions in AU and velocities in AU/hr)
- I convert G to the appropriate units.
- I use a timestep of 0.001
- I run till time = 2000.

- Integrating using leapfrog method.
- Using numpy's longdouble for better precision.
- Calculating acceleration using following code (brute particle to particle method).

               # Distances between the two bodies.
                dx = other.x - body.x
                dy = other.y - body.y
                dz = other.z - body.z

                r = math.sqrt(dx * dx + dy * dy + dz * dz)

                # Gravitational force between particles.
                f = -G / (r *r * r) * other.m

                body.ax += f * dx
                body.ay += f * dy
                body.az += f * dz

t = 0

SUN, m=1.00002250368 x=-0.00379862800765 y=0.00681683726634 z=0.0029805622517
MERCURY, m=1.66054877697e-07 x=-0.0671159633132 y=-0.40333325859 z=-0.209555180056
VENUS, m=2.44832880699e-06 x=0.719400450883 y=0.0713573963636 z=-0.0137387980893
EARTH, m=3.04135120774e-06 x=-0.170126299491 y=0.895975454466 z=0.388429986171
MARS, m=3.22781112204e-07 x=-1.32389966792 y=-0.811358049279 z=-0.336669172693

t = 2000

SUN, m=1.00002250368 x=-0.00449593006396 y=0.00665448086074 z=0.00293125410067 
MERCURY, m=1.66054877697e-07 x=2.20774164708 y=-1.9084499026 z=-1.24941970561 
VENUS, m=2.44832880699e-06 x=1.52493838842 y=2.00993141858 z=0.80753519142 
EARTH, m=3.04135120774e-06 x=-1.93442316424 y=1.25283838345 z=0.543155283391 
MARS, m=3.22781112204e-07 x=-0.886516277866 y=-1.83209935391 z=-0.81665499933

Questions.
- If I plot it for a longer time duration, the bodies are moving further and further away from the Sun. Am I doing something wrong? Is it because Python cannot handle the right precision, my timestep is not small enough, leapfrog not good enough?

If you could give me some direction or hints, that would be great!

 

[Leonarte]

I have implemented this very same problem (a N-body simulation, which ended up as the solar system simulation).
By the time, I faced the same issue of the planets moving far away from the Sun, and I was not using numpy (just "raw" python).

The problem is not the precision, but the method used to calculate position and/or speed. I started too using this basic method you are using. How I solved it? I changed to Runge-Kutta 4th order method, found out the Differential Equations to use and all worked with great precision.

I think this is the way you should follow. However, there are other methods to solve DE that you could use, but RK4 does the job :)

 

[TurtleMeister]

A few years ago I had some fun playing around with n-body simulations using the leapfrog method. I also used data from Horizons to seed my simulation. I first did it in Java and later in VB.Net. Your methods for calculating dx, dy, dz, and r are the same as mine. But I don't know what you're doing after that. I think you should calculate f using the mass of both bodies then calculate it for each of the three planes. Here's the process I used:

Calculate force:
f = G * body.m * other.m / r^2

Calculate the force in each plane:
f.x = f * dx / r
f.y = f * dy / r
f.z = f * dz / r

Calculate body acceleration in each plane:
body.ax += f.x / body.m
body.ay += f.y / body.m
body.az += f.z / body.m

Calculate body velocity in each plane:
body.vx += Ts * body.ax
body.vy += Ts * body.ay
body.vz += Ts * body.az

Calculate new position for body:
body.rx += Ts * body.vx
body.ry += Ts * body.vy
body.rz += Ts * body.vz

 

[iLoveBigData]

@Leonarte I'll implement RK4 and report back. Yeah plenty of ODEs. The library uses IAS15 15th order integrator but that might take some time to implement.

@TurtleMeister for each step

# Calculate drift phase of symplectic

        for b in self.bodies:
            b.x += 0.5 * dt * b.vx
            b.y += 0.5 * dt * b.vy
            b.z += 0.5 * dt * b.vz

        time += dt / 2

# Calculate accelerations (using code posted in earlier reply..)
# Calculate kick-drift phase of symplectic integration.

        for b in self.bodies:
            b.vx += dt * b.ax
            b.vy += dt * b.ay
            b.vz += dt * b.az

            b.x += 0.5 * dt * b.vx
            b.y += 0.5 * dt * b.vy
            b.z += 0.5 * dt * b.vz

        time += dt / 2

I originally had something similar (using Euler integration) but after integrating for a while the results weren't looking that convincing on my end.

 

[m4r35n357]

Hey guys!

Glad to see I'm not the only one doing this stuff ;)

I'm off out tonight so I'll have to be brief for now, but I've got an (unadvertised) N-body simulator with Horizons solar system data and a symplectic integrator here. Data is in the ic directory. There is a Python version on GitHub too but it's a bit out of date now, I might merge it into the Vala project at some point.

As far as I can tell it is 100% working. I don't scale anything and my time step is 10,000s (about 3 hours!).

I'll check back in tomorrow (UK).

 

[m4r35n357]

Glad to see I'm not the only one doing this stuff ;)

I've just remembered, found and uploaded another partially written n-body sim I started a couple of years ago. This one uses three.js. It works, but I don't have enough useful initial conditions to do much ATM. Demo http://m4r35n357.github.io/nbody/, github here.

 

[D H]

Leapfrog is needed when the number of particles is very, very large. (Ridiculously large. Think of simulating molecules in a fluid, or stars in a galaxy. In both cases the number of objects is so very large that it is not possible to simulate each and every one of them.) Using leapfrog to propagate the main objects in a star system is highly suboptimal.

There's always a tradeoff between accuracy, stability, and efficiency. You cannot have all three. You are lucky if you can have even two. Leapfrog is stable but it is incredibly inefficient (you have to make very small steps) and it is incredibly lousy at long term accuracy. RK4 is considerably more efficient than leapfrog, is somewhat accurate for shortish (a few tens or hundreds of orbits) periods of time, but it is unstable. Higher order techniques such as IAS15 offer benefits, but they can also be problematic.

For example, I found that 16th order Gauss-Jackson is considerably less accurate and less stable than is 8th order Gauss-Jackson if one uses the standard IEEE double precision arithmetic as implemented in microcode. Using higher precision arithmetic means no hardware support, and that in turn means the propagation is incredibly slow.

 

[m4r35n357]

Leapfrog is needed when the number of particles is very, very large. (Ridiculously large. Think of simulating molecules in a fluid, or stars in a galaxy. In both cases the number of objects is so very large that it is not possible to simulate each and every one of them.) Using leapfrog to propagate the main objects in a star system is highly suboptimal.

There's always a tradeoff between accuracy, stability, and efficiency. You cannot have all three. You are lucky if you can have even two. Leapfrog is stable but it is incredibly inefficient (you have to make very small steps) and it is incredibly lousy at long term accuracy. RK4 is considerably more efficient than leapfrog, is somewhat accurate for shortish (a few tens or hundreds of orbits) periods of time, but it is unstable. Higher order techniques such as IAS15 offer benefits, but they can also be problematic.

Hmm, I was under the impression that leapfrog is a second order symplectic scheme, so it should have excellent long term stability. As I understand it you need small steps only because it is not amenable to variable step size. Are either of these assumptions incorrect?

All my simulators can use a composed second order (which I believe is effectively leapfrog) symplectic base (configurable between order 2 and 10), or RK4. I agree that RK4 is more efficient as it does not need to evaluate derivatives.