Numerical aspects of time elements

Numerical results have shown that the use of time elements with time transformations provides increased accuracy in the numerical solution of gravitational systems.To gain additional accuracy improvements, it appears that the time and the time element should be calculated from quantities that have been adjusted so as to satisfy the energy integral exactly.We also have found that by reducing the growth of the time element to being nearly linear rather than quadratic causes an increase in the magnitude of the local truncation error in the solution but with a decrease in the rate of growth of the truncation error.     

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

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

Analysis of stellar occultation data: Effects of photon noise and initial conditions

A new inversion technique for obtaining temperature, pressure, and number density profiles of a planetary atmosphere from an occultation light curve is described. This technique employs an improved boundary condition to begin the numerical inversion and permits the computation of errors in the profiles caused by photon noise in the light curve. We present our assumptions about the atmosphere, optics, and noise and develop the equations for temperature, pressure, and number density and their associated errors. By inverting in equal increments of altitude, Δh, rather than in equal increments of time, Δt, the inversion need not be halted at the first negative point on the light curve as required by previous methods. The importance of the boundary condition is stressed, and a new initial condition is given. Numerical results are presented for the special case of inversion of a noisy isothermal light curve. From these results, simple relations are developed which can be used to predict the noise quality of an occultation. It is found that fractional errors in temperature profiles are comparable to those of pressure and number density profiles. An example of the inversion method is shown. Finally, we discuss the validity of our assumptions. In an appendix we demonstrate that minimum fractional errors in scale height determined from the inversion are comparable to those from an isothermal fit to a noisy isothermal light curve.     

Behaviour of the SMF method for the numerical integration of satellite orbits

The SMF algorithms were recently developed by the authors as a multistep generalization of the ScheifeleG-functions one-step method. Like the last, the proposed codes integrate harmonic oscillations without truncation error and the perturbing parameter appears as a factor of that error when integrating perturbed oscillations. Therefore they seemed to be convenient for the accurate integration of orbital problems after the application of linearizing transformations, such as KS or BF. In this paper we present several numerical experiments concerning the propagation of Earth satellite orbits, that illustrate the performance of the the SMF method. In general, it provides greater accuracy than the usual standard algorithms for similar computational cost.     

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.     

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.     

Extreme sensitivity of the YORP effect to small-scale topography

Radiation recoil (YORP) torques are shown to be extremely sensitive to small-scale surface topography, using numerical simulations. Starting from a set of “base objects” representative of the near-Earth object population, random realizations of three types of small-scale topography are added: Gaussian surface fluctuations, craters, and boulders. For each, the expected relative errors in the spin and obliquity components of the YORP torque caused by the observationally unresolved small-scale topography are computed. Gaussian power, at angular scales below an observational limit, produces expected errors of order 100% if observations constrain the surface to a spherical harmonic order l?10. For errors under 10%, the surface must be constrained to at least l=20. A single crater with diameter roughly half the object's mean radius, placed at random locations, results in expected errors of several tens of percent. The errors scale with crater diameter D as D² for D>0.3 and as D³ for D<0.3 mean radii. Objects that are identical except for the location of a single large crater can differ by factors of several in YORP torque, while being photometrically indistinguishable at the level of hundredths of a magnitude. Boulders placed randomly on identical base objects create torque errors roughly 3 times larger than do craters of the same diameter, with errors scaling as the square of the boulder diameter. A single boulder comparable to Yoshinodai on 25143 Itokawa, moved by as little as twice its own diameter, can alter the magnitude of the torque by factors of several, and change the sign of its spin component at all obliquities. Most of the total torque error produced by multiple unresolved craters is contributed by the handful of largest craters; but both large and small boulders contribute comparably to the total boulder-induced error. A YORP torque prediction derived from groundbased data can be expected to be in error by of order 100% due to unresolved topography. Small surface changes caused by slow spin-up or spin-down may have significant stochastic effects on the spin evolution of small bodies. For rotation periods between roughly 2 and 10 h, these unpredictable changes may reverse the sign of the YORP torque. Objects in this spin regime may random-walk up and down in spin rate before the rubble-pile limit is exceeded and fissioning or loss of surface objects occurs. Similar behavior may be expected at rotation rates approaching the limiting values for tensile-strength dominated objects.     

Tidal torques on infrequently colliding particle disks in binary systems and the truncation of the asteroid belt

We study numerically and analytically the conditions leading to the truncation, at the 2:1 resonance, of a disk of infrequently colliding particles surrounding the primary of a binary system. We focus on systems with small mass ratios, q, such as the Sun-Jupiter system with q = 10^?3. Previous studies showed that if collisions are frequent with respect to the orbital period, truncation 3nly occurs if the Reynolds number is greater than q^?2. This corresponds to particle eccentricity, e, less than of order q for a particle disk of optical depth unity. In thepresent case collisions are less frequent than $\text{q}{}^{\mathrm{<mtext>?2</mtext><mtext>3</mtext>}}$ orbital periods (the period of the forced eccentricity at the 2:1 resonance), and truncation occurs and (Kirkwood) gaps are produced only if e is less than some critical value which we estimate to be of order $\text{q}{}^{\mathrm{<mtext>5</mtext><mtext>9</mtext>}}$, or ～0.02 for the Sun-Jupiter case. We mention several means whereby the eccentricities may have been subsequently increased to their observed values.     

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.     

Elements of spin motion

For use in numerical studies of rotational motion, a set of elements is introduced for the torque-free rotational motion of a rigid body around its barycenter. The elements are defined as the initial values of a modification of the Andoyer canonical variables. A computational procedure is obtained for determining these elements from the combination of the spin angular momentum vector and a triad defining the orientation of the rigid body. A numerical experiment shows that the errors of transformation between the elements and variables are sufficiently small. The errors increase linearly with time for some elements and quadratically for some others.     

Practical considerations in the numerical solution of theN-body problem by Taylor series methods

It is shown that in the numerical integration ofN-body problems, as much importance should be given to considerations of the computer programming language to be used as to questions of the accummulation of round-off and truncation error, the stability of the method chosen and the problem being treated. By careful programming processing time may be cut by a factor of 2 or 3 which is an important consideration in extended numerical investigations. The relative usefulness of differing strategies for determining the step size is discussed and in addition the usefulness is shown of treatingN-body problems by a Taylor series method.     

An optimal (m+3)[m+4] Runge Kutta algorithm

Recurrent power series methods are particularly applicable to problems in celestial mechanics since the Taylor coefficients may be expressed by recurrence relations. However, as the number of Taylor coefficients increases as is often necessary because of accuracy requirements, the computing time grows prohibitively large. In order to avoid this unfavorable situation, Dr E. Fehlberg introduced in 1960 Runge-Kutta methods that use the firstm Taylor coefficients obtained by recursive relations, or some other technique.Optimalm-fold Runge-Kutta methods are introduced. Embedded methods of order (m+3)[m+4] and (m+4)[m+5] are presented which have coefficients that produce minimum local truncation errors for the higher order pair of solutions of the method, as well as providing a near maximum absolute stability region. It is emphasized that the methods are formulated such that the higher order pair of solutions is to be utilized. These optimal methods are compared to the existingm-fold methods for several test problems. The numerical comparisons show that the optimal methods are more efficient. It is stressed that these optimal methods are particularly efficient whenm is small.     

The H N method for Milne problem with polarization

The H _N method, employed for studies in neutron transport theory, is used to establish numerical results basic to the vector equation describing the transfer of polarized light in a Rayleigh scattering atmosphere with true absorption. The method has been applied to the classical Milne problem. The exit distribution is defined as a series in powers of the zenith observation angle. The numerical results are computed and compared with exact values obtained using the exit distribution in terms of the H-matrix. The numerical results are in good agreement with previously published findings.     

Bivariate least squares linear regression: Towards a unified analytic formalism. I. Functional models

Concerning bivariate least squares linear regression, the classical approach pursued for functional models in earlier attempts (York, 1966, York, 1969) is reviewed using a new formalism in terms of deviation (matrix) traces which, for unweighted data, reduce to usual quantities leaving aside an unessential (but dimensional) multiplicative factor. Within the framework of classical error models, the dependent variable relates to the independent variable according to the usual additive model. The classes of linear models considered are regression lines in the general case of correlated errors in X and in Y for weighted data, and in the opposite limiting situations of (i) uncorrelated errors in X and in Y, and (ii) completely correlated errors in X and in Y. The special case of (C) generalized orthogonal regression is considered in detail together with well known subcases, namely: (Y) errors in X negligible (ideally null) with respect to errors in Y; (X) errors in Y negligible (ideally null) with respect to errors in X; (O) genuine orthogonal regression; (R) reduced major-axis regression. In the limit of unweighted data, the results determined for functional models are compared with their counterparts related to extreme structural models i.e. the instrumental scatter is negligible (ideally null) with respect to the intrinsic scatter (Isobe et al., 1990, Feigelson and Babu, 1992). While regression line slope and intercept estimators for functional and structural models necessarily coincide, the contrary holds for related variance estimators even if the residuals obey a Gaussian distribution, with the exception of Y models. An example of astronomical application is considered, concerning the [O/H]–[Fe/H] empirical relations deduced from five samples related to different stars and/or different methods of oxygen abundance determination. For selected samples and assigned methods, different regression models yield consistent results within the errors (?σ) for both heteroscedastic and homoscedastic data. Conversely, samples related to different methods produce discrepant results, due to the presence of (still undetected) systematic errors, which implies no definitive statement can be made at present. A comparison is also made between different expressions of regression line slope and intercept variance estimators, where fractional discrepancies are found to be not exceeding a few percent, which grows up to about 20% in the presence of large dispersion data. An extension of the formalism to structural models is left to a forthcoming paper.     

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

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

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.