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

for the position and velocity of an individual particle, where the index indicates the values for time and for the time after one more time step has been taken: . The acceleration induced on a particle by the gravitational forces of all other particles is indicated by . 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!

平成18年10月10日