title: Efficiently compute n-body gravitation in python
body: <p>I am trying to compute the accelerations due to gravity for an n-body problem in 3-space (I'm using symplectic Euler).</p>
  
  <p>I have position and velocity vectors for each time step, and am using the below (working) code to calculate accelerations and update velocity and position. Note that the accelerations are vectors in 3-space, not just magnitudes.</p>
  
  <p>I would like to know if there's a more efficient way to compute this with numpy to avoid the loops.</p>
  
  <pre><code>def accelerations(positions, masses):
      '''Params:
      - positions: numpy array of size (n,3)
      - masses: numpy array of size (n,)
      Returns:
      - accelerations: numpy of size (n,3), the acceleration vectors in 3-space
      '''
      n_bodies = len(masses)
      accelerations = numpy.zeros([n_bodies,3]) # n_bodies * (x,y,z)
  
      # vectors from mass(i) to mass(j)
      D = numpy.zeros([n_bodies,n_bodies,3]) # n_bodies * n_bodies * (x,y,z)
      for i, j in itertools.product(range(n_bodies), range(n_bodies)):
          D[i][j] = positions[j]-positions[i]
  
      # Acceleration due to gravitational force between each pair of bodies
      A = numpy.zeros((n_bodies, n_bodies,3))
      for i, j in itertools.product(range(n_bodies), range(n_bodies)):
          if numpy.linalg.norm(D[i][j]) &gt; epsilon:
              A[i][j] = gravitational_constant * masses[j] * D[i][j] \
              / numpy.linalg.norm(D[i][j])**3
  
      # Calculate net acceleration of each body (vectors in 3-space)
      accelerations = numpy.sum(A, axis=1) # sum of accel vectors for each body of shape (n_bodies,3)
  
      return accelerations
  </code></pre>