Previous ToC Up Next

2. The Body Class

2.1. Code Listing

Alice: It would be good to print out that whole Body class, to get an idea of what is happening on the ground level.

Bob: My pleasure! Here it all is. 

 class Body
 
   attr_accessor :acc, :jerk, 
                 :pred_pos, :pred_vel,
                 :time, :next_time
 
   def autonomous_step(ba, dt_param)
     take_one_step(ba, @next_time, dt_param)
   end
 
   def forced_step(ba, t, dt_param)
     take_one_step(ba, t, dt_param)
   end
 
   def take_one_step(ba, t, dt_param)
     ba.each do |b|
       b.predict_step(t)
     end
     correct_step(ba, t, dt_param)
   end
 
   def predict_step(t)
     if t > @next_time
       STDERR.print "predict_step: t = ", t, " > @next_time = "
       STDERR.print @next_time, "\n"
       exit
     end
     dt = t - @time
     @pred_pos = @pos + @vel*dt + @acc*(dt*dt/2.0) + @jerk*(dt*dt*dt/6.0)
     @pred_vel = @vel + @acc*dt + @jerk*(dt*dt/2.0)
   end
 
   def correct_step(ba, t, dt_param)
     dt = t - @time
     new_acc, new_jerk = get_acc_and_jerk(ba)
     new_vel = @vel + (@acc + new_acc)*(dt/2.0) +        # first compute new_vel
                       (@jerk - new_jerk)*(dt*dt/12.0)   # since new_vel is used
     new_pos = @pos + (@vel + new_vel)*(dt/2.0) +        # to compute new_pos
                       (@acc - new_acc)*(dt*dt/12.0)
     @pos = new_pos
     @vel = new_vel
     @acc = new_acc
     @jerk = new_jerk
     @pred_pos = @pos
     @pred_vel = @vel
     @time = t
     @next_time = @time + collision_time_scale(ba) * dt_param
   end  
 
   def collision_time_scale(body_array)
     time_scale_sq = VERY_LARGE_NUMBER
     body_array.each do |b|
       unless b == self
         r = b.pred_pos - @pred_pos                                           
         v = b.pred_vel - @pred_vel                                           
         r2 = r*r
         v2 = v*v
         estimate_sq = r2 / v2            # [distance]^2/[velocity]^2 = [time]^2
         if time_scale_sq > estimate_sq
           time_scale_sq = estimate_sq
         end
         a = (@mass + b.mass)/r2
         estimate_sq = sqrt(r2)/a         # [distance]/[acceleration] = [time]^2
         if time_scale_sq > estimate_sq
           time_scale_sq = estimate_sq
         end
       end
     end
     sqrt(time_scale_sq)
   end
 
   def get_acc_and_jerk(body_array)
     a = j = @pos*0                         # null vectors of the correct length
     body_array.each do |b|
       unless b == self
         r = b.pred_pos - @pred_pos
         r2 = r*r
         r3 = r2*sqrt(r2)
         v = b.pred_vel - @pred_vel
         a += r*(b.mass/r3)
         j += (v-r*(3*(r*v)/r2))*(b.mass/r3)
       end
     end
     [a, j]
   end    
 
   def ekin                         # kinetic energy
     0.5*@mass*(@vel*@vel)
   end
 
   def epot(body_array)             # potential energy
     p = 0
     body_array.each do |b|
       unless b == self
         r = b.pos - @pos
         p += -@mass*b.mass/sqrt(r*r)
       end
     end
     p
   end
 
 end

2.2. A Familiar Part

Bob: The second half of the Body class is familiar, but with some subtle differences. Take the collision_time_scale for example. There are only two lines that are not exactly the same as we had before, line four and five of the body of the definition. In the shared time step case, we had: 

         r = b.pos - @pos                                                     
         v = b.vel - @vel

where in the individual time step case, we have instead: 

         r = b.pred_pos - @pred_pos                                           
         v = b.pred_vel - @pred_vel

This is exactly what is needed here: in order to push a particle forward, we predict the position at it latest sales date time @next_time, and we ask all other particles to predict where they will be at that very same time. We can then compute the differences in position and velocity at that time, and estimate the new time step.

Similarly, the computation of acceleration and jerk are also almost the same as they were before. Since we are here working only with the Hermite scheme, I thought it would be simpler to combine the previous methods acc and jerk into one combined method acc_and_jerk. And here, too, I have switched to predicted positions for both position and velocity.

Finally, the methods ekin and epot are completely unchanged.

Alice: Shouldn't they be based on predicted positions and velocity too?

Bob: No, since I invoke these diagnostics methods only after synchronizing the whole system, in which case we can use the old form of the methods verbatim. We'll come to that. First let's look at the top of the Body class, where the most important differences reside.

2.3. Predicting and Correcting

Alice: Ah, there at the very top is our friend autonomous_step, which invited us to come here in the first place.

Bob: Yes. This method takes as a first argument the reference to the array of all bodies in the N-body system. It needs that as a handle to reach all other particles. The second parameter is the factor with which we multiply the time scale to obtain the new time step. That part is the same as before, the only difference being that each particle now computes its own time scale; the parameter dt_param is a constant, an overall scaling factor that the user can set at the beginning of the run.

Alice: And all that autonomous_step does, is take a step from the present time time to the @next_time.

Bob: Exactly. And you see in the next definition why I have added that extra time argument in the middle: when we need to synchronice all particles, we can no longer give them the luxury to step autonomously. Instead of letting them choose their own next time, we need to force them to all converge at a fixed target time t, the middle argument of the method forced_time_step. The only difference is that in the autonomous case, a particle can determine itself where it wants to force itself to go to, so to speak.

Alice: So to speak. Trying to explain code in words is always a tricky business, but with the code in hand, I see what you mean. Okay, let's move to where the real work is being done, in the method take_one_step. For starters, you ask all other particles to predict their position to time t.

Bob: All other particles, yes, and also the particle itself that is doing the asking. Note that I have written ba.each, without taking exception to self. So in one stroke I push all particles virtually to the target time t, self as well as others.

Alice: And then you correct something, I presume, in correct_step. What is it that you correct?

Bob: It is in correct_step that we finally encounter the full Hermite algorithm. Let's compare this with the shared time step case. There we sent the Hermite algorithm to individual particles using the following method on the NBody level: 

   def hermite
     calc(" @old_pos = @pos ")
     calc(" @old_vel = @vel ")
     calc(" @old_acc = acc(ba) ")
     calc(" @old_jerk = jerk(ba) ")
     calc(" @pos += @vel*dt + @old_acc*(dt*dt/2.0) + @old_jerk*(dt*dt*dt/6.0) ")
     calc(" @vel += @old_acc*dt + @old_jerk*(dt*dt/2.0) ")
     calc(" @vel = @old_vel + (@old_acc + acc(ba))*(dt/2.0) +
                       (@old_jerk - jerk(ba))*(dt*dt/12.0) ")
     calc(" @pos = @old_pos + (@old_vel + vel)*(dt/2.0) +
                       (@old_acc - acc(ba))*(dt*dt/12.0) ")
   end

Lines five and six of the body of this method do exactly what is now done in predict_step, which I will show here again: 

   def predict_step(t)
     if t > @next_time
       STDERR.print "predict_step: t = ", t, " > @next_time = "
       STDERR.print @next_time, "\n"
       exit
     end
     dt = t - @time
     @pred_pos = @pos + @vel*dt + @acc*(dt*dt/2.0) + @jerk*(dt*dt*dt/6.0)
     @pred_vel = @vel + @acc*dt + @jerk*(dt*dt/2.0)
   end

And the next two lines in the old code do exactly what is now done in correct_step in the new code: 

   def correct_step(ba, t, dt_param)
     dt = t - @time
     new_acc, new_jerk = get_acc_and_jerk(ba)
     new_vel = @vel + (@acc + new_acc)*(dt/2.0) +        # first compute new_vel
                       (@jerk - new_jerk)*(dt*dt/12.0)   # since new_vel is used
     new_pos = @pos + (@vel + new_vel)*(dt/2.0) +        # to compute new_pos
                       (@acc - new_acc)*(dt*dt/12.0)
     @pos = new_pos
     @vel = new_vel
     @acc = new_acc
     @jerk = new_jerk
     @pred_pos = @pos
     @pred_vel = @vel
     @time = t
     @next_time = @time + collision_time_scale(ba) * dt_param
   end

Alice: I see. It is beginning to make sense for me now. And I am also starting to remember how we derived the Hermite algorithm in the first place. In the old code, two steps were involved. Starting from a given position, we made a tentative move to a new position, in order to calculate the new acceleration and jerk. Then we combined that result with the knowledge of the old acceleration and jerk, and together we were able to reconstruct the higher derivatives snap and crackle at the start of the step. Finally, we built a Taylor series using the acceleration, jerk, snap and crackle at the starting time, to make a more accurate determination of the ending position.

Bob: Yes, that is exactly what happens. And what you called here a tentative step is what I call predict_step, while what you called a more accurate determination I call correct_step.

2.4. Shifting Perspective

Alice: That all makes a lot of sense. And the next four lines, following the actual correction part of correct_step are the same as the first four lines in our old hermite method. The only difference is one of notation: in the nbody_sh1.rb we compared the old values with the current values, while here in nbody_ind1.rb we compare the current values with the new values.

Bob: Yes, I hadn't quite realized that I had implicitly shifted my perspective by one time step. Why did I do that? Ah, of course. Interesting! In the past we were dealing with shared time steps. So as soon as we were ready to move, we could put ourselves one step ahead, so to speak, in our imagination, while looking back on the previous values as the old values.

However, in the individual time step case, most of the time we are moving forward in time only in a virtual way: when particles predict their position, in order to help another particle find its way, they don't really move themselves. Since they stay at their old place, it makes more sense to called the old place the current place, current for each particle that is, and to call the predicted time the new time.

Alice: So your shift in perspective makes sense. To predict, you have to look ahead. But if you're walking in lock step, you move with the whole group, and all you can do is look back.

Bob: You may be stretching the analogy here, but if it helps you remember the notation, fine! Now the remaining two lines of correct_step update the two time variables that belong to our particle: the new current time @time, for which @pos and @vel are actually calculated, is t, and the new @next_time is found by calculating the next time step, as we have already seen.
Previous ToC Up Next