Previous | ToC | Up | Next |

** 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

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:

** 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:

** 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 :

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:

** 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

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:

** 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:

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:

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:

** 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!

** 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-01I'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

** 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

** 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 |