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.     

A momentum-conserving N-body scheme with individual time steps

We present an N-body code called Taichi for galactic dynamics and controlled numerical experiments. The code includes two high-order hierarchical multipole expansion methods: the Barnes-Hut (BH) tree and the fast multipole method (FMM). For the time integration, the code can use either a conventional adaptive KDK or a Hamiltonian splitting integrator. The combination of FMM and the Hamiltonian splitting integrator leads to a momentum-conserving N-body scheme with individual time steps. We find Taichi performs well in the typical applications in galactic dynamics. In the isolated and interacting galaxies tests, the momentum conserving scheme produces the same result as a conventional BH tree code. But for similar force accuracies, FMM significantly speeds up the simulations compared to the monopole BH tree. In the cold collapse test, we find the inner structure after relaxation can be sensitive to the force accuracies. Taichi is ready to incorporate special treatment of close encounters thanks to the Hamiltonian splitting integrator, suitable for studying dynamics around central massive bodies.     

辛积分器中沿迹误差的一种补偿方法

辛积分器严格描述了一摄动Ｈａｍｉｌｔｏｎ系统的流，因而导致天体轨道的沿迹误差随时间呈线性增长趋势。本文利用这一特点，提出了一种对其沿迹误差进行估算的数值方法，从而达到了对数值结果进行沿迹误差补偿的目的，数值结果证实了此方法在较大积分步长和较长积分时间的数值计算中是有效的。     

Lie group variational integrators for the full body problem in orbital mechanics

Equations of motion, referred to as full body models, are developed to describe the dynamics of rigid bodies acting under their mutual gravitational potential. Continuous equations of motion and discrete equations of motion are derived using Hamilton’s principle. These equations are expressed in an inertial frame and in relative coordinates. The discrete equations of motion, referred to as a Lie group variational integrator, provide a geometrically exact and numerically efficient computational method for simulating full body dynamics in orbital mechanics; they are symplectic and momentum preserving, and they exhibit good energy behavior for exponentially long time periods. They are also efficient in only requiring a single evaluation of the gravity forces and moments per time step. The Lie group variational integrator also preserves the group structure without the use of local charts, reprojection, or constraints. Computational results are given for the dynamics of two rigid dumbbell bodies acting under their mutual gravity; these computational results demonstrate the superiority of the Lie group variational integrator compared with integrators that are not symplectic or do not preserve the Lie group structure.     

Numerical experiments on the efficiency of second-order mixed-variable symplectic integrators for N-body problems

We discuss the efficiency of the so-called mixed-variable symplectic integrators for N-body problems. By performing numerical experiments, we first show that the evolution of the mean error in action-like variables is strongly dependent on the initial configuration of the system. Then we study the effect of changing the stepsize when dealing with problems including close encounters between a particle and a planet. Considering a previous study of the slow encounter between comet P/Oterma and Jupiter, we show that the overall orbital patterns can be reproduced, but this depends on the chosen value of the maximum integration stepsize. Moreover the Jacobi constant in a restricted three-body problem is not conserved anymore when the stepsize is changed frequently: over a 10⁵ year time span, to keep a relative error in this integral of motion of the same order as that given by a Bulirsch-Stoer integrator requires a very small integration stepsize and much more computing time. However, an integration of a sample including 10⁴ particles close to Neptune shows that the distributions of the variation of the elements over one orbital period of the particles obtained by the Bulirsch-Stoer integrator and the symplectic integrator up to a certain integration stepsize are rather similar. Therefore, mixed-variable symplectic integrators are efficient either for N-body problems which do not include close encounters or for statistical investigations on a big sample of particles.     

几类辛方法的数值稳定性研究

主要对一阶隐式Euler辛方法M1、二阶隐式Euler中点辛方法M2、一阶显辛Euler方法M3和二阶leapfrog显辛积分器M4共4种辛方法及一些组合算法进行了通常意义下的线性稳定性分析．针对线性哈密顿系统,理论上找到每个数值方法的稳定区,然后用数值方法检验其正确性．对于哈密顿函数为实对称二次型的情况,为了理论推导便利,特推荐采用相似变换将二次型的矩阵对角化来研究辛方法的线性稳定性．当哈密顿分解为一个主要部分和一个小摄动次要部分且二者皆可积时,无论是线性系统还是非线性系统,这种主次分解与哈密顿具有动势能分解相比,明显扩大了辛方法的稳定步长范围．     

A symplectic integrator for the symmetry reduced and regularised planar 3-body problem with vanishing angular momentum

We construct an explicit reversible symplectic integrator for the planar 3-body problem with zero angular momentum. We start with a Hamiltonian of the planar 3-body problem that is globally regularised and fully symmetry reduced. This Hamiltonian is a sum of 10 polynomials each of which can be integrated exactly, and hence a symplectic integrator is constructed. The performance of the integrator is examined with three numerical examples: The figure eight, the Pythagorean orbit, and a periodic collision orbit.     

Extended Schubart Averaging

The purpose of this paper is the presentation of an integrator for the average motion of an asteroid in mean motion commensurability with Jupiter. The program is valid for any (p+q)/p mean motion commensurability (except whenq=0) and uses a double precision version of DE (Shampine and Gordon 1975) as propagator. The averaged equations of motion of the asteroid are evaluated in a non-singular way for any value of the eccentricities and the inclinations and the orbit of Jupiter is described by the most important terms in Longstop 1B (Nobiliet al. 1989). This integrator can be considered as an extension of the well known Schubart Averaging (Schubart 1978) in which Jupiter is moving on a fixed ellipse.     

A pseudo third order symplectic integrator

The symplectic integrator has been regarded as one of the optimal tools for research on qualitative secular evolution of Hamiltonian systems in solar system dynamics. An integrable and separate Hamiltonian system H = H₀ + Σ_i=1^N ε_iH_i (ε_i ≪ 1) forms a pseudo third order symplectic integrator, whose accuracy is approximately equal to that of the first order corrector of the Wisdom-Holman second order symplectic integrator or that of the Forest-Ruth fourth order symplectic integrator. In addition, the symplectic algorithm with force gradients is also suited to the treatment of the Hamiltonian system H = H₀(q,p) + εH₁(q), with accuracy better than that of the original symplectic integrator but not superior to that of the corresponding pseudo higher order symplectic integrator.     

Explicit Symplectic Integrator for Rotating Satellites

An explicit symplectic integrator is constructed for the problem of a rotating planetary satellite on a Keplerian orbit. The spin vector is fixed perpendicularly to the orbital plane. The integrator is constructed according to the Wisdom-Holman approach: the Hamiltonian is separated in two parts so that one of them is multiplied by a small parameter. The parameter depends on the satellite’s shape or the eccentricity of its orbit. The leading part of the Hamiltonian for small eccentricity orbits is similar to the simple pendulum and hence integrable; the perturbation does not depend on angular momentum which implies a trivial ‘kick’ solution. In spite of the necessity to evaluate elliptic function at each step, the explicit symplectic integrator proves to be quite efficient. This revised version was published online in July 2006 with corrections to the Cover Date.     

The long-term evolution and initial size of comets 46P/Wirtanen and 67P/Churyumov–Gerasimenko

We present a new method to study the long-term evolution of cometary nuclei in order to estimate their original size, and we consider the case of comets 46P/Wirtanen (hereafter 46P) and 67P/Churyumov–Gerasimenko (hereafter 67P). We calculate the past evolution of the orbital elements of both comets over 100 000 yr using a Bulirsch–Stoer integrator and over 450 000 yr using a Radau integrator, and we incorporate a realistic model of the erosion of their nucleus. Their long-term orbital evolution is prominently chaotic, resulting from several close encounters with planets, and this result is independent of the choice of the integrator and of the presence or not of non-gravitational forces. The dynamical lifetime of comet 46P is estimated at ∼133 000 yr and that of comet 67P at ∼105 000 yr. Our erosion model assumes a spherical nucleus composed of a macroscopic mixture of two thermally decoupled components, dust and pure water ice. Erosion strongly depends upon the active fraction and the density of the nucleus. It mainly takes place at heliocentric distances <4 au and lasts for only ∼7 per cent of the lifetime. Assuming a density of 300 kg m⁻³ and an average active fraction over time of 10 per cent, we find an initial radius of ∼1.3 km for 46P and ∼2.8 km for 67P. Upper limit are obtained assuming a density of 100 kg m⁻³ and an active fraction of 100 per cent, and amounts to 21 km for 46P and 25 km for 67P. Erosion acts as a rejuvenating process of the surface so that exposed materials on the surface may only contain very little quantities of primordial materials. However, materials located just under it (a few centimetres to metres) may still be much less evolved. We will apply this method to several other comets in the future.     

Explicit Symplectic Algorithms For Time‐Transformed Hamiltonians

By Hamiltonian manipulation we demonstrate the existence of separable time‐transformed Hamiltonians in the extended phase‐space. Due to separability explicit symplectic methods are available for the solution of the equations of motion. If the simple leapfrog integrator is used, in case of two‐body motion, the method produces an exact Keplerian ellipse in which only the time‐coordinate has an error. Numerical tests show that even the rectilinear N‐body problem is feasible using only the leapfrog integrator. In practical terms the method cannot compete with regularized codes, but may provide new directions for studies of symplectic N‐body integration. This revised version was published online in July 2006 with corrections to the Cover Date.     

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.     

QYMSYM: A GPU-accelerated hybrid symplectic integrator that permits close encounters

We describe a parallel hybrid symplectic integrator for planetary system integration that runs on a graphics processing unit (GPU). The integrator identifies close approaches between particles and switches from symplectic to Hermite algorithms for particles that require higher resolution integrations. The integrator is approximately as accurate as other hybrid symplectic integrators but is GPU accelerated.     

Satellite attitude dynamics and estimation with the implicit midpoint method

We describe the application of the implicit midpoint integrator to the problem of attitude dynamics for low-altitude satellites without the use of quaternions. Initially, we consider the satellite to rotate without external torques applied to it. We compare the numerical solution with the exact solution in terms of Jacobi’s elliptic functions. Then, we include the gravity-gradient torque, where the implicit midpoint integrator proves to be a fast, simple and accurate method. Higher-order versions of the implicit midpoint scheme are compared to Gauss–Legendre Runge–Kutta methods in terms of accuracy and processing time. Finally, we investigate the performance of a parameter-adaptive Kalman filter based on the implicit midpoint integrator for the determination of the principal moments of inertia through observations.     

Numerical integration for the real time production of fundamental ephemerides over a wide time span

A simplified model of the solar system has been developed along with an integration method, enabling to compute planetary and lunar ephemerides to an accuracy better than 1 and 2 milliarcsecs, respectively. On current personal computers, the integration procedure (SOLEX) is fast enough that by using a relatively small ( 20 Kbytes/Cy) database of starting conditions, any epoch in the time interval (up to ±100 Cy) covered by the database can be reached by the integrator in a few seconds. This makes the algorithm convenient for the direct computation of high precision ephemerides over a time span of several millennia.     

Explicit Adaptive Symplectic (Easy) Integrators: A Scaling Invariant Generalisation of the Levi-Civitaand KS Regularisations

We present a generalisation of the Levi-Civita and Kustaanheimo-Stiefel regularisation. This allows the use of more general time rescalings. In particular, it is possible to find a regularisation which removes the singularity of the equations and preserves scaling invariance. In addition, these equations can, in certain cases, be integrated with explicit symplectic Runge-Kutta-Nyström methods. The combination of both techniques gives an explicit adaptive symplectic (EASY) integrator. We apply those methods to some perturbations of the Kepler problem and illustrate, by means of some numerical examples, when scaling invariant regularisations are more efficient that the LC/KS regularisation.     

天体轨道长期数值积分的误差估计方法

数值积分方法是进行天体力学研究的重要工具, 尤其对于行星历表的研究工作而言. 由于在使用数值方法计算天体轨道时, 最终误差通常是难以预知的, 所以在面对精度要求较高或者积分时间较长的工作时具体积分方案的设计---尤其是当使用定步长方法时的步长选择---需要十分谨慎, 因为这将意味着是否能在时间成本可以被接受的范围内使解的精度达到要求. 因此, 在使用数值方法解决实际问题时如何快速寻找效率与精度之间的最佳平衡点是每一个数值积分方法的设计者与使用者都会面临的难题. 为解决这一问题, 在定步长条件下对数值积分方法的舍入误差概率分布函数以及截断误差积累量对步长的依赖关系和随时间的增长关系进行了深入研究. 基于所得结论, 提出了一种仅需较少的数值实验资料即可对选择任意时间步长积分至任意积分时刻时的舍入误差概率分布函数与截断误差积累量进行准确估计的方法, 并使用Adams-Cowell方法对该误差估计方法在圆周期轨道条件下进行了验证. 该误差估计方法在未来有望用于不同数值算法的性能对比研究, 同时也可以对数值积分方法求解实际轨道问题时的决策工作带来重要帮助.     

The Long-term Error Estimation Method for the Numerical Integrations of Celestial Orbits

Numerical methods have become a very important type of tool for celestial mechanics, especially in the study of planetary ephemerides. The errors generated during the computation are hard to know beforehand when applying a certain numerical integrator to solve a certain orbit. In that case, it is not easy to design a certain integrator for a certain celestial case when the requirement of accuracy were extremely high or the time-span of the integration were extremely large. Especially when a fixed-step method is applied, the caution and effort it takes would always be tremendous in finding a suitable time-step, because it is about whether the accuracy and time-cost of the final result are acceptable. Thus, finding the best balance between efficiency and accuracy with the least time cost appeared to be a major obstruction in the face of both numerical integrator designers and their users. To solve this problem, we investigate the variation pattern of truncation error and the pattern of rounding error distributions with time-step and time-span of the integration. According to those patterns, we promote an error estimation method that could predict the distribution of rounding errors and the total truncation errors with any time-step at any time-spot with little experimental cost, and test it with the Adams-Cowell method in the calculation of circular periodic orbits. This error estimation method is expected to be applied to the comparison of the performance of different numerical integrators, and also it can be of great help for finding the best solution to certain cases of complex celestial orbits calculations.     

Comparing maps to symplectic integrators in a galactic type Hamiltonian

We obtain thex - p _xPoincare phase plane for a two dimensional, resonant, galactic type Hamiltonian using conventional numerical integration, a second order symplectic integrator and a map based on the averaged Hamiltonian. It is found that all three methods give good results, for small values of the perturbation parameter, while the symplectic integrator does a better job than the mapping, for large perturbations. The dynamical spectra are used to distinguish between regular and chaotic motion.