This code gives two versions of an N-body integrator that are functional. The first one uses the approach of the mutable methods that does the minimum number of distance calculations, but at the cost of using updated on Vector. The second version does twice as many distance calculations, but doesn't have to rebuild data structures nearly as much…

NBodyFunctional.scala

object NBodyFunctional {

  def initBodies(numBodies: Int): Vector[ImmutableBody] = {
    Vector.tabulate(numBodies) { i =>
      ImmutableBody(Vect3(i,0,0), Vect3(0,math.sqrt(1.0/i),0),
        if(i==0) 1 else 1e-10)
    }
  }

  def forSim(bodies: Vector[ImmutableBody], steps: Int, dt: Double): Vector[ImmutableBody] = {
    if(steps <= 0) bodies else {
      def calcAccel(accel: Vector[Vect3], i: Int, j: Int): Vector[Vect3] = {
        if(i >= bodies.length) accel
        else if(j >= bodies.length) calcAccel(accel, i+1, i+2)
        else {
          val pi = bodies(i)
          val pj = bodies(j)
          val d = pi.p - pj.p
          val dist = math.sqrt(d.x*d.x+d.y*d.y+d.z*d.z)
          val a = accel.updated(i, accel(i) - d*(pj.mass/(dist*dist*dist)))
            .updated(j, accel(j) + d*(pi.mass/(dist*dist*dist)))
          calcAccel(a, i, j+1)
        }
      }
      val accel = calcAccel(Vector.fill(bodies.length)(Vect3(0,0,0)), 0, 1)
      forSim((bodies,accel).zipped.map((b, a) => b.step(a, dt)), steps-1, dt)
    }
  }

  def forSim2(bodies: Vector[ImmutableBody], steps: Int, dt: Double): Vector[ImmutableBody] = {
    if(steps <= 0) bodies else {
      val accel = for(i <- 0 until bodies.length) yield {
        val pi = bodies(i)
        (0 until bodies.length).foldLeft(Vect3(0,0,0)) { (a, j) =>
          if(j == i) a else {
            val pj = bodies(j)
            val d = pi.p - pj.p
            val dist = math.sqrt(d.x*d.x+d.y*d.y+d.z*d.z)
            a - d*(pj.mass/(dist*dist*dist))
          }
        }
      }
      forSim2((bodies,accel).zipped.map((b, a) => b.step(a, dt)), steps-1, dt)
    }
  }
}