4.2 A Simple Leapfrog for the 2-Body Problem

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 $a_{i+1}$ is calculated, based upon the new position $r_{i+1}$. 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 $10^{10}$ 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 $z$ axis (the third component in three dimensions) are not zero, as they should be. Since we started moving in the $\{x, y\}$ plane, with no force components acting perpendicular to the plane, it is clear that $r_z = v_z = 0$ 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 $\{x, y\}$ values of the position we might not have noticed that there was anything wrong! What must have happened is that $a_z$, 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 $x$ and $y$ 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>