| Previous | ToC | Up | Next |
Bob: Now that we're pretty sure that we can integrate the two-body
problem with our new code, how about trying our hand at a three-body
system? After all, that was the reason I wrote this program, to go
beyond two bodies.
Alice: We will have to choose some initial conditions. Rather than picking
something at random, how about trying to integrate a figure-eight
configuration?
Bob: You mean this new `classical' solution to the equal-mass three-body
system that was discovered recently by Montgomery and Chenciner? I
remember how surprised I was when I first read about that. I had
always assumed that the classic celestial mechanicians had found all
the interesting solutions already centuries ago, you know, Legendre,
Lagrange, Laplace, and perhaps a few others Le's and La's. I was
thrilled to see such an elegant new solution.
Alice: Me too. As soon as I read about it, I tried it out for myself,
and by golly, the configuration was stable: three stars chasing each other
on a figure eight, without the system falling apart. You could even perturb
the initial conditions by a fraction of a percent, and preserve stability.
Bob: This all goes to show that we really should implement some form of
graphics soon. I'd love to see how the new code will handle that
configuration.
Alice: I agree. But at least for now, we can use the figure-eight triple
to test the code. I stored the initial conditions somewhere. Ah,
here they are. I'll put them in a file named figure8.in:
Bob: Let's predict what will happen. The third particle has velocity
components that are negative, both in the x and y direction. This
means that it moves to the left and downward. That would suggest that
it will in due time reach the position of the second particle, the
particle that starts off with a large negative x value. The second
particle then has no choice but to replace the first particle, and the
first one will in turn wind up in the center. So I predict that an
integration of 2.1088 time units will produce the following output:
Alice: That seems reasonable. This will be a nice test. Let's try it!
Alice: I'd say! And you were right about the order of the particles.
The first particle did indeed reach very close to the center, and all
other particles reached your predicted position to within a fraction
of a tenth of a percent: in all cases the first three digits are
correct. Congratulations!
Bob: Thanks! So let's move on.
Alice: Well, before we move on, let us revisit the choice you made of
installing double links between the particles in an N-body system and
their parents. In the Nbody class you have an array of bodies
@body that gives access to the members of the Body class,
while in the Body class there is an instance variable @nb
that gives each particle access to its parent.
Bob: Yes, I remember that you were not happy with that choice.
Alice: The problem is, I can see the danger of modifying, say, the
downward link and forgetting to make a corresponding modification in
the upward link. Let's play with some alternatives.
Bob: Be my guest! I'll keep a copy of the last, debugged, version,
and it is fun to try out some alternative implementations. After all,
Ruby invites this type of prototyping. Here, why don't you take the
key board.
Alice: The simplest way
to do away with the need for a backward pointer @nb would be
to give the Body method acc an extra parameter. Let me take out
the declaration of @nb in the attr_accessor of the
Body class, and let me replace your Body method
by a new version, that takes the array of bodies as an explicit parameter:
Of course, we have to do the same thing for the potential energy calculation,
where we replace
by the equivalent version:
Bob: Okay, and let us put the code into a separate file rknbody3.rb.
Life gets a bit simpler for the Nbody code: when it creates its daughters,
it no longer has to tell them who their parent is. Instead of the old version:
we now will use the simpler and more natural version:
And similarly we can leave out the corresponding line in the
simple_read method: instead of
we now have only:
Alice: Good! I was a bit worried about that repetition. It would have
been an easy mistake to make to modify the initialize method
in one way, and to either forget to modify the simple_read method
or worse, modify it in a different way. Your original approach violated
the DRY principle.
Bob: I wasn't aware of violating anything, let alone a principle I
hadn't even heard of. What is the DRY principle?
Alice: Don't Repeat Yourself. The idea is that you should try to avoid
stating the same information in more than one place, exactly because
it is so easy to update information in one place, and to forget to update
it in one or more other places.
Bob: But in many cases it could never hurt you. In C++, for
example, you normally declare all your functions, and specify the
types of the arguments and the return values, in different places
from where you actually define those functions. Now if you make a
mistake, and provide conflicting information, the compiler will give
you an error message. So there is nothing that can go wrong there.
Alice: Even if there is little danger for errors, it is annoying to
have to provide the same information in two places, especially if you
work with header files, where the information may even reside in one
or more different files. Such requirements go against any notion of
rapid prototyping.
While violating the DRY principle may not always lead to likely
errors, the consequences are equally bad if it discourages free
experimentation. In that way you are likely to miss a more optimal
solution. Miss out on more optimal solutions often enough in a
large software project, and you get an end product that can be very
far from optimal, through the aggregate effect of all the suboptimal
decisions you have made.
Bob: You sound like a manager. I must say that I don't like to
think in terms of principles, and I tend do develop an allergy against
people who follow principles blindly. But I must admit that this
particular idea makes a lot of sense. Having to repeat information in
different places in C++ is certainly one of the things that I like
least about that language. In fact, I've developed an allergy against
that feature already quite a while ago.
Alice: I'm glad to hear that your allergy against C++ is stronger
than your allergy against me, if I parse you correctly.
Bob: So far, certainly. Let's see how things develop.
Alice: Given that C++ is less likely to change in any fundamental
way soon, I guess it's up to us to see how our collaboration develops
further. As for me, I think our different attitudes have been
balanced quite well, so far.
Bob: I agree. And I'm glad you take my recalcitrance with a sense of
humor, since I'm not likely to become `principled' any time soon!
Alice: This discussion about the DRY principle started when I told
you I was glad that the initialize method no longer had to
repeat what the simple_read method also did, namely the
initialization of the upward links from particles to N-body system,
links that you have now removed. However, I still see some repetition
left in the the initialization method, of stuff that is already taken
care of in simple_read. Here is your latest version of
initialize:
If you compare this with the input method simple_read:
you see that the number of particles N and the time are being
used here only. After reading in N, simple_read extends
the array of bodies to contain N elements, and in the process it
reads in the values for each body. And after reading in the time,
it assigns that value directly to @time. It seems to me that
we can simplify the initialize method to just:
Bob: That is a lot shorter, and yes, that would work fine in this case.
The reason that I included the two arguments to initialize is that I
had been thinking about a situation where you may want to construct an
initial model, say a Plummer model or a King model. In that case, you
probably do want to set up an array of N bodies, at a specified
time, before you internally assign values to the masses, positions and
velocities of those bodies.
Alice: Ah, yes, in that case these extra arguments would come in handy,
I agree. But so far you have often made it clear to me that we should
be demand-driven, and that we should not build in extra options for
possible future use, unless we actually plan to use those pretty soon.
Bob: I agree. Yes, it would be fine with me to take your simpler
version, at least for now. We can easily extend initialize back to
the longer form later on, when needed.
Alice: great, four lines saved and one line simplified. What else
is there left to be done?
5. A Single-Links Version
5.1. A Figure-8 Triple
3
0
1
0.9700436 -0.24308753
0.466203685 0.43236573
1
-0.9700436 0.24308753
0.466203685 0.43236573
1
0 0
-0.93240737 -0.86473146
The third particle starts in in the center of the figure eight, in the
exact center of the coordinate system. The other two bodies move
symmetrically with respect to each other, each one third of a period
displaced. Ah, my notes tell me that I found the total period of
revolution to be about 6.3264 time units. Let us integrate our system
for 1/3 of that time: since all particles have the same mass, after 1/3
of the period they should have changed places. That means an
integration for a total of 2.1088 time units.
5.2. Switching Places
3
0
1
0 0
-0.93241 -0.86473
1
0.97004 -0.24309
0.46620 0.43237
1
-0.97004 0.24309
0.46620 0.43237
I have only given five digits, since you gave the total time duration
to that accuracy, and therefore there is no guarantee that we will halt
the calculation exactly after 1/3 of an orbit; in fact we are bound to
either overshoot or undershoot, making an error in the fifth or sixth
significant digit in all coordinates.
|gravity> ruby rknbody2b_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
Bob: Not quite as good as I had hoped, but I certainly came close.
5.3. Single Links in Body
5.4. Single Links in Nbody
def initialize(n=0, time = 0)
@time = time
@body = []
for i in 0...n
@body[i] = Body.new
end
end
5.5. The DRY Principle
5.6. Simplifying Further
def initialize(n=0, time = 0)
@time = time
@body = []
for i in 0...n
@body[i] = Body.new
end
end
| Previous | ToC | Up | Next |