3.1 Choosing an Algorithm

Carol:

Now that we have the equations of motion for the relative position of one particle with respect to the other, I guess we need an algorithm to integrate these equations.

Alice:

Indeed, and there is a large choice! If you pick up any book on numerical methods, you will see that you can select from a variety of lower-order and higher-order integrators, and for each one there are additional choices as to the precise structure of the algorithm.

Bob:

What is the order of an algorithm?

Alice:

It signifies the rate of convergence. Since no algorithm with a finite time step size is perfect, they all make numerical errors. In a fourth-order algorithm, for example, this error scales as the fourth power of the time step - hence the name fourth order.

Carol:

If that is the case, why not take a tenth order or even a twentieth order algorithm. By only slightly reducing the time step, we would read machine accuracy, of order $10^{-15}$ for the usual double precision (8 byte, i.e. 64 bit) representation of floating point numbers.

Alice:

The drawback of using high-order integrators is two-fold: first, they are far more complex to code; and secondly, they do not allow arbitrarily large time steps, since their region of convergence is limited. As a consequence, there is an optimal order for each type of problem. When you want to integrate a relatively well-behaved system, such as the motion of the planets in the solar system, a twelfth-order integrator may well be optimal. Since all planets follow well-separated orbits, there will be no sudden surprises there. But when you integrate a star cluster, where some of the stars can come arbitrarily close to each other, experience shows that very high order integrators lose their edge. In practice, fourth-order integrators turn out to be optimal for the job.

Bob:

How about starting with the lowest-order integrator we can think off? A zeroth-order integrator would make no sense, since the error would remain constant, independent of the time step size. So the simplest one must be a first-order integrator.

Alice:

Indeed. And the simplest version of a first-order integrator is called the forward Euler integrator.

Bob:

Was Euler so forward-looking, or is there also a backward Euler algorithm?

Alice:

There is indeed. In the forward version, at each time step you simply take a step tangential to the orbit you are on. After that, at the next step, the new value of the acceleration forces you to slightly change direction, and again you move for a time step $dt$ in a straight line in that direction. Your approximate orbit is thus constructed out of a number of straight line segments, where each one has the proper direction at the beginning of the segment, but the wrong one at the end.

Bob:

And the backward Euler algorithm must have the right direction at the end of a time step, and the wrong one at the beginning. Let's see. That seems much harder to construct. How do you know at the beginning of a time step in what direction to move so that you come out with the right direction tangential to a correct orbit at that point?

Alice:

You can only do that through iteration. You guess a direction, and then you correct for the mistake you find yourself making, so that your second iteration is much more accurate, in fact first-order accurate. Given this extra complexity, I suggest that we start with the forward Euler algorithm.

Carol:

Can't we do both, i.e. make half the mistakes of each of the two, while trying to strike the right balance between forward and backward Euler?

Alice:

Aha! That is a good way to construct better algorithms, which then become second-order accurate, because you have canceled the first-order errors. Examples are second-order Runge Kutta, and leapfrog. We'll soon come to that, but for now let's keep it simple, and stay with first order. Here is the mathematical notation:

$${\bf r}_{i+1} = {\bf r}_i + {\bf v}_i dt$$ (3.1)

$${\bf v}_{i+1} = {\bf v}_i + {\bf a}_i dt$$ (3.2)

for the position ${\bf r}$ and velocity ${\bf v}$ of an individual particle, where the index $i$ indicates the values for time $t_i$ and $i+1$ for the time $t_{i+1}$ after one more time step has been taken: $dt = t_{i+1} - t_i$. The acceleration induced on a particle by the gravitational forces of all other particles is indicated by ${\bf a}$. So, all we have to do now is to code it up. By the way, let's rename the file. Rather than a generic name nbody.C, let's call it forward_euler1.C.

Bob:

Alice, go for it!