./ch05.html

ToC

5.1. Equations of Motion

Alice: That was great, to be able to do such rapid prototyping in a language we hardly knew. I can see the advantages of an interpreted, dynamically typed computer language. If we had tried to do this in C++, it would have taken quite a bit more time, and we would have had to write many more lines of code.

Bob: Yes. Defining a class, getting it to behave, providing I/O, and chaining it together by piping data from one program to another, all that is a nontrivial beginning. This is encouraging! Let's move on, to see how much we have to add before we can let the integrator integrate.

Alice: Not much, I think, at least if we start with the simplest integrator possible, forward Euler, for the simplest N-body system, namely to solve the 2-body problem.

Bob: The other day I was asked by a friend outside astronomy why we call it the N-body problem. What was the problem with those bodies, she asked. I had to think for a moment, because we are just in the habit to talk about it that way.

Alice: We talk about solving the equations of motion; normally when you solve something, you solve a problem. I guess the terminology comes from seeing a differential equation as posing a problem, that you then solve, either analytically or numerically. Newton's equations of motion for N bodies form a system of N differential equations, and the challenge to obtain solutions is called the N-body problem.

I must admit, I don't like the term "N-body" at all. Why N, and not some other symbol? I would much rather use the term "many body", as they do in solid state physics, say. Don't you agree that the "gravitational many-body problem" sounds much more elegant than the "gravitational N-body problem"?

Bob: I really don't care, as long as I can solve it. Okay, let's start with the 2-body problem, which is really a 1-body problem of course, since you only have to solve the relative motion between the two particles.

Alice: You say "of course", but if this is going to be a presentation for students, we'd better be more explicit. We should start with Newton's equations of motion for a gravitational many-body system:

Here and are the mass and position vector of particle , and is the gravitational constant. To bring out the inverse square nature of gravity, we can define , with , and unit vector . The gravitational acceleration on particle then becomes:

5.2. Relative Motion

For a 2-body system, everything simplifies a lot. Instead of dealing with position for the first particle and for the second particle, we can write the above system of equations of motion as a single equation instead, in terms of the relative position, defined as:

This can be visualized as a vector pointing from particle 1 to particle 2, in other words with its head at the position of particle 2, and its tail at the position of particle 1.

Introducing , we then get:

We can choose physical units for mass, length, and time in such a way that . With the shorter notation , we then get

You often see an even more abbreviated `dot' notion for time derivatives. If we place a dot on top of a variable for each time derivative we take, we wind up with

Bob: End of lecture. Yes, let's start solving that last equation. Here is how the forward Euler integration scheme works:

for the position , the velocity , and the acceleration of our single particle, that describe aspects of the relative motion between two particles. The index indicates the values for time and for the time after one more time step has been taken: .

Of course, we have to provide intial conditions for the position and velocity vectors, before we can solve the equations, and we will do that as soon as we have a working integrator. For now, all we have to do is to code up these equations.

Alice: That was an even shorter lecture. Yes, let's do it.

5.3. Modularity

Bob: Shall we type the code in the same file test.rb where we put the Body class? We can still use the Body format for our `relative' particle, as long as we remember that the mass of that particle corresponds to the sum of the masses of the original particle.

Alice: We can certainly use the Body class, but I suggest that we put the definition of the Body class in a file body.rb, and only the actual integrator in our file test.rb.

Bob: A modular approach, I take it?

Alice: Sure, whenever I can get away with it!

Bob: It may not be a bad idea, in this case. And to show you how modular I can be, I will even call my first attempt body0.rb. When we create other versions, we can give them higher numbers, and when we are really happy, we can call the final version body.rb.

Alice: That sort-of goes in a modular direction. It all depends on what you do with that stack of versions. But I appreciate your attempt!

Bob: Coming back to your suggestion, yes, we can put the Body class definition in a separate file. This is similar to what you do in C and C++ with an include file -- only here it is easier to do so.

Alice I remember well those constructs in C where you had to write things like #ifndef this and #ifndef that before you could be sure that it was safe to include a file without including it more than once? You're saying that Ruby makes life easier?

Bob: Yes, in Ruby you can use a construct called require "filename": it only includes the file if it hasn't been included yet, directly or indirectly.

Here is the file body0.rb:

 class Body
 
   attr_accessor :mass, :pos, :vel
 
   def initialize(mass = 0, pos = [0,0,0], vel = [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

It is just as we left it, but without my one-liner I/O hacks. Now give me some time to figure out how to implement the forward Euler idea . . . .

5.4. The First Integrator

. . . Here it is, the new version of test.rb. As you can see, it starts with requiring that body0.rb gets included at the top.

 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

Alice: Short and simple, but it looks simple only because we are used to coding up integrators. For the students, it would be good to spell out, just one time, what the equations above look like, in the simplest case where we are working with vectors in two dimensions. To be specific in our notation, we can use subscript x for the first component of each vector and subscript y for the second component, as follows:

The last three equations in the earlier sections then become six equations, one for each vector component.

There are two equations for the calculation of the relative acceleration components and , each of which can be calculated in terms of the relative position components and , since :

Two equations tell us how to find the position vector components at the new time step: we take the components at the previous time, and we add the increment in position provided by the product of the corresponding velocity vector component and the time step size:

Similarly, we update the velocity vector components, by incrementing those using the acceleration, which is the time derivative of the velocity:

If we would do this calculation in three-dimensional space, all vectors would acquire a third component: , and so on. Each pair of equations above would then be replaced by a triple of equations, because an extra equation would appear for the corresponding component.

To conclude the story, the forward Euler algorithm for the two-body problem can be summarized as follows. Given the relative position and the relative velocity between the two particles at a given time , first calculate the relative acceleration in terms of the relative position at . Then compute the new values for the relative positions at , in terms of the relative positions and velocities at . Similarly, compute the new values for the relative velocities at , in terms of the relative velocities and accelerations at .

Bob: And after you have once said that in words, it becomes obvious how much more compact and efficient equations are in component notation. And then it is clear that vector notation is even more efficient than component notation.

5.5. Writing Clean Code

Alice: But now I want to understand how exactly you coded these six equations in the file test.rb above.

Bob: The line include Math tells Ruby to make the Math module visible to the program. This is an example of what is called a mixin.

However, this is merely a convenience. We could have left that out, but then the square root sqrt would not have been recognized; we would have had to write Math.sqrt instead. It is just more convenient to mixin Math right at the beginning.

Alice: Then you define the size of the time step and the number of steps you want to take. It is a good idea, to introduce them as variables at the top of your program, instead of hard coding them below inside the loops. It will be much easier to change the values only once at the top, rather than having to inspect the whole program to see where the values occur.

Bob: You would be surprised how many scientists of the hard-coding type you can still find!

Alice: That is because most scientists have never been taught how to include levels of data abstraction in their program, right from the beginning in the conceptual view of the architecture of a program. One of the main principles of data abstraction is: no constants except 0 and 1 are allowed in the body of the code. Everything else should be given a symbolic name.

Bob: I'm in principle against principles, but I like the content of what you just said. As for those many scientists, at least they all use subroutines.

Alice: Not all! There are still legacy codes being used widely that jump around through many pages of programs using do loops between blocks of code that are effectively a poor man's subroutine. One problem is that scientists view subroutines merely as a convenient way to save time: if you use the same piece of code in three places in your program, you save time by writing it only once as a subroutine, and calling it from three places. But that is the least important reason to use functions or subroutines in a program.

A much more important reason is the issue of code maintenance. If the same piece of code occurs in more than one place, it becomes essentially impossible to maintain the code in a consistent way. Change something in one place in someone's legacy code, and most likely you don't even know that it would have to be changed in the other place to. Bugs appear, for no apparent reason, while you are sure that you did the right thing; you just didn't know about the other place which has now become incompatible. I bet that literally tens of person-millennia have gone down the drain this way, during the last fifty years, chasing those bugs.

Bob: You mean more than 10,000 person-years?

Alice: I wouldn't be surprised. I would estimate there to be at least a hundred thousand programmers in the world. Most of them have struggled for a total of a month in their life, at the very least, chasing bugs that originated in codes which they tried to update, only to find out that there were unexpected and typically undocumented side effects. That already makes a person-millennium worth of effort. And I think that is a vast underestimate.

Bob: Impressive. I never thought about it quite that way.

Alice: So the name of the game is modularity, and I'm glad to see you giving the students the right example, in the third and fourth line of your integrator.

Bob: This type of what you call modularity at least I agree with. It would have been easy to write 100.times in the seventh line instead of ns.times, and change that number by hand later. I thought about doing that. But when I saw that dt occurs not once but twice in the body of the do loop, I realized that it would be all to easy to change the numerical value of one of them, and not the other.

Alice: With the result of having the position and velocity stepping forward in time at different speeds -- quite a nightmare, when you want to debug it. Most likely you will start your debugging on the assumption that you made a mistake in the physics of gravitational interaction, or in the mathematical equations, or in the way you solve them numerically, or just a typo in the way you coded it all. The last thing you suspect would be that you would update the position with time step 0.01 and the velocity with time step 0.001, say.

Bob: I know it all too well, from past experience. That's why I've wizened up.

5.6. Where the Work is Done

Alice: Talking about the do loop, I presume it is being traversed ns times, almost as if you read Ruby like English. How can that possibly work?

Bob: The period means that times is a method associated with the class of which the variable ns points to an instance. Since ns has been initialized with an integer, it now has become an instance of class Integer. And conveniently, this class has a method times built in that takes the value of the instance of the class, here the value of ns, and iterates the following block of code ns times.

As often is the case in Ruby, the explanation of a construction sounds far more complicated than the way you use it. As you already remarked, ns.times reads like English. Another example of the principle of least surprise.

Alice: Within the do loop, you first compute by introducing a variable r2, initialize it to zero, and then accumulating the result of squaring the value of each component of the vector . What you need for the acceleration is the 3/2 power, so you compute that in the next line.

Finally, you solve the three pairs of equations I wrote above. The array method each_index presumably does what it says, it executes the next block of code once for each possible value of the index of the array?

Bob: Yes. In the case of a two-dimensional array a, the two components are a[0] and a[1]; remember that Ruby arrays start at zero by default. For such an array, the block of code following a.each_index will be executed once for k=0, and once for k=1.

Alice: And if we would have used three-dimensional arrays for positions and velocities and accelerations, the block would be executed also for k=2. How neat to see an integrator without any need to remind the computer ad nauseam to go through an inner loop for k is 1 to 3, or something like that -- just as we saw earlier when we wrote the I/O routines.

Bob: But we can do even better. Granted, we have avoided the use of an explicit variable NDIM for the number of dimensions, but there is still the lingering smell of components hanging in the air, in terms of the use of k in the last two blocks of code. I wanted to get the code working quickly, so I didn't think about it to much, but I have an idea of how to get rid of even k.