Previous ToC Up Next

7. A Final Version

7.1. Clarity

Bob: It was fun to play with so many different versions, but I'm beginning to get a little confused as to which version did what. Maybe we should pick just one version, and use that while we add more features and move toward specific applications.

Alice: I agree. And since we decided not to worry, for now at least, about performance, let us concentrate on clarity of expression. I must say, I like the last one best, in file rknbody6.rb, where everything fits on one line. However, the version where we gave the Body class explicit helper variables, in file rknbody4.rb, was even shorter.

Bob: Yeah, I did not particularly like the idea of giving this poor Body class all possible helper variables for all possible application. Even though we discussed more clever ways to do this, frankly, I don't care too much what we will choose in the end. I liked your suggestion to send a command string to be evaluated dynamical, thereby generating the proper helper variables, mostly because of its novelty.

Alice: And it goes with the spirit of the times: just-in-time-delivery! But what I like best about this latest method is that it obeys the DRY principle: we are not repeating ourselves.

Bob: Apart from the fact that you repeatedly bring up particular principles.

Alice: I'll repeatedly ignore that. Apart from that point, when clarity is really the criterion, I am not sure whether the last version is really clearer. Let us compare the forward Euler algorithm in both cases, in rknbody4.rb: 

   def forward(dt)
     @body.each{|b| b.old_acc = b.acc(@body)}
     @body.each{|b| b.pos += b.vel*dt}                                        
     @body.each{|b| b.vel += b.old_acc*dt}
   end

and in rknbody6.rb 

   def forward(dt)
     calc(dt," @old_acc = acc(ba) ")
     calc(dt," @pos += @vel*dt ")
     calc(dt," @vel += @old_acc*dt ")
   end

Bob: The last one is clearly shorter.

Alice: That I don't mind so much. I'm just not happy with the fact that it is not clear, for a casual reader, what that variable dt is doing there, as the first two arguments of calc, and the appearance of ba is also a mystery; there is no indication here that ba stands for `body array' and will get its value from @body, somewhere else. In contrast, the earlier version has nothing hidden: the @body.each{|b| . . . } construct is vanilla flavor for someone familiar with Ruby.

7.2. Brevity

Bob: Perhaps we can improve the calc method of Nbody further. How about redefining it in such a way that we can leave out the first argument altogether?

Alice: Ah, that's a good idea. It is also a logical next step, after introducing the shortcut notion of sending a string in the first place. Once we do something that is somewhat dirty and not so self-explanatory, we might as well go all the way.

Bob: I suppose that we would have to introduce an extra instance variable @dt for Nbody. Otherwise it will not be possible to remove the first argument dt from the current calc. In fact, that would make the definition even shorter. So it would look very clean, like this: 

   def forward
     calc(" @old_acc = acc(ba) ")
     calc(" @pos += @vel*dt ")
     calc(" @vel += @old_acc*dt ")
   end

Note that I have created yet another version of our code, in file rknbody7.rb.

Alice: That is short and sweet, indeed. And we have to modify calc on the Nbody level from what we had before: 

   def calc(y,s)
     @body.each{|b| b.calc(@body,y,s)}
   end

to a version with only one parameter, namely the command string: 

   def calc(s)
     @body.each{|b| b.calc(@body, @dt, s)}
   end

Since the other two parameters to the calc method of Body are now both instance variables, their value is known here.

Bob: What else do we have to change? We have to set the value of the new variable @dt in the evolve method, and we have to leave out the dt argument when invoking the integration methods. Two small changes, which leaves us with evolve looking like this: 

   def evolve(integration_method, dt, dt_dia, dt_out, dt_end)
     nsteps = 0
     e_init
     write_diagnostics(nsteps)
     @dt = dt
     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)
       @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

7.3. Correctness

Alice: And the only other changes, besides the change in the forward Euler algorithm we already saw, are the simplified readings of the three other integrators. They now become: 

   def leapfrog
     calc(" @vel += acc(ba)*0.5*dt ")
     calc(" @pos += @vel*dt ")
     calc(" @vel += acc(ba)*0.5*dt ")
   end


   def rk2
     calc(" @old_pos = @pos ")
     calc(" @half_vel = @vel + acc(ba)*0.5*dt ")
     calc(" @pos += @vel*0.5*dt ")
     calc(" @vel += acc(ba)*dt ")
     calc(" @pos = @old_pos + @half_vel*dt ")
   end


   def rk4
     calc(" @old_pos = @pos ")
     calc(" @a0 = acc(ba) ")
     calc(" @pos = @old_pos + @vel*0.5*dt + @a0*0.125*dt*dt ")
     calc(" @a1 = acc(ba) ")
     calc(" @pos = @old_pos + @vel*dt + @a1*0.5*dt*dt ")
     calc(" @a2 = acc(ba) ")
     calc(" @pos = @old_pos + @vel*dt + (@a0+@a1*2)*(1/6.0)*dt*dt ")
     calc(" @vel += (@a0+@a1*4+@a2)*(1/6.0)*dt ")
   end

Bob: Time to check whether the new code does the same thing as all the older ones: 

 |gravity> ruby rknbody7a_driver.rb < figure8.in
 dt = 0.001
 dt_dia = 2.1088
 dt_out = 2.1088
 dt_end = 2.1088
 method = rk4
 at time t = 0, after 0 steps :
   E_kin = 1.21 , E_pot =  -2.5 , E_tot = -1.29
              E_tot - E_init = 0
   (E_tot - E_init) / E_init = -0
 at time t = 2.109, after 2109 steps :
   E_kin = 1.21 , E_pot =  -2.5 , E_tot = -1.29
              E_tot - E_init = -2e-15
   (E_tot - E_init) / E_init = 1.55e-15
 3
   2.1089999999998787e+00
   1.0000000000000000e+00
  -1.6047303546488470e-04 -1.9320664965417420e-04
  -9.3227640249930266e-01 -8.6473492670753516e-01
   1.0000000000000000e+00
   9.7020367429337440e-01 -2.4296620300772800e-01
   4.6595057278750124e-01  4.3244644507801255e-01
   1.0000000000000000e+00
  -9.7004320125790211e-01  2.4315940965738195e-01
   4.6632582971180025e-01  4.3228848162952316e-01
Alice: And so it does. Good! I think we can stick with this version as our work horse, for a while at least.

7.4. More Information

Bob: There is one improvement I would like to add. Remember when we were debugging our code and we were analyzing a single time step, trying to figure out whether the position and velocity were updated correctly?

Alice: Yes, that was a simple situation, in the case of the forward Euler algorithm, and we could analytically calculate the acceleration.

Bob: Of course, in general we can't be so lucky, and I would prefer to have an easy way to ask for an output of the acceleration as well. Ah, come to think of it, we have this pretty print output version, that we wrote long ago, but that never used since. Here it is: 

   def pp                                # pretty print
     print "     N = ", @body.size, "\n"
     print "  time = ", @time, "\n"
     @body.each{|b| b.pp}
   end

All it does is invoke the pretty print version for each particle: 

   def pp               # pretty print
     print to_s
   end

which does nothing else but print the Body instance in its standard string form: 

   def to_s
     "  mass = " + @mass.to_s + "\n" +
     "   pos = " + @pos.join(", ") + "\n" +
     "   vel = " + @vel.join(", ") + "\n"
   end

Well, let us just add an extra line for the acceleration. Since pretty print became pp, a pretty print with extras should be called a ppx: 

   def ppx(body_array)              # pretty print, with extra information (acc)
     STDERR.print to_s + "   acc = " + acc(body_array).join(", ") + "\n"
   end

I will put this version in file rknbody8.rb.

Alice: That makes a lot of sense. You choose the standard error output stream STDERR because you don't want this extra debugging information to be mixed up with the particle output snapshot, which is written on the default standard output stream, and which can be piped into another program.

And of course, you have to give ppx a parameter, namely the array of all the particles in the N-body system, otherwise our particle cannot compute the acceleration that it receives from all other particles. Let's see. This means that ppx on the Body level must be invoked from the Nbody as: 

   def ppx                          # pretty print, with extra information (acc)
     print "     N = ", @body.size, "\n"
     print "  time = ", @time, "\n"
     @body.each{|b| b.ppx(@body)}
   end

7.5. An Initial Snapshot Output

Bob: Indeed. Now what else is there left to do? We don't want to have this extra information all the time, since it would clutter up the output. Let us introduce a special flag, x_flag, a boolean variable that will be set to `true' if the extra output is desired, and set to `false' when we don't need it.

This means that the method that writes the diagnostics will now get x_flag as a second parameter, and a few extra lines at the end, to invoke ppx if the flag is set to be `true': 

   def write_diagnostics(nsteps, x_flag)
     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
     if x_flag
       STDERR.print "  for debugging purposes, here is the internal data ",
                    "representation:\n"
       ppx
     end
   end

Alice: Now that we are adding features, I would like to have the option to echo the initial snapshot, the data file that is read in before any integration step is taken. I know that I can always read that information myself from the input file, but sometimes it is nice to have it right in front of you, together with the other data,

Bob: Yes, especially when are debugging and you want to take a single time step. Okay, let us introduce a second flag init_out. If the value of this flag is `true', we will require an initial output, in the form of a snapshot on the standard output stream. If the value is `false', we skip the initial output.

We can implement the effects of both flags, x_flag and init_out, by passing them as additional parameters to evolve. The modifications to the body of this method are then very minor. We only have to change three lines. First, the two calls to write_diagnostics acquire x_flag as extra parameter, as we have already seen. Second, we add an extra line 

     simple_print if init_out

just before we start the while loop. The new version of the evolve method thus becomes: 

   def evolve(integration_method, dt, dt_dia, dt_out, dt_end, init_out, x_flag)
     nsteps = 0
     e_init
     write_diagnostics(nsteps, x_flag)
     @dt = dt
     t_dia = dt_dia - 0.5*dt
     t_out = dt_out - 0.5*dt
     t_end = dt_end - 0.5*dt
 
     simple_print if init_out                                                 
 
     while @time < t_end
       send(integration_method)                                               
       @time += dt
       nsteps += 1
       if @time >= t_dia
         write_diagnostics(nsteps, x_flag)
         t_dia += dt_dia
       end
       if @time >= t_out
         simple_print
         t_out += dt_out
       end
     end
   end

Alice: An alternative would have been to make both flags into instance variables for the Nbody class.

Bob: Yes, that would perhaps look more tidy, giving us fewer arguments to pass around. Either way, I don't care very much: most class definitions involve a trade off between the number of instance variables and the number of arguments being passed around. I'm happy just to keep them as arguments for now.

7.6. A New Driver

Alice: Let's make the necessary changes in the driver file as well: 

 require "rknbody8.rb"
 
 include Math
 
 dt = 0.001           # time step
 dt_dia = 2.1088      # diagnostics printing interval
 dt_out = 2.1088      # output interval
 dt_end = 2.1088      # duration of the integration
 init_out = false     # initial output requested ?
 x_flag = false       # extra diagnostics requested ?
 ##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",
       "init_out = ", init_out, "\n",
       "x_flag = ", x_flag, "\n",
       "method = ", method, "\n"
 
 nb = Nbody.new
 nb.simple_read
 nb.evolve(method, dt, dt_dia, dt_out, dt_end, init_out, x_flag)

7.7. A Final Test

Bob: All very straightforward indeed: at the end the two extra parameters for evolve, and above the introduction of the two flags in analogy with the other parameters. Let's test it out, to see whether we still get the same results for our figure out orbit: 

 |gravity> ruby rknbody8a_driver.rb < figure8.in
 dt = 0.001
 dt_dia = 2.1088
 dt_out = 2.1088
 dt_end = 2.1088
 init_out = false
 x_flag = false
 method = rk4
 at time t = 0, after 0 steps :
   E_kin = 1.21 , E_pot =  -2.5 , E_tot = -1.29
              E_tot - E_init = 0
   (E_tot - E_init) / E_init = -0
 at time t = 2.109, after 2109 steps :
   E_kin = 1.21 , E_pot =  -2.5 , E_tot = -1.29
              E_tot - E_init = -2e-15
   (E_tot - E_init) / E_init = 1.55e-15
 3
   2.1089999999998787e+00
   1.0000000000000000e+00
  -1.6047303546488470e-04 -1.9320664965417420e-04
  -9.3227640249930266e-01 -8.6473492670753516e-01
   1.0000000000000000e+00
   9.7020367429337440e-01 -2.4296620300772800e-01
   4.6595057278750124e-01  4.3244644507801255e-01
   1.0000000000000000e+00
  -9.7004320125790211e-01  2.4315940965738195e-01
   4.6632582971180025e-01  4.3228848162952316e-01
Alice: All is well, clearly. Good! So now we have a version of the code that is both lean and easy to read on the level of the integrators, and has extra debugging options. Progress!
Previous ToC Up Next