# Solar system simulation errors

• I
Hi, as a side project I am making a solar system simulator in Python, but I am getting extreme inaccuracies.

As the image shows - I am comparing the simulation to data obtained from NASA's Horizon.

I have written a a basic Verlet integrator in the Python class posten below:
Python:
import numpy as np

class Body:
"""A celestial body class, with all initial values in SI units """

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

# Gravitational parameter
self.GM = gm

# 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

# Radius of the body (m)

# Storage of last integration results
self.last_coords = np.zeros(shape=(2, 3), dtype=float)

# Integration counter
self.i = 0

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

# Deltas squared
delta_x2 = delta_x ** 2
delta_y2 = delta_y ** 2
delta_z2 = delta_z ** 2

# Acceleration in the x-direction (m/s^2)
ax = self.GM / (delta_x2 + delta_y2 + delta_z2) * \
delta_x / np.sqrt(delta_x2 + delta_y2 + delta_z2)

# Acceleration in the y-direction (m/s^2)
ay = self.GM / (delta_x2 + delta_y2 + delta_z2) * \
delta_y / np.sqrt(delta_x2 + delta_y2 + delta_z2)

# Acceleration in the z-direction (ms/s^2)
az = self.GM / (delta_x2 + delta_y2 + delta_z2) * \
delta_z / np.sqrt(delta_x2 + delta_y2 + delta_z2)

return np.array([ax, ay, az])

def step(self, dt, targets):
"""Verlet integrator"""

# dt squared
dt2 = dt ** 2
x0 = np.array([self.x0, self.y0, self.z0])
v0 = np.array([self.vx0, self.vy0, self.vz0])

if self.i == 0:

self.x = x0[0]
self.y = x0[1]
self.z = x0[2]

self.last_coords[0][:] = x0

elif self.i == 1:

a = np.zeros(shape=3, dtype=float)
for o in targets:
a += o.compute_acceleration(x0[0], x0[1], x0[2])

x1 = x0 + v0 * dt + 0.5 * dt2 * a

self.x = x1[0]
self.y = x1[1]
self.z = x1[2]

self.last_coords[1][:] = x1

elif (self.i % 2) == 0:
xnminusone = self.last_coords[0][:]
xn = self.last_coords[1][:]

a = np.zeros(shape=3, dtype=float)

for o in targets:
a += o.compute_acceleration(xn[0], xn[1], xn[2])

xnplusone = 2 * xn - xnminusone + dt2 * a

self.x = xnplusone[0]
self.y = xnplusone[1]
self.z = xnplusone[2]

self.last_coords[0][:] = xnplusone
elif (self.i % 2) != 0:
xnminusone = self.last_coords[1][:]
xn = self.last_coords[0][:]

a = np.zeros(shape=3, dtype=float)

for o in targets:
a += o.compute_acceleration(xn[0], xn[1], xn[2])

xnplusone = 2 * xn - xnminusone + dt2 * a

self.x = xnplusone[0]
self.y = xnplusone[1]
self.z = xnplusone[2]

self.last_coords[1][:] = xnplusone

self.i += 1
And then calling the Body.step() method for each planet individially - but I am unsure whether this is the correct approach. Any advice would be welcome.

Related Classical Physics News on Phys.org
mfb
Mentor
The main loop is missing.
What is the time step?

Are you sure the initial values are accurate and used correctly?

The main loop is missing.
What is the time step?

Are you sure the initial values are accurate and used correctly?
The initial values are from the same dataset as the diagnostics. The main loop is given below:
Python:
def gen_coords():
if not os.path.exists('trajectories/'):
os.makedirs('trajectories/')
# Generate coordinates and save to file

# NumPy arrays for coordinates
trajectory_earth = np.zeros(shape=(n, 3))
trajectory_mars = np.zeros(shape=(n, 3))
trajectory_venus = np.zeros(shape=(n, 3))
trajectory_jupiter = np.zeros(shape=(n, 3))
trajectory_mercury = np.zeros(shape=(n, 3))
trajectory_sun = np.zeros(shape=(n, 3))
trajectory_saturn = np.zeros(shape=(n, 3))
trajectory_uranus = np.zeros(shape=(n, 3))
trajectory_neptune = np.zeros(shape=(n, 3))
trajectory_pluto = np.zeros(shape=(n, 3))

for i in range(n):
trajectory_earth[i][0] = earth.x
trajectory_earth[i][1] = earth.y
trajectory_earth[i][2] = earth.z

trajectory_mars[i][0] = mars.x
trajectory_mars[i][1] = mars.y
trajectory_mars[i][2] = mars.z

trajectory_venus[i][0] = venus.x
trajectory_venus[i][1] = venus.y
trajectory_venus[i][2] = venus.z

trajectory_jupiter[i][0] = jupiter.x
trajectory_jupiter[i][1] = jupiter.y
trajectory_jupiter[i][2] = jupiter.z

trajectory_mercury[i][0] = mercury.x
trajectory_mercury[i][1] = mercury.y
trajectory_mercury[i][2] = mercury.z

trajectory_sun[i][0] = sun.x
trajectory_sun[i][1] = sun.y
trajectory_sun[i][2] = sun.z

trajectory_saturn[i][0] = saturn.x
trajectory_saturn[i][1] = saturn.y
trajectory_saturn[i][2] = saturn.z

trajectory_uranus[i][0] = uranus.x
trajectory_uranus[i][1] = uranus.y
trajectory_uranus[i][2] = uranus.z

trajectory_neptune[i][0] = neptune.x
trajectory_neptune[i][1] = neptune.y
trajectory_neptune[i][2] = neptune.z

trajectory_pluto[i][0] = pluto.x
trajectory_pluto[i][1] = pluto.y
trajectory_pluto[i][2] = pluto.z

earth.step(dt, earth_array)
mars.step(dt, mars_array)
venus.step(dt, venus_array)
jupiter.step(dt, jupiter_array)
mercury.step(dt, mercury_array)
sun.step(dt, sun_array)
saturn.step(dt, saturn_array)
uranus.step(dt, uranus_array)
neptune.step(dt, neptune_array)
pluto.step(dt, pluto_array)

np.savetxt("trajectories/earth.gz", trajectory_earth, delimiter=',')
np.savetxt("trajectories/mars.gz", trajectory_mars, delimiter=',')
np.savetxt("trajectories/venus.gz", trajectory_venus, delimiter=',')
np.savetxt("trajectories/jupiter.gz", trajectory_jupiter, delimiter=',')
np.savetxt("trajectories/mercury.gz", trajectory_mercury, delimiter=',')
np.savetxt("trajectories/sun.gz", trajectory_sun, delimiter=',')
np.savetxt("trajectories/saturn.gz", trajectory_saturn, delimiter=',')
np.savetxt("trajectories/uranus.gz", trajectory_uranus, delimiter=',')
np.savetxt("trajectories/neptune.gz", trajectory_neptune, delimiter=',')
np.savetxt("trajectories/pluto.gz", trajectory_pluto, delimiter=',')

print('Trajectories generated and saved to file')
The image is with a timestep of dt = 86400 seconds (24 hrs) but I have tested down to dt = 86400/256 and there was no visible improvement.

mfb
Mentor
Diagnostic is from actual positions at different times? It looks like Mars is too slow and Earth is too fast. That rules out a wrong GM as source.
The initial values are from the same dataset as the diagnostics.
That doesn't answer my question, as the speed doesn't enter the diagnostics data if I understand that dataset correctly.

The initial values for the simulations, x, y, z position and velocities in all three dimensions, are the first entries in the two diagnostic data sets. The simulator does not estimate the velocities, even though the diagnostic data set does include velocities. The two diagnostic data sets where I got the initial values, covers the same time period. In the image I posted the timesteps in the two data sets were identical - 24 hrs. I can only assume that the data provided by Nasa is accurate.

CWatters
Homework Helper
Gold Member
I noted that the diameter of the orbit is smaller in both cases. Was wondering if you had correct value for G or wrong calculation for the gravitational force or something like that? I'm not a programmer (Python or otherwise) but are you using gm for the gravitational force or GMm/r2 or ??

Sorry if that's way off beam.

Edit: Actually I'm wrong. The simulated orbit for earth is bigger but the simulated orbit for mars is smaller.

I noted that the diameter of the orbit is smaller in both cases. Was wondering if you had correct value for G or wrong calculation for the gravitational force or something like that? I'm not a programmer (Python or otherwise) but are you using gm for the gravitational force or GMm/r2 or ??

Sorry if that's way off beam.

Edit: Actually I'm wrong. The simulated orbit for earth is bigger but the simulated orbit for mars is smaller.
I am using GM, I rechecked the values and they seem correct. Thanks for the suggestion :)

Hi there. It's not clear to me how you're computing the force on each body. Remember that you have to compute the net force on a planet due to all other planets and the sun at each instant in time. My advice would be create a matrix whose elements are the force vectors on each pairs of objects. For instance the (1,2) element would be (vector) force on 1st object due to second object. By Newton's 3rd law you can immediately conclude that the matrix should be anti-symmetric. And if you want to the net force on any body (say ith body) you just sum up all elements in ith row. Make sense? The force matrix needs to be updated every timestep.

It would be more useful to present your code structure rather than the code itself because it would be easier for us to read and critique. The structure of your code should look something like this -

1) The Main program where you will run the simulation (i.e. the time integration)
2) A planet class whose each instance defines a particular body (i.e. a mass a name to label the planet and the dynamic info (position, velocity) at each instant of time). You should have a series of methods in this class to compute things like the K.E. (if required), the acceleration (given the net force on the object),
3) A Solar system class whose single instance consists of a number of instances of the planet class. This class defines your system.

I have written a a basic Verlet integrator [...]
In order to test if the problem comes from the integrator or not, I suggest to replace it with another simple method (e.g. Euler) and then to check if the error remains.

CWatters
Homework Helper
Gold Member
Thinking aloud.. if the problem was the integrator would you expect to get a spiral rather than a closed eclipse?

if the problem was the integrator would you expect to get a spiral rather than a closed eclipse?
I don't think that speculations about the expected output of integrators with unknown errors are a useful strategy to solve this problem.

CWatters
Homework Helper
Gold Member
No I agree, trying a different integrator as you suggested is the way to go.

Hi all.
Before I tried Verlet, I had written a RK4 integrator where I got the same results.
I have rewritten the program and (still) get the same results. Would it be helpful for me to post the code for the new integrator (still Verlet) here, or should I post it on the programming forum?

mfb
Mentor
I still think it is a problem with the initial values. It would help if you could post them.

Do the other planets have similar issues?

Before I tried Verlet, I had written a RK4 integrator where I got the same results.
That makes it very unlikely that the problem results from the integrator. The remaining code also seems to be OK. The error is most likely somewhere else. Maybe in the initialisation or in the starting values.

I still think it is a problem with the initial values. It would help if you could post them.

Do the other planets have similar issues?
Hi,
Below are the starting values, each line is one body and the order is
1. Sun
2. Mercury
3. Venus
4. Earth
5. Mars
6. Jupiter
7. Saturn
8. Uranus
9. Neptune
10. Pluto
The format is x, y, z, vx, vy, vz and the units are meters and meters/second.
Code:
X,                      Y,                      Z,                     VX,                     VY,                     VZ,
0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00,0.000000000000000000e+00
-5.800688680354694366e+10,-5.263914309864824295e+09,4.888657867308663368e+09,-5.109148337453860222e+03,-4.468780042797738133e+04,-3.301153657093841957e+03
1.086277585922677917e+11,-1.111676131126495171e+10,-6.466689718936344147e+09,4.254227849581348892e+03,3.644856710598556674e+04,2.240072398777410854e+02
4.686576603567236328e+10,-1.452026741849145813e+11,-7.969451561190932989e+07,2.881751478690566000e+04,1.070411872270270032e+04,-5.981058904023539924e+01
-1.069678494192601013e+11,2.178946590300487976e+11,7.127760024113416672e+09,-1.987968885764122024e+04,-6.890821978535064773e+03,2.692590721284270785e+02
-7.380730089508210449e+11,-3.454859091579780884e+11,1.787694413697813034e+10,6.314240278907134780e+03,-9.470825232410625176e+03,-1.339925757010264817e+02
-1.304888817871955109e+11,-1.499135374963784912e+12,3.118376082520866013e+10,1.003472377895323007e+04,8.797144891037893331e+02,-4.072343316799688182e+02
2.697672386110554199e+12,1.267068661324523926e+12,-3.030033252663791275e+10,-2.004016234153325286e+03,7.595358087569742565e+03,-6.229132445572105325e-01
4.267155615470120117e+12,-1.365384675364279053e+12,-7.030556256294619751e+10,2.561894377946528948e+03,6.958596862656445410e+03,-2.054695549871400999e+02
1.532853186821589111e+12,-4.749999433589528320e+12,6.476653007042766571e+10,6.248392678614260149e+03,2.299716536952594197e+03,-1.656088331201209940e+03
The other planets do have similar issues.

Thanks!

mfb
Mentor
The speed of Earth (30.74 km/s) is 3% above the average orbital speed, while the distance (152.6 million km) is beyond the aphelion distance. Something is wrong with your initial values.