Previous ToC Up Next

8. Formula Translating

8.1. A Body with Vectors

Bob: Okay, here is the next step. We first have to tell our Body class that the positions and velocities are no longer arrays, but vectors. Let me create a new file vbody1.rb for this modified class. I will call the modified class Body as well. If we are happy with it, we can discard the previous Body class definition.

 require "vector3.rb"
 
 class Body
 
   attr_accessor :mass, :pos, :vel
 
   def initialize(mass = 0, pos = Vector[0,0,0], vel = Vector[0,0,0])
     @mass, @pos, @vel = mass, pos, vel
   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}
     @vel = gets.split.map{|x| x.to_f}
   end
 
 end

Here it is. I have added one line and modified one line, that's all! The line on top require "vector3.rb" I have added so as to make Ruby aware of our new class Vector; as we have seen, the Ruby interpreter effectively replaces this line by the file vector3.rb. Now the only thing I changed is the initialization line for pos and vel. They are now bona fide vectors, right when they see the light of day.

The next step is to adapt our integrator, so that it can handle particles with vectors for their positions and velocities. Remember, that was the whole reason for the exercise: to simplify the notation in our forward Euler integrator. You were unhappy with the |k| notation, and that started us off on this long trek.

Alice: But we learned a lot on our journey. Yes, I remember now. I was complaining about component notation, and then you decided you can give me real vector notation, without components. Seems like a while ago!

8.2. Shrinking Code

Bob: But now we're getting close. Let me rename the old file euler1.rb, and then add vectors to euler2.rb. No, let me be more careful. Let me call it euler2a.rb, for alpha version of euler2.rb. If and when it works, we can rename it euler2.rb.

Alice: You're getting cautious! Can you remind me what euler1.rb looked like?

Bob: Here it is, or was; or will be -- but not for much longer:

 require "body0.rb"
 
 include Math
 
 dt = 0.01            # time step size
 ns = 100             # number of time steps
 
 b = Body.new
 b.simple_read
 
 ns.times do
   r2 = 0
   b.pos.each {|p| r2 += p*p}
   r3 = r2 * sqrt(r2)
   acc = b.pos.map { |x| -b.mass * x/r3 }
   b.pos.each_index { |k| b.pos[k] += b.vel[k] * dt }
   b.vel.each_index { |k| b.vel[k] += acc[k] * dt }
 end
 
 b.pp

All we have to do is to change the first require line to require vbody1.rb to be loaded instead of body.rb, and then we are free to purge the do loop from the obnoxious k component notation. Are you ready for this?

Alice: Sure thing!

Bob: Here it is. One line shorter in total, with most of the lines shorter in length as well. I've put it into euler2a.rb for now:

 require "vbody1.rb"                                                          #1
 
 include Math
 
 dt = 0.01            # time step size
 ns = 100             # number of time steps
 
 b = Body.new
 b.simple_read
 
 ns.times do
   r2 = b.pos * b.pos
   r3 = r2 * sqrt(r2)
   acc = b.pos * (-b.mass / r3)
   b.pos += b.vel * dt
   b.vel += acc * dt
 end
 
 b.pp

8.3. FORTRAN

Alice: One of these days you'll shorten it so much that you can hand me a one-line integrator, just like you produced one-line read and write functions, earlier!

Bob: Don't count on it; I think a five-line integrator is good enough already. What a difference with euler1.rb! Instead of asking b.pos to add the square of each of its components to r2, after setting r2 to zero, we simply order r2 to be the inner product of b.pos with itself, in the first line of the do loop.

The next line is unchanged. But then the fun really starts: the acceleration vector is now directly given as the product of the velocity vector with the terms mass over distance cubed.

Alice: Just as we saw it in formula form earlier:

Bob: So you see, now we can translate our formulas directly into computer code!

Alice: Perhaps Ruby should be called FORmula TRANslator.

Bob: I'm afraid that name has already been taken, quite a few decades ago . . .

Alice: And the next two lines clearly mean:

also just like we wrote it before. What a relief, to write a computer code as if you write equations directly down into the software!

Bob: Indeed. I've seen C++ doing that, what they call overloading operators, and you can get a similar notation, but with far more work.

Alice: And with a complex class definition and declaration structure that is very hard to play with and to change at will. In contrast, here with Ruby we can really do what they always advertise as rapid prototyping.

Bob: But before congratulating ourselves too much, let us first see whether this now runs. I've grown a bit more cautious after those surprises we got earlier. Ruby is certainly powerful, but as we have seen, we'd better know what we are doing, in order to use all that power wisely.

Alice: Well said, Bob! Let's run the new version of the integrator with the same values we had before, which I see you have still there in the file euler2a.rb: a hundred time steps of 0.01 time units each.

Bob: Here goes:

 |gravity> ruby euler2a.rb < euler.in
 euler2a.rb:12:in `*': can't convert Array into Integer (TypeError)
        from euler2a.rb:12
        from euler2a.rb:11:in `times'
        from euler2a.rb:11

8.4. An Old Friend

Alice: Where did we see that error message before?

Bob: I know, it is an old friend, isn't it? I had hoped we had banned it forever. But at least we have an idea now where it could come from. I bet that there is still a vector somewhere that somehow thinks it is an array . . .

Alice: Which line is line 12 in euler2a.rb

Bob: it's the very first line in our do loop. Just a few lines above, b is freshly minted, so perhaps b.simple_read is the bad influence that turns some vector into an array? Could be. But before jumping to conclusions, let me sprinkle in a few class statements.

Alice: You are getting cautious!

Bob: Yeah, but I don't want to waste more time, this time around. Let me do it in irb directly, so that we can get everything echoed.

    |gravity> irb --prompt short_prompt -r vbody1.rb
    001:0> include Math
    Object
    002:0> dt = 0.01
    0.01
    003:0> ns = 100
    100
    004:0> b = Body.new
    #<Body:0x400cda40 @vel=[0, 0, 0], @pos=[0, 0, 0], @mass=0>
    005:0> b.class
    Body
    006:0> b.pos.class
    Vector
    007:0> b.simple_read
And now it hangs. Hello, computer!? What are you waiting for?

Alice: You haven't given it the euler.in input file.

Bob: Ah, of course. I can see that there are drawbacks of using an interactive interpreter if you're testing something that has dependencies with other things. Well, the file is short enough to type in by hand. Here goes:

    007:0> b.simple_read
    1
    1 0
    0 0.5
    [0.0, 0.5]
    008:0> b.pos.class
    Array
Alice: Our old friend, the vector-turned-array bug! So you were right, after all, suspecting that simple_read was the culprit.

Bob: What do you mean `after all'? I was almost certain, but now I know for sure.

8.5. Magic

Alice: Right you are! Now what is wrong with simple_read?

Bob: Let's inspect. Here is the definition:

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

And, yes, it is obvious! @pos receives its value from a gets.split.map operation, in other words first gets is invoked to read a line, then the string method split splits the string that gets just returned, in the form of an array, and finally map does something on all the elements of that array, returning a new array with the results. And that last part is where the trouble arises: map returns an array that is assigned to @pos.

Alice: And poor @pos, Ruby chameleon that it is, instantly degrades into an array, forgetting its whole rich vectorness. Of course. Well, when you get some more experience, debugging becomes a breeze, doesn't it?

Bob: Ah, but then you start writing more complicated programs and then you run into more subtle bugs. It's like an arms race.

Alice: We'll see. For now I'm just glad that we tracked this one down quickly. So what do we do about it?

Bob: Simple. We just turn the right-hand side of the offending statement into a vector, rather than an array.

Alice: But how do you do that?

Bob: Ah, when I browsed in the manual, I saw some neat trick. Just a moment. Here it is. Are you ready for this? This is good job for irb. Neat and simple.

    |gravity> irb --prompt short_prompt -r vector3.rb
    001:0> a = [1, 2, 3]
    [1, 2, 3]
    002:0> a.class
    Array
    003:0> b = Vector[*a]
    [1, 2, 3]
    004:0> b.class
    Vector
Alice: Huh? How did that work?? What type of magic did you just invoke, putting a star in front of an array, then wrapping it in square brackets, and pulling a Vector spell on it?

Bob: Another convenient Ruby trick. It takes a while to learn all these tricks, but they sure are useful. Here is what happened. The * notation in front of an array effectively dissolves the outer [] brackets, turning an array into a comma separated list of its elements.

Alice: Ah, and then you feed that list into the initializer of a new Vector object. Very clever, if not devious!

Bob: So I'll just clothe the right hand side of the @pos and @vel assignments in Vector[* ] garb. Let me call this file euler2b.rb with b for beta version. The only difference is that in euler2a.rb I started at the top with

 require "vbody1.rb"                                                          

while in euler2b.rb I started at the top with

 require "vbody2.rb"                                                          

instead. I think we're getting closer now. Here it is, in vbody2.rb:

   def simple_read
     @mass = gets.to_f
     @pos = Vector[*gets.split.map{|x| x.to_f}]
     @vel = Vector[*gets.split.map{|x| x.to_f}]
   end

And let me try it:

 |gravity> ruby euler2b.rb < euler.in
   mass = 1.0
    pos = 0.443229564341409, 0.384215558686636
    vel = -1.29436531289392, 0.0258120006669036
Alice: Congratulations!

8.6. A Matter of Taste

Bob: Thanks. So now we have a FORmula TRANslation device that actually works!

Alice: But it is ugly.

Bob: I beg your pardon?

Alice: It is ugly.

Bob: How so?

Alice: I don't like that [*...] notation. Looks like a hack.

Bob: I like hacks.

Alice: I know you do. But can't we come up with something cleaner?

Bob: Like what?

Alice: Well, we have seen that to_s will happily turn something into a string; we even gave our Body class its own to_s method, so that it could show its face in public. I bet there is a method called to_a that turns objects into arrays.

Bob: There is. And it does just that.

Alice: Can't we then write our own to_v method which turns something into a vector?

Bob: Now that's a thought. Sure. That would be fun! And I just know how to do that! I only have to add a few lines to the vector3.rb file. I think this is really what we want. To celebrate, let me call the file vector.rb now.

Here, it is short, let me show the whole file:

 class Vector < Array
   def +(a)
     sum = Vector.new
     self.each_index{|k| sum[k] = self[k]+a[k]}
     sum
   end
   def *(a)
     if a.class == Vector              # inner product
       product = 0
       self.each_index{|k| product += self[k]*a[k]}
     else
       product = Vector.new           # scalar product
       self.each_index{|k| product[k] = self[k]*a}
     end
     product
   end
 end
 
 class Array
   def to_v
     Vector[*self]
   end
 end

8.7. Voila

Alice: That's much better. I don't mind the [*...] per se, it is only that the line with gets.split.map followed by a block of code just got too crowded for comfort. But your definition of to_v is clean and simple. So let me check: by writing class Array and so on, you are just extending the Ruby-defined standard array definition, right?

Bob: Right. Not only can you extend your own classes by giving new bits and pieces like this, as we have seen before, you can do the same thing to the predefined classes. In general, almost nothing in Ruby is off limits; you can typically treat the built-in stuff as if you had just written it all yourself.

Alice: That's really neat. So now can we get rid of the clutter at the bottom of the vbody2.rb file?

Bob: Sure thing. Here is our simple_read, no clutter, and with the new to_v at the end. And let me call this one vbody.rb, to celebrate that now we know what we are doing:

   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

Alice: Very nice. And presumably it runs?

Bob: I bet you it will! And to show my confidence, I will call the file that includes vbody.rb now euler2.rb, since I think we will be beyond alpha and beta testing once we see this.

Alice: You do sound confident. Let's see.

 |gravity> ruby euler2.rb < euler.in
   mass = 1.0
    pos = 0.443229564341409, 0.384215558686636
    vel = -1.29436531289392, 0.0258120006669036
Bob: Voila!

Alice: Very good.
Previous ToC Up Next