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.     

Parallel Construction of Nonlinear Force-Free Fields

A numerical approach to calculating nonlinear force-free fields is presented. The approach is similar to Sakurai (1981) being a current-field iteration scheme using the integral solution to Ampere's law (the Biot–Savart law). However, the method of solution presented here is simpler than Sakurai's approach, in that the field is directly constructed on a grid without the intermediate solution of a large system of nonlinear equations. The method also permits straightforward implementation on parallel computers. Results of applying the method to a number test cases, including boundary conditions with substantial currents, are presented.     

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

Fast Procedure Solving Universal Kepler's Equation

We developed a procedure to solve a modification of the standard form of the universal Kepler’s equation, which is expressed as a nondimensional equation with respect to a nondimensional variable. After reducing the domain of the variable and the argument by using the symmetry and the periodicity of the equation, the method first separates the case where the solution is so small that it is given an inverted series. Second, it separates the cases where the elliptic, parabolic, or hyperbolic standard forms of Kepler’s equation are suitable. Here the separation is done by judging whether detouring these nonuniversal equations will cause a 1-bit loss of information to their nonuniversal solutions or not. Then the nonuniversal equations are solved by the author’s procedures to solve the elliptic Kepler’s equation (Fukushima, 1997a), Barker’s equation (Fukushima, 1998), and the hyperbolic Kepler’s equation (Fukushima, 1997b), respectively. And their nonuniversal solutions are transformed back to the solution of the universal equation. For the rest of the case, we obtain an approximate solution by solving roughly the approximated cubic equation as we did in solving Barker’s equation. Then the correction to the approximate solution is obtained by Halley’s method precisely. There the special function appeared in the universal equation is rewritten into a combination of similar special functions of small arguments, so that they are efficiently evaluated by their Taylor series. Numerical measurements showed that, in the case of Intel Pentium II processor, the new method is 10–25 times as fast as Shepperd’s method (Shepperd, 1985) and 7–13 times as fast as the standard Newton method. This revised version was published online in July 2006 with corrections to the Cover Date.     

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.     

Review of planetary and satellite theories

Planetary and satellite theories have been historically and are presently intimately related to the available computing capabilities, the accuracy of observational data, and the requirements of the astronomical community. Thus, the development of computers made it possible to replace planetary and lunar general theories with numerical integrations, or special perturbation methods. In turn, the availability of inexpensive small computers and high-speed computers with inexpensive memory stimulated the requirement to change from numerical integration back to general theories, or representative ephemerides, where the ephemerides could be calculated for a given date rather than using a table look-up process. In parallel with this progression, the observational accuracy has improved such that general theories cannot presently achieve the accuracy of the observations, and, in turn, it appears that in some cases the models and methods of numerical integration also need to be improved for the accuracies of the observations. Planetary and lunar theories were originally developed to be able to predict phenomena, and provide what are now considered low accuracy ephemerides of the bodies. This proceeded to the requirement for high accuracy ephemerides, and the progression of accuracy improvement has led to the discoveries of the variable rotation of the Earth, several planets, and a satellite. By means of mapping techniques, it is now possible to integrate a model of the motion of the entire solar system back for the history of the solar system. The challenges for the future are: Can general planetary and lunar theories with an acceptable number of terms achieve the accuracies of observations? How can numerical integrations more accurately represent the true motions of the solar system? Can regularly available observations be improved in accuracy? What are the meanings and interpretations of stability and chaos with respect to the motions of the bodies of our solar system? There has been a parallel progress and development of problems in dealing with the motions of artificial satellites. The large number of bodies of various sizes in the limited space around the Earth, subject to the additional forces of drag, radiation pressure, and Earth zonal and tesseral forces, require more accurate theories, improved observational accuracies, and improved prediction capabilities, so that potential collisions may be avoided. This must be accomplished by efficient use of computer capabilities.     

The use of integrals in numerical integrations of theN-body problem

The numerical integration of systems of differential equations that possess integrals is often approached by using the integrals to reduce the number of degrees of freedom or by using the integrals as a partial check on the resulting solution, retaining the original number of degrees of freedom.Another use of the integrals is presented here. If the integrals have not been used to reduce the system, the solution of a numerical integration may be constrained to remain on the integral surfaces by a method that applies corrections to the solution at each integration step. The corrections are determined by using linearized forms of the integrals in a least-squares procedure.The results of an application of the method to numerical integrations of a gravitational system of 25-bodies are given. It is shown that by using the method to satisfy exactly the integrals of energy, angular momentum, and center of mass, a solution is obtained that is more accurate while using less time of calculation than if the integrals are not satisfied exactly. The relative accuracy is ascertained by forward and backward integrations of both the corrected and uncorrected solutions and by comparison with more accurate integrations using reduced step-sizes.     

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 Mars 1:2 resonant population

An excess of around 400 asteroids in the distribution of the semimajor axes of the asteroids is identified by means of numerical integrations as generated by a population of approximately 1000 asteroids evolving inside the exterior resonance 1:2 with Mars. Approximately 200 asteroids are librating around the asymmetric libration centers and their evolution in a time-scale of 1 million years appears stable but with a strong influence of Mars' eccentricity. The biggest Mars 1:2 resonant asteroid is (142) Polana.     

A new family of multistep numerical integration methods based on Hermite interpolation

In this paper, a new family of explicit and implicit multistep methods is presented both for the error-controlled and uncontrolled modes. The main concept is to replace the Newton interpolation with the Hermite interpolation, where the Hermite polynomial is fitted to the function values and its derivatives. This idea is very useful in the numerical solution of problems (e.g., orbit propagation problem) where higher-order derivatives can easily be computed. In addition to the theoretical concept, the stability regions of the proposed methods are determined. The new methods are more stable than the well-known multistep numerical integrators (i.e., Adams–Bashforth and Adams–Bashforth–Moulton) in the explicit, implicit, and predictor–corrector forms. Using the second-order derivatives gives smaller error constants in the proposed method. The new integrators are numerically tested for a few examples, and the solutions are compared with those of the well-known multistep methods. Moreover, the CPU time and absolute integration error are compared in the satellite orbit propagation problem using various integration methods. The CHAMP mission, i.e., a German small-satellite mission for geoscientific and atmospheric research and applications, is considered as a case study for comparing the achievable accuracy of the proposed method with the existing method for solving the two-body problem.     

GPU accelerated manifold correction method for spinning compact binaries

The graphics processing unit (GPU) acceleration of the manifold correction algorithm based on the compute unified device architecture (CUDA) technology is designed to simulate the dynamic evolution of the Post-Newtonian (PN) Hamiltonian formulation of spinning compact binaries. The feasibility and the efficiency of parallel computation on GPU have been confirmed by various numerical experiments. The numerical comparisons show that the accuracy on GPU execution of manifold corrections method has a good agreement with the execution of codes on merely central processing unit (CPU-based) method. The acceleration ability when the codes are implemented on GPU can increase enormously through the use of shared memory and register optimization techniques without additional hardware costs, implying that the speedup is nearly 13 times as compared with the codes executed on CPU for phase space scan (including \(314 \times 314\) orbits). In addition, GPU-accelerated manifold correction method is used to numerically study how dynamics are affected by the spin-induced quadrupole–monopole interaction for black hole binary system.     

Non-canonical perturbations in symplectic integration

The inclusion of non-canonical perturbations in symplectic integration schemes has been discussed. A rigorous derivation of an analog for theWisdom–Holman (1991) method, such that velocity dependent forces can be included, has been outlined. This is done both by using the δfunction formalism and also by means of formal Hamiltonization. Application to the relativistic corrections in Solar System integrations is discussed as an example. Numerical experiments confirm the usefulness of the method. This revised version was published online in July 2006 with corrections to the Cover Date.     

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.     

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.     

Propeller effect during magnetocentrifugal plasma acceleration

We study the effect of magnetic-field axial asymmetry on the magnetocentrifugal acceleration of plasma when it flows in a source’s rotating magnetosphere (propeller effect). For an axisymmetric steady plasma flow, the first corrections to the energy that arise when the source rotates slowly are proportional to Ω⁴, suggesting a highly inefficient plasma acceleration. Magnetic-field axial asymmetry is shown to substantially modify the acceleration. The first corrections arise even in the first order in Ω. The plasma acceleration turns out to be considerably more efficient in a nonaxisymmetric magnetic field.     

A Method Solving Kepler’s Equation for Hyperbolic Case

We developed a method to solve Keplers equation for the hyperboliccase. The solution interval is separated into three regions; F 1, F 1, F 1. For the region F is large, we transformed the variablefrom F to exp F and applied the Newton method to the transformed equation.For the region F is small, we expanded the equation with respect to F andapplied the Newton method to the approximated equations. For the middleregion, we adopted a discretization of the Newton method and used the Taylorseries expansion for the evaluation of transcendental functions (Fukushima,1997). Numerical measurements showed that, in the case of Intel Pentiumprocessor, the new method is more than 3.7 times as fast as the combinationof a fourth order correction formula (Fukushima, 1997) and Roysstarting procedure (Roy, 1988) and others.     

A mapping method for the gravitational few-body problem with dissipation

Recently a new class of numerical integration methods — mixed variable symplectic integrators — has been introduced for studying long-term evolution in the conservative gravitational few-body problem. These integrators are an order of magnitude faster than conventional ODE integration methods. Here we present a simple modification of this method to include small non-gravitational forces. The new scheme provides a similar advantage of computational speed for a larger class of problems in Solar System dynamics.     

Secular Dynamics of Asteroids in the Inner Solar System

In the last three years we have carried out numerical and semi-analytical studies on the secular dynamical mechanisms in the region (semimajor axis a < 2 AU) where the NEA orbits evolve. Our numerical integrations (over a time span of a few Myr) have shown that: (i) the linear secular resonances with both the inner and the outer planets may play an important role in the dynamical evolution of NEAs; (ii) the apsidal secular resonance with Mars could provide an important dynamical transport mechanism by which asteroids in the Mars-crossing region eventually achieve Earth-crossing orbits; (iii) in this region, due to the interaction with the terrestrial planets, the Kozai resonance can occur at small inclinations, with the argument of perihelion ω librating around 0° or 180°, providing a temporary protection mechanism against close approaches to the planets. The location of the linear secular resonances in this zone has also been obtained by an automatic procedure using a semi-numerical method valid for all values of the inclinations and eccentricities of the small bodies, and also in the case of libration of the argument of perihelion. A map of the secular resonances in the (a, i) plane shows — in agreement with the numerical integrations — that all the resonances with the terrestrial and giant planets are present, and also that some of them overlap. Thus the way is now open to fully take into account secular resonances in modelling the dynamical evolution of NEAs. This revised version was published online in July 2006 with corrections to the Cover Date.     

Evolution of Dust Shells and Jets in the Inner Coma of Comet C/1995 O1 (Hale-Bopp)

Comet Hale-Bopp was observed with the 80 cm reflector + CCD at the Haute-Provence observatory (OHP) and with the 62 cm reflector + CCD at the Saint-Véran observatory (Queyras, France). The morphology of the shells was followed from their first appearence on 1997 Jan. 30, until their disappearance on May 9. These shells spread from the nucleus region with a velocity in agreement with a nuclear rotation period of about 11.33 hours. We report also a short and bright dust ejection on May 8. CN images show a long spiral jet in the tailward side invisible on continuum images. The circumnuclear structures have been followed at Saint-Véran from Apr. 5 to Apr. 11, 1997 with a high spatial resolution (200 km/pixel). We have followed the emergence of a recurrent linear polar jet. Measurements of its expansion show a constant acceleration of material with typical expansion velocity of 1 km/s. The CCD frames show the interconnection between spiral jets and the successive shells. This revised version was published online in July 2006 with corrections to the Cover Date.