The foregoing method is not self-starting, so a separate initial routine is used to start the process. From the initial conditions , the respective Taylor series derivatives are formed by successive differentiation of equation (5.15), analytically. Having calculated the force and derivatives for each particle, the divided differences to third order and the time-steps are determined, and the main program is started with these quantities as if they had been obtained by the main program itself.
Since each particle is assigned its own time-step, the integration cycle itself begins by determining the next particle i to be advanced, which is the one with the minimum time-step. The coordinates and velocities of this particle are then determined to order by using equations (5.16) and (5.17), whereupon the current force may be obtained by direct summation. At this stage, the four past time epochs are updated to be consistent with the definition that denotes the most recent force evaluation. New differences are now formed, including . Together with the new , these correction terms are combined to improve the current coordinates and velocities to highest order, whereupon the primary coordinates are initialized by .
New time-steps are assigned initially for all particles and at the end of each integration cycle for particle i. The cycle is then repeated on the particle which now has the minimum time-step.
Figure: Action trajectories for 8 galaxies in the LG and the LN (solid lines), including N-body orbits (dotted lines). Cf. Figure 5.1.
In order to be able to compare the N-body orbits with the the AVP solution, the physical coordinates of the NBODY1 trajectories needs to be transformed into comoving coordinates. Since each positions in NBODY1 is given at time t, the expansion parameter at that time has to be determined. This was done by solving (3.33) for each time t, using the Runge-Kutta method. The code for this subroutine is given in Appendix B.