Numerical integration of dynamical systems with Lie series: Relativistic acceleration and non-gravitational forces
aa r X i v : . [ a s t r o - ph . E P ] A ug Noname manuscript No. (will be inserted by the editor)
Numerical integration of dynamical systems with Lie series relativistic acceleration and non-gravitational forces
D. Bancelin · D. Hestro ff er · W. Thuillot
Received: date / Accepted: date
Abstract
The integration of the equations of motion in gravitational dynamical systems—either in our Solar System or for extra-solar planetary systems— being non integrable in theglobal case, is usually performed by means of numerical integration. Among the di ff erentnumerical techniques available for solving ordinary di ff erential equations, the numerical in-tegration using Lie series has shown some advantages. In its original form [Hanslmeier and Dvorak,1984], it was limited to the N − body problem where only gravitational interactions are takeninto account.We present in this paper a generalisation of the method by deriving an expression of theLie-terms when other major forces are considered. As a matter of fact, previous studies hadbeen made but only for objects moving under gravitational attraction. If other perturbationsare added, the Lie integrator has to be re-built. In the present work we consider two casesinvolving position and position-velocity dependent perturbations: relativistic accelerationin the framework of General Relativity and a simplified force for the Yarkovsky e ff ect. Ageneral iteration procedure is applied to derive the Lie series to any order and precision. Wethen give an application to the integration of the equation of motions for typical Near-Earthobjects and planet Mercury. Keywords
Numerical Methods · Celestial Mechanics · Lie Series · General Relativity · Yarkovsky e ff ect D. BancelinIMCCE, Paris Observatory, CNRS, UPMC77, Av. Denfert-Rochereau 75014 Paris FranceTel.: +
33 1 4051 2271Fax: +
33 1 4633 2834E-mail: [email protected]. Hestro ff erTel.: +
33 1 4051 2260E-mail: [email protected]. ThuillotTel.: +
33 1 4051 2262E-mail: [email protected] D. Bancelin et al.
The integration of the equations of motion in gravitational dynamical systems, since theyare not integrable in closed form with more than 2 massive bodies involved, is usually per-formed by means of numerical integration. This can be done over times scales of a fewdays or a century for establishing ephemerides, or longer periods of several million yearsto establish stability properties [P´al and S¨uli, 2007]. It concerns bodies in our own SolarSystem as well as extra-solar planetary systems. Several algorithms have been used or de-veloped in particular in celestial mechanics to solve ordinary di ff erential equations (here-after ODEs) such as Runge-Kutta method, Everhart’s integrator, Bulirsch-Stoer, Adams, LieSeries, etc. Numerical integration by means of Lie series [Hanslmeier and Dvorak, 1984],which is based on generating the Taylor expansion of an ODEs solution, has shown someinterest. It is a high speed symplectic integrator adapted to the case of the N − body problem.The construction of Lie-series integrator was first investigated by Hanslmeier and Dvorak[1984], Delva [1985] and Lichtenegger [1984]. It has also been developed for studying thestability and dynamical evolution of Near-Earth Objects (NEOs) in the inner Solar System[Dvorak and Pilat-Lohinger, 1999, P´al and S¨uli, 2007], or the stability of extra-solar systems[Schwarz et al., 2005]. It has been developed to also integrate the case of varying masses[Dvorak and Lichtenegger, 1983] and the damped oscillator [Dvorak and Hanslmeier, 1983].In the previous works, a recurrence formula for the Lie-terms has been developed butonly for masses moving under their mutual gravitational attraction. As a consequence, itcannot be used as such for masses (e.g. asteroids and comets) for which relativistic or non-gravitational forces are non-negligible.When dealing with comets, NEOs, highly eccentric orbits, and extra-solar systems, otherforces and acceleration can come into play. Objects close to their star with significant ec-centricity will show a precession of their perihelion from the relativistic acceleration. Suchperihelion precession will a ff ect the transit time determination in extra-solar planets quests[P´al and Kocsis, 2008]. It can also a ff ect general studies of long-term stability as it has beenput into evidence in our Solar System [Laskar, 2009], since secular resonances can eventu-ally be avoided or induced from this additional frequency. It will also a ff ect the ephemeridesof comets and NEOs because these can have small perihelion distance and large eccentrici-ties. In addition to the relativistic acceleration, NEOs can also be a ff ected by the Yarkovskye ff ect (Vokrouhlick´y [1998], Bottke et al. [2002]) . This e ff ect is caused by a recoil force ofanisotropically emitted thermal radiation and is of importance for small bodies ( < ∼
20 km indiameter) close to the Sun, for much smaller bodies the solar radiation pressure will havea larger influence. The Yarkovsky e ff ect can be divided into a seasonal and diurnal e ff ect,depending on the orientation of the spin axis. It is the diurnal force that will have a ma-jor e ff ect on the orbit of an asteroid, by yielding a secular variation of the semi-major axis[Bottke et al., 2002].Although the use of Lie-series for integrating the N − body problem has been studiedby Gr¨obner [1967] it was not successfully implemented until the major breakthrough fromHanslmeier and Dvorak [1984] who provided a recurrence formula for practical use of theseries expansion. As a matter of fact, the Lie integrator has to be redesigned when otherforces are taken into account.After a brief introduction of the Lie operator and use of Lie series for the numericalintegration of the N − body problem (see section 2), we present in section 3 a generalisedexpression of the Lie Series integration by providing the recurrence formula applicable tothe case of a relativistic acceleration and a velocity-dependent force (applicable to the caseof the Yarkovsky e ff ect). We give some numerical tests and an application to the integration umerical integration of dynamical systems with Lie series 3 of the equation of motions for NEOs and planet Mercury in section 4. We then conclude anddiscuss future possible developments in section 5. D over a manifold of dimension N is defined as: D = N X k = θ k ∂ · ∂ z k (1)where θ = ( θ ( z ), ..., θ N ( z )) is a holomorphic function over a domain D in the z -space, thatis,it can be expanded in a converging power series, and z = (z , ..., z N ).Applying n − times this operator on a function f( z ), holomorphic over the same domain D , we have: D f = N X k = θ k ∂ f ∂ z k ; D n f = D ( D n − f ) (2)We remind here some useful properties of linearity and Leibniz-rule of Lie operator such as: ( D n ( f + g ) = D n f + D n gD n ( f · g ) = P nk = (cid:16) nk (cid:17) D k f D n − k g (3)The Lie-series are then defined as: L ( z , t ) ≡ ∞ X µ = t µ µ ! D µ f ( z ) µ ∈ N (4)which is converging over D [Gr¨obner, 1967]. And from its similarity with the expansion ofthe exponential function, we can write symbolically: L ( z , t ) ≡ e tD f ( z ) (5)This expression can be used to solve a system of ODEs such as:˙ z k = θ k ( z ) (6)The solution is z k = e tD ( ξ k ), where ξ k are the initial conditions, and D = P Nk = θ k ( ξ ) ∂ · ∂ξ k .Similarly, the approximate solution of ˙z = θ ( z ) at time t + ∆ t is: z ( t + ∆ t ) = L ( z , ∆ t ) = e ∆ tD z ( t ) (7) D. Bancelin et al. N − body problemWe can apply the Lie integration method to solve the equations of motion of the N − bodyproblem. Since we will generalise this problem by including additional forces, we remindhere the general steps of the procedure for the purely gravitational N − body case as given inHanslmeier and Dvorak [1984].According to the law of attraction, the Newtonian acceleration resulting from the forcesacting on a particle ν is: ¨ x ν = G N X µ = ,µ , ν M µ ( x µ − x ν ) k x µ − x ν k (8)where G is the gravitational constant and x ν the barycentric position vector of particle ν withmass M ν . By introducing a new variable v , the 2 nd order system (8) can be transform into asystem of 1 st order di ff erential equations.˙ x ν = v ν (9)˙ v ν = G N X µ = ,µ , ν M µ ( x µ − x ν ) k x µ − x ν k = GM ⋆ N X µ = ,µ , ν m µ x νµ ρ − νµ (10)such as ρ νµ = k x µ − x ν k and x νµ = x µ − x ν .Here M ⋆ is a conversion factor to express the mass M µ of the perturbing body in the unit ofstar mass. Its numeric value is equal to the central star’s mass of the planetary system. Thusm µ is the mass of the perturbing body in the unit of star mass (or Solar mass if GM ⋆ = k with k representing the Gauss constant ).The Lie operator for the Newtonian gravitational N − body system has been given byEggl and Dvorak [2010] and can be expressed as : D = N X ν = v ν ∂ · ∂ x ν + GM ⋆ N X µ = ,µ , ν m µ x νµ ρ − νµ ∂ · ∂ v ν (11)According to Eq. (1) and (6) the solutions of Eq. (9) and (10) are given by the series expan-sion: x ν ( τ ) = e τ D x ν (0) = ∞ X n = τ n D n n ! x ν (0) (12) v ν ( τ ) = e τ D v ν (0) = ∞ X n = τ n D n n ! v ν (0) (13)where τ is the current step-size and is defined by: τ = t j − t j − For this method to be of practical use, a recurrence formula has to be given to deriveany order of the Lie operator. This was given in Eggl and Dvorak [2010] for the N − bodyproblem: D n x ν = GM ⋆ N X µ = ,µ , ν m µ n − X k = n − k ! D k Φ νµ D n − − k x νµ (14) G has to be converted first in AU kg − day − umerical integration of dynamical systems with Lie series 5 where Φ νµ = ρ − νµ .The evolution of Φ νµ is given by: D n Φ νµ = ρ − νµ n − X k = a n , k + D n − − k Φ νµ D k Λ νµ (15)the a i , j coe ffi cients being defined for n ≥ a n , n = − a n , k = a n − , k − + a n − , k for k > a n , = a n − , − Λ νµ is defined by: Λ νµ = x νµ · v νµ and its evolution is ruled by: D n Λ νµ = n X k = nk ! D n − k x νµ D k v νµ (16)where D n v νµ = D n + x νµ (17)Now if additionnal forces or accelerations act on the particule ν , we have to re-write thesystem (9), (10), the Lie operator (11), and derive a new recurrence formula to obtain the n th derivative D n of the Lie operator. This section is dedicated to a generalised expression of the Lie operator for position-dependentand velocity-dependent forces acting on a particle ν . Let H ν be the contribution of all theaccelerations derived from those forces acting on body ν .As in the previous section, the second order ODE is split into a system of six first-orderdi ff erential equations: ˙ x ν = v ν (18)˙ v ν = H ν (19)The Lie operator becomes: D = N X µ = " v µ · ∂ · ∂ x µ + H µ · ∂ · ∂ v µ (20)Applying D to x , the construction of the Lie terms begins as: D x ν = x ν D x ν = v ν D x ν = H ν D x ν = D H ν D. Bancelin et al.
From which we can deduce the recurrence formula for D n x ν and D n v ν as following: ( D n x ν = D n − H ν n ≥ D n v ν = D n + x ν (21)The next step is to express Eq. (21) according to the expression of the forces which meansthat we have to find a recurrence formula for D n H ν . This will be done in the follow-ing sections for three kinds of accelerations: gravitational ( γγγ G ), relativistic ( γγγ R ) and non-gravitational accelerations ( γγγ Y ). Thus H = γγγ G + γγγ R + γγγ Y .3.1 Lie terms for the gravitational accelerationThe gravitational acceletation γγγ G /ν derived from the gravitational forces of N bodies actingon a body ν is: γγγ G /ν = GM ⋆ N X µ = ,µ , ν m µ Φ νµ, x νµ (22)Here we have introduced the new variable Φ νµ, p = ρ − p νµ as it will be useful in the other sections.From the first applications of the Lie operator on γγγ G /ν : D γγγ G /ν = GM ⋆ N X µ = ,µ , ν m µ (cid:0) Φ νµ, D x νµ + D Φ νµ, x νµ (cid:1) (23) D γγγ G /ν = GM ⋆ N X µ = ,µ , ν m µ (cid:16) Φ νµ, D x νµ + D Φ νµ, x νµ + D Φ νµ, x νµ (cid:17) (24)we can deduce the evolution of D n γγγ G /ν : D n γγγ G /ν = GM ⋆ N X µ = ,µ , ν m µ n X k = nk ! D k Φ νµ, D n − k x νµ (25)The evolution of the D n Φ νµ, p is ruled by: D n Φ νµ, p = ρ − νµ n − X k = a n , k + D n − − k Φ νµ, p D k Λ νµ (26)with a n , n = − pa n , k = a n − , k − + a n − , k for k > a n , = a n − , − n x νµ and D n Λ νµ are given respectively by Eq. (21) and (16). umerical integration of dynamical systems with Lie series 7 ff man) equationsis very slow even using modern computers and are not suitable for long time integration.Following Beutler [2005], the expression of the relativistic acceleration γγγ R /ν used to generatethe Lie-terms is a lighter version (provided that the mass m ν of the massless body is negliglecompared to the mass of the central star expressed as: γγγ R /ν = GM ⋆ c Φ µν, h(cid:16) GM ⋆ Φ µν, − ˙r (cid:17) r + r · ˙r ) ˙r i (27)where c represents the speed of light and µ the Sun.Here, r and ˙r are the heliocentric position and velocity of the body ν and are related with thebarycentric coordinates by: r = x ν − x µ = x µν ˙r = v ν − v µ = v µν r = ρ µν Equation (27) can be written as follows: γγγ R /ν = GM ⋆ c Φ µν, ( γγγ + γγγ ) (28)where ( γγγ = x µν (cid:16) GM ⋆ Φ µν, − v µν (cid:17) γγγ = v µν Λ µν (29)In the same way, the first derivatives of D n γγγ R /ν lead to the algorithm: D n γγγ R /ν = GM ⋆ c n X k = nk ! D k Φ µν, D n − k ( γγγ + γγγ ) (30)We will now express the evolution of the D n γγγ i by watching the evolution of their firstderivatives. Using the Lie properties expressed in Eq. (3) there is no di ffi culty to find thisevolution. This approach leads to the recurrence formulas: D n γγγ = P nk = (cid:16) nk (cid:17) D n − k x µν (cid:16) GM ⋆ D k Φ µν, − P kk ′ = (cid:16) kk ′ (cid:17) D k ′ v µν D k − k ′ v µν (cid:17) D n γγγ = P nk = (cid:16) nk (cid:17) D k v µν D n − k Λ µν (31)We remind that the evolution of the D n x ν and D n v ν are given by Eq. (21) and the evolu-tion of the D n Φ νµ, p and D n Λ νµ are given by Eq. (26) and Eq. (16). D. Bancelin et al. ≈ ff ects the orbital velocity, the transversal component causes achange of semimajor axis. Thus the Yarkovsky e ff ect is the strongest non-gravitational forceacting on small bodies when comparing with the intensity of Poynting-Robertson and solarradiation pressure e ff ect.This force can be decomposed into two e ff ects [Vokrouhlick´y et al., 2000]: – the diurnal acceleration that depends on the rotation frequency of the body around itsspin axis, is maximum when the spin vector is perpendicular to the orbital plane (itmeans that obliquity is equal to 0 ◦ or 180 ◦ ) and causes a drift in semimajor axis (positiveor negative depending on the value of the spin obliquity). – the seasonal acceleration that depends on the mean motion frequency of the objectaround the Sun, acts when the spin is in-plane (with obliquity equals to 90 ◦ ) and causesa steady decrease of the semimajor axis.Many tests have been done to see the impact of Yarkovsky e ff ect on the orbital mo-tion of asteroids and potentially hazardous asteroids (PHAs) (Giorgini et al. [2008], Chesley[2006]). Those authors show that the Yarkovsky force cannot be neglected as it has an ef-fect on the post-close encounter orbit of PHAs. Although some physical parameters can beunknown, it is possible to take this e ff ect into account without any assumptions regardingphysical characteristics of asteroids.Following Marsden [1976], the non-gravitational acceleration acting on comets can bedecomposed into three components: radial, transverse and normal. For asteroids, the mainnon-gravitational perturbation due to the Yarkovsky e ff ect is caused by the transverse com-ponent. So, we can simply express the acceleration like: γγγ Y = A g ( r ) T (32)where T = ( r ˙r − ( ˙r · r ) r / r ) / h is the transverse unit vector (in the direction of motion). Here( r , ˙r ) are the heliocentric position and velocity of the asteroid and h = k r ∧ ˙r k is the norm ofthe angular momentum of the asteroid. The coe ffi cient A is a non gravitational parameterdepending on the body. This parameter is generally accurately determined with three appari-tions of a comet, or with optical and radar observations for asteroids. For instance, the valueof A for an asteroid is ≈ − AU / day and ≈ − AU / day for comets.Finally, g(r) is a function depending on the heliocentric distance of the object. For Near-Earth Objects, we can simply express this function as: g(r) = (cid:16) A . U . r (cid:17) Using the barycentric coordinates: ( T = (cid:16) ρ µν x µν − ( v µν · x µν ) x µν ρ − µν (cid:17) / hh = k x µν ∧ v µν k To express the Lie terms for the Yarkovsky acceleration, it is better to write γγγ Y like: umerical integration of dynamical systems with Lie series 9 γγγ Y = A Φ µν, J ( T − T ) (33)with: T = Φ µν, − v µν T = Λ µν Φ µν, x µν J = h − Φ µν, = g ( r )As in the previous sections, it is easy to find the recurrence formula for the evolution of γγγ Y : D n γγγ Y = A n X k = nk ! D n − k ( T − T ) k X k ′ = kk ′ ! D k ′ J D k − k ′ Φ µν, (34)Let J = ( h · h ) − / . Thus, D J = − ( h · h ) − / D h · h . If Γ = D h · h then D J = − J Γ The construction of D n J goes as: D J = J ( − Γ D J − J D Γ ) D J = J (cid:16) − Γ D J − D Γ D J − J D Γ (cid:17) and we can deduce that: D n J = J n − X k = a n , k + D n − − k J D k Γ (35)where the coe ffi cients a i , j are defined, for n ≥
1, by: a n , n = − a n , k = a n − , k − + a n − , k for k > a n , = a n − , − D ( x ∧ v ) = D x ∧ v + x ∧ D v Thus D n h = P nk = (cid:16) nk (cid:17) D k x µν ∧ D n − k v µν D n Γ = P nk = (cid:16) nk (cid:17) D k + h · D n − k h (36)Finally, the expression of D n T and D n T can readily be found: D n T = n X k = nk ! D k Φ µν, − , D n − k v µν (37)and D n T = n X k = nk ! D n − k x µν k X k ′ = kk ′ ! D k ′ Λ µν D k − k ′ Φ µν, (38)The evolution of D n x ν and D n v ν are again given by Eq. ( 21). In this section, we provide some numerical tests of validation for the redesigned Lie series.Previous CPU and accuracy tests have already been provided in Eggl and Dvorak [2010]wherein the authors compared two packages of integrators among them Radau 15, Bulirsch-Stoer and Lie. As a complement to the previous CPU tests, we tested the time computationfor those integrators when considering massless bodies simultaneously integrated. Then, wecompare the close encounter results between the Earth and asteroid (99942) Apophis withthe Radau and Bulirsh-Stoer integrators. Finally, we propose some tests for the relativis-tic algorithm (towards computation of the precession of Mercury’s perihelion) and for theYarkovsky e ff ect algorithm (towards computation of secular drift of semi-major axis forsome asteroids).4.1 CPU testsAs a complement to the previous CPU tests already done [Eggl and Dvorak, 2010], we pro-pose in this part to investigate the integration time for a complete Solar System (Sun toNeptune, including Moon) and some varying number of massless bodies . Those minorbodies are taken from a synthetic population of NEOs moving only under gravitation. Thepurpose of this test is to compare the CPU time e ffi ciency of Lie and Radau integrators whenintegrating simultaneously massless bodies. For this test, the internal accuracy for both in-tegrators is set to 10 − . The ODE solved for Radau is a Nclass = Fig. 1
CPU time for Lie and Radau integrators. The abscissa represents the number of minor bodies simul-taneously integrated.
We represented in Fig. 1 the CPU time for both integrators. The abscissa represents thenumber of massless bodies integrated over 100 years. One can see that the time computationfor Lie increases faster than Radau as a function of the number of massless bodies simulate-nously integrated. As a consequence, if Lie is faster for a number less than 50 minor bodies,up to this number, Radau becomes more e ffi cient. More investigations have to be done toexplain this behaviour e.g. requiered work-steps for both integrators, use of the L2 cacheinstead of L1 cache, etc, .... So, as noted in Eggl and Dvorak [2010] the use of Lie integratordepends on the problem to be solved. umerical integration of dynamical systems with Lie series 11 Table 1
Date and distance of the 2029-close approach of asteroid Apophis computed with three numericalintegrators: Radau, Bulirsch and Lie.Integrator Date of close approach Distance of close approach (AU)Lie 2029 04 13.90716 0.000253446941Bulirsh-Stoer 2029 04 13.90714 0.000253446071Radau 2029 04 13.90725 0.000253444524
Table 1 shows the results for the 2029-close approach of asteroid Apophis computedwith three di ff erent integrators: Lie, Bulirsh and Radau. We can see that the value of thedistance and date of close approach computed with Lie integrator is in a good agreementwith the value computed with Radau and Burlirsh-Stoer integrator. This test shows that nowLie integrator can handle deep close approach determination with a good precision, compa-rable to other integrators. It could thus be used for close-approach analysis and detection ofimpacts with Earth.4.3 Validation test for General Relativity algorithmFirst of all, we tested our general relativity algorithm by integrating the equation of motionof asteroid (1566) Icarus. This asteroid is a good example for this test in as much as, itseccentricity is ∼ ff ect. We calcu-late the relativistic acceleration with Lie series and then we compare our results with Radauintegrator. For both integrators, the numerical precision is set to 10 − . Figure 2 representsthe value of this acceleration (between ∼ − and ∼ − ) calculated with Lie and Radauand shows that those value are very close. The absolute di ff erence between those valuesvalidates the Lie algorithm as this di ff erence lies between ∼ − and ∼ − .We also tested our algorithm by calculating the perihelion precession of Mercury andcomparing with the expected value. The General Relativity predicts a secular precession ofMercury’s perihelion which expression is [Balogh and Giampieri, 2002]: d ω dt = π GM ⊙ c a (1 − e ) ≈ ′′ / revolution (39)where ω , a and e are respectively, the perihelion, semimajor axis and eccentricity of Mer-cury, M ⊙ the mass of the Sun. Fig. 2 relativistic acceleration calculated with Lie ( + ) and Radau Integrator (—) and the absolute di ff erence ∆γγγ (- - -) between the value of γγγ calculated with Radau and Lie integrators. The y axis is on a logarithmicscale. Figure 3 represents the relativistic e ff ect on orbital parameter ω for 5 Mercury orbital pe-riods. This graphic shows that Mercury’s perihelion position shifts of about 0.1 arcsec foreach Mercury year. Moreover, we also compute the absolute di ff erence ∆ ( ∆ω ) between ∆ω calculated with Eq. (39) and with the Lie integrator. This di ff erence is about 1.10 − arc-sec / revolution what shows the Lie value is in a good agreement with the analytical result.4.4 Validation test for the Yarkovsky e ff ect algorithmIn this section, we test our Yarkovsky e ff ect algorithm by calculating the secular drift ( da / dt )of the semimajor axis caused by this e ff ect. From the knowledge of da / dt , we can deduce avalue for the constant component A in Eq. (32). As a matter of fact, the expression of thesemimajor axis variations due to Yarkovsky e ff ect along the transverse component, given byGauss equations, is: dadt = p (1 − e ) n (1 + e cos θ ) k γγγ Y k with θ the true anomaly. Thus, we deduce that the secular drift of the semimajor axis perrevolution of the object is: ∆ a = π ak (1 − e ) A (40)We propagated the motion of nine candidate asteroids for Yarkovsky perturbation overfive Keplerian revolutions ( ∆ T ) around the Sun. The asteroids considered here are taken umerical integration of dynamical systems with Lie series 13 Fig. 3 relativistic e ff ect on Mercury’s perihelion calculated with Lie integrator in arcsec ( - ). ∆ ( ∆ω ) representsthe absolute di ff erence between the theory and Lie value in arcsec ( + ). We also represented the perihelionprecession per revolution ( • ). Table 2
Secular drift of semimajor axis per revolution calculated with Lie and Radau integrator.
Asteroids Period A da / dt (Chesley) da / dt (Radau) da / dt (Lie)(day) (AU / day ) km / rev km / rev km / revGolevka 1442.31 -1.47 × − -0.3785 -0.3616 -0.3604Apollo 651.17 -4.67 × − -0.0642 -0.0612 -0.0611Ra-Shalom 277.26 -1.20 × − -0.0807 -0.0781 -0.0782Bacchus 408.90 -2.22 × − -0.1778 -0.1675 -0.1675YORP 368.43 -5.49 × − -0.3801 -0.3595 -0.3594Hathor 283.32 -2.35 × − -0.1622 -0.1405 -0.1405Cerberus 409.94 -1.46 × − -0.1313 -0.1324 -0.1323Geographos 507.70 -2.68 × − -0.0246 -0.0242 -0.0242Toro 583.92 -1.13 × − -0.0125 -0.0118 -0.0117 from Chesley et al. [2008], with the correspondent value of the da / dt given. We comparedthe mean value of da / dt calculated with Lie integrator, with Radau integrator to the oneprovided in Chesley et al. [2008]. The results listed in table 2 show that the value of da / dt found with Lie integrator is in good agreement with Radau’s one. Moreover, those valuesare not far from Chesley’s one ( ∼
4% or 5% higher or lower). This di ff erence can come fromthe fact that the Yarkovsky model used in our simulation may be di ff erent than Chesley’sone. We presented in this paper a redesigned Lie integrator including relativistic and a simpli-fied Yarkovsky forces. In addition to the CPU and accuracy tests already performed, thenumerical tests presented here show the great capacities of this extended Lie integrator thatmake it a useful and complete integrator. This integrator can be useful when dealing withlong-time integration and it can also handle simultaneous integration of massless bodies aswell as Radau integrator. Besides, the integration of PHAs and the study of close encountersand impact probabilities can also be done now thanks to the relativistic acceleration andYarkovsky e ff ect algorithms.For the future, the Lie integrator can be extended to the problem of the integration ofcomets’s motion. The Yarkovsky e ff ect algorithm is already a first approach for this study inthat, the transverse vector is already implemented. Besides, one can consider to generalisethe Lie integrator for any position and velocity dependent forces. References
A. Balogh and G. Giampieri. Mercury: the planet and its orbit.
Reports on Progress inPhysics , 65:529–560, April 2002.G. Beutler.
Methods of celestial mechanics. Vol. I: Physical, mathematical, and numericalprinciples . Beutler, G., 2005.W. F. Bottke, Jr., D. Vokrouhlick´y, D. P. Rubincam, and M. Broz.
The E ff ect of YarkovskyThermal Forces on the Dynamical Evolution of Asteroids and Meteoroids . 2002.S. R. Chesley. Potential impact detection for Near-Earth asteroids: the case of 99942Apophis (2004 MN 4 ). In L. Daniela, M. Sylvio Ferraz, & F. J. Angel, editor, Aster-oids, Comets, Meteors , volume 229 of
IAU Symposium , pages 215–228, 2006.S. R. Chesley, D. Vokrouhlick´y, S. J. Ostro, L. A. M. Benner, J.-L. Margot, R. L. Matson,M. C. Nolan, and M. K. Shepard. Direct Estimation of Yarkovsky Accelerations on Near-Earth Asteroids.
LPI Contributions , 1405:8330– + , 2008.M. Delva. A Lie integrator program and test for the elliptic restricted three body problem. A & AS , 60:277–284, May 1985.R. Dvorak and A. Hanslmeier. Numerical Integration with Lie-Series. In S. Ferraz-Mello& P. E. Nacozy, editor, The Motion of Planets and Natural and Artifical Satellites , pages65– + , 1983.R. Dvorak and H. Lichtenegger. On the two-body problem with variable masses. InS. Ferraz-Mello & P. E. Nacozy, editor, The Motion of Planets and Natural and Artifi-cal Satellites , pages 11–17, 1983.R. Dvorak and E. Pilat-Lohinger. On the dynamical evolution of the Atens and the Apollos.
Planet. Space Sci. , 47:665–677, May 1999.S. Eggl and R. Dvorak. An Introduction to Common Numerical Integration Codes Usedin Dynamical Astronomy. In J. Souchay & R. Dvorak, editor,
Lecture Notes in Physics,Berlin Springer Verlag , volume 790 of
Lecture Notes in Physics, Berlin Springer Verlag ,pages 431–480, March 2010.J. D. Giorgini, L. A. M. Benner, S. J. Ostro, M. C. Nolan, and M. W. Busch. Predicting theEarth encounters of (99942) Apophis.
Icarus , 193:1–19, January 2008.W. Gr¨obner. Die Lie-Reihen und ihre Anwendungen.
VEB Deutscher Verlag , 1967.A. Hanslmeier and R. Dvorak. Numerical Integration with Lie Series. A & A , 132:203– + ,March 1984. umerical integration of dynamical systems with Lie series 15 J. Laskar. Large Scale Solar Sytem Simulations. In
AAS / Division of Dynamical AstronomyMeeting , volume 40 of
AAS / Division of Dynamical Astronomy Meeting , May 2009.H. Lichtenegger. The dynamics of bodies with variable masses.
Celestial Mechanics , 34:357–368, December 1984.B. G. Marsden. Nongravitational forces on comets.
NASA Special Publication , 393:465–488, 1976.A. P´al and B. Kocsis. Periastron precession measurements in transiting extrasolar planetarysystems at the level of general relativity.
MNRAS , 389:191–198, September 2008.A. P´al and ´A. S¨uli. Solving linearized equations of the N-body problem using theLie-integration method.