Solar system simulation errors

  • #1
30
2
Hi, as a side project I am making a solar system simulator in Python, but I am getting extreme inaccuracies.
diagnostic.png

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)
        self.radius = radius

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

Answers and Replies

  • #2
34,681
10,825
The main loop is missing.
What is the time step?

Are you sure the initial values are accurate and used correctly?
 
  • #3
30
2
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.
 
  • #4
34,681
10,825
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.
 
  • #5
30
2
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.
 
  • #6
CWatters
Science Advisor
Homework Helper
Gold Member
10,529
2,297
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.
 
  • #7
30
2
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 :)
 
  • #8
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.
 
  • Like
Likes madsmh
  • #9
1,896
332
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.
 
  • #10
CWatters
Science Advisor
Homework Helper
Gold Member
10,529
2,297
Thinking aloud.. if the problem was the integrator would you expect to get a spiral rather than a closed eclipse?
 
  • #11
1,896
332
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.
 
  • #12
CWatters
Science Advisor
Homework Helper
Gold Member
10,529
2,297
No I agree, trying a different integrator as you suggested is the way to go.
 
  • #13
30
2
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?
 
  • #14
34,681
10,825
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?
 
  • #15
1,896
332
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.
 
  • #16
30
2
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!
 
  • #17
34,681
10,825
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.
 
  • Like
Likes madsmh
  • #18
30
2
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.
Errh.. I had accidentally chosen a another coordinate center than the one I intended. It simulates beautifully now!
So thank you very much, all of you, for your help.
orbit.png
 
  • Like
Likes CWatters and mfb

Related Threads on Solar system simulation errors

Replies
15
Views
547
Replies
18
Views
3K
  • Last Post
Replies
6
Views
2K
Replies
2
Views
2K
  • Last Post
Replies
10
Views
699
  • Last Post
Replies
3
Views
2K
  • Last Post
Replies
16
Views
3K
  • Last Post
Replies
1
Views
2K
  • Last Post
Replies
10
Views
11K
Replies
9
Views
2K
Top