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.     

Symplectic integrator for general near-integrable Hamiltonian system

In this paper, following the idea of constructing the mixed symplectic integrator (MSI) for a separable Hamiltonian system, we give a low order mixed symplectic integrator for an inseparable, but nearly integrable, Hamiltonian system, Although the difference schemes of the integrators are implicit, they not only have a small truncation error but, due to near integrability, also a faster convergence rate of iterative solution than ordinary implicit integrators, Moreover, these second order integrators are time-reversible.     

辛方法的校正公式

1996年Wisdom等提出了对辛方法进行校正的概念和实践，现在继续对辛校正进行详尽讨论和数值比较，尤其对哈密顿函数可分解为一个主要部分和多个次要部分的一般情形，用Lie级数推导任意阶的各种辛算法的一次和二次辛校正公式并对一些算法给出具体的辛校正公式。又以日、木、土三体问题为模型进行数值实验，结果表明一次辛校正能提高精度，改善数值稳定性。计算效率也比较高，因而值得推荐使用，辛方法通常用大步长数值积分，这时二次辛校正并没有显著提高结果的精度，却大大增加了计算时间，不应予以推荐。     

Deflecting NEOs in Route of Collision with the Earth

We consider a small sample of known near Earth objects (NEOs), both asteroids and comets, with low minimum orbital intersection distance (MOID). Through a simple numerical procedure we generate slightly different orbits from this sample in such a way that these bodies will collide with the Earth at a specific epoch. Then we study the required change in orbital velocity (along track Δv) in order to deflect these NEOs at different epochs before the impact event. The orbital evolution of these NEOs is performed through a full N-body numerical integrator. A comparison with analytical estimates is also performed in selected cases. Interesting features in the Δv/time before impact plots are found; as a prominent result, we find that close approaches to the Earth before the epoch of the impact can make the overall deflection easier.     

A stability study of asteroid families near the 3:1 and 5:2 resonance with Jupiter

By using theD-criterion Lindblad (1992) has identified 14 asteroid families from a sample of 4100 numbered asteroids with proper elements from Milani and Kneevi (1990). Taxonomic types and other physical properties for a significant number of objects in five of the families show strong homogeneity within each family, further strengthening their internal relationship.To test the hypothesis of a common origin in, e.g., a catastrophic collision event, we have set out to integrate the orbits of the members of the Maria, Dora and Oppavia-Gefion families over some 10⁶ years. The mean distance for the Maria family is close to the 3:1 mean-motion resonance with Jupiter, while the other two families lie close to the 5:2 resonance.We used a simplified solar system model which included the perturbations by Jupiter and Saturn only and implemented Everhart's variable stepsize integrator RA15. All close encounters between the family members (within 0.1 AU) were recorded as well. Preliminary results from integrations over 4×10⁵ years are presented here.The statistics of close encounters show pronounced peaks for several members within each family, while for others no significant levels above the background of random encounters or even very low frequencies were found. This indicates a subclustering within the families. Quite a lot of very close (<0.005 AU) mutual encounters are found, which suggest that, at least for the larger members in a family, the mutual gravitational interactions could be of some importance for the real orbital evolutions.The encounter statistics between the Dora and Oppavia family members suggest a possible interrelationship between this two groups.     

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.     

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.     

Symplectic integrators with adaptive timestep applied to spinning compact binaries

This paper deals mainly with the application of the mixed leapfrog symplectic integrators with adaptive timestep to a conservative post-Newtonian Hamiltonian formulation with canonical spins for spinning compact binaries. The adaptive timestep depends on the two body separation r and the magnitude of the spins. Various numerical tests including a chaotic high-eccentricity orbit show that the fixed step symplectic integrators lost drastically the good long term behaviour in the test cases with large eccentricity, the adaptive timestep integrator is always superior to the constant step in the integral precision.     

A note on the algorithm of symplectic integrators

For a Hamiltonian that can be separated into N+1(N\geq 2) integrable parts, four algorithms can be built for a symplectic integrator. This research compares these algorithms for the first and second order integrators. We found that they have similar local truncation errors represented by error Hamiltonian but rather different numerical stability. When the computation of the main part of the Hamiltonian, H ₀, is not expensive, we recommend to use S ^* type algorithm, which cuts the calculation of the H ₀ system into several small time steps as Malhotra(1991) did. As to the order of the N+1 parts in one step calculation, we found that from the large to small would get a slower error accumulation. This revised version was published online in July 2006 with corrections to the Cover Date.     

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.     

New adaptive time step symplectic integrator:an application to the elliptic restricted three-body problem

The time-transformed leapfrog scheme of Mikkola Aarseth was specifically designed for a second-order differential equation with two individually separable forms of positions and velocities.It can have good numerical accuracy for extremely close two-body encounters in gravitating few-body systems with large mass ratios,but the non-time-transformed one does not work well.Following this idea,we develop a new explicit symplectic integrator with an adaptive time step that can be applied to a time-dependent Hamiltonian.Our method relies on a time step function having two distinct but equivalent forms and on the inclusion of two pairs of new canonical conjugate variables in the extended phase space.In addition,the Hamiltonian must be modified to be a new time-transformed Hamiltonian with three integrable parts.When this method is applied to the elliptic restricted three-body problem,its numerical precision is explicitly higher by several orders of magnitude than the nonadaptive one's,and its numerical stability is also better.In particular,it can eliminate the overestimation of Lyapunov exponents and suppress the spurious rapid growth of fast Lyapunov indicators for high-eccentricity orbits of a massless third body.The present technique will be useful for conservative systems including N-body problems in the Jacobian coordinates in the the field of solar system dynamics,and nonconservative systems such as a time-dependent barred galaxy model in a rotating coordinate system.     

Simulations for the Original Distribution of KBOs

Using the N-body dynamical model that includes the sun, the 8 planets, Pluto, UB313 and massless particles, we simulate the orbital evolution of 551 Kuiper Belt Objects (KBOs) with known parameters. The initial conditions of the simulations are the currently observed orbital parameters. The integration backtracks from now to -10×10⁸ yr. The results show that about 10×10⁸ years ago, more than 1/3 of the presently observed KBOs resided in the region of the present Kuiper main belt, a few were located inside the Neptune orbit, and the rest were beyond 50AU; and that about 4.5×10⁸ years ago, all the objects in the Kuiper main belt exhibited a rather good normal distribution, without so many objects concentrated in the Neptune's 3:2 resonance region, as at present time.     

N-body integrators with individual time steps from Hierarchical splitting

We review the implementation of individual particle time-stepping for N-body dynamics. We present a class of integrators derived from second order Hamiltonian splitting. In contrast to the usual implementation of individual time-stepping, these integrators are momentum conserving and show excellent energy conservation in conjunction with a symmetrized time step criterion. We use an explicit but approximate formula for the time symmetrization that is compatible with the use of individual time steps. No iterative scheme is necessary. We implement these ideas in the HUAYNO¹ code and present tests of the integrators and show that the presented integration schemes shows good energy conservation, with little or no systematic drift, while conserving momentum and angular momentum to machine precision for long term integrations.     

The orbital evolution of NEA 30825 1900 TG1

The orbital evolution of the near-Earth asteroid (NEA) 30825 1990 TG1 has been studied by numerical integration of the equations of its motion over the 100 000-year time interval with allowance for perturbations from eight major planets and Pluto, and the variations in its osculating orbit over this time interval were determined. The numerical integrations were performed using two methods: the Bulirsch-Stoer method and the Everhart method. The comparative analysis of the two resulting orbital evolutions of motion is presented for the time interval examined. The evolution of the asteroid motion is qualitatively the same for both variants, but the rate of evolution of the orbital elements is different. Our research confirms the known fact that the application of different integrators to the study of the long-term evolution of the NEA orbit may lead to different evolution tracks.     

Improved Leap—Frog Symplectic Integrators for Orbits of Small Eccentricity in the Perturbed Kepler Problem

I have improved the precision of the leap–frog symplectic integrators for perturbed Kepler problems at small eccentricities, without significant loss of CPU time. The integration scheme proposed is competitive, in some situations, with the so-called mixed variable integrators.     

Evolution of the obliquities of the giant planets in encounters during migration

Tsiganis et al. [Tsiganis, K., Gomes, R., Morbidelli, A., Levison, H.F., 2005. Nature 435, 459-461] have proposed that the current orbital architecture of the outer Solar System could have been established if it was initially compact and Jupiter and Saturn crossed the 2:1 orbital resonance by divergent migration. The crossing led to close encounters among the giant planets, but the orbital eccentricities and inclinations were damped to their current values by interactions with planetesimals. Brunini [Brunini, A., 2006. Nature 440, 1163-1165] has presented widely publicized numerical results showing that the close encounters led to the current obliquities of the giant planets. We present a simple analytic argument which shows that the change in the spin direction of a planet relative to an inertial frame during an encounter between the planets is very small and that the change in the obliquity (which is measured from the orbit normal) is due to the change in the orbital inclination. Since the inclinations are damped by planetesimal interactions on timescales much shorter than the timescales on which the spins precess due to the torques from the Sun, especially for Uranus and Neptune, the obliquities should return to small values if they are small before the encounters. We have performed simulations using the symplectic integrator SyMBA, modified to include spin evolution due to the torques from the Sun and mutual planetary interactions. Our numerical results are consistent with the analytic argument for no significant remnant obliquities.     

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 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.     

Rotation rate and velocity dispersion of planetary ring particles with size distribution II. Numerical simulation for gravitating particles

We examine rotation rates of gravitating particles in low optical depth rings, on the basis of the evolution equation of particle rotational energy derived by Ohtsuki [Ohtsuki, K., 2006. Rotation rate and velocity dispersion of planetary ring particles with size distribution. I. Formulation and analytic calculation. Icarus 183, 373-383]. We obtain the rates of evolution of particle rotation rate and velocity dispersion, using three-body orbital integration that takes into account distribution of random velocities and rotation rates. The obtained stirring and friction rates are used to calculate the evolution of velocity dispersion and rotation rate for particles in one- and two-size component rings as well as those with a narrow size distribution, and agreement with N-body simulation is confirmed. Then, we perform calculations to examine equilibrium rotation rates and velocity dispersion of gravitating ring particles with a broad size distribution, from 1 cm up to 10 m. We find that small particles spin rapidly with 〈ω²〉^1/2/Ω?¹⁰²-¹⁰³, where ω and Ω are the particle rotation rate and its orbital angular frequency, respectively, while the largest particles spin slowly, with 〈ω²〉^1/2/Ω?1. The vertical scale height of rapidly rotating small particles is much larger than that of slowly rotating large particles. Thus, rotational states of ring particles have vertical heterogeneity, which should be taken into account in modeling thermal infrared emission from Saturn's rings.