Previous ToC Up Next

11. The Hermite Algorithm

11.1. A Self-Starting Fourth-Order Scheme

Alice: I'm quite impressed with the large collection of integrators we have built up. But perhaps it is time to move on, and adapt all of them to the real N-body problem, now that we have their two-body form?

Bob: Yes, but there is one more integrator that I wanted to add to our pile. I had mentioned the Hermite algorithm, when we started to look at multiple-step methods, as an alternative to get fourth-order accuracy, without using multiple steps and without using half steps. Actually, coding it up was even simpler than I had expected. You could argue that it is the simplest self-starting Fourth-Order Scheme, once you allow the direct calculation of the jerk, in addition to that of the acceleration.

Alice: Can you show me first the idea behind the algorithm?

Bob: The trick is to differentiate Newton's gravitational equations of motion

which I have written here for the general N-body problem. Differentiation with respect to time gives us the jerk directly:

where .

Note that the jerk has one very convenient property. Although the expression above looks quite a bit more complicated than Newton's original equations, they can still be evaluated through one pass over the whole -body system. This is no longer true for higher derivatives. For example, we can obtain the fourth derivative of the position of particle , the snap, by differentiating the above equation one more time:

where , and this is the expression that thickens the plot. Unlike the and expressions, that are given by the initial conditions, has to be calculated from the positions and velocities. However, this calculation does not only involve the pairwise attraction of particle on particle , but in fact all pairwise attractions of all particles on each other! This follows immediately when we write out what the shorthand implies:

When we substitute this back into Eq. (88), we see that we have to do a double pass over the -body system, summing over both indices and in order to compute a single fourth derivative for the position of particle .

11.2. A Higher-Order Leapfrog Look-Alike

Alice: I see. So it is quite natural to differentiate Newton's equations of motion just once, and then to use both the acceleration and the jerk, obtained directly from Eqs. (86) and (87). Can you show me what the expressions for the position and velocity look like, in terms of those two quantities?

Bob: Here they are:

Alice: That looks like a direct generalization of the leapfrog! I had not expected such simplicity and elegance.

Bob: Yes, I think you can look at the Hermite scheme as a generalization of the leapfrog, in some sense. If you neglect the terms that are quadratic in , you get exactly the leapfrog back.

Alice: Yes, since in that case you can write the terms linear in as effectively half-step quantities, up to higher powers in :

since the quadratic terms are eliminated by the clever choice of half-step, which characterizes the leapfrog.

In general, you might have expected two extra terms, one of order and one of order , for the position, as well as for the velocity. Instead there is just one term, in each case. And that one term is written with a misleading factor of , but of course there is at least one more factor lurking in the factor , which to first order in is nothing else than ; or for that matter. And similarly, there is a factor hidden in .

The question is, can we understand why the remaining terms are what they are, and especially why they can be written in such a simple way?

Bob: I thought you would ask that, and yes, this time I've done my homework, anticipating your questions. Here is the simplest derivation that I could come up with.

Since we are aiming at a fourth-order scheme, all we have to do is to expand the position and velocity up to fourth-order terms in the time step , while we can drop one factor of for the acceleration, and two factors for the jerk, to the same order of consistency.

Alice: You mean that a term proportional to, say, in the jerk would be derived from a term proportional to in the velocity or from a term proportional to in the position, both of which are beyond our approximation.

Bob: Exactly. Here it is, in equation form:

Alice: I'll buy that.

11.3. A Derivation

Bob: We can now eliminate snap and crackle at time , expressing them in terms of the acceleration and jerk at times and , using the last two lines of Eq. (92). This gives us for the snap

Substituting this back into the last line of Eq. (92), we get for the crackle:

Substituting these expressions for the snap and crackle in the second line of Eq. (92) we find:

Indeed, we have recovered the second line of Eq. (90), and thereby explained the mysterious factor in the last term.

The only thing left to do is to check the expression for the position, the first line in Eq. (90). Let us bring everything to the left-hand side there, except the acceleration terms. In other words, let us split off the leapfrog expression, and see what is left over:

To within our order of accuracy, this is indeed what we were trying to prove, since the right-hand term of the first line in Eq. (90) can be written as

This proves the desired result:

11.4. Implementation

Alice: Great! Can you show me how you have implemented this scheme?

Bob: It was surprisingly simple. I opened a file called hbody.rb, with yoshida6.rb as my starting point. And all I had to do was to add a method jerk to compute the jerk, and then a method hermite to do the actual Hermite integration.

The jerk method, in fact, is rather similar to the acc method: 

   def acc
     r2 = @pos*@pos
     r3 = r2*sqrt(r2)
     @pos*(-@mass/r3)
   end

apart from the last line: 

   def jerk
     r2 = @pos*@pos
     r3 = r2*sqrt(r2)
     (@vel+@pos*(-3*(@pos*@vel)/r2))*(-@mass/r3)
   end

Alice: Before you show me the integrator, let me look at the equations again, which you have derived, to get an idea as to what to expect. Here they are:

Now I remember what bothered me about these expressions. On the left-hand side you express the positions and velocities at time . However, on the right-hand side you rely also on the acceleration and jerk at time , which you can only compute after you have determined the position and velocity at time . In addition, you also have the future velocity at the right-hand side as well.

In other words, these are implicit equations for the position and the velocity. How do you deal with them?

Bob: Through iteration. I first determine trial values for the position and velocity, simply by expanding both as a Taylor series, using only what we know at time , which are the position, velocity, acceleration and jerk. In fact, this is a type of predictor-corrector method. If I express those trial values with a subscript p, I start with:

Then, with those trial values, I determine the trial values for the acceleration and jerk. Using all those trial values as proxies for the actual values at time , I can solve the equations for the corrector step, where I indicate the final values for the position and velocity at time with a subscript s :

Alice: Does that really work? Let's see whether your procedure is actally fourth-order accurate. In Eq. (101), you leave out terms of order in the velocity, and terms of order in the position, but those predictor values for position and velocity are only used to compute acceleration and jerk, terms that are always multiplied with factors of , so the resulting errors are of order . Good!

However, I am worried about that velocity term in the first line on the right-hand side of Eq. (102). According to Eq. (101), your predicted velocity has an error of order , so by multiplying that with a single factor , you wind up with an error of order , which is certainly not acceptable.

Bob: Huh, you're right. How can that be? Is that really what I implemented? That cannot be correct. I'd better check my code.

Here it is, the hermite method that is supposed to do the job: 

   def hermite(dt)
     old_pos = @pos
     old_vel = @vel
     old_acc = acc
     old_jerk = jerk
     @pos += @vel*dt + old_acc*(dt*dt/2.0) + old_jerk*(dt*dt*dt/6.0)
     @vel += old_acc*dt + old_jerk*(dt*dt/2.0)
     @vel = old_vel + (old_acc + acc)*(dt/2.0) + (old_jerk - jerk)*(dt*dt/12.0)
     @pos = old_pos + (old_vel + vel)*(dt/2.0) + (old_acc - acc)*(dt*dt/12.0)
   end

Ah, I see what I did! Of course, and now I remember also why I did it that way. I was wrong in what I wrote above. I should have written Eq. (102) as follows:

Alice: I see. By switching the order of computation for the corrected form of the position and velocity, you are able to use the corrected version of the velocity, rather than the predicted version, in determining the corrected version of the position. So all is well!

Bob: Yes, indeed. And it is all coming back now: at first I computed and in the normal order, and I got a pretty bad energy behavior, until I realized that I should switch the order of computations. How easy to make a quick mistake and a quick fix, and then to forget all about it!

Alice: It shows once again the need for detailed documentation, not only for the mythical user, but at least as much for the actual producer of the code!

11.5. Testing

Bob: To test hermite, let us start with a comparison run, taking our old rk method: 

 |gravity> ruby integrator_driver2j.rb < euler.in
 dt = 0.01
 dt_dia = 0.1
 dt_out = 0.1
 dt_end = 0.1
 method = rk4
 at time t = 0, after 0 steps :
   E_kin = 0.125 , E_pot =  -1 , E_tot = -0.875
              E_tot - E_init = 0
   (E_tot - E_init) / E_init =-0
 at time t = 0.1, after 10 steps :
   E_kin = 0.129 , E_pot =  -1 , E_tot = -0.875
              E_tot - E_init = 1.79e-12
   (E_tot - E_init) / E_init =-2.04e-12
   1.0000000000000000e+00
   9.9499478009063858e-01  4.9916426216739009e-02
  -1.0020902861389222e-01  4.9748796005932194e-01
I'll give the Hermite scheme the same task: 

 |gravity> ruby integrator_driver5a.rb < euler.in
 dt = 0.01
 dt_dia = 0.1
 dt_out = 0.1
 dt_end = 0.1
 method = hermite
 at time t = 0, after 0 steps :
   E_kin = 0.125 , E_pot =  -1 , E_tot = -0.875
              E_tot - E_init = 0
   (E_tot - E_init) / E_init =-0
 at time t = 0.1, after 10 steps :
   E_kin = 0.129 , E_pot =  -1 , E_tot = -0.875
              E_tot - E_init = 5.31e-13
   (E_tot - E_init) / E_init =-6.07e-13
   1.0000000000000000e+00
   9.9499478009151798e-01  4.9916426220332356e-02
  -1.0020902857150518e-01  4.9748796006319129e-01
Alice: Not bad! Can you give it a somewhat larger time step, to check whether you really have a fourth-order scheme?

Bob: I'll increase the time step by a factor two: 

 |gravity> ruby integrator_driver5b.rb < euler.in
 dt = 0.02
 dt_dia = 0.1
 dt_out = 0.1
 dt_end = 0.1
 method = hermite
 at time t = 0, after 0 steps :
   E_kin = 0.125 , E_pot =  -1 , E_tot = -0.875
              E_tot - E_init = 0
   (E_tot - E_init) / E_init =-0
 at time t = 0.1, after 5 steps :
   E_kin = 0.129 , E_pot =  -1 , E_tot = -0.875
              E_tot - E_init = 7.93e-12
   (E_tot - E_init) / E_init =-9.07e-12
   1.0000000000000000e+00
   9.9499478011948561e-01  4.9916426283208984e-02
  -1.0020902812740490e-01  4.9748796010457508e-01
Alice: Close to a factor of 16 degradation, as expected. So we have a working fourth-order Hermite implementation.

Bob: And all that by just adding two short methods to an existing code, without any need for any other modifications! Ruby really invites you to a game of playful prototyping and exploring new algorithms. I guess I've become converted to using Ruby, by now.

Alice: I'm glad to hear that, but perhaps it is time to come back to our original project, to provide some toy models for your students, to do some actual N-body integrations.

Bob: Yes, I'm afraid I got a bit too much carried away with all these integrators. But it was so much fun! I've never in my life written an eighth-order integrator. Somehow, working in Fortran or C++ made such a thing seem like a tremendous project. But somehow, while I was playing around, one thing led to another, and before I knew it, I had no less than two very different eighth-order integrators.

You're right, though, let's get back to our original plan. While we have reached considerable sophistication on the level of algorithms, we are still stuck with . High time to go to large values.
Previous ToC Up Next