6.4 Implementing Hermite

Let us return to Alice, Bob, and Carol, who are about to implement the Hermite scheme. Since Bob was most eager to test out this new algorithm, he got his turn behind the computer.

Bob:

Let's give this new Hermite scheme a stress test. I'll take leapfrog2a.C, the one with the off-set in initial velocity of $10^{-4}$, call it hermite1a.C for now, and simply substitute the leapfrog equations 4.4, 4.5 by the Hermite equations 6.1, 6.2.

Carol:

Why not call the code hermite1.C?

Bob:

Even though the update seems almost trivial, I don't have the audacity to believe we'll get everything right the first time around. When it all works, we will rename the code to hermite1.C. So here is the new version:

Alice:

Ah, I see now that you were wise giving the code an preliminary name. Notice where you have added the jerk calculation. You replaced

            for (int k = 0; k < 3; k++){
                a[i][k] += m * rji[k] / r3;
                a[j][k] -= m * rji[k] / r3;
            }

by:

            for (int k = 0; k < 3; k++){
                a[i][k] += m * rji[k] / r3;
                a[j][k] -= m * rji[k] / r3;
                jk[i][k] += m * (vji[k] - 3 * rv * rji[k]) / r3;
                jk[j][k] -= m * (vji[k] - 3 * rv * rji[k]) / r3;
            }

in addition to some extra changes, such as introducing and calculating the variable vji[] for the relative velocities for particles $i$ and $j$. All of this is fine, as far as I can see. The statements are correctly written and reflect the Hermite, but there is one problem. In the last two lines above you use the relative velocities $vij[]$, but at this point in the program they have not been assigned yet.

Bob:

Ah, you are right. That happens just a few lines below, with this `clever' trick of using the magnitude of the centrifugal acceleration to determine the magnitude of the velocities. Thanks! Okay, so I'll make another version, hermite1b.C, which has two initial loops, one for the accelerations, followed by the assignment of velocities, and then followed by the second loop which computes the jerks.

Carol:

Maybe we can just run the program using our friend /dev/null to see how the energy behaves, before start making pictures. I'm really curious to see whether we have now reached fourth-order accuracy.

Bob:

Easy to test:

|gravity> g++ -o hermite1b hermite1b.C
|gravity> hermite1b > /dev/null
Please provide a value for the time step
0.01
and for the duration of the run
100
Initial total energy E_in = -0.866025
Final total energy E_out = -0.975389
absolute energy error: E_out - E_in = -0.109364
relative energy error: (E_out - E_in) / E_in = 0.126282
|gravity> !!
hermite1b > /dev/null
Please provide a value for the time step
0.001
and for the duration of the run
100
Initial total energy E_in = -0.866025
Final total energy E_out = -0.866586
absolute energy error: E_out - E_in = -0.000560173
relative energy error: (E_out - E_in) / E_in = 0.000646832
|gravity> !!
hermite1b > /dev/null
Please provide a value for the time step
0.0001
and for the duration of the run
100
Initial total energy E_in = -0.866025
Final total energy E_out = -0.866026
absolute energy error: E_out - E_in = -5.52839e-07
relative energy error: (E_out - E_in) / E_in = 6.38364e-07
|gravity> !!
hermite1b > /dev/null
Please provide a value for the time step
0.00001
and for the duration of the run
100
Initial total energy E_in = -0.866025
Final total energy E_out = -0.866025
absolute energy error: E_out - E_in = -5.53711e-10
relative energy error: (E_out - E_in) / E_in = 6.39371e-10
|gravity>

Carol:

What is this? We seem to have constructed a third-order algorithm!

Alice:

Hmmm. It seems that way. Each refinement of a factor ten in step size gives a reduction of the error of a factor 1000. But this was not supposed to happen

Bob:

This is strange indeed. Look, at the end, there they are, the real Hermite expressions, what can be wrong with such simple statements??

        for (int i = 0; i < n; i++){
            for (int k = 0; k < 3; k++){
                r[i][k] = old_r[i][k] + (old_v[i][k] + v[i][k])*dt/2
                                      + (old_a[i][k] - a[i][k])*dt*dt/12;
                v[i][k] = old_v[i][k] + (old_a[i][k] + a[i][k])*dt/2
                                      + (old_j[i][k] - jk[i][k])*dt*dt/12;
            }
        }

Carol:

Yes, they are just what we derived before, as Eqs. 6.1, 6.2.

Alice:

Could we have made a similar mistake as we did just before, when the expressions were absolutely correct, but the order of execution was wrong?

Bob:

Well ...hmmm ...aha, of course, you are right! Look, that is exactly the problem. In the first line, where the position is updated, we indicate that we are using both the old velocity values and the new velocity values. However, the new ones have not been computed yet -- that happens only on the next line! And the solution is obvious: fortunately, the second line updating the velocity does not have this problem, since it uses only accelerations and jerks, all of which have already been computed, both for the old and the new values. So we can simply swap the lines, and this should now work:

        for (int i = 0; i < n; i++){
            for (int k = 0; k < 3; k++){
                v[i][k] = old_v[i][k] + (old_a[i][k] + a[i][k])*dt/2
                                      + (old_j[i][k] - jk[i][k])*dt*dt/12;
                r[i][k] = old_r[i][k] + (old_v[i][k] + v[i][k])*dt/2
                                      + (old_a[i][k] - a[i][k])*dt*dt/12;
            }
        }