# $N$-body simulation of the solar system [closed]

I'm trying to write a Solar System simulator in Python, with the Body class given below, and I have used initial values obtained from Nasa's Horizon page. However I get inaccurate predictions from the simulator - after one simulated day the predictions are off by an order of at least 10^6 m no matter how small a timestep I use, when compared to the data from Nasa. In the step method I use RK4 and add the accelerations from each interacting body to the second derivative. I can't figure out what's wrong and would appreciate some help :)

Thanks in advance.

EDIT: I forgot to mention another area of concern: as seen in the image, the orbit of Mars passes awfully close to earth - which I don't think is realistic. class Body:

def __init__(self, name, x0, y0, z0, vx0, vy0, vz0, mass, radius):

# Constants of nature
# Universal constant of gravitation
self.G = 6.67408e-11

# Name of the body (string)
self.name = name

# Initial position of the body (m)
self.x0 = x0
self.y0 = y0
self.z0 = z0

# Position (m). Set to initial value.
self.x = self.x0
self.y = self.y0
self.z = self.z0

# Initial velocity of the body (m/s)
self.vx0 = vx0
self.vy0 = vy0
self.vz0 = vz0

# Velocity (m/s). Set to initial value.
self.vx = self.vx0
self.vy = self.vy0
self.vz = self.vz0

# Mass of the body (kg)
self.M = mass

# Radius of the body (m)
self.radius = radius

def compute_acceleration(self, x, y, z):
"""Computes the gravitational acceleration due to self at position (x, y) (m)"""
# Deltas
delta_x = self.x - x
delta_y = self.y - y
delta_z = self.z - z

# Acceleration in the x-direction (m/s^2)
ax = self.G * self.M / (delta_x ** 2 + delta_y ** 2 + delta_z ** 2) * \
delta_x / np.sqrt(delta_x ** 2 + delta_y ** 2 + delta_z ** 2)

# Acceleration in the y-direction (m/s^2)
ay = self.G * self.M / (delta_x ** 2 + delta_y ** 2 + delta_z ** 2) * \
delta_y / np.sqrt(delta_x ** 2 + delta_y ** 2 + delta_z ** 2)

# Acceleration in the z-direction (ms/s^2)
az = self.G * self.M / (delta_x ** 2 + delta_y ** 2 + delta_z ** 2) * \
delta_z / np.sqrt(delta_x ** 2 + delta_y ** 2 + delta_z ** 2)

return ax, ay, az

def step(self, dt, targets):
"""4th order Runge-Kutta integration"""

# Acceleration due to targets (NumPy array)
a = np.zeros(shape=(1, 3))
for o in targets:
a = a + np.array(o.compute_acceleration(self.x, self.y, self.z))

k1x = self.vx
k1y = self.vy
k1z = self.vz

k1vx = a
k1vy = a
k1vz = a

k2x = self.vx + dt / 2 * k1vx
k2y = self.vy + dt / 2 * k1vy
k2z = self.vz + dt / 2 * k1vz

# Acceleration due to targets (NumPy array)
a = np.zeros(shape=(1, 3))
for o in targets:
a += np.array(o.compute_acceleration(self.x + dt / 2 * k1x,
self.y + dt / 2 * k1y,
self.z + dt / 2 * k1z))

k2vx = a
k2vy = a
k2vz = a

k3x = self.vx + dt / 2 * k2vx
k3y = self.vy + dt / 2 * k2vy
k3z = self.vz + dt / 2 * k2vz

# Acceleration due to targets (NumPy array)
a = np.zeros(shape=(1, 3))
for o in targets:
a += np.array(o.compute_acceleration(self.x + dt / 2 * k2x,
self.y + dt / 2 * k2y,
self.z + dt / 2 * k2z))

k3vx = a
k3vy = a
k3vz = a

k4x = self.vx + dt * k3vx
k4y = self.vy + dt * k3vy
k4z = self.vz + dt * k3vz

# Acceleration due to targets (NumPy array)
a = np.zeros(shape=(1, 3))
for o in targets:
a += np.array(o.compute_acceleration(self.x + dt * k3x,
self.y + dt * k3y,
self.z + dt * k3z))

k4vx = a
k4vy = a
k4vz = a

# Update position
self.x = self.x + dt / 6 * (k1x + 2 * k2x + 2 * k3x + k4x)
self.y = self.y + dt / 6 * (k1y + 2 * k2y + 2 * k3y + k4y)
self.z = self.z + dt / 6 * (k1z + 2 * k2z + 2 * k3z + k4z)

# Update velocity
self.vx = self.vx + dt / 6 * (k1vx + 2 * k2vx + 2 * k3vx + k4vx)
self.vy = self.vy + dt / 6 * (k1vy + 2 * k2vy + 2 * k3vy + k4vy)
self.vz = self.vz + dt / 6 * (k1vz + 2 * k2vz + 2 * k3vz + k4vz)


## closed as off-topic by David Hammen, sammy gerbil, Qmechanic♦Jul 12 '17 at 20:29

• This question does not appear to be about physics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

• I quickly calculate the Earth travels about 900 Gm per day relative to the sun. So 1 Mm is about a 1 ppm error in this number. You've only specified (and it looks like we only know) G to 6 sig figs, so I don't think you should be surprised by this magnitude of error. Maybe astronomers have some better way of predicting celestial motions without relying on the uncertain knowledge of G? – The Photon Jul 12 '17 at 18:09
• Don't use $G$ and $M$. Instead, use the product $GM$, which are known to much higher precision. You can find somewhat recent values for these gravitational parameters in en.wikipedia.org/wiki/Standard_gravitational_parameter . – David Hammen Jul 12 '17 at 18:20
• I'm voting to close this question as off-topic because while computational physics is on topic, we are not a programming site. If your question is about implementing computational code - in particular, if it's about writing, compiling, debugging or optimizing code, or about a specific language or library - then it is off topic. – David Hammen Jul 12 '17 at 18:37
• This is off topic here. – Wrichik Basu Jul 12 '17 at 19:10
• BTW, there are a few issues with your code. Eg, in the compute_acceleration method you needlessly compute delta_x ** 2 + delta_y ** 2 + delta_z ** 2 6 times, and its square root 3 times (which is a relatively slow computation). There's not much point to using Numpy here. It can certainly improve the speed of numerical calculations when working with arrays, but you aren't really harnessing its power with this code. But these are topics for Code Review, not Physics. – PM 2Ring Jul 13 '17 at 12:08