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