1. An N-Body Code

1.1. First Attempt

 require "vector.rb"
 
 class Body
 
   attr_accessor :mass, :pos, :vel, :nb
 
   def initialize(mass = 0, pos = Vector[0,0,0], vel = Vector[0,0,0])
     @mass, @pos, @vel = mass, pos, vel
   end
 
   def acc
     a = @pos*0                    # null vector of the correct length        
     @nb.body.each do |b|
       unless b == self                                                       
         r = b.pos - @pos                                                     
         r2 = r*r
         r3 = r2*sqrt(r2)
         a += r*(b.mass/r3)                                                   
       end
     end
     a
   end    
 
   def ekin                        # kinetic energy
     0.5*(@vel*@vel)
   end
 
   def epot                        # potential energy
     p = 0
     @nb.body.each do |b|
       unless b == self
         r = b.pos - @pos
         p += -@mass*b.mass/sqrt(r*r)
       end
     end
     p
   end
 
   def to_s
     "  mass = " + @mass.to_s + "\n" +
     "   pos = " + @pos.join(", ") + "\n" +
     "   vel = " + @vel.join(", ") + "\n"
   end
 
   def pp               # pretty print
     print to_s
   end
 
   def simple_print
     printf("%24.16e\n", @mass)
     @pos.each{|x| printf("%24.16e", x)}; print "\n"
     @vel.each{|x| printf("%24.16e", x)}; print "\n"
   end
 
   def simple_read
     @mass = gets.to_f
     @pos = gets.split.map{|x| x.to_f}.to_v
     @vel = gets.split.map{|x| x.to_f}.to_v
   end
 
 end
 
 class Nbody
 
   attr_accessor :time, :body
 
   def initialize(n=0, time = 0)
     @time = time
     @body = []
     for i in 0...n
       @body[i] = Body.new
       @body[i].nb = self                                           
     end
   end
 
   def evolve(integration_method, dt, dt_dia, dt_out, dt_end)
     nsteps = 0
     e_init
     write_diagnostics(nsteps)
 
     t_dia = dt_dia - 0.5*dt
     t_out = dt_out - 0.5*dt
     t_end = dt_end - 0.5*dt
 
     while @time < t_end
       send(integration_method,dt)
       @time += dt
       nsteps += 1
       if @time >= t_dia
         write_diagnostics(nsteps)
         t_dia += dt_dia
       end
       if @time >= t_out
         simple_print
         t_out += dt_out
       end
     end
   end
 
   def forward(dt)
     old_acc = []
     @body.each_index{|i| old_acc[i] = @body[i].acc}
     @body.each{|b| b.pos += b.vel*dt}
     @body.each_index{|i| @body[i].vel += old_acc[i]*dt}
   end
 
   def leapfrog(dt)
     @body.each{|b| b.vel += b.acc*0.5*dt}
     @body.each{|b| b.pos += b.vel*dt}
     @body.each{|b| b.vel += b.acc*0.5*dt}
   end
 
   def rk2(dt)
     old_pos = []
     @body.each_index{|i| old_pos[i] = @body[i].pos}
     half_vel = []
     @body.each_index{|i| half_vel[i] = @body[i].vel + @body[i].acc*0.5*dt}
     @body.each{|b| b.pos += b.vel*0.5*dt}
     @body.each{|b| b.vel += b.acc*dt}
     @body.each_index{|i| @body[i].pos = old_pos[i] + half_vel[i]*dt}
   end
 
   def rk4(dt)
     old_pos = []
     @body.each_index{|i| old_pos[i] = @body[i].pos}
     a0 = []
     @body.each_index{|i| a0[i] = @body[i].acc}
     @body.each_index{|i|
       @body[i].pos = old_pos[i] + @body[i].vel*0.5*dt + a0[i]*0.125*dt*dt}
     a1 = []
     @body.each_index{|i| a1[i] = @body[i].acc}
     @body.each_index{|i|
       @body[i].pos = old_pos[i] + @body[i].vel*dt + a1[i]*0.5*dt*dt}
     a2 = []
     @body.each_index{|i| a2[i] = @body[i].acc}
     @body.each_index{|i|
       @body[i].pos = old_pos[i] + @body[i].vel*dt +
                      (a0[i]+a1[i]*2)*(1/6.0)*dt*dt}
     @body.each_index{|i| @body[i].vel += (a0[i]+a1[i]*4+a2[i])*(1/6.0)*dt}
   end
 
   def ekin                        # kinetic energy
     e = 0
     @body.each{|b| e += b.ekin}
     e
   end
 
   def epot                        # potential energy
     e = 0
     @body.each{|b| e += b.epot}
     e/2                           # pairwise potentials were counted twice
   end
 
   def e_init                      # initial total energy
     @e0 = ekin + epot
   end
 
   def write_diagnostics(nsteps)
     etot = ekin + epot
     STDERR.print <<END
 at time t = #{sprintf("%g", time)}, after #{nsteps} steps :
   E_kin = #{sprintf("%.3g", ekin)} ,\
  E_pot =  #{sprintf("%.3g", epot)} ,\
  E_tot = #{sprintf("%.3g", etot)}
              E_tot - E_init = #{sprintf("%.3g", etot - @e0)}
   (E_tot - E_init) / E_init = #{sprintf("%.3g", (etot - @e0)/@e0 )}
 END
   end
 
   def pp                                # pretty print
     print "     N = ", @body.size, "\n"
     print "  time = ", @time, "\n"
     @body.each{|b| b.pp}
   end
 
   def simple_print
     print @body.size, "\n"
     printf("%24.16e\n", @time)
     @body.each{|b| b.simple_print}
   end
 
   def simple_read
     n = gets.to_i
     @time = gets.to_f
     for i in 0...n
       @body[i] = Body.new
       @body[i].nb = self
       @body[i].simple_read
     end
   end
 
 end

1.2. Driver

 require "rknbody1.rb"
 
 include Math
 
 dt = 0.0001          # time step
 dt_dia = 10          # diagnostics printing interval
 dt_out = 10          # output interval
 dt_end = 10          # duration of the integration
 ##method = "forward"   # integration method
 ##method = "leapfrog"   # integration method
 ##method = "rk2"   # integration method
 method = "rk4"   # integration method
 
 STDERR.print "dt = ", dt, "\n",
       "dt_dia = ", dt_dia, "\n",
       "dt_out = ", dt_out, "\n",
       "dt_end = ", dt_end, "\n",
       "method = ", method, "\n"
 
 nb = Nbody.new
 nb.simple_read
 nb.evolve(method, dt, dt_dia, dt_out, dt_end)

Alice: The last three lines are almost exactly the same as the last three lines in the driver for our pseudo-one-body integrator:

 b = Body.new                                                    
 b.simple_read                                                   
 b.evolve(method, dt, dt_dia, dt_out, dt_end)                    

apart for three n's and one N. The rest is exactly the same, with of course an extra n in the file name that is being required on top. A very minimal change.

Bob: Yes, I preferred to stay as close as I could to our initial 2-body relative coordinates integrator.

Alice: I'm very curious to see how you integrated the combined equations of motion for an arbitrary number of particles, but let me follow the logic of the driver first, to feel my way around in the program flow. You start by creating a whole N-body system, and then reading in the data. And besides the old class Body, you have introduced a class Nbody, to contain the data for that whole system.

Bob: Yes. The initializer has two parameters, as you can see:

   def initialize(n=0, time = 0)
     @time = time
     @body = []
     for i in 0...n
       @body[i] = Body.new
       @body[i].nb = self                                           
     end
   end

By default an N-body system is created empty, containing 0 particles, at starting at time 0, but if you like, you can create it with several particles right from the start. In our case, the driver creates an empty default system, leaving it up to the simple_read input function to create the necessary particle slots, depending on how many are needed to store the input data.

Alice: Since all particles share the same time step in this code, it makes sense to make the shared simulation time an instance variable @time for the Nbody class. And @body must be an array of Body instances, since if you specify n particles to be present, you fill the array by creating n new bodies, from @body[0] through @body[n-1]. But what about the for loop? Ah, yes, I remember now: in Ruby 0...n means that you exclude n, while 0..n does indeed include n.

Bob: A handy short notation, but like every short notation, potentially confusing if you're not yet used to it.

1.3. Input

       @body[i].nb = self                                           

Bob: I have given the Body class an extra instance variable @nb, which for each body instance will contain the address of the parent, an instance of the Nbody class. In that way, each Body daughter is doubly linked to her Nbody parent. The parent can call the appropriate daughter, by selecting her from the @body[] array, and the daughter can call the parent directly through her own @nb variable. Remember that the expression self gives the address of the Nbody instance itself, which then gets handed down to each body.

Alice: Hmm. In general, I am quite wary of doubly linked list. It is all too easy to change the link in one direction and to forget to change the link in the other direction, or to change it in the wrong way.

Bob: Typical one-off errors that often happen in C and C++ are less likely to occur in Ruby, because so much of the bookkeeping is handled behind the scenes, as long as you don't confuse 0...n and 0..n. But I see your point, and perhaps we should change that, later on. For now, let's just go through the code, and then we can decide whether it will be easy to unlink the backward pointers from daughter to parent.

My motivation to provide backward links was to give each daughter the possibility to communicate with her siblings. If one daughter wants to compute her acceleration, she would need to find the positions of all other daughters, and the simplest way to do that, I thought, was to give her a way to ask her parents how to find all the others.

Alice: Fine for now. In the simple_read program you also provide those backward links for each particle, after which you invoke the simple_read method for that particle:

   def simple_read
     @mass = gets.to_f
     @pos = gets.split.map{|x| x.to_f}.to_v
     @vel = gets.split.map{|x| x.to_f}.to_v
   end

The format you have chosen is to start with the number of particles and the time, followed by the data for each particle.

Bob: Yes. It seemed safer to tell the input routine how many particles to expect, rather than to let it read in everything to the end of the file. In some cases you might want to store more than one snapshot, for example, in one file.

In fact, my code normally will output a series of snapshots, one after each dt_out interval, just as we did it for the single pseudo-body case. This will make it possible to restart a run: you can later sequentially read in a number of those snapshots, each with a nb.simple_read statement in the driver, selecting the proper one by checking the time variable specified.

For example, when you invoke the code with:

  dt_out = 5
  dt_end = 10
  nb = Nbody.new
  nb.simple_read
  nb.evolve(method, dt, dt_dia, dt_out, dt_end)

  nb = Nbody.new
  nb.simple_read
  nb.simple_read
  nb.evolve(method, dt, dt_dia, dt_out, dt_end)

  nb = Nbody.new
  nb.simple_read
  while (nb.time < 10)
    nb.simple_read
  end