N-body Solar System questions

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).

Python:
               # 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
Code:
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
Code:
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!

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 :)

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

@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
Code:
# 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.

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).

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 here, github here.

D H
Staff Emeritus