Previous ToC Up Next

4. Partitioned Runge-Kutta Algorithms

Bob: You promised a better approach to solving second-order differential equations, using Runge-Kutta schemes. What did you call such an algorithm again?

Alice: It is called a partitioned Runge-Kutta algorithm. The idea is to combine the force calculations in different ways for the position and for the velocity. The word `partitioned' here means that separate the treatment of from the treatment of . We already saw an example at the end of our previous discussion, where we had found a scheme that was almost, but not quite, a leapfrog scheme. If we would have tinkered with that scheme, we could have turned it into a leapfrog, but it would then no longer be a vector generalization of a Runge-Kutta scheme.

Bob: So you're saying that we have a lot more freedom, when we allow separate ways to update position and velocity, after first calculation a number of force evaluations.

Alice: Exactly. And we have already done this, for our fourth-order integrator, the one we plucked from Abramowitz and Stegun.

Bob: Does this mean that we can write our good old leapfrog as a partitioned Runge-Kutta scheme? That would be interesting! I have always thought about Runge-Kutta methods and the leapfrog scheme as two completely different animals, pardon the pun. Do you think that the leapfrog can be view as a type of Runge-Kutta algorithm?

Alice: I'm not sure. One reason to do this systematic landscape exploration is to find the answers to that type of question! And I'm sure we'll find out soon. By exploring all possible schemes with up to two new force calculations per step, we're bound to encounter the leapfrog, if indeed it is a citizen of the Runge-Kutta world.

So let us return to our special second-order differential equation

Let us first gather some useful expressions, starting with the two first-order equations

As before, we expand the position and the velocity of the orbit in Taylor series:

and when we differentiate the set of differential equations several times, we obtain the following equations:

The last three lines can be derived in the same way as the second line, by fully writing out the differentiations, using the chain rule. This derivation is completely analogous to what we did for our first-order differential equation.

4.1. One Force Evaluation per Step

Bob: I guess we will forget about force recycling, at least for now.

Alice: Yes. To keep things simple, let us look at a single integration step. But we have another choice to make.

In the case of a first-order differential equation, at the start of our integration we can only evaluate the right-hand side at time zero, at the beginning of the integration time step. If we simply follow that example, we start with:

This leads to the following dimensionally correct expressions:

We can write the expression for the position, substituting , as

We have to compare this with the Taylor series

Using Eqs. (XXXXX) we can write this as

Comparing Eqs. (149) and (151), we find:

and

We cannot satisfy the third order term in , so as far as the expression for the position is concerned, we can make our scheme second-order accurate.

Looking now at the velocity, we have

We have to compare this with the Taylor series

Using Eqs. (XXXXX) we can write this as

Comparing Eqs. (154) and (156), we find from the terms that are first order in :

Going to second order in would require that

which cannot be true for general and . Even though we can construct a second-order algorithm for the position, we can only find a first-order algorithm for the velocity.

Bob: And I presume that there is no point in using a higher-order algorithm for the position than for the velocity, since the overall order of the integration scheme must be the lowest order of that of the components. Hmmm. Is that so?

Alice: Yes, that is correct. For the very first step, it is possible in this case to find a new position that is second-order accurate. But as soon as we take the second step, we use the velocity that we arrived at in the first step, which is only first-order accurate. The same is true for each subsequent step: we always use the velocity value from the previous step.

Bob: But each time we multiply the velocity with . So even though the velocity is first-order, the product of the velocity with must be second-order correct, leading to an error term that is third-order in .

Alice: Yes, that is formally correct. However, \dots

We have to conclude that our approach only leads to a first-order correct algorithm, which is of course the forward Euler algorithm:

%\subsubsubsection{A Delayed Force Evaluation}

{\bf 4.1.1.2. A Delayed Force Evaluation}

Given the special form of our second-order differential equation, it is not necessary to start with a force evaluation at time zero. The first equation in the set

allows us to make a first-order prediction of the position, as:

which in turn allows us to postpone the first force evaluation to this non-zero time:

Note that this trick is not possible for a general force term that would depend on velocity as well. In that case, the last equation would read

which would mean that we need an initial force evaluation at time zero, before we can perform a subsequent force evaluation at time .

4.2. xxx

Let us exploit this extra freedom, for our special differential equation, by repeating our previous analysis for a delayed force evaluation. Our first force evaluation can now take place at time where is a free parameter. Using the linear extrapolation of the position, as sketched above, we obtain:

This leads to the following dimensionally correct expressions:

If we expand to first order in , we obtain:

Eq. (165) can thus be written as:

Comparing this with

EXPAND

let us first consider the terms, which leads to the requirement:

which leads to

Similarly, we can use the Taylor series for :

EXPAND

Considering the term, we have:

which leads to

Considering the terms, we have:

This implies:

In this way, all three free parameters are fixed, by requiring the algorithm to be second-order in in both position and velocity.

Could it be that we are in luck, and that this fixed solution can give us expressions for and that are also third-order correct in ? Let us start with the position equation. This would require:

Alas, since we are forced to use , the coefficient on the right-hand side is while the one on the left-hand side is . We have no freedom left, so this equation has no solutions for a general function and a general initial velocity .

%\subsubsubsection{Verlet-St\"ormer-Delambre Scheme}

{\bf 4.1.1.3. Verlet-St\"ormer-Delambre Scheme}

Using one evaluation of the right-hand side of the differential equation, we have thus arrived at the following second-order integration scheme:

Even though the equation for the velocity looks first-order, it is actually second-order accurate, through the clever choice of time at which the right-hand side of the differential equation is evaluated, namely in between the times at which and are determined.

In fact, this scheme is nothing else than the good old leapfrog algorithm, also known as the Verlet-St\"ormer-Delambre scheme, as we will show now.

Define

and

where

Our new scheme can then be written as

Alternatively, we can express the first of these equations in terms of

We have thus derived at expressions that show the leapfrog nature of the algorithm most clearly:

This representation shows clearly that the equations are fully time symmetric. This fact was rather hidden in the original formulation, eqs. (177). However, if we explicitly take a step forward and then another step backward, using eqs. (177), we can recover the time symmetry inherent in these equations, as follows. Let us denote the resulting position and velocity by , which are obtained from by taking a time step with size . Our task is to show that actually coincide with , not only to second order, as would be guaranteed in any second order scheme, but in fact to all orders in .

Finally, we can also write

and

We thus find

We can now shift our zero point in time by half a time step, to arrive at the more convenient notation:

Comparing this with our starting point, eqs. (177):

we have arrived at the alternative formulation:

The second formulation seems rather different, in that it requires two force calculations. Note, however, that the position at which the second force calculation takes place is exactly the same position at which the first force calculation for the next step will take place. Therefore, the second force calculation of each step can be recycled as the first force calculation of the next step. Effectively, we thus use only one force calculation per step. This trick is known in the literature as FSAL, short for First-Same-As-Last. We will come back to this point below.

%\subsubsubsection{An Historical Note}

{\bf 4.1.1.4. An Historical Note}

Almost everywhere in the literature, Runge-Kutta methods are assumed to start with : letting the first evaluation of the right-hand side of the differential equation take place at the very beginning of the step. This is necessary in the general case, but not for the special case of a second-order differential equation where there is no velocity dependence in the force term. The only place we have found so far in the literature, which mentions the possibility of starting with the force evaluation already at a later time is a paragraph in (1925), the original paper introducing what is now known as the algorithms.

In his section 2, p. 7, near the bottom, he remarks that, to be consistent, we should allow the freedom to write a general expression of the type we have done above in Eq. (164). He then adds that he decided against considering this extra freedom, for two reasons, both pragmatic, the first related to speed of execution of the algorithms, the second related to speed of derivation of the expressions fot the algorithms. Here are his arguments.

First of all, we often know already the force evaluation at the beginning of the step, from the last stage of the calculation of the previous step (at least approximately; and using even earlier force calculations, we can further improve the accuracy, without having to perform new force evaluations). Secondly, he adds, starting from such a general expression has led him to such unwieldy expressions that he was more or less forced to put in his equivalent to our Eq. (164).

Of course, current availability of algebraic manipulation programs have now invalidated his second argument. Curiously, all text books seem to propagate the simplifying assumption without questioning what the basis for this assumption may have been.

4.3. Two Force Evaluations per Step

%\subsubsubsection{General Form}

{\bf 4.1.2.1. General Form}

If we allow two evaluations of the right-hand side of the differential equation, we can work with the following general expression that is dimensionally correct

which leads to the following expressions for position and velocity steps:

Substituting the values, these equations expand into

So far, everything is still the exact prescription, given in the original algorithmic scheme. If we now expand the expressions in powers of , we get:

This should be equal to the Taylor series:

EXPAND

This leads to the following conditions:

These can be solved in terms of , as follows:

%\subsubsubsection{Examples}

{\bf 4.1.2.2. Examples}

If we take the standard assumption , we get:

This produces exactly what (1925) gives this as his simplest algorithm:

Henrici (1962) also lists this algorithm, refering back to (1925). However, Henrici's expressions contain a typo: he lists the last coefficient as . does list the term correctly, as .

If we try the other obvious choice , we find that some of the coefficients diverge: . With two force evaluations, it seems not to be possible to let the first one start right in the middle.

There is a natural choices that leads to relatively simple expressions for the coefficients: . Here the first force evaluation takes place after one third of the duration of the time step. In this case we get:

This leads to the following equations:

A complementary choise is , for which the first force evaluation takes place after two third of the duration of the time step. In this case we get:

This leads to the following equations:

There seem to be no other sets of simple coefficients. We might be tempted to try, say, , but in that case we get the much more complicated looking set:

This leads to the following equations:

xxx
Previous ToC Up Next