Recent progress in the theory and application of symplectic integrators

In this paper various aspect of symplectic integrators are reviewed. Symplectic integrators are numerical integration methods for Hamiltonian systems which are designed to conserve the symplectic structure exactly as the original flow. There are explicit symplectic schemes for systems of the formH=T(p)+V(q), and implicit schemes for general Hamiltonian systems. As a general property, symplectic integrators conserve the energy quite well and therefore an artificial damping (excitation) caused by the accumulation of the local truncation error cannot occur. Symplectic integrators have been applied to the Kepler problem, the motion of minor bodies in the solar system and the long-term evolution of outer planets.     

综合评价线性多步积分公式轨道积分性能的两项新指标

概述了制约线性多步积分公式轨道积分状态的多种因素 ;提出了综合评价线性多步积分公式积分性能的两项新指标。建议在对数值计算有较高精度要求的科研项目中 ,应将构造并选择适合研究项目的线性多步积分公式以及高效的积分方式列为课题前期工作的重要部分     

Variational and symplectic integrators for satellite relative orbit propagation including drag

Orbit propagation algorithms for satellite relative motion relying on Runge–Kutta integrators are non-symplectic—a situation that leads to incorrect global behavior and degraded accuracy. Thus, attempts have been made to apply symplectic methods to integrate satellite relative motion. However, so far all these symplectic propagation schemes have not taken into account the effect of atmospheric drag. In this paper, drag-generalized symplectic and variational algorithms for satellite relative orbit propagation are developed in different reference frames, and numerical simulations with and without the effect of atmospheric drag are presented. It is also shown that high-order versions of the newly-developed variational and symplectic propagators are more accurate and are significantly faster than Runge–Kutta-based integrators, even in the presence of atmospheric drag.     

适用于人造卫星轨道数值积分的线性多步法的研究

总结了多年来构造适合卫星轨道计算的线性多步积分方法的研究进展;介绍了在构造积分公式的过程中选取伪根的三条原则。结合具体算例,对推荐的适用于不同类型卫星轨道的五线性多积分公式的性能与对称文坛及科威尔方法作了详细地比对和评述。     

Bandlimited implicit Runge–Kutta integration for Astrodynamics

We describe a new method for numerical integration, dubbed bandlimited collocation implicit Runge–Kutta (BLC-IRK), and compare its efficiency in propagating orbits to existing techniques commonly used in Astrodynamics. The BLC-IRK scheme uses generalized Gaussian quadratures for bandlimited functions. This new method allows us to use significantly fewer force function evaluations than explicit Runge–Kutta schemes. In particular, we use a low-fidelity force model for most of the iterations, thus minimizing the number of high-fidelity force model evaluations. We also investigate the dense output capability of the new scheme, quantifying its accuracy for Earth orbits. We demonstrate that this numerical integration technique is faster than explicit methods of Dormand and Prince 5(4) and 8(7), Runge–Kutta–Fehlberg 7(8), and approaches the efficiency of the 8th-order Gauss–Jackson multistep method. We anticipate a significant acceleration of the scheme in a multiprocessor environment.     

Parallel/Vector Integration Methods for Dynamical Astronomy

This paper reviews three recent works on the numerical methods to integrate ordinary differential equations (ODE), which are specially designed for parallel, vector, and/or multi-processor-unit(PU) computers. The first is the Picard-Chebyshev method (Fukushima, 1997a). It obtains a global solution of ODE in the form of Chebyshev polynomial of large (> 1000) degree by applying the Picard iteration repeatedly. The iteration converges for smooth problems and/or perturbed dynamics. The method runs around 100-1000 times faster in the vector mode than in the scalar mode of a certain computer with vector processors (Fukushima, 1997b). The second is a parallelization of a symplectic integrator (Saha et al., 1997). It regards the implicit midpoint rules covering thousands of timesteps as large-scale nonlinear equations and solves them by the fixed-point iteration. The method is applicable to Hamiltonian systems and is expected to lead an acceleration factor of around 50 in parallel computers with more than 1000 PUs. The last is a parallelization of the extrapolation method (Ito and Fukushima, 1997). It performs trial integrations in parallel. Also the trial integrations are further accelerated by balancing computational load among PUs by the technique of folding. The method is all-purpose and achieves an acceleration factor of around 3.5 by using several PUs. Finally, we give a perspective on the parallelization of some implicit integrators which require multiple corrections in solving implicit formulas like the implicit Hermitian integrators (Makino and Aarseth, 1992), (Hut et al., 1995) or the implicit symmetric multistep methods (Fukushima, 1998), (Fukushima, 1999). This revised version was published online in July 2006 with corrections to the Cover Date.     

Explicit adaptive symplectic integrators for solving Hamiltonian systems

We consider Sundman and Poincaré transformations for the long-time numerical integration of Hamiltonian systems whose evolution occurs at different time scales. The transformed systems are numerically integrated using explicit symplectic methods. The schemes we consider are explicit symplectic methods with adaptive time steps and they generalise other methods from the literature, while exhibiting a high performance. The Sundman transformation can also be used on non-Hamiltonian systems while the Poincaré transformation can be used, in some cases, with more efficient symplectic integrators. The performance of both transformations with different symplectic methods is analysed on several numerical examples.     

Numerical integration methods for orbital motion

The present report compares Runge-Kutta, multistep and extrapolation methods for the numerical integration of ordinary differential equations and assesses their usefulness for orbit computations of solar system bodies or artificial satellites. The scope of earlier studies is extended by including various methods that have been developed only recently. Several performance tests reveal that modern single- and multistep methods can be similarly efficient over a wide range of eccentricities. Multistep methods are still preferable, however, for ephemeris predictions with a large number of dense output points.     

Satellite onboard orbit propagation using Deprit’s radial intermediary

Short-term satellite onboard orbit propagation is required when GPS position measurements are unavailable due to an obstruction or a malfunction. In this paper, it is shown that natural intermediary orbits of the main problem provide a useful alternative for the implementation of short-term onboard orbit propagators instead of direct numerical integration. Among these intermediaries, Deprit’s radial intermediary (DRI), obtained by the elimination of the parallax transformation, shows clear merits in terms of computational efficiency and accuracy. Indeed, this proposed analytical solution is free from elliptic integrals, as opposed to other intermediaries, thus speeding the evaluation of corresponding expressions. The only remaining equation to be solved by iterations is the Kepler equation, which in most of cases does not impact the total computation time. A comprehensive performance evaluation using Monte-Carlo simulations is performed for various orbital inclinations, showing that the analytical solution based on DRI outperforms a Dormand–Prince fixed-step Runge–Kutta integrator as the inclination grows.     

Boundary-value problem formulations for computing invariant manifolds and connecting orbits in the circular restricted three body problem

We demonstrate the remarkable effectiveness of boundary value formulations coupled to numerical continuation for the computation of stable and unstable manifolds in systems of ordinary differential equations. Specifically, we consider the circular restricted three-body problem (CR3BP), which models the motion of a satellite in an Earth–Moon-like system. The CR3BP has many well-known families of periodic orbits, such as the planar Lyapunov orbits and the non-planar vertical and halo orbits. We compute the unstable manifolds of selected vertical and halo orbits, which in several cases leads to the detection of heteroclinic connections from such a periodic orbit to invariant tori. Subsequent continuation of these connecting orbits with a suitable end point condition and allowing the energy level to vary leads to the further detection of apparent homoclinic connections from the base periodic orbit to itself, or the detection of heteroclinic connections from the base periodic orbit to other periodic orbits. Some of these connecting orbits are of potential interest in space mission design.     

Multistep methods of numerical integration using back-corrections

A new class of linear multistep methods is proposed for the solution of the equations of motion of certain dynamical systems encountered in celestial mechanics and astrodynamics. These methods are distinguished from the classical predictor-corrector methods in that they permit back-corrections of the solution to be made. As the integration advances in time, the numerical solution is corrected or improved at certain points in the past. The enhanced numerical stability of these methods allows the meaningful application of high-order algorithms. Consequently, stepsizes larger than those attainable with the classical methods may be adopted and thus greater over-all efficiency may be realized. The application of these methods to the problem of determining the orbit of an artificial satellite is accomplished and the results are compared with those obtained using classical methods.     

Two multiderivate multistep (MDMS) integrators

A preliminary survey of multiderivative multistep integrators is carried out. It is found that all of them are much more accurate than the classical linear multistep methods, but most of them have poor stability. After parameter adjustment, two of them (called MDMS I and MDMS II by us) are competitive with or superior to the classical methods in some aspects, such as accuracy and stability. MDMS I behaves especially well in all the cases which have been studied.     

On a time-symmetric Hermite integrator for planetary N-body simulation

We describe a P(EC) ⁿ Hermite scheme for planetary N -body simulation. The fourth-order implicit Hermite scheme is a time-symmetric integrator that has no secular energy error for the integration of periodic orbits with time-symmetric time-steps. In general N -body problems, however, this advantage is of little practical significance, since it is difficult to achieve time-symmetry with individual variable time-steps. However, we can easily enjoy the benefit of the time-symmetric Hermite integrator in planetary N -body systems, where all bodies spend most of the time on nearly circular orbits. These orbits are integrated with almost constant time-steps even if we adopt the individual time-step scheme. The P(EC) ⁿ Hermite scheme and almost constant time-steps reduce the integration error greatly. For example, the energy error of the P(EC)² Hermite scheme is two orders of magnitude smaller than that of the standard PEC Hermite scheme in the case of an N  = 100,  m  = 10²⁵ g planetesimal system with the rms eccentricity 〈 e ²〉^1/2 ≲0.03.     

Sixth- and eighth-order Hermite integrator for N-body simulations

We present sixth- and eighth-order Hermite integrators for astrophysical N-body simulations, which use the derivatives of accelerations up to second-order (snap) and third-order (crackle). These schemes do not require previous values for the corrector, and require only one previous value to construct the predictor. Thus, they are fairly easy to implement. The additional cost of the calculation of the higher-order derivatives is not very high. Even for the eighth-order scheme, the number of floating-point operations for force calculation is only about two times larger than that for traditional fourth-order Hermite scheme. The sixth-order scheme is better than the traditional fourth-order scheme for most cases. When the required accuracy is very high, the eighth-order one is the best. These high-order schemes have several practical advantages. For example, they allow a larger number of particles to be integrated in parallel than the fourth-order scheme does, resulting in higher execution efficiency in both general-purpose parallel computers and GRAPE systems.     

Symplectic Integrators: Rotations and Roundoff Errors

We investigate the numerical implementation of a symplectic integrator combined with a rotation (as in the case of an elongated rotating primary). We show that a straightforward implementation of the rotation as a matrix multiplication destroys the conservative property of the global integrator, due to roundoff errors. According to Blank et al. (1997), there exists a KAM-like theorem for twist maps, where the angle of rotation is a function of the radius. This theorem proves the existence of invariant tori which confine the orbit and prevent shifts in radius. We replace the rotation by a twist map or a combination of shears that display the same kind of behaviour and show that we are able not only to recover the conservative properties of the rotation, but also make it more efficient in term of computing time. Next we test the shear combination together with symplectic integrator of order 2, 4, and 6 on a Keplerian orbit. The resulting integrator is conservative down to the roundoff errors. No linear drift of the energy remains, only a divergence as the square root of the number of iterations is to be seen, as in a random walk. We finally test the three symplectic integrators on a real case problem of the orbit of a satellite around an elongated irregular fast rotating primary. We compare these integrators to the well-known general purpose, self-adaptative Bulirsch–Stoer integrator. The sixth order symplectic integrator is more accurate and faster than the Bulirsch–Stoer integrator. The second- and fourth- order integrators are faster, but of interest only when extreme speed is mandatory. This revised version was published online in August 2006 with corrections to the Cover Date.     

The equations of relative motion in the orbital reference frame

The analysis of relative motion of two spacecraft in Earth-bound orbits is usually carried out on the basis of simplifying assumptions. In particular, the reference spacecraft is assumed to follow a circular orbit, in which case the equations of relative motion are governed by the well-known Hill–Clohessy–Wiltshire equations. Circular motion is not, however, a solution when the Earth’s flattening is accounted for, except for equatorial orbits, where in any case the acceleration term is not Newtonian. Several attempts have been made to account for the \(J_2\) effects, either by ingeniously taking advantage of their differential effects, or by cleverly introducing ad-hoc terms in the equations of motion on the basis of geometrical analysis of the \(J_2\) perturbing effects. Analysis of relative motion about an unperturbed elliptical orbit is the next step in complexity. Relative motion about a \(J_2\)-perturbed elliptic reference trajectory is clearly a challenging problem, which has received little attention. All these problems are based on either the Hill–Clohessy–Wiltshire equations for circular reference motion, or the de Vries/Tschauner–Hempel equations for elliptical reference motion, which are both approximate versions of the exact equations of relative motion. The main difference between the exact and approximate forms of these equations consists in the expression for the angular velocity and the angular acceleration of the rotating reference frame with respect to an inertial reference frame. The rotating reference frame is invariably taken as the local orbital frame, i.e., the RTN frame generated by the radial, the transverse, and the normal directions along the primary spacecraft orbit. Some authors have tried to account for the non-constant nature of the angular velocity vector, but have limited their correction to a mean motion value consistent with the \(J_2\) perturbation terms. However, the angular velocity vector is also affected in direction, which causes precession of the node and the argument of perigee, i.e., of the entire orbital plane. Here we provide a derivation of the exact equations of relative motion by expressing the angular velocity of the RTN frame in terms of the state vector of the reference spacecraft. As such, these equations are completely general, in the sense that the orbit of the reference spacecraft need only be known through its ephemeris, and therefore subject to any force field whatever. It is also shown that these equations reduce to either the Hill–Clohessy–Wiltshire, or the Tschauner–Hempel equations, depending on the level of approximation. The explicit form of the equations of relative motion with respect to a \(J_2\)-perturbed reference orbit is also introduced.     

Constructive-analytical solution of the problem of the secular evolution of polar satellite orbits

The well-known twice-averaged Hill problem is considered by taking into account the oblateness of the central body. This problem has several integrable cases that have been studied qualitatively by many scientists, beginning with M.L. Lidov and Y. Kozai. However, no rigorous analytical solution can be obtained in these cases due to the complexity of the integrals. This paper is devoted to studying the case where the equatorial plane of the central body coincides with the plane of its orbital motion relative to the perturbing body, while the satellite itself moves in a polar orbit. A more detailed qualitative study is performed, and an approximate constructive-analytical solution of the evolution system in the form of explicit time dependences of the eccentricity and pericenter argument of the satellite orbit is proposed. The methodical accuracy for the polar orbits of lunar satellites has been estimated by comparison with the numerical solution of the system.     

Orbital maneuver optimization using time-explicit power series

Orbital maneuver transfer time optimization is traditionally accomplished using direct numerical sampling to find the mission design with the lowest delta-v requirements. The availability of explicit time series solutions to the Lambert orbit determination problem allows for the total delta-v of a series of orbital maneuvers to be expressed as an algebraic function of only the individual transfer times. The delta-v function is then minimized for a series of maneuvers by finding the optimal transfer times for each orbital arc. Results are shown for the classical example of the Hohmann transfer, a noncoplanar transfer as well as an interplanetary fly-by mission to the asteroids Pallas and Juno.     

Orbit covariance propagation via quadratic-order state transition matrix in curvilinear coordinates

In this paper, an analytical second-order state transition matrix (STM) for relative motion in curvilinear coordinates is presented and applied to the problem of orbit uncertainty propagation in nearly circular orbits (eccentricity smaller than 0.1). The matrix is obtained by linearization around a second-order analytical approximation of the relative motion recently proposed by one of the authors and can be seen as a second-order extension of the curvilinear Clohessy–Wiltshire (C–W) solution. The accuracy of the uncertainty propagation is assessed by comparison with numerical results based on Monte Carlo propagation of a high-fidelity model including geopotential and third-body perturbations. Results show that the proposed STM can greatly improve the accuracy of the predicted relative state: the average error is found to be at least one order of magnitude smaller compared to the curvilinear C–W solution. In addition, the effect of environmental perturbations on the uncertainty propagation is shown to be negligible up to several revolutions in the geostationary region and for a few revolutions in low Earth orbit in the worst case.     

Sequential solution to Kepler’s equation

Seven sequential starter values for solving Kepler’s equation are proposed for fast orbit propagation. The proposed methods have constant complexity (not iterative), do not require pre-computed data, and can be implemented in just a few lines of code. The resulting sequential orbit propagation techniques can be done at different levels of accuracy and speed, depending essentially on the value of orbit eccentricity. Accuracy and algorithmic complexity are evaluated for all the proposed approaches and compared with several existing single-point techniques to solve Kepler’s equation. The new methods obtain improved accuracy at lower computational cost as compared to the best existing methods.