In the study of the LG, the number of particles is well within the limit of the PP method, but as seen in the previous section, this method is rather slow. Apart from the softening parameter, are there any other ways to speed things up?

Aa-63 proposed a method known as the individual time-step scheme. The idea is
that each particle is treated separately in the integration. The only coupling occurs at
the time when the force on a particle is calculated. Also, the assignment of a time-step
to each particle avoids the calculation of the force on every particle when only the
recalculation of the force on a particular particle is required. This is especially
important during close encounters when the force on a particle is changing very quickly
and, hence, necessitates recalculation at much higher frequency than for the rest of the
system. An improvement to this scheme has been made, called ``double individual''
time-steps. One then divide the force acting on a particle into two parts: one slowly
varying part due to the ``distant'' particles, and a more rapid varying part due to
particles in the immediate neighbourhood. This modification to the individual time-step
scheme has not been adapted into Sverre Aarseth's N-body code `NBODY1`, which I will
make use of, and will therefore not be discussed any further.

The basic idea of the method used in `NBODY1`, is that one calculates the force
acting on a particle at a given time *t* by extrapolating from the forces, and its
derivatives, calculated at earlier times. The particles will be assigned with two types
of coordinates: a primary and a secondary, and . The secondary,
evaluated at *t*, will be extrapolated from the primary, evaluated at . The same
procedure is used for the velocities. The method of extrapolation used in this case is
the Divided Difference Polynomial method.

Let us now take a closer look at this scheme Aa-85,AC. Mathematically, the
N-body gravitational problem used in `NBODY1` involves the solution of *N*
second-order differential equations,

The summation is over all particles of mass and coordinates .
This equation is only valid when the system in consideration is purely gravitational;
i.e. that all other forces than gravitation can be ignored. To be able to compare
`NBODY1` to `AVP`, we need to modify this picture somewhat. In AVP, the
cosmological constant play a major role for certain cosmological models, and has
therefore a substantial effect on the orbits of galaxies in these models. This means
that we have to introduce to `NBODY1`. From the Poisson equation for a
pressureless universe with a non-negligible cosmological constant, equation
(3.17), we see that the modified version of equation (5.13) has to
be

The cosmological constant , or rather , can be looked upon as a cosmic repulsive force, which counter the attractive gravitational force. Notice that equation (5.14) is the physical version of the Lagrangian in equation (3.20) solved in the AVP. The difference between AVP and N-body is therefore the way that these equations of motion are solved.

To avoid singularities as the mutual separation , thereby making
the numerical solution well behaved, and to reduce the computing time (see previous
section), one can introduce a softening parameter *c*. Equation (5.14)
then becomes

where .

To evolve the positions and the velocities for a particle from to *t*, a Taylor
series up to the fourth order in the force derivative is used:

where is the *k* derivative of the force with
respect to time, evaluated at time . The force is determined by extrapolating it
from the values of at four successive past epochs , and ,
where is the most recent. The fourth-order force polynomial is defined as

The three first divided differences are defined as (in compact notation)

where and the square brackets refer to the appropriate
time intervals. The term is defined similarly by and
. The force derivatives in terms of the divided differences
are obtained by derivation of equation (5.18) at successive times and by setting
. It should be noted that at , the fourth difference can not be determined
due to it needing the knowledge of the force at the following time *t*. Thus, the
contribution from is included as an correcting factor at the end of each
integration step. This so-called ``semi-iteration'' gives increased accuracy at little
extra cost and no extra memory.

Mon Jul 5 02:59:28 MET DST 1999