1. Limited time only! Sign up for a free 30min personal tutor trial with Chegg Tutors
    Dismiss Notice
Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

I Solar system simulation errors

  1. Jul 14, 2017 #1
    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:
    Code (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.
  2. jcsd
  3. Jul 14, 2017 #2


    User Avatar
    2017 Award

    Staff: Mentor

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

    Are you sure the initial values are accurate and used correctly?
  4. Jul 14, 2017 #3
    The initial values are from the same dataset as the diagnostics. The main loop is given below:
    Code (Python):
    def gen_coords():
        if not os.path.exists('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.
  5. Jul 14, 2017 #4


    User Avatar
    2017 Award

    Staff: 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.
    That doesn't answer my question, as the speed doesn't enter the diagnostics data if I understand that dataset correctly.
  6. Jul 14, 2017 #5
    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.
  7. Jul 14, 2017 #6


    User Avatar
    Science Advisor
    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.
  8. Jul 14, 2017 #7
    I am using GM, I rechecked the values and they seem correct. Thanks for the suggestion :)
  9. Jul 14, 2017 #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.
  10. Jul 14, 2017 #9
    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.
  11. Jul 15, 2017 #10


    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    Thinking aloud.. if the problem was the integrator would you expect to get a spiral rather than a closed eclipse?
  12. Jul 15, 2017 #11
    I don't think that speculations about the expected output of integrators with unknown errors are a useful strategy to solve this problem.
  13. Jul 15, 2017 #12


    User Avatar
    Science Advisor
    Homework Helper
    Gold Member

    No I agree, trying a different integrator as you suggested is the way to go.
  14. Jul 15, 2017 #13
    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?
  15. Jul 15, 2017 #14


    User Avatar
    2017 Award

    Staff: 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?
  16. Jul 15, 2017 #15
    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.
  17. Jul 15, 2017 #16
    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 (Text):
                        X,                      Y,                      Z,                     VX,                     VY,                     VZ,
    The other planets do have similar issues.

  18. Jul 15, 2017 #17


    User Avatar
    2017 Award

    Staff: 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.
  19. Jul 15, 2017 #18
    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.
Share this great discussion with others via Reddit, Google+, Twitter, or Facebook

Have something to add?
Draft saved Draft deleted