title: What is wrong with my gravity simulation?
body: <p>As per advice given to me in <a href="https://stackoverflow.com/a/15111343/924313">this answer</a>, I have implemented a Runge-Kutta integrator in my gravity simulator.</p>
  
  <p>However, after I simulate one year of the solar system, the positions are still off by cca 110 000 kilometers, which isn't acceptable.</p>
  
  <p>My initial data was provided by NASA's HORIZONS system. Through it, I obtained position and velocity vectors of the planets, Pluto, the Moon, Deimos and Phobos at a specific point in time.</p>
  
  <p>These vectors were 3D, however, some people told me that I could ignore the third dimension as the planets aligned themselves in a plate around the Sun, and so I did. I merely copied the x-y coordinates into my files.</p>
  
  <p>This is the code of my improved update method:</p>
  
  <pre><code>"""
  Measurement units:
  
  [time] = s
  [distance] = m
  [mass] = kg
  [velocity] = ms^-1
  [acceleration] = ms^-2
  """
  
  class Uni:
      def Fg(self, b1, b2):
          """Returns the gravitational force acting between two bodies as a Vector2."""
  
          a = abs(b1.position.x - b2.position.x) #Distance on the x axis
          b = abs(b1.position.y - b2.position.y) #Distance on the y axis
  
          r = math.sqrt(a*a + b*b)
  
          fg = (self.G * b1.m * b2.m) / pow(r, 2)
  
          return Vector2(a/r * fg, b/r * fg)
  
      #After this is ran, all bodies have the correct accelerations:
      def updateAccel(self):
          #For every combination of two bodies (b1 and b2) out of all bodies:
          for b1, b2 in combinations(self.bodies.values(), 2):
              fg = self.Fg(b1, b2) #Calculate the gravitational force between them
  
              #Add this force to the current force vector of the body:
              if b1.position.x &gt; b2.position.x:
                  b1.force.x -= fg.x
                  b2.force.x += fg.x
              else:
                  b1.force.x += fg.x
                  b2.force.x -= fg.x
  
  
              if b1.position.y &gt; b2.position.y:
                  b1.force.y -= fg.y
                  b2.force.y += fg.y
              else:
                  b1.force.y += fg.y
                  b2.force.y -= fg.y
  
          #For body (b) in all bodies (self.bodies.itervalues()):
          for b in self.bodies.itervalues():
              b.acceleration.x = b.force.x/b.m
              b.acceleration.y = b.force.y/b.m
              b.force.null() #Reset the force as it's not needed anymore.
  
      def RK4(self, dt, stage):
          #For body (b) in all bodies (self.bodies.itervalues()):
          for b in self.bodies.itervalues():
              rd = b.rk4data #rk4data is an object where the integrator stores its intermediate data
  
              if stage == 1:
                  rd.px[0] = b.position.x
                  rd.py[0] = b.position.y
                  rd.vx[0] = b.velocity.x
                  rd.vy[0] = b.velocity.y
                  rd.ax[0] = b.acceleration.x
                  rd.ay[0] = b.acceleration.y
  
              if stage == 2:
                  rd.px[1] = rd.px[0] + 0.5*rd.vx[0]*dt
                  rd.py[1] = rd.py[0] + 0.5*rd.vy[0]*dt
                  rd.vx[1] = rd.vx[0] + 0.5*rd.ax[0]*dt
                  rd.vy[1] = rd.vy[0] + 0.5*rd.ay[0]*dt
                  rd.ax[1] = b.acceleration.x
                  rd.ay[1] = b.acceleration.y
  
              if stage == 3:
                  rd.px[2] = rd.px[0] + 0.5*rd.vx[1]*dt
                  rd.py[2] = rd.py[0] + 0.5*rd.vy[1]*dt
                  rd.vx[2] = rd.vx[0] + 0.5*rd.ax[1]*dt
                  rd.vy[2] = rd.vy[0] + 0.5*rd.ay[1]*dt
                  rd.ax[2] = b.acceleration.x
                  rd.ay[2] = b.acceleration.y
  
              if stage == 4:
                  rd.px[3] = rd.px[0] + rd.vx[2]*dt
                  rd.py[3] = rd.py[0] + rd.vy[2]*dt
                  rd.vx[3] = rd.vx[0] + rd.ax[2]*dt
                  rd.vy[3] = rd.vy[0] + rd.ay[2]*dt
                  rd.ax[3] = b.acceleration.x
                  rd.ay[3] = b.acceleration.y
  
              b.position.x = rd.px[stage-1]
              b.position.y = rd.py[stage-1]
  
      def update (self, dt):
          """Pushes the uni 'dt' seconds forward in time."""
  
          #Repeat four times:
          for i in range(1, 5, 1):
              self.updateAccel() #Calculate the current acceleration of all bodies
              self.RK4(dt, i) #ith Runge-Kutta step
  
          #Set the results of the Runge-Kutta algorithm to the bodies:
          for b in self.bodies.itervalues():
              rd = b.rk4data
              b.position.x = b.rk4data.px[0] + (dt/6.0)*(rd.vx[0] + 2*rd.vx[1] + 2*rd.vx[2] + rd.vx[3]) #original_x + delta_x
              b.position.y = b.rk4data.py[0] + (dt/6.0)*(rd.vy[0] + 2*rd.vy[1] + 2*rd.vy[2] + rd.vy[3])
  
              b.velocity.x = b.rk4data.vx[0] + (dt/6.0)*(rd.ax[0] + 2*rd.ax[1] + 2*rd.ax[2] + rd.ax[3])
              b.velocity.y = b.rk4data.vy[0] + (dt/6.0)*(rd.ay[0] + 2*rd.ay[1] + 2*rd.ay[2] + rd.ay[3])
  
          self.time += dt #Internal time variable
  </code></pre>
  
  <p>The algorithm is as follows:</p>
  
  <ol>
  <li>Update the accelerations of all bodies in the system</li>
  <li>RK4(first step)</li>
  <li>goto 1</li>
  <li>RK4(second)</li>
  <li>goto 1</li>
  <li>RK4(third)</li>
  <li>goto 1</li>
  <li>RK4(fourth)</li>
  </ol>
  
  <p>Did I mess something up with my RK4 implementation? Or did I just start with corrupted data (too few important bodies and ignoring the 3rd dimension)?</p>
  
  <p>How can this be fixed?</p>
  
  <hr />
  
  <p>Explanation of my data etc...</p>
  
  <p>All of my coordinates are relative to the Sun (i.e. the Sun is at (0, 0)).</p>
  
  <pre><code>./my_simulator 1yr
  
  Earth position: (-1.47589927462e+11, 18668756050.4)
  
  HORIZONS (NASA):
  
  Earth position: (-1.474760457316177E+11,  1.900200786726017E+10)
  </code></pre>
  
  <p>I got the <code>110 000 km</code> error by subtracting the Earth's x coordinate given by NASA from the one predicted by my simulator.</p>
  
  <pre><code>relative error = (my_x_coordinate - nasa_x_coordinate) / nasa_x_coordinate * 100
                 = (-1.47589927462e+11 + 1.474760457316177E+11) / -1.474760457316177E+11 * 100
                 = 0.077%
  </code></pre>
  
  <p>The relative error seems miniscule, but that's simply because Earth is really far away from the Sun both in my simulation and in NASA's. The distance is <em>still</em> huge and renders my simulator useless.</p>