Abstract
The Earth orbits the Sun in a sort of cosmic shooting gallery, subject to impacts from comets and asteroids. It is only fairly recently that we have come to appreciate that these impacts by asteroids and comets (often called Near Earth Objects, or NEO) pose a significant hazard to life and property. Although the probability of the Earth being struck by a large NEO is extremely small, the consequences of such a collision are so catastrophic that it is adviseable to assess the nature of the threat and prepare to deal with it. One of the major issues in determining whether an asteroid can be dangerous for the Earth is given by the uncertainties in the determination of its position and velocity. Important studies arose from the previous issue, which dealt with the problem of getting accurate uncertainty estimates of the state of orbiting objects by means of several approaches, among which statistical theory has showed important results [1,2,3]. Recently, several tools and techniques have been developed for the robust prediction of Earth close encounters and for the identification of possible impacts of NEO with the Earth [4,5,6]. However, these methods might suffer from being either not sufficiently accurate when relying on simplifications (e.g., linear approximations) or computationally intensive when based on several integration runs (e.g., the Monte Carlo approach). Moreover, the standard integration schemes are affected by numerical integration errors, which might unacceptably accumulate during the integration, especially when long term integrations are performed. The resulting inaccuracies might strongly affect the validity of the results, precluding the use of such integrators. The necessity of solving these problems brought about a strong interest in self validated integration methods [7], which are based on the use of interval analysis. Interval analysis was originally formalized by Moore in 1966 [8]. The main idea beneath this theory is the substitution of real numbers with intervals of real numbers; consequently, interval arithmetic and analysis are developed in order to operate on the set of interval numbers in place of the classical analysis of real numbers. This turned out to be an effective tool for error and uncertainty propagation, as both the numerical errors and the uncertainties can be bounded by means of intervals, which are rigorously propagated in the computation by operating on them using interval analysis. Unfortunately, the naïve application of interval analysis might result in an unacceptable overestimation of the solution set of an ordinary differential equation (ODE). The reasons for such an overestimation are the so-called dependency problem and wrapping effect [9]. The dependency problem is related to the fact that, when a given variable occurs more than once in an interval computation, it is treated as a different variable at each occurrence; thus, for an interval number X = [a,b], even the operation X-X coincides with the operation X-Y, with Y equal to X but independent from it; i.e., X-X = [a-b,b-a], which includes the true solution X-X=0, but overestimates it. The wrapping effect is related to the fact that the naïve interval analysis computes interval enclosures of usually non-interval shaped sets. Focusing on the integration of ODE, the flow of nonlinear equations might turn out to deform and stretch an initial box of initial conditions. Thus, the exact solution set is non-interval shaped and non-convex in general. The Taylor models (TM) introduced in [10] and [11] solve both the dependency problem and the wrapping effect. The TM approach combines high-order multivariate polynomial techniques and the interval technique for verification. In particular, it represents the multivariate functional dependence of an arbitrary function by a high order multivariate Taylor polynomial P and a remainder interval I. The n-th order Taylor polynomial P captures the bulk of the functional dependency. Because the manipulation of those polynomials can be performed by operations on the coefficients where the minor errors due to their floating point nature are moved into the remainder bound, the major source of interval overestimation is removed. Thus, the overestimation only occurs in the remainder interval, the size of which scales with order n of the width of the domain. When applied to the verified integration of ODE, the relationships between the state vector at a generic time t and the initial conditions are expressed in terms of a Taylor model (P, I) and a tight enclosure for the action of the differential equations on an extended region is then provided. The TM-based integrator implemented in COSY VI [12] has been already succesfully exploited for the long-term rigourous integration of asteroids motion [13]. In this paper, an improved version of the integrator that exploits dynamic domain decomposition, automatic step size control, and a flow operator based on Lie derivatives [14], is applied to the more challenging task of Apophis close approach rigorous integration. The paper is organized as follows. The mathematical models adopted to describe Apophis dynamics and to evaluate the planetary ephemerides are presented first. Then, an introduction to Taylor models is given, focused on their application to the verified integration of ODE. The last part of the paper illustrates the results of the rigorous integration of Apophis motion. Conclusions end the paper.