高斯型龙格库塔积分器

提出并建立了一种新型的常微分方程数值解法－高斯型隐式Runge-Kutta算法,该算法具有κ级2κ－1阶的代数精确度。数值结果表明,高期型隐式Runge-Kutta算法比传统显示Runge-Kutta算法精度高,稳定性好。     

The photogravitational Hill problem with oblateness: equilibrium points and Lyapunov families

We introduce a new version of Hill’s problem that incorporates the effects of radiation of the primary and oblateness of the secondary and study the basic dynamical features of this new model-problem. This formulation is more appropriate for some astronomical applications as an approximation to the corresponding restricted three-body problem. We use iterative methods for deriving approximate expressions of the equilibrium point locations and study their stability properties by using a linear stability analysis. All equilibrium points are unstable. We also employ singular perturbations methods for obtaining approximate expressions of the Lyapunov families emanating from equilibrium points, in both coplanar and spatial case, and numerical techniques for their continuation.     

稳定化方法与后稳定化方法的比较

对BaumgaLrte的稳定化和Chin的后稳定化进行了详尽讨论与数值比较．用经典数值方法并结合这两种稳定化方式都能提高数值精度和改善数值稳定性．在最佳稳定参数下稳定化精度一般不等价于后稳定化．两者精度优劣并无常定．考虑到Baumgarte的稳定化使得数值积分的右函数更复杂和增加计算耗费,尤其是存在稳定参数最佳选取的麻烦,故推荐后稳定化投入实算．但值得注意的是用后稳定化与没有经过稳定化处理的经典积分器来比不宜扩大积分步长．     

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.     

Secular variations in the restricted problem: A numerical technique for studying stability of equilibrium solutions

In this paper we develop and implement an explicit numerical technique for studying stability of equilibrium solutions concerning the secular variations in the restricted problem. In our implementation, we use fourth-order expansions for the secular terms, but the method itself is independent of the particular order used as upper limit in the required expansions.     

关于数值求解天体运动方程的几个问题

本文讨论三个问题：1.在采用各种非辛（Symplectic）的数值积分器积分天体运动方程时，截断误差将引起人为的能量耗散，这一问题是不能用简单地在相应的力模型中加进一个人为的阻力因子而得以解决的，被歪曲的能量（或数值轨道）必须在积分过程的每一步用能量关系来进行校正，此即能量控制方法．2.当摄动加速度涉及到坐标轴的旋转时，如何在各种积分器中采用能量控制方法．3.对于大偏心率轨道，用数值方法求解相应运动方程时，积分步长必须随运动天体与中心天体之间的距离变化而改变，显然，这对所有积分器都是不方便的，特别是多步积分器．本义给出了一种步长均匀化的处理，可以使上述大偏心率轨道积分问题按定步长计算．     

Constraints and numerical integration

In an autonomous Hamiltonian system, one constraint always exists, namely the energy integral or a constant magnitude of the 4-velocity in relativistic dynamics. The constraint should bring better numerical stability if it can be kept at every step of the numerical integration. In Newtonian mechanics, the order of the equations of motion can not be reduced by use of the constraint in most cases, because its kinetic energy is usually a quadratic form of the elliptic type and one would meet difficulty when trying the order reducing. However, the metric in general relativity is hyperbolic. In particular, when the spacetime bears some symmetries there exists a global transformation so that at least one element in its main diagonal vanishes. As a result, the constraint can be solved for a certain velocity or a momentum without any difficulty, and so reducing the order. Similarly, this technique can also be applied to the evolution of the Mixmaster universe. It is shown that this technique can raise the precision and improve the numerical stability dramatically even when a classical integrator is used, although it might not keep the symplectic structure of the system.     

利用正交变换方法计算协方差分析的统计量

证明用Givens—Gentleman正交变换给出的加仅最小二乘解与统计定轨理论求得的解一致．采用正交变换方法计算其它的一些重要的统计量,如考察协方差矩阵、摄动矩阵等,并计算这些量随时间的传播．这种算法的优点是通过降低法方程的条件数提高计算的稳定性,同时可以方便地对不同的参数组合情况求解而不需多次解算法方程．     

Expansion of the planetary disturbing function

Some methods are described for the expansion of the disturbing function in planetary theory. One method uses the classical binomial expansion theorem or a successive approximation process derived from it. Another method is a direct application of the Laplace series expansions. For both methods it is proposed to first prepare the series to be manipulated by a scaling operation. These methods can be applied either in a literal or in a numerical form, or any combination of both, but they are especially designed for use on a large scale digital computer with standard Poisson series programs. No usage is made of Newcomb operators or derivatives of Laplace coefficients.     

Bandlimited implicit Runge–Kutta integration for Astrodynamics

We describe a new method for numerical integration, dubbed bandlimited collocation implicit Runge–Kutta (BLC-IRK), and compare its efficiency in propagating orbits to existing techniques commonly used in Astrodynamics. The BLC-IRK scheme uses generalized Gaussian quadratures for bandlimited functions. This new method allows us to use significantly fewer force function evaluations than explicit Runge–Kutta schemes. In particular, we use a low-fidelity force model for most of the iterations, thus minimizing the number of high-fidelity force model evaluations. We also investigate the dense output capability of the new scheme, quantifying its accuracy for Earth orbits. We demonstrate that this numerical integration technique is faster than explicit methods of Dormand and Prince 5(4) and 8(7), Runge–Kutta–Fehlberg 7(8), and approaches the efficiency of the 8th-order Gauss–Jackson multistep method. We anticipate a significant acceleration of the scheme in a multiprocessor environment.     

Detection of Ordered and Chaotic Motion Using the Dynamical Spectra

Two simple and efficient numerical methods to explore the phase space structure are presented, based on the properties of the "dynamical spectra". 1) We calculate a "spectral distance" D of the dynamical spectra for two different initial deviation vectors. D → 0 in the case of chaotic orbits, while D → const ≠ 0 in the case of ordered orbits. This method is by orders of magnitude faster than the method of the Lyapunov Characteristic Number (LCN). 2) We define a sensitive indicator called ROTOR (ROtational TOri Recongnizer) for 2D maps. The ROTOR remains zero in time on a rotational torus, while it tends to infinity at a rate ∝ N = number of iterations, in any case other than a rotational torus. We use this method to locate the last KAM torus of an island of stability, as well as the most important cantori causing stickiness near it. This revised version was published online in July 2006 with corrections to the Cover Date.     

Linear stability for some symmetric periodic simultaneous binary collision orbits in the four-body problem

We apply the analytic-numerical method of Roberts to determine the linear stability of time-reversible periodic simultaneous binary collision orbits in the symmetric collinear four-body problem with masses 1, m, m, 1, and also in a symmetric planar four-body problem with equal masses. In both problems, the assumed symmetries reduce the determination of linear stability to the numerical computation of a single real number. For the collinear problem, this verifies the earlier numerical results of Sweatman for linear stability with respect to collinear and symmetric perturbations.     

Error propagation in the numerical solutions of the differential equations of orbital mechanics

The relationship between the eigen values of the linearized differential equations of orbital mechanics and the stability characteristics of numerical methods is presented. It is shown that the Cowell, Encke, and Encke formulation with an independent variable related to the eccentric anomaly all have a real positive eigen value when linearized about the initial conditions. The real positive eigen value causes an amplification of the error of the solution when used in conjunction with a numerical integration method. In contrast an element formulation has zero eigen values and is numerically stable.     

Numerical Instruments for the Analysis of Secular Dynamics of Exoplanetary Systems

A review is given of modern numerical methods for the analysis of resonant and chaotic dynamics: calculation of the Lyapunov characteristic exponents, the MEGNO method, and the maximum eccentricity method. These methods are used to construct stability diagrams for the planetary systems γ Cep, HD 196885, and HD 41004. The diagrams are analyzed to determine the most probable values taken by the orbital parameters of the exoplanets and obtain estimates for the Lyapunov time of their orbital dynamics. The stability diagrams constructed using the different methods are compared to analyze their effectiveness in the study of secular dynamics of exoplanetary systems.     

约束条件和数值积分

自治的哈密顿系统存在约束条件，例如能量积分或广义相对论中的4速度大小为常数，它能否在数值积分过程中始终满足将直接影响数值稳定性．在牛顿力学中哈密顿系统的动能一般为椭圆型，直接运用约束条件对方程进行降阶存在开平方判断正负号的困难，导致应用高精度的经典数值积分器时能量存在耗散．然而相对论力学的度规为双曲型，利用约束条件有可能实行方程降阶．在时空具有一定对称性的情况下，能够找到整个时空的一个全局变换使变换后的度规的主对角线某一元素为零，于是从约束方程中不需开平方能够解出某一动量，顺利实现运动方程的降阶．相对论力学中另一个可以降阶的模型是Mixmaster宇宙模型．数值实验表明将经典算法用于降阶后的运动方程能够严格地满足约束，但不一定能保持辛结构。     

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.     

Smoothing methods comparison for CMB E-and B-mode separation

The anisotropies of the B-mode polarization in the cosmic microwave background radiation play a crucial role in the study of the very early Universe. However, in real observations, a mixture of the Emode and B-mode can be caused by partial sky surveys, which must be separated before being applied to a cosmological explanation. The separation method developed by Smith(2006) has been widely adopted,where the edge of the top-hat mask should be smoothed to avoid numerical errors. In this paper, we compare three different smoothing methods and investigate leakage residuals of the E-B mixture. We find that, if less information loss is needed and a smaller region is smoothed in the analysis, the sin- and cos-smoothing methods are better. However, if we need a cleanly constructed B-mode map, the larger region around the mask edge should be smoothed. In this case, the Gaussian-smoothing method becomes much better. In addition, we find that the leakage caused by numerical errors in the Gaussian-smoothing method is mostly concentrated in two bands, which is quite easy to reduce for further E-B separations.     

含辐射和扁率的圆型限制性三体问题的轨道稳定性研究

基于最小二乘法原理的速度因子方法是保流形结构算法中效率最高、稳定性最好、应用最广的方法.利用速度因子方法讨论了主星为辐射源,伴星为扁球的平面圆型限制性三体问题的稳定性问题.数值研究表明:(1)仅考虑扁状摄动项时,系统混沌运动的轨道数量会增多;(2)仅考虑辐射项时,系统有序运动的轨道数量会增多;(3)同时存在辐射和扁状摄动时,辐射占主导作用,系统有序运动的几率会增加.     

Algorithmic regularization of the few-body problem

Using a modified leapfrog method as a basic mapping, we produce a new numerical integrator for the stellar dynamical few-body problem. We do not use coordinate transformation and the differential equations are not regularized, but the leapfrog algorithm gives regular results even for collision orbits. For this reason, application of extrapolation methods gives high precision. We compare the new integrator with several others and find it promising. Especially interesting is its efficiency for some potentials that differ from the Newtonian one at small distances.