Previous | ToC | Up | Next |

** 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:

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.

** 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 . . . .

. . . 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.

** 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.

** 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`.

** Alice**: That would be even better. But shall we first check whether
everything works as advertised?

Previous | ToC | Up | Next |