Previous ToC Up Next

6. Returning to Simplicity

6.1. Extra Body Variables

Bob: Now that we've begun to modify my original N-body code, I would like to write a third version, where I keep local copies of the auxiliary variables, such as old_acc, half_vel, and the like, as instance variables of the Body class. That should make the notation within the integrator methods a lot simpler. I don't like the notion of proliferating instance variables, though.

Alice: I don't think that would be such a bad thing. One could argue that in object-oriented programming, it is appropriate to modify the object if you use it for a different purpose. And especially in Ruby, such a modification will be easy. You can always add a few lines in a new class definition, and as you know, those new lines will be directly added to the existing definitions.

Bob: Yeah, well, I still prefer to be parsimonious. But I'm curious to see how much simpler the integrators will become. The first step will be to add those auxiliary variables to the Body class. In the first version of the code, I had: 

   attr_accessor :mass, :pos, :vel, :nb

In the alternative version, we left out the backward link: 

   attr_accessor :mass, :pos, :vel

And now let me put in the extra variables: 

   attr_accessor :mass, :pos, :vel, :old_pos, :half_vel, :a0, :a1, :a2

Ah, this is the only modification that we have to make to the Body class, so you're right, it is not as bad as I thought. All other changes happen only in the Nbody class!

6.2. Alternatives

Alice: But is this really necessary? I thought that in Ruby you can take any class definition, whether you defined the class yourself or whether you get it from a library, and extend that class by adding something like 

  class Body
    . . .
  end
Bob: Yes, that is true. Indeed, we have used that procedure to define the method to_v for the Array class. However, in our code we have two class definitions directly following each other in a single file. The first one gives the Body class, and the second one gives the Nbody class. I suppose I could have given three definitions: first a standard one for Body, then the particular extension in the form of a second, additional definition for Body as you just indicated, and finally the Nbody class definition. But since everything is bundled in the same file anyway, that seemed to me to be a rather unnecessary maneuver.

Alice: Well, you know I like modularity. And you yourself were just complaining about the fact that a modification in the Nbody class forced you to intrude in the earlier definition of the Body class. At the very least, writing three definitions that way will make it more clear what is going on.

In fact, I would probably prefer to put the standard Body class definition in one file, and the Nbody class definition in another file. In that case, you can add the necessary Body extensions to the second file, together with the Nbody class definition that triggered the extensions.

Bob: I don't like to create a plethora of small files. But I see your point. And come to think of it, there must be other solutions.

Alice: Well, the best thing would be if you could modify the definition of the Body class from within Nbody, but I'm pretty sure that would violate some Ruby rules.

Bob: Actually, I'm not so sure. In Ruby you can do almost anything. The trick is to find out how. Hmmm. I have read something about a Ruby class called Binding, and a method binding, both of which seem to be related to creating the possibility to do just the sort of thing you brought up.

Alice: What do you know. Well, not now, I suggest.

Bob: I agree.

6.3. Forward

Alice: Let's start again with the forward Euler method, listing now all three versions in order: your original one; the one we got by including an explicit parameter in the acc call; and the new version you're writing now with extra variables in the Body class: 

   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 forward(dt)
     old_acc = []
     @body.each_index{|i| old_acc[i] = @body[i].acc(@body)}
     @body.each{|b| b.pos += b.vel*dt}
     @body.each_index{|i| @body[i].vel += old_acc[i]*dt}
   end


   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

6.4. Clean Code

Bob: I must say, it is very gratifying to see how much cleaner this last version looks. Let me rewrite the other three integration methods as well, and call this file rknbody4.rb. Here they all are: 

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


   def rk2(dt)
     @body.each{|b| b.old_pos = b.pos}
     @body.each{|b| b.half_vel = b.vel + b.acc(@body)*0.5*dt}
     @body.each{|b| b.pos += b.vel*0.5*dt}
     @body.each{|b| b.vel += b.acc(@body)*dt}
     @body.each{|b| b.pos = b.old_pos + b.half_vel*dt}
   end


   def rk4(dt)
     @body.each{|b| b.old_pos = b.pos}
     @body.each{|b| b.a0 = b.acc(@body)}
     @body.each{|b| b.pos = b.old_pos + b.vel*0.5*dt + b.a0*0.125*dt*dt}
     @body.each{|b| b.a1 = b.acc(@body)}
     @body.each{|b| b.pos = b.old_pos + b.vel*dt + b.a1*0.5*dt*dt}
     @body.each{|b| b.a2 = b.acc(@body)}
     @body.each{|b| b.pos = b.old_pos + b.vel*dt + (b.a0+b.a1*2)*(1/6.0)*dt*dt}
     @body.each{|b| b.vel += (b.a0+b.a1*4+b.a2)*(1/6.0)*dt}
   end

Alice: A lot easier to read. A great improvement over the previous two versions, though not quite as clean as the two-body version. Let's put up the rk4 method for the old two-body code: 

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

Bob: The difference is that in our latest version we still have to indicate which body b it is that gets the instructions, hence the "b." in front of each pos, vel, etc., call.

6.5. Sending a String

Alice: Well, I have an idea. Perhaps we can get rid of the "b." in front of each physical variable, after all.

Bob: How? I mean, you have to loop over each particle! I can't see how you can get rid of that.

Alice: You can't get rid of that, I agree, but you don't have to repeat the fact that you are looping many times in one line, as we are doing now. I think we can express the idea of a loop just once in each line.

Bob: I still don't see how you can do that. In the old two-body code, we could get away with writing pos += vel*dt because there was only one pos and one vel, but here we have little choice but writing b.pos += b.vel*dt, at a minimum. And I already chose the shortest name I could think of, b for the body whose pos gets updated!

Alice: My idea is to let the Nbody class give a command to the Body class, specifying directly to do pos += vel*dt for particle b, without repeating the b presence separately for pos and vel.

Bob: That would be nice, yes, but it still looks impossible to do, since it would involve something sending that whole line to the Body class!

Alice: Exactly!

Bob: I beg your pardon?

Alice: You got it! Let us send that line to the Body class! We can ask the integration methods in Nbody to specify what needs to be done, but instead of making the actual command calls, these methods can write the commands into a string, and then pass that string down to the Body class.

Bob: Hey, that is a great idea! I would never have thought about that. Can you really do that? Well, of course you can. I guess I'm still thinking too much in Fortran and C terms. But I must say, I'm beginning to get really worried about speed; this may slow down execution of the code even more. And Ruby is already slow to begin with.

Alice: Let us postpone efficiency discussions for later. I'm pretty hopeful that we can find ways to speed things up later. In the worst case we can switch back to a more traditional form in our final production code. If we can make the prototyping process easier and more transparent, that would already be a big gain. I care less about elegance of the final version than I do about the intermediate test versions.

Bob: Well, okay, but we shouldn't wait too long in checking the speed. I would feel more comfortable if we had some hard-nosed proof that a Ruby-based N-body code can really compete with codes written in C or Fortran. Meanwhile, I must admit, I do like the flexibility of Ruby. What a trick, to implement message passing between the Nbody system scheduler and the individual bodies, literally by passing messages in the form of a string!

Alice: Yes, Ruby invites a different way of thinking. Even if you could find a way to do such a trick in C++, in some convolved way, the point is that it would not occur to you to do so, since the language does not invite that way of thinking.

Bob: Quite convolved: for one thing, you would need to find a way to implement a form of dynamic loading! C++ is a compiled language, not an interpreted language like Ruby.

Alice: Ah, yes, of course, I'm already beginning to forget what it means to work with compile cycles in prototyping code!

6.6. Wishful Thinking

Bob: How did you get this idea, of using message passing?

Alice: When I browsed through some Ruby code on the web, I came across something similar as what I just suggested, and then I realized that that would be a very natural way of passing instructions between classes.

Bob: Let's try it! I'll create yet another file, rknbody5.rb, and see how we can implement your idea.

Alice: This time, I suggest we start on the level of the Nbody class. In a wishful thinking sort of way, let us assume that we can tell a Body named b to calculate something, by issuing a command b.calc(s), where s is a string that will get executed by b. The forward Euler method would then look like this: 

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

Bob: And indeed, with no mention of b anymore within the string that is passed as the third argument of calc.

Alice: The first two arguments are needed, because otherwise the Body class will not know what to do with the string: it knows about its own instance variables such as @pos, and it knows about the method acc, but not about the argument to acc, which has to be specified explicitly.

Bob: I see. What does ba stand for?

Alice: Body array. When calc executes the call, it will replace ba by its first argument, @body. Similarly, it will replace dt by dt. I'll show you in a moment. I hope it will all work. And I'm pretty sure it will.

6.7. Implementation

Bob: You used the calc method in the middle line of forward as well, even though you could have used the simpler statement 

     @body.each{|b| b.pos += b.vel*dt}

which we used before, in the previous version.

Alice: Yes, but my intention was to not break the symmetry between the lines. By treating all of them in the same way, your eye can be guided to what is different on each line, forgetting the left half of each line, which is the same in all cases. Here, let me write the other three methods as well, and then it will become more clear how the actual strings that contain the commands will stand out: 

   def leapfrog(dt)
     @body.each{|b| b.calc(@body,dt," @vel += acc(ba)*0.5*dt ")}
     @body.each{|b| b.calc(@body,dt," @pos += @vel*dt ")}
     @body.each{|b| b.calc(@body,dt," @vel += acc(ba)*0.5*dt ")}
   end


   def rk2(dt)
     @body.each{|b| b.calc(@body,dt," @old_pos = @pos ")}
     @body.each{|b| b.calc(@body,dt," @half_vel = @vel + acc(ba)*0.5*dt ")}
     @body.each{|b| b.calc(@body,dt," @pos += @vel*0.5*dt ")}
     @body.each{|b| b.calc(@body,dt," @vel += acc(ba)*dt ")}
     @body.each{|b| b.calc(@body,dt," @pos = @old_pos + @half_vel*dt ")}
   end


   def rk4(dt)
     @body.each{|b| b.calc(@body,dt," @old_pos = @pos ")}
     @body.each{|b| b.calc(@body,dt," @a0 = acc(ba) ")}
     @body.each{|b| b.calc(@body,dt," @pos = @old_pos + 
                                             @vel*0.5*dt + @a0*0.125*dt*dt ")}
     @body.each{|b| b.calc(@body,dt," @a1 = acc(ba) ")}
     @body.each{|b| b.calc(@body,dt," @pos = @old_pos + 
                                             @vel*dt + @a1*0.5*dt*dt ")}
     @body.each{|b| b.calc(@body,dt," @a2 = acc(ba) ")}
     @body.each{|b| b.calc(@body,dt," @pos = @old_pos +
                                         @vel*dt + (@a0+@a1*2)*(1/6.0)*dt*dt ")}
     @body.each{|b| b.calc(@body,dt," @vel += (@a0+@a1*4+@a2)*(1/6.0)*dt ")}
   end

Bob: It's an improvement over the first two versions, which used an index i. In the previous version, we could leave out that obnoxious i index at the expense of introducing extra variables on the Body level. But now we can have our cake and eat it: no more i and no more extra variables with Body, I presume.

Alice: Correct! The content of the string effectively will declare the extra variables for us when the string is evaluated in the Body class. Here is how I would write the calc method for the Body class: 

   def calc(body_array, time_step, s)
     ba  = body_array
     dt = time_step
     eval(s)
   end

Bob: Simplicity itself. I see now what you meant, when you described the way the two parameters ba and dt were going to be substituted in an actual call to calc.

Alice: Almost too simple to be true, but I think this is all correct.

6.8. Indirect String Sending

Bob: Shall we move on? Or do you have a suggestion for further improvements?

Alice: I must say that the integration methods still look too cluttered for my taste. They miss the simple elegance and brevity of expression of our previous version. For one thing, in that version we did not have to break any statement up over two lines. What bothers me especially is that for most statements, more than half of the line gets repeated exactly.

Bob: I wonder whether we can do something about that.

Alice: I think we can. So far, we have introduced a calc function on the Body level. How about introducing a second calc function on the Nbody level? Let's create one more file, rknbody6.rb, in which we give the Nbody class the following extra method: 

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

Bob: Brevity indeed. I see what you mean. Forward Euler then becomes, instead of 

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

which we just wrote, quite a bit shorter as: 

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

Alice: Exactly. And the following three methods become: 

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


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


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

6.9. The Same, Yet Different

Bob: I like that: every statement now fits on one line, we don't have to modify or extend the Body class when we introduce a new integration scheme, and we don't mention the variable @body anymore in each line. I'm satisfied: this combines all good things in one. And I must say, I'm growing fond of working with Ruby.

Alice: So do I. Even though I knew that Ruby was a well designed language, I was a bit skeptical at first about how much that would really buy us. But as we already have seen, it buys us quite a bit in terms of clarity of expression. And when writing complicated programs and packages, something we will start doing soon, clarity of expression is more important than anything else. Nothing else will allow you to maintain an overview over the whole situation.

However, just one point: you mentioned that the Body class did not need to be extended. But that is not true. We have extended the Body class by adding a calc method to i.

Bob: I meant that we don't have to make a different and separate extension to the Body class, each time we introduce a new integrator in Nbody. We only modify Body once, by adding calc, and the same calc will do different things for different integrators.

Alice: Which means that in practice, dynamically, we are augmenting the Body class, each time we add an integrator to the Nbody class.

Bob: That is true, but it is invisible as far as the code is concern. If the Body class would be written in a separate file, that file would not have to be changed, upon adding an integrator.

Alice: But I thought you didn't like to cut up our file into, what did you call it again, ah yes, a plethora of small files.

Bob: Very funny. I'm not really suggesting to split up the file, I just tried to make a clear point even clearer. But of course we are both right: I am right to say that the written definition of the Body class remains unchanged, and you are right if you say that the actual definition, as it is dynamically changed during execution of the code by the interpreted, does change.

Alice: Right you are! And right I am. Okay, let's move on.

6.10. Testing

Bob: With all this bickering, we haven't tested yet any of our latest three versions. For now, let's just try the figure-8 triple. I'll run the first singly-linked code version again, to make sure we got the right output. 

 |gravity> ruby rknbody3b_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
And then our version with the extra variables added by hand to the Body class: 

 |gravity> ruby rknbody4a_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
Good! Now the version that is sending a string from Nbody to Body: 

 |gravity> ruby rknbody5a_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
Also perfect. Finally our last version with the two calc methods: 

 |gravity> ruby rknbody6a_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
Great! It all works.

Alice: That's very nice. As we already said, almost too good to believe.
Previous ToC Up Next