2.1 Numerical integration

For the computation of highly accurate asteroid proper elements the numerical procedures involved, in particular the numerical integration, must also be as accurate and reliable as possible. For the purpose of our project, however, the choice of the integrator and the tuning of the parameters controlling the performance had to satisfy also another, practically opposite requirement, namely computational speed. This because we were to integrate thousands of orbits, computing not only the differential equations of motion, but also the corresponding variational equations to estimate maximum Lyapunov Characteristic Exponents (LCE); and this for a time span long enough to enable an accurate and efficient averaging out of both the short and the long-periodic effects. For all these reasons we decided to make use of the ORBIT9 integrator (version 9c), developed by one of us (AM), employing different dynamical models in different regions of the asteroid main belt. For the integrations we present in this paper (of the orbits of asteroids located in the outer main belt) we used only four outer planets as perturbing bodies, applying the ``barycentric correction'' to the initial conditions in order to account for the major part of the indirect effect of inner planets. The ORBIT9 propagator-interpolator is described elsewhere (e.g. Milani and Nobili [1988]); it uses as starter a symplectic single step method (implicit Runge-Kutta-Gauss), and a multi-step predictor performs most of the propagation; the latter has a constant stepsize h and uses only one evaluation of the right hand side of the equations of motion for each step. The states at the discrete steps are not interpolated, because the filter input interval is automatically selected to be a multiple of the step h. The latter is selected in an optimum way by using the estimates of the truncation error from Milani and Nobili [1988]. However, we have improved the algorithm for this stepsize selection by replacing the value of e in the corresponding equations by the value of (e² + sinI²)^1/2, which is more predictive of the long term behavior of the eccentricity. In this way a suitable stepsize is computed for each asteroid. Then the input catalog is sorted by the stepsize and integrated in batches with similar stepsizes, by using the minimum of the recommended stepsizes. We propagated the orbits of all the asteroids in the input catalog from the present state until 2 My in the past. Then proper elements were computed as discussed in Sections 2.3 and 2.4, and tested for stability in time as discussed there. For most asteroids the results turned out to be satisfactory; for the minority for which there were indications of instability and/or strong chaotic motion we performed a further integration from the initial state to 10 Myr in the future. The final proper elements catalog has been produced by joining the output files and replacing the 2 My results with the 10 My ones when available.