Our friends decided to write their leapfrog scheme using the second way of implementing the algorithm, centering all physical quantities on integer time step times.

The first leapfrog code that Bob wrote, `leapfrog1.C`, is similar
to the last forward Euler code that they produced, at the end of the
previous chapter, `forward_euler3.C`. There are really only two
differences. One is a computation of the acceleration before we get
into the integration loop, and the other is the difference in the
integration steps themselves. Here is the code:

Let us start with the integration loop `for (double t = 0; ...`.
A direction implementation of Eqs. 4.4, 4.5
would have give for the third and fourth line of the loop:

r[k] += v[k] * dt + 0.5 * a[k] * dt * dt; v[k] += 0.5 * a[k] * dt;

thus updating the position completely, and leaving the second half of the velocity increment until after the new acceleration is calculated, based upon the new position . However, they realized that switching the order of updating position and velocity would simplify the above two lines to:

v[k] += 0.5 * a[k] * dt; r[k] += v[k] * dt;

since the last term in the position update has now been taken care off already by the velocity update. In professional codes you will often find such short cuts, sometimes explained by a comment, more often not (we will address the question of code comments in the next part of this book).

The other change with respect to `forward_euler3.C` stems from
the fact that we need to know the acceleration at the beginning of the
loop, in order to let the velocity begin to step forward, *and* we
need to calculate the acceleration at the end of the loop, in order to
complete the velocity step. The simplest solution would be to
calculate the acceleration twice, at the beginning and at the end of
the loop. However, this would be a terrible waste of computer time.
Calculating accelerations is the most expensive part of an N-body code,
in fact it is the one place where all the computational muscle is
needed, since it scales with the square of the number of particles.
Therefore, it is important to see whether we can avoid such excess.

The solution is simple: rather than calculating the acceleration at
the beginning of the loop, as we did in `forward_euler3.C`, we
calculate it only at the end of the loop. In that way, the next round
through the loop will conveniently start with the old value of the
acceleration, while the new value will become available at the end of
that round.

The only problem, and a potential source for run time bugs, is that the
loop is no longer complete in and by itself. Unlike the case for the
`forward_euler3.C` code, we cannot simply enter the loop at the
beginning, since the accelerations `a[k]` would be undefined, and
have arbitrary values. This is the reason for the three lines
preceding the integration loop, where we calculate the initial
acceleration between the two particles:

double r2 = r[0]*r[0] + r[1]*r[1] + r[2]*r[2]; for (int k = 0; k < 3; k++) a[k] = - r[k] / (r2 * sqrt(r2));

Just for curiosity, and also to show what type of nonsense can be
produced so easily by forgetting to `spin the wheels' once before
getting into an integration loop, let us take out these three lines,
and call the bug version `leapfrog1a.C`. Notice yet another
common source of bugs. Deleting or commenting out the three lines
above will give an error message from the compiler, with a complaint
along the lines of `leapfrog1a.C:34: `r2' undeclared`. Indeed,
after deleting the three lines, we have to put the declaration of the
variable `r2` back into the loop, something that was not necessary
in `leapfrog1.C`, where `r2` had already been declared.

Here are the first few output steps, on one particular computer
(with a bug present, all bets are off with respect to reproducibility
on another computer; and even on the same computer we can get
different results at different times, if our program addresses memory
that has random values in it). Rather than redirecting the output to
an output file, we will use the `head` command, the counterpart of
the `tail` command. Typing `head -5` will give us only the
first five lines of output.

|gravity> g++ -o leapfrog1a leapfrog1a.C |gravity> leapfrog1a | head -5 Please provide a value for the time step 0.01 Initial total energy E_in = -0.875 1.00012 0.0102179 0.00021032 -0.00412328 0.510817 0.0105144 1.00003 0.0153255 0.000315454 -0.0141193 0.510689 0.0105118 0.999835 0.0204317 0.000420556 -0.0241158 0.510511 0.0105081 0.999544 0.0255357 0.000525616 -0.034114 0.510281 0.0105034 0.999153 0.0306373 0.000630624 -0.0441151 0.51 0.0104976 |gravity>

The problem is that the output does not look so unreasonable. If we would have gotten values like or larger somewhere, we would immediately know that we had made a blunder somewhere. But one dead giveaway is the fact that the values of position and velocity along the axis (the third component in three dimensions) are not zero, as they should be. Since we started moving in the plane, with no force components acting perpendicular to the plane, it is clear that should hold for all times.

The moral of the story is that it is always a good idea to inspect a
few lines of output visually - simply by plotting the
values of the position we might not have noticed that there was
anything wrong! What must have happened is that , represented by
the third component `a[2]` in the acceleration array, was non-zero
at the start of the program. And since we did not calculate it before
using, that value first caused the velocity the deviate out of the
plane, and subsequently also the position.

Here are the correct values, different also for the and components of position and velocity, but not dramatically so:

|gravity> g++ -o leapfrog1 leapfrog1.C |gravity> leapfrog1 | head -5 Please provide a value for the time step 0.01 Initial total energy E_in = -0.875 0.9998 0.0099995 0 -0.0200019 0.4999 0 0.99955 0.014998 0 -0.0300059 0.499775 0 0.9992 0.019995 0 -0.0400138 0.4996 0 0.99875 0.02499 0 -0.0500266 0.499374 0 0.998199 0.0299825 0 -0.0600457 0.499098 0 |gravity>

平成18年10月10日