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.     

Recent developments of integrating the gravitational problem ofN-bodies

This paper discusses the formulation and the numerical integration of large systems of differential equations occurring in the gravitational problem ofn-bodies.Different forms of the pertinent differential equations of motion are presented, and various regularizing and smoothing transformations are compared. Details regarding the effectiveness and the efficiency of the Kustaanheimo-Stiefel and of other methods are discussed. In particular, a method is described in which some of the phase variables are treated in the regularized system and others in the ordinary system. This mixed method of numerical regularization offers some advantages.Several numerical integration techniques are compared. A high order Runge-Kutta method, Steffensen's method, and a finite difference method are investigated, especially with regard to their adaptability to regularization.The role of integrals and integral invariants is displayed in controlling the accuracy of the numerical integration.Numerical results are described with 5, 25 and 500 bodies participating. These examples compare the various integration techniques, several regularization methods and different logics in treating binaries.     

Numerical analysis of the symmetric methods

Aimed at the initial value problem of the particular second-order ordinary differential equations,y =f(x, y), the symmetric methods (Quinlan and Tremaine, 1990) and our methods (Xu and Zhang, 1994) have been compared in detail by integrating the artificial earth satellite orbits in this paper. In the end, we point out clearly that the integral accuracy of numerical integration of the satellite orbits by applying our methods is obviously higher than that by applying the same order formula of the symmetric methods when the integration time-interval is not greater than 12000 periods.     

Expression of Cassini's third law for Callisto, and theory of its rotation

The rotation of the main natural satellites of the Solar System is widely assumed to be synchronous, because this corresponds to an equilibrium state. In the case of the Moon, 3 laws have been formulated by Cassini, assuming a spin-orbit resonance and a 1:1 nodal resonance. The recent gravitational data collected by the spacecrafts Galileo (in the jovian system) and Cassini (in the saturnian system) allows us to study the rotation of other natural satellites, and to check the universality of Cassini's laws. This paper deals with the rotation of the Galilean satellites of Jupiter J-4 Callisto. In this study we use both analytical (like Lie transforms) and numerical methods (numerical detection of chaos, numerical integration, frequency analysis) to first check the reliability of Cassini Laws for Callisto, and then to give a first theory of its rotation, Callisto's being considered as a rigid body. We first show that the Third Cassini Law (i.e. the nodal resonance), is not satisfied in every reference frame, in particular in the most natural one (i.e. the J2000 jovian equator). The difference of the nodes presents a chaotic-like behavior, that we prove to be just a geometrical illusion. Moreover, we give a mathematical condition ruling the choice of an inertial reference frame in which the Third Cassini Law is fulfilled. Secondly, we give a theory of Callisto's rotation in the International Celestial Reference Frame (ICRF). We highlight a small motion (i.e. <200 m) of its rotation axis about its body figure, a 11.86-yr periodicity in Callisto's length-of-day, and the proximity of a resonance that forces 182-yr librations in Callisto's obliquity.     

Treatment of close approaches in the numerical integration of the gravitational problem ofN bodies

One of the main difficulties encountered in the numerical integration of the gravitationaln-body problem is associated with close approaches. The singularities of the differential equations of motion result in losses of accuracy and in considerable increase in computer time when any of the distances between the participating bodies decreases below a certain value. This value is larger than the distance when tidal effects become important, consequently,numerical problems are encounteredbefore the physical picture is changed. Elimination of these singularities by transformations is known as the process of regularization. This paper discusses such transformations and describes in considerable detail the numerical approaches to more accurate and faster integration. The basic ideas of smoothing and regularization are explained and applications are given.     

Enhancement of on-board forecasting accuracy of GEO SC COG motion using the compensative transversal acceleration

A method is suggested for enhancing the on-board forecasting accuracy of the COG motion of a GEO SC with a long time of independent operation. The suggested method consists of introducing so-called compensative transversal acceleration (CTA), along with zonal harmonics into the right sides of the differential equations of SC motion among other disturbances due to the Earth’s gravitational field eccentricity. The CTA compensates the integral effect of the sectoral and tesseral harmonics; its value is constant for a specified point of GEO SC location (standing point) and is calculated on the Earth from numerical integration of differential equations of motion taking into account the complete set of gravitational field harmonics. The CTA value is transmitted on-board of an SC as program command data. The method is implemented in algorithms of on-board forecasting of Electro-L SC motion and can be used to enhance the on-board forecasting accuracy of the COG motion of GEO SCs with a long time of independent operation.     

Explicit algorithmic regularization in the few-body problem for velocity-dependent perturbations

A new algorithm is presented for the numerical integration of second-order ordinary differential equations with perturbations that depend on the first derivative of the dependent variables with respect to the independent variable; it is especially designed for few-body problems with velocity-dependent perturbations. The algorithm can be used within extrapolation methods for enhanced accuracy, and it is fully explicit, which makes it a competitive alternative to standard discretization methods.     

A Time-Transformed Leapfrog Scheme

We present a time-transformed leapfrog scheme combined with the extrapolation method to construct an integrator for orbits in N-body systems with large mass ratios. The basic idea can be used to transform any second-order differential equation into a form which may allow more efficient numerical integration. When applied to gravitating few-body systems this formulation permits extremely close two-body encounters to be considered without significant loss of accuracy. The new scheme has been implemented in a direct N-body code for simulations of super-massive binaries in galactic nuclei. In this context relativistic effects may also be included.     

Studies in the application of recurrence relations to special perturbation methods

The numerical integration of the differential equations describing dynamical systems has been shown in previous papers of this series to be most effectively accomplished by an explicit Taylor series method.In this paper we show that one explicit Taylor series method, developed earlier in this series and which appears to possess a high degree of versatility, yields considerable gains in efficiency over classical single-step and multi-step methods. (In this context efficiency is a measure of the time taken to carry out a calculation of a specific accuracy).For a given accuracy criterion governing the local truncation error (LTE) it is found that the Taylor series method is generallytwice as fast as the classical multi-step method and up totwenty times faster than the classical single-step method.     

The symplectic character of periodic orbits in the three dipole problem

In this work we reveal for the first time that in the three dipole problem only asymmetric periodic orbits exist.For these periodic orbits — planar and three dimensional — of a charged particle moving under the influence of the electromagnetic field of the three dipoles we give their symplectic relations using the Hamiltonian formulation which is related to the symplectic matrix. Also we study the properties of the symplectic matrix and we give the relations there are among the variations of a periodic solution. These relations have been used to check the accuracy of numerical integration of equations of first order variations.     

The haloing effect of the third integral

Whether Contopoulos's galactic system is separable (unlikely) or not (likely), the fact is that there exists a vicinity of the equilibrium in which numerical integration of high accuracy cannot separate the system from its image through Birkhoff's normalization of high order. To all practical purposes, stellar dynamics is then justified in pretending that the model is, in that region, structured by a so-called third integral.     

天文动力学方程数值积分中的一种有效变步法

利用积分曲线的曲率控制步长的技巧，使天文动力学方程数值解法的精度和速度有较大提高，这种方法适用于天体精密定轨以及一些精度要求高的常微分方程初值问题的数值积分。     

A multiple shooting approach for the numerical treatment of stellar structure and evolution

We present a new numerical method for solving the system of partial differential equations describing the structure and evolution of a spherically symmetric star. As usual, we employ the transversal method of lines in order to split the equations into a coupled spatial and temporal part. The novel features of the algorithm are the following: (a) Instead of using the Lagrangian picture we formulate the system of partial differential equations in the Eulerian picture. (b) We reformulate the equations of stellar structure as a multipoint boundary-value problem. By means of this reformulation the rather clumsy iterative matching procedure of stellar atmosphere and interior is avoided. (c) The multipoint boundary-value problem is solved by the multiple shooting method. This approach not only ensures a high accuracy of the stellar models calculated at each time step but also allows the free boundaries inside the star due to different energy transport mechanisms to be located exactly. (d) The time derivatives involved in the stellar-structure equations are discretized implicitly to second order accuracy. Moreover, at each time step, the chemical abundances are determined by using a sophisticated update procedure. In this way, a high accuracy is achieved with respect to the integration in time. The algorithm has turned out to be exceedingly reliable and numerically accurate. This is shown by the evolution of a 1 M_⊙ star up to the hydrogen-shell burning phase. In this example, the virial theorem, the law of mass conservation, and the law of energy conservation is fulfilled to a hitherto unattainable degree of accuracy. Since the multiple shooting method, which is at the heart of our approach, is a perfect example of a parallel algorithm, the computational speed of the algorithm might be substantially improved provided easy-to-program, high-performance parallel computers with sufficiently many processors become available in the near future.     

Integral surfaces with space-time coordinates,in the gravitational field of a rotating system

A new integration theory is formulated for dynamical systems with two degrees of freedom, in the gravitational field of a rotating system. Four integrals of motion may be determined from complete solutions of a system of three first-order, partial differential equations in three independent variables. The solutions of this system define two integral surfaces with space-time coordinates. These surfaces represent two independent solutions of a second-order kinematic system to which the original fourth-order system has been reduced. An integral curve may be represented as the locus of intersection points of the integral surfaces. The new theory is the theoretical basis for a method of analytic continuation of periodic orbits of the circular restricted problem.     

Non-Existence of the Modified First Integral by Symplectic Integration Methods II: Kepler Problem

Symplectic integration methods conserve the Hamiltonian quite well because of the existence of the modified Hamiltonian as a formal conserved quantity. For a first integral of a given Hamiltonian system, the modified first integral is defined to be a formal first integral for the modified Hamiltonian. It is shown that the Runge-Lenz vector of the Kepler problem is not well conserved by symplectic methods, and that the corresponding modified first integral does not exist. This conclusion is given for a one-parameter family of symplectic methods including the symplectic Euler method and the Störmer/Verlet method.     

Integration error over very long time spans

The main limit to the time span of a numerical integration of the planetary orbits is no longer set by the availability of computer resources, but rather by the accumulation of the integration error. By the latter we mean the difference between the computed orbit and the dynamical behaviour of the real physical system, whatever the causes. The analysis of these causes requires an interdisciplinary effort: there are physical model and parameters errors, algorithm and discretisation errors, rounding off errors and reliability problems in the computer hardware and system software, as well as instabilities in the dynamical system. We list all the sources of integration error we are aware of and discuss their relevance in determining the present limit to the time span of a meaningful integration of the orbit of the planets. At present this limit is of the order of 10⁸ years for the outer planets. We discuss in more detail the truncation error of multistep algorithms (when applied to eccentric orbits), the coefficient error, the method of Encke and the associated coordinate change error, the procedures used to test the numerical integration software and their limitations. Many problems remain open, including the one of a realistic statistical model of the rounding off error; at present, the latter can only be described by a semiempirical model based upon the simpleN ² formula (N=number of steps, =machine accuracy), with an unknown numerical coefficient which is determined only a posteriori.     

On error bounds and initialization in satellite orbit theories

The order of magnitude of the error is investigated for a first-order von Zeipel theory of satellite orbits in an axisymmetric force field, i.e., first-order long period and short-period effects are included along with second order secular rates. The treatment is valid for zero eccentricity and/or inclination. In the case where initial position and velocity vectors are known, the in-track position error over time intervals of order 1/J ₂ is kept at 0(J ₂ ²), like the other position errors and velocity errors, by calibration of the mean motion with the aid of the energy integral. The results are specifically applicable to accuracy comparisons of the Brouwer orbit prediction method with numerical integration. A modified calibration is presented for the general asymmetric force field which includes tesseral harmonics.     

A velocity scaling method with least-squares correction of several constraints

A new scale transformation to the integrated velocity vector is designed to monitor the accumulation of numerical errors in several integrals of motion. The scale factor is derived from the least-squares correction that minimizes the sum of the squares of the errors of these integrals. In order to preserve an invariant, we employ the velocity scaling method for rigorously satisfying the constraint. When adjusting many constants, the new scheme like other existing methods is valid to typically reduce the integration errors below those of an uncorrected integrator. Via integral invariant relations, the new method is also able to treat slowly-varying quantities, such as the Keplerian energy and the Laplace vector, for a perturbed Keplerian problem or each of multiple bodies in the solar system dynamics. Consequently it does nearly agree with the rigorous dual scaling method in the sense of drastically improving the integration accuracy. As one of its advantages, the implementation of the new method is significantly easier than that of other methods. In particular, the method can be simply applied to a complicated dynamical system with some constraints.     

A note on obtaining an approximate value for the period of a periodic solution of\ddot x = - V\prime (x)

The first integral of yields an integral for the period of a periodic solution, if such exists. In general, this integral cannot be evaluated. By means of an approximate solution along with the minimization of a mean-square error, one can obtain an approximate value for the period in terms of the amplitude of the motion. The calculated period agrees very well with the period obtained by means of numerical integration for the case of orbit-orbit resonance involving the motion of two satellites of a planet.The same method is applied to obtain an alternative derivation of the first Krylov-Bogoliuboff averaging method in non-linear mechanics.     

Stabilization by manipulation of the Hamiltonian

Numerical integration of unstable differential equations should be avoided since a numerical error during thenth step produces erroneous initial values for the next step and thus deteriorates the subsequent integration in an unstable manner. A method is offered to stabilize the equations of motion corresponding to a given HamiltonianH by transformingH into a new HamiltonianH ^* which is equivalent to the Hamiltonian of a harmonic oscillator. In contrast to other methods of stabilization the realm of canonical mechanics is thus not abandoned. Perturbations are discussed and as examples the Keplerian motion and the motion of a gyroscope are presented.