title: Why is my n-body simulation not producing orbits (using Earth and Venus around the Sun)?
body: <p>I am trying to make an n-body simulation in python. I have chosen to test the algorithm on a 3-body system, consisting of Earth, Venus, and the Sun, over a duration of 10 years using a time-step of one day. My algorithm will run to completion, but the output does not seem correct; I am expecting Earth and Venus to have nearly circular paths around the Sun; the actual output is shown below:</p>
  
  <p><a href="https://i.stack.imgur.com/RS6TL.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/RS6TL.png" alt="This output is wrong"></a></p>
  
  <p>I am hoping someone can show me either where my implementation is wrong or where my logic is wrong. The implementation is below:</p>
  
  <pre><code>import numpy as np
  import matplotlib.pyplot as plt
  from mpl_toolkits.mplot3d import Axes3D # view 3d
  
  class CelestialBody():
  
      def __init__(self, identifier, facecolor, mass, position, velocity=None, acceleration=None, notes=None):
          """ """
          self.identifier = identifier
          self.facecolor = facecolor
          self.mass = mass
          self.position = np.array(position, dtype=float)
          self.velocity = self.autocorrect_parameter(velocity)
          self.acceleration = self.autocorrect_parameter(acceleration)
          self._previous_positions = [self.position]
          self._previous_velocities = [self.velocity]
          self._previous_accelerations = [self.acceleration]
          self.notes = notes
  
      @property
      def previous_positions(self):
          return np.array(self._previous_positions)
  
      @property
      def previous_velocities(self):
          return np.array(self._previous_velocities)
  
      @property
      def previous_accelerations(self):
          return np.array(self._previous_accelerations)
  
      def autocorrect_parameter(self, args):
          """ """
          if args is None:
              return np.zeros(self.position.shape, dtype=float)
          else:
              return np.array(args, dtype=float)
  </code></pre>
  
  <p>From the class defined above, we can instantiate <code>Earth</code>, <code>Venus</code>, and the <code>Sun</code>. Mass is given in units of kg, positions in meters, velocity in m/s, acceleration in m/s^2. The Sun is initially located at the origin; the initial position of Earth is completely in the x-direction (with its velocity directed in the y-direction, tangential ...); and the initial position of Venus is completely in the y-direction (with is velocity directed in the x-direction). The values used are <a href="https://www.wolframalpha.com/input/?i=venus%20orbital%20speed" rel="nofollow noreferrer">obtained and converted via wolframalpha</a>.</p>
  
  <pre><code>Sun = CelestialBody(identifier='Sun', facecolor='yellow', mass=1.988e30, position=[0, 0, 0])
  Earth = CelestialBody(identifier='Earth', facecolor='blue', mass=5.972e24, position=[1.496e11, 0, 0], velocity=[0, 29766.5, 0])
  Venus = CelestialBody(identifier='Venus', facecolor='darkorange', mass=4.867e24, position=[0, 1.082e11, 0], velocity=[34995.75, 0, 0])
  </code></pre>
  
  <p>Below is the method to run the simulation. First, a for-loop is used to continually sum the acceleration caused by the force of the <code>j-th</code> mass onto the <code>i-th</code> mass over the distance that separates them (where <code>i =/= j</code>). From the net acceleration, velocity is obtained by multiplying the acceleration with the time-step; the position is calculated by multiplying the velocity with the time-step. This happens for each time-step in the outer loop; for each time-step, the newly calculated position, velocity, and acceleration are appended into their respective attributes of the class <code>CelestialBody</code>. </p>
  
  <pre><code>class NBodySimulation():
  
      def __init__(self, bodies, gravitational_constant=6.67e-11, t=0):
          """ """
          self.bodies = bodies
          self.gravitational_constant = gravitational_constant
          self.t = t
          self.collision_times = []
  
      def update_accelerations(self):
          """ """
          for ith_body, ibody in enumerate(self.bodies):
              ibody.acceleration *= 0
              for jth_body, jbody in enumerate(self.bodies):
                  if ith_body != jth_body:
                      # dpos = ibody.position - jbody.position
                      dpos = jbody.position - ibody.position
                      r = np.sum(dpos**2)
                      acc = -1 * self.gravitational_constant * jbody.mass * dpos / np.sqrt(r**3)
                      ibody.acceleration += acc
              ibody._previous_accelerations.append(ibody.acceleration)
              self.bodies[ith_body] = ibody
  
      def update_velocities_and_positions(self, dt):
          """ """
          for ith_body, ibody in enumerate(self.bodies):
              ibody.velocity += ibody.acceleration * dt
              ibody.position += ibody.velocity * dt
              ibody._previous_velocities.append(ibody.velocity)
              ibody._previous_positions.append(ibody.position)
              self.bodies[ith_body] = ibody
  
      def run_simulation(self, dt, duration):
          """ """
          nsteps = int(duration / dt)
          for ith_step in range(nsteps):
              self.update_accelerations()
              self.update_velocities_and_positions(dt)
              self.t += dt
  
  
      def view_positions(self):
          """ """
          if self.bodies[0].position.size == 2:
              fig, ax = plt.subplots()
              ndim = 2
          elif self.bodies[0].position.size == 3:
              fig = plt.figure()
              ax = fig.add_subplot(1, 1, 1, projection='3d')
              ndim = 3
          npos = len(self.bodies[0].previous_positions)
          for ibody in self.bodies:
              if ndim == 2:
                  ax.plot(ibody.previous_positions[0, :], ibody.previous_positions[1, :], color=ibody.facecolor, label=ibody.identifier, alpha=0.5)
              else:
                  ax.plot(ibody.previous_positions[0, :], ibody.previous_positions[1, :], ibody.previous_positions[2, :], color=ibody.facecolor, label=ibody.identifier, alpha=0.5)
          # ax.view_init(azim=60, elev=45)
          fig.legend(mode='expand', loc='lower center',     ncol=len(self.bodies))
          plt.show()
          plt.close(fig)
  </code></pre>
  
  <p>To run the code, </p>
  
  <pre><code>bodies = [Sun, Earth, Venus]
  dt = 86400 # seconds in one day
  duration = 315360000 # seconds in 10 years
  NB = NBodySimulation(bodies)
  NB.run_simulation(dt, duration)
  NB.view_positions()
  </code></pre>
  
  <p>Is my output incorrect because of faulty implementation or faulty logic? Any help would be appreciated.</p>