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.     

Jacobi Dynamics And The N-Body Problem With Variable Masses

An appropriate generalization of the Jacobi equation of motion for the polar moment of inertia I is considered in order to study the N-body problem with variable masses. Two coupled ordinary differential equations governing the evolution of I and the total energy E are obtained. A regularization scheme for this system of differential equations is provided. We compute some illustrative numerical examples, and discuss an average method for obtaining approximate analytical solutions to this pair of equations. For a particular law of mass loss we also obtain exact analytical solutions. The application of these ideas to other kind of perturbed gravitational N-body systems involving drag forces or a different type of mass variation is also considered. This revised version was published online in July 2006 with corrections to the Cover Date.     

Limits of predictability of gravitational systems

The concept of finite predictability of gravitational many-body systems is related to the non-deterministic nature of celestial mechanics and of dynamics, in general. The basic, fundamental reasons for the uncertainty of predictions are as follows: (1) the initial conditions are known only approximately since they are obtained either from observations or from approximate computations; (2) the equations of motion given by a selected model describe the actual system only approximately; (3) the physical constants of the dynamical system have error limits; (4) the differential equations of motion are non-integrable and numerical integration methods must be used for solution, generating errors in the final result at every integration step.In addition to these reasons, mostly depending on our techniques, there are some more fundamental reasons depending on the nature of the dynamical system investigated. These are the appearance of regions of instability, non-integrability and chaotic motion.Details, effects and controls of these regions for finite predictability are discussed for various dynamical systems of importance in celestial mechanics with special emphasis on planetary systems.     

Regularization of Motion Equations with L-Transformation and Numerical Integration of the Regular Equations

The sets of L-matrices of the second, fourth and eighth orders are constructed axiomatically. The defining relations are taken from the regularization of motion equations for Keplerian problem. In particular, the Levi-Civita matrix and KS-matrix are L-matrices of second and fourth order, respectively. A theorem on the ranks of L-transformations of different orders is proved. The notion of L-similarity transformation is introduced, certain sets of L-matrices are constructed, and their classification is given. An application of fourth order L-matrices for N-body problem regularization is given. A method of correction for regular coordinates in the Runge–Kutta–Fehlberg integration method for regular motion equations of a perturbed two-body problem is suggested. Comparison is given for the results of numerical integration in the problem of defining the orbit of a satellite, with and without the above correction method. The comparison is carried out with respect to the number of calls to the subroutine evaluating the perturbational accelerations vector. The results of integration using the correction turn out to be in a favorable position.     

Special perturbations employing osculating reference states

The concept of employing osculating reference position and velocity vectors in the numerical integration of the equations of motion of a satellite is examined. The choice of the reference point is shown to have a significant effect upon numerical efficiency and the class of trajectories described by the differential equations of motion. For example, when the position and velocity vectors on the osculating orbit at a fixed reference time are chosen, a universal formulation is yielded. For elliptical orbits, however, this formulation is unattractive for numerical integration purposes due to Poisson terms (mixed secular) appearing in the equations of motion. Other choices for the reference point eliminate this problem but usually at the expense of universality. A number of these formulations, including a universal one, are considered here. Comparisons of the numerical characteristics of these techniques with those of the Encke method are presented.     

Application of the integral variation method to satellite orbit prediction

The Integral Variation (IV) method is a technique to generate an approximate solution to initial value problems involving systems of first-order ordinary differential equations. The technique makes use of generalized Fourier expansions in terms of shifted orthogonal polynomials. The IV method is briefly described and then applied to the problem of near Earth satellite orbit prediction. In particular, we will solve the Lagrange planetary equations including the first three zonal harmonics and drag. This is a highly nonlinear system of six coupled first-order differential equations. Comparison with direct numerical integration shows that the IV method indeed provides accurate analytical approximations to the orbit prediction problem.Advanced Systems Studies; Bldg. 254EElectro-Optical Systems Laboratory; Bldg. 201.     

Numerical integration of dynamical systems with Lie series

The integration of the equations of motion in gravitational dynamical systems—either in our Solar System or for extra-solar planetary systems—being non integrable in the global case, is usually performed by means of numerical integration. Among the different numerical techniques available for solving ordinary differential equations, the numerical integration using Lie series has shown some advantages. In its original form (Hanslmeier and Dvorak, Astron Astrophys 132, 203 1984), it was limited to the N-body problem where only gravitational interactions are taken into account. We present in this paper a generalisation of the method by deriving an expression of the Lie terms when other major forces are considered. As a matter of fact, previous studies have been done but only for objects moving under gravitational attraction. If other perturbations are added, the Lie integrator has to be re-built. In the present work we consider two cases involving position and position-velocity dependent perturbations: relativistic acceleration in the framework of General Relativity and a simplified force for the Yarkovsky effect. A general iteration procedure is applied to derive the Lie series to any order and precision. We then give an application to the integration of the equation of motions for typical Near-Earth objects and planet Mercury.     

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.     

Accurate computation of highly eccentric satellite orbits

For computing highly eccentric (e0.9) Earth satellite orbits with special perturbation methods, a comparison is made between different schemes, namely the direct integration of the equations of motion in Cartesian coordinates, changes of the independent variable, use of a time element, stabilization and use of regular elements. A one-step and a multi-step integration are also compared.It is shown that stabilization and regularization procedures are very helpful for non or smoothly perturbed orbits. In practical cases for space research where all perturbations are considered, these procedures are no longer so efficient. The recommended method in these cases is a multi-step integration of the Cartesian coordinates with a change of the independent variable defining an analytical step size regulation. However, the use of a time element and a stabilization procedure for the equations of motion improves the accuracy, except when a small step size is chosen.     

Extrapolation of symplectic Integrators

We build high order numerical methods for solving differential equations by applying extrapolation techniques to a Symplectic Integrator of order 2n. We show that, in general, the qualitative properties are preserved at least up to order 4n+1. This new procedure produces much more efficient methods than those obtained using the Yoshida composition technique. This revised version was published online in July 2006 with corrections to the Cover Date.     

Direct integration methods of theN-body problem

A fourth-order polynomial method for the integration ofN-body systems is described in detail together with the computational algorithm. Most particles are treated efficiently by an individual time-step scheme but the calculation of close encounters and persistent binary orbits is rather time-consuming and is best performed by special techniques. A discussion is given of the Kustaanheimo-Stiefel regularization procedure which is used to integrate dominant two-body encounters as well as close binaries. Suitable decision-making parameters are introduced and a simple method is developed for regularizing an arbitrary number of simultaneous two-body encounters.     

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.     

Regularization of Equations of Motion of Celestial Bodies. 2. Correction of the KS-Coordinate System in Numerical Integration

The idea of using various L-matrices in numerical integration of the regular equations, which describe the motion of small bodies of the Solar System, is developed. The problem of the optimal position of the radius vector and velocity at numerical integration in the KS-coordinate system is posed. The solution of this problem, which reduces the number of calculations of the vector of perturbing accelerations, is given. The transformation providing this optimal solution is suggested, and the results of numerical integration are given.     

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.     

A slow-down treatment for close binaries

We present a new method for fast numerical integration of close binaries inN-body systems. The basic idea is to slow down the motion of the binary artificially, which makes a faster numerical integration possible but still maintains correct treatment of secular and long-period effects on the motion. We discuss the general principle, with application to close binaries inN-body codes and in the chain regularization.     

Coronal Holes and Icosahedral Symmetry

The study of the expansion of the solar wind out of a system of coronal holes is continued. To this end, we consider the numerical integration of partial differential equations for problems with icosahedral symmetry, in general. First, employing Weyl theory, orbifold coordinates are introduced. Second, the icosahedral coordinates are discussed in detail. Third, following an analysis of the properties of these coordinates and the derivation of a few expressions useful for grid construction, various alternatives for the distribution of lattice points required for numerical integration are considered. A comparison of these numerical grids motivates the choice of a specific grid optimized for the numerical integration carried out in the accompanying paper by Kalish et al.(2002). This revised version was published online in July 2006 with corrections to the Cover Date.     

Numerical stabilization of all laws of conservation in the many body problem

When a system of differential equations admits a first integral (e.q. the law of energy), the value of that integral may be used as a check during the numerical integration. Often this check is satisfied with poor accuracy since the existence of the first integral is unknown to the computer. The aim of the paper is to show how such a first integral can be satisfied with better accuracy and in a stabilized manner by adding an appropriate control term to the differential system. The accuracy of the numerical integration is thereby improved. The basic idea is applied to the problem ofN bodies.Presented at the Conference on Celestial Mechanics, Oberwolfach, Germany, August 27–September 2, 1972.     

轨道改进中计算状态转移矩阵的分析方法

对当今人卫轨道改进问题,由于力学模型的复杂,精密星历和状态转移矩阵的计算均采用数值方法,这就需要积分两组常微分方程．本文针对状态转移矩阵在定轨中的作用,对定轨弧段不太长的情况,给出了状态转移矩阵的一种分析算法,从而避免数值求解两组常微分方程的问题,并以实际算例证实了这种算法的有效性     

Radiative Transfer in Inhomogeneous Atmospheres. I

This series of papers is devoted to multiple scattering of light in plane parallel, inhomogeneous atmospheres. The approach proposed here is based on Ambartsumyan's method of adding layers. The main purpose is to show that one can avoid difficulties with solving various boundary value problems in the theory of radiative transfer, including some standard problems, by reducing them to initial value problems. In this paper the simplest one dimensional problem of diffuse reflection and transmission of radiation in inhomogeneous atmospheres with finite optical thicknesses is considered as an example. This approach essentially involves first determining the reflection and transmission coefficients of the atmosphere, which, as is known, are a solution of the Cauchy problem for a system of nonlinear differential equations. In particular, it is shown that this system can be replaced with a system of linear equations by introducing auxiliary functions P and S. After the reflectivity and transmissivity of the atmosphere are determined, the radiation field in it is found directly without solving any new equations. We note that this approach can be used to obtain the required intensities simultaneously for a family of atmospheres with different optical thicknesses. Two special cases of the functional dependence of the scattering coefficient on the optical thickness, for which the solutions of the corresponding equations can be expressed in terms of elementary functions, are examined in detail. Some numerical calculations are presented and interpreted physically to illustrate specific features of radiative transport in inhomogeneous atmospheres.     

A theory of integration for the extended Delaunay method

In this paper a method for the integration of the equations of the extended Delaunay method is proposed. It is based on the equations of the characteristic curves associated with the partial differential equation of Delaunay-Poincaré. The use of the method of characteristics changes the partial differential equation for higher order approximations into a system of ordinary differential equations. The independent variable of the equations of the characteristics is used instead of the angular variables of the Jacobian methods and the averaging principle of Hori is applied to solve the equations for higher orders. It is well known that Jacobian methods applied to resonant problems generally lead to the singularity of Poincaré. In the ideal resonance problem, this singularity appears when higher order approximations of the librational motion are considered. The singularity of Poincaré is non-essential and is caused by the choice of the critical arguments as integration variables. The use of the independent variable of the equation of the characteristics in the place of the critical angles eliminates the singularity of Poincaré.