这是用户在 2025-4-14 21:26 为 https://app.immersivetranslate.com/pdf-pro/4dfae176-5263-4b17-9849-cab4648ca527/?isTrial=true 保存的双语快照页面，由 沉浸式翻译 提供双语支持。了解如何保存？

Secular evolution of the solar system over 10 million years
太阳系一千万年来的长期演化

J. Laskar  J. 拉斯卡Service de Mécanique Céleste du Bureau des Longitudes, Equipe de recherche associée au CNRS, 77 Avenue Denfert-Rochereau, F-75014 Paris, France

Received August 4, accepted September 21, 1987
1987年8月4日收到，9月21日接受

Abstract  抽象的

Summary. We have computed the differential system giving the secular evolution of the 8 main planets of the solar system up to the order 2 with respect to the masses and degree 5 in ec − 5 in ec − 5inec-5 \mathrm{in} \mathrm{ec-} centricity and inclination including lunar and relativistic contributions (Laskar, 1985, 1986). This secular system is numerically integrated over 30 million years. A modified Fourier analysis is performed to obtain a solution for the secular evolution of the orbits on a quasi-periodical form. Comparison with Bretagnon’s ephemeris VSOP82 allows to derive uncertainties for the determination of the main frequencies of the secular system. These uncertainties range from 0.001 arcsec / year 0.001 arcsec / year 0.001arcsec//year0.001 \mathrm{arcsec} / \mathrm{year} to 0.02 arcsec / year 0.02 arcsec / year 0.02arcsec//year0.02 \mathrm{arcsec} / \mathrm{year}. Comparisons are made with the results of long term numerical integrations of Applegate et al. (1986) and Carpino et al., (1986) and with the analytical theory of Bretagnon (1974, 1984). The solutions of the outer solar system appear to be more stable than the solutions of the inner solar system where lots of quasiresonances in the long period frequencies prevent a good convergence of the solutions.
摘要。我们计算了微分系统，给出了太阳系八大行星在质量、 5 in ec − 5 in ec − 5inec-5 \mathrm{in} \mathrm{ec-} 度中心性和倾角方面的长期演化，包括月球和相对论的贡献，最高可达二阶（Laskar，1985，1986）。该长期系统在 3000 万年内进行了数值积分。采用改进的傅里叶分析，获得了准周期形式的轨道长期演化的解。与 Bretagnon 的星历表 VSOP82 进行比较，可以得出确定长期系统主要频率的不确定性。这些不确定性的范围从 0.001 arcsec / year 0.001 arcsec / year 0.001arcsec//year0.001 \mathrm{arcsec} / \mathrm{year} 到 0.02 arcsec / year 0.02 arcsec / year 0.02arcsec//year0.02 \mathrm{arcsec} / \mathrm{year} 。并与 Applegate 等人（1986）和 Carpino 等人（1986）的长期数值积分结果以及 Bretagnon（1974，1984）的解析理论进行了比较。外太阳系的解似乎比内太阳系的解更稳定，因为内太阳系的长周期频率中的大量准共振阻碍了解的良好收敛。

Key words: Celestial mechanics - planetary theories - solar system - secular perturbations
关键词：天体力学 - 行星理论 - 太阳系 - 长期扰动

1. Introduction  1. 简介

This paper is the third one on the construction of a solution for the secular evolution of the solar system which was undertaken at Bureau des Longitudes since 1984, following the earlier works of Bretagnon (1974) and Duriez (1979, 1982). Our first paper (Laskar, 1985) describes the methods used for the construction of the secular system at the second order with respect to the masses and up to degree 5 in the variables eccentricity and inclination. In the second paper (Laskar, 1986), a numerical integration of the whole system over 10,000 years allowed us to derive the secular terms in the form of polynomial series of the times which are then directly comparable to the secular terms of classical Le Verrier type solutions like VSOP82 (Bretagnon, 1982). The use of these secular terms allowed us to extend the length of validity of the classical ephemeris solution over 10,000 years (Bretagnon et al., 1985; Bretagnon and Simon, 1986) and to derive new formulae for the precessional quantities up to t 10 t 10 t^(10)t^{10} (Laskar, 1986). These two papers will be extensively quoted in the present work, and we shall refer to them respectively as ( P 1 P 1 P1P 1 ) and ( P 2 P 2 P2P 2 ). For all the definitions of the variables, the reader is supposed to refer to these previous papers.
本文是经度局自 1984 年以来开展的第三篇关于太阳系长期演化解构建的研究论文，此前，Bretagnon（1974）和 Duriez（1979、1982）也参与了相关研究。我们的第一篇论文（Laskar，1985）描述了构建二阶长期系统的方法，该系统涉及质量以及偏心率和倾角等变量，最高阶为 5 阶。在第二篇论文（Laskar，1986）中，我们对整个系统进行了一万年的数值积分，从而导出了时间多项式级数形式的长期项，这些长期项可以直接与经典勒威耶型解（如 VSOP82）的长期项进行比较（Bretagnon，1982）。这些长期项的使用使我们能够将经典星历表解的有效期延长至 10,000 年以上（Bretagnon 等人，1985 年；Bretagnon 和 Simon，1986 年），并推导出最高至 t 10 t 10 t^(10)t^{10} 的岁差量的新公式（Laskar，1986 年）。这两篇论文将在本书中被广泛引用，我们分别将其称为( P 1 P 1 P1P 1 )和( P 2 P 2 P2P 2 )。所有变量的定义，请读者参考这些先前的论文。
In the present work, the same secular system is integrated over 30 million years ( -10 to +20 millions). This integration is made possible by the use of a supercomputer CRAY-1. The results are then Fourier analysed in order to derive a quasi-periodic solution which is then comparable to previous works on analytical theories (Bretagnon, 1974 and Duriez, 1979) and to the recent numerical integrations of the outer planets made over such an extended period with a dedicated computer (Applegate et al., 1986) or a CRAY computer (Carpino et al., 1986).
在本研究中，我们将同一长期系统在 3000 万年（-1000 万到+2000 万年）内进行积分。该积分借助超级计算机 CRAY-1 实现。然后对结果进行傅里叶分析，以得出准周期解，该解可与先前关于解析理论的研究（Bretagnon，1974 年和 Duriez，1979 年）以及近期使用专用计算机（Applegate 等人，1986 年）或 CRAY 计算机（Carpino 等人，1986 年）对如此长周期的外行星进行的数值积分进行比较。

2. Description of the secular system
2. 世俗体系描述

The organization of the secular system showed itself to be critical as far as its computation and numerical integration are concerned. From the beginning, we started with a different attitude than in the usual analytical works dealing with Poisson quasi-periodical series: the computations are made in a very extensive way, without neglecting any of the terms which we could possibly encounter (of course, up to a given order with respect to the masses and to a given degree in the variables eccentricity and inclination) ( P 1 P 1 P1P 1 ). This led to a huge number of terms (153824 in the secular system) but these terms can be handled in a much easier way because the structure of the series is then made very simple and remains always the same.
就其计算和数值积分而言，世俗系统的组织至关重要。从一开始，我们就以一种不同于处理泊松准周期序列的常规分析工作的态度开始：计算以非常广泛的方式进行，不会忽略我们可能遇到的任何项（当然，在质量方面达到给定的阶数，在变量偏心率和倾角方面达到给定的阶数）（ P 1 P 1 P1P 1 ）。这导致了大量的项（世俗系统中有153824个），但这些项可以以更容易的方式处理，因为序列的结构变得非常简单并且始终保持不变。
The secular system is a polynomial system with imaginary coefficients and can be expressed in the form
世俗系统是一个具有虚系数的多项式系统，可以表示为以下形式
α ˙ = − 1 ( Φ 1 α + Φ 3 ( α , α ¯ ) + Φ 5 ( α , α ¯ ) ) α ˙ = − 1 Φ 1 α + Φ 3 ( α , α ¯ ) + Φ 5 ( α , α ¯ ) alpha^(˙)=sqrt(-1)(Phi_(1)alpha+Phi_(3)(alpha,( bar(alpha)))+Phi_(5)(alpha,( bar(alpha))))\dot{\alpha}=\sqrt{-1}\left(\Phi_{1} \alpha+\Phi_{3}(\alpha, \bar{\alpha})+\Phi_{5}(\alpha, \bar{\alpha})\right)
where α = ( z ¯ 1 , z ¯ 2 , … , z ¯ 8 , ζ ¯ 1 , ζ ¯ 2 , … , ζ ¯ 8 ) α = z ¯ 1 , z ¯ 2 , … , z ¯ 8 , ζ ¯ 1 , ζ ¯ 2 , … , ζ ¯ 8 alpha=( bar(z)_(1), bar(z)_(2),dots, bar(z)_(8), bar(zeta)_(1), bar(zeta)_(2),dots, bar(zeta)_(8))\alpha=\left(\bar{z}_{1}, \bar{z}_{2}, \ldots, \bar{z}_{8}, \bar{\zeta}_{1}, \bar{\zeta}_{2}, \ldots, \bar{\zeta}_{8}\right)  其中 α = ( z ¯ 1 , z ¯ 2 , … , z ¯ 8 , ζ ¯ 1 , ζ ¯ 2 , … , ζ ¯ 8 ) α = z ¯ 1 , z ¯ 2 , … , z ¯ 8 , ζ ¯ 1 , ζ ¯ 2 , … , ζ ¯ 8 alpha=( bar(z)_(1), bar(z)_(2),dots, bar(z)_(8), bar(zeta)_(1), bar(zeta)_(2),dots, bar(zeta)_(8))\alpha=\left(\bar{z}_{1}, \bar{z}_{2}, \ldots, \bar{z}_{8}, \bar{\zeta}_{1}, \bar{\zeta}_{2}, \ldots, \bar{\zeta}_{8}\right)
and where Φ 1 Φ 1 Phi_(1)\Phi_{1} is a real matrix with constant coefficients, Φ 3 Φ 3 Phi_(3)\Phi_{3} gathers all the terms of degree 3 and Φ 5 Φ 5 Phi_(5)\Phi_{5} the terms of degree 5 .
其中 Φ 1 Φ 1 Phi_(1)\Phi_{1} 是具有常数系数的实数矩阵， Φ 3 Φ 3 Phi_(3)\Phi_{3} 收集了所有 3 次项， Φ 5 Φ 5 Phi_(5)\Phi_{5} 收集了所有 5 次项。
All the terms are generated according to ( P 1 P 1 P1P 1 ) and obey the following rules:
所有词条均依据（ P 1 P 1 P1P 1 ）生成，并遵循以下规则：
  • The series in z ¯ z ¯ bar(z)\bar{z} are even in ζ , ζ ¯ ζ , ζ ¯ zeta, bar(zeta)\zeta, \bar{\zeta} while the series in ζ ζ zeta\zeta are odd in ζ , ζ ¯ ζ , ζ ¯ zeta, bar(zeta)\zeta, \bar{\zeta}.
    z ¯ z ¯ bar(z)\bar{z} 中的系列在 ζ , ζ ¯ ζ , ζ ¯ zeta, bar(zeta)\zeta, \bar{\zeta} 中为偶数，而 ζ ζ zeta\zeta 中的系列在 ζ , ζ ¯ ζ , ζ ¯ zeta, bar(zeta)\zeta, \bar{\zeta} 中为奇数。
  • For any given term T T TT, the characteristic of the monomial c M ( T ) c M ( T ) c_(M)(T)c_{M}(T) is equal to + 1 ( c M ( T ) = δ ( α ¯ ) − δ ( α ) + 1 c M ( T ) = δ ( α ¯ ) − δ ( α ) +1(c_(M)(T)=delta(( bar(alpha)))-delta(alpha):}+1\left(c_{M}(T)=\delta(\bar{\alpha})-\delta(\alpha)\right. where δ ( α ) δ ( α ) delta(alpha)\delta(\alpha) is the degree in α α alpha\alpha ). This last rule is the result of d’Alembert’s relations in the expansions of the disturbing function (See the Appendix).
    对于任何给定的项 T T TT ，单项式 c M ( T ) c M ( T ) c_(M)(T)c_{M}(T) 的特征等于 + 1 ( c M ( T ) = δ ( α ¯ ) − δ ( α ) + 1 c M ( T ) = δ ( α ¯ ) − δ ( α ) +1(c_(M)(T)=delta(( bar(alpha)))-delta(alpha):}+1\left(c_{M}(T)=\delta(\bar{\alpha})-\delta(\alpha)\right. ，其中 δ ( α ) δ ( α ) delta(alpha)\delta(\alpha) 是 α α alpha\alpha 的次数。最后这条规则是达朗贝尔在扰动函数展开式中的关系的结果（参见附录）。
These two simple rules allow to generate all the possible monomials that will actually appear in the computations. Up to the order 2 with respect to the masses, a monomial can only involve the elements of 1,2 or 3 planets. Up to the degree 5 , the number of monomials even or odd in inclination (or in ζ ζ zeta\zeta ) is given
这两个简单的规则可以生成所有可能在计算中出现的单项式。在质量阶数不超过2的单项式中，只能包含1、2或3个行星的元素。在阶数不超过5的单项式中，倾角（或 ζ ζ zeta\zeta ）为偶数或奇数的单项式的数量如下：
Table 1. Number of monomials involving 1, 2, or 3 planets in the secular system
表 1. 长期系统中涉及 1、2 或 3 个行星的单项式数量

涉及行星的数量
No. of planets
involved
No. of planets involved| No. of planets | | :--- | | involved |

每个方程的组合数
No. of combinations
in each equations
No. of combinations in each equations| No. of combinations | | :--- | | in each equations |
 Monomials  单项式 Even in inclination  即使倾斜
"No. of planets involved" "No. of combinations in each equations" Monomials Even in inclination| No. of planets <br> involved | No. of combinations <br> in each equations | Monomials | Even in inclination | | :--- | :--- | :--- | :---: |
Odd in inclination  倾向奇怪
in Table 1. We can construct a standard list of monomials for each case of Table 1 . This represents only 2 ( 10 + 102 + 318 ) = 2 ( 10 + 102 + 318 ) = 2(10+102+318)=2(10+102+318)= 860 typical monomials. All the other monomials are obtained by changing the name of the planets involved in the different list of monomials. We have 8 planets, 28 couples of planets and 56 groups of three planets. Each of these combinations will generate a different list of monomials, using the corresponding typical list. For any of the 16 equations in (1), the left member contains 8 × 10 + 28 × 102 + 21 × 318 = 9614 8 × 10 + 28 × 102 + 21 × 318 = 9614 8xx10+28 xx102+21 xx318=96148 \times 10+28 \times 102+21 \times 318=9614 monomials (in the monomial involving 3 planets, one of the planets is the planet corresponding to the given equation) which are stored in three different arrays T 1 ( 10 × 8 ) , T 2 ( 102 × 28 ) , T 3 ( 318 × 21 ) T 1 ( 10 × 8 ) , T 2 ( 102 × 28 ) , T 3 ( 318 × 21 ) T1(10 xx8),T2(102 xx28),T3(318 xx21)T 1(10 \times 8), T 2(102 \times 28), T 3(318 \times 21). In each array, the first index is the monomial index and the second one the combination of planets index. The total system (1) contains then 16 × 9614 = 153824 16 × 9614 = 153824 16 xx9614=15382416 \times 9614=153824 monomials which are then all generated by the 860 typical monomials and the different combinations of planets. The whole system is stored in three real arrays T T 1 ( 10 , 8 , 16 ) , T T 2 ( 102 , 28 , 16 ) T T 1 ( 10 , 8 , 16 ) , T T 2 ( 102 , 28 , 16 ) TT1(10,8,16),TT2(102,28,16)T T 1(10,8,16), T T 2(102,28,16) and T T 3 ( 318 , 21 , 16 ) T T 3 ( 318 , 21 , 16 ) TT3(318,21,16)T T 3(318,21,16). Apart from the typical monomials, there is no need to keep any information about the exponents in a given term: everything is contained in the location of the coefficient in the arrays.
在表 1 中。我们可以为表 1 中的每种情况构建一个标准单项式列表。这仅代表 2 ( 10 + 102 + 318 ) = 2 ( 10 + 102 + 318 ) = 2(10+102+318)=2(10+102+318)= 860 个典型单项式。所有其他单项式都通过更改不同单项式列表中涉及的行星的名称获得。我们有 8 颗行星、28 对行星和 56 组三颗行星。每种组合都会使用相应的典型列表生成不同的单项式列表。对于 (1) 中的 16 个方程中的任何一个，左边的成员包含 8 × 10 + 28 × 102 + 21 × 318 = 9614 8 × 10 + 28 × 102 + 21 × 318 = 9614 8xx10+28 xx102+21 xx318=96148 \times 10+28 \times 102+21 \times 318=9614 个单项式（在涉及 3 颗行星的单项式中，其中一个行星是与给定方程相对应的行星），它们存储在三个不同的数组 T 1 ( 10 × 8 ) , T 2 ( 102 × 28 ) , T 3 ( 318 × 21 ) T 1 ( 10 × 8 ) , T 2 ( 102 × 28 ) , T 3 ( 318 × 21 ) T1(10 xx8),T2(102 xx28),T3(318 xx21)T 1(10 \times 8), T 2(102 \times 28), T 3(318 \times 21) 中。在每个数组中，第一个索引是单项式索引，第二个索引是行星组合索引。整个系统 (1) 包含 16 × 9614 = 153824 16 × 9614 = 153824 16 xx9614=15382416 \times 9614=153824 个单项式，它们全部由 860 个典型单项式和不同的行星组合生成。整个系统存储在三个实数数组 T T 1 ( 10 , 8 , 16 ) , T T 2 ( 102 , 28 , 16 ) T T 1 ( 10 , 8 , 16 ) , T T 2 ( 102 , 28 , 16 ) TT1(10,8,16),TT2(102,28,16)T T 1(10,8,16), T T 2(102,28,16) 和 T T 3 ( 318 , 21 , 16 ) T T 3 ( 318 , 21 , 16 ) TT3(318,21,16)T T 3(318,21,16) 中。除了典型的单项式外，无需保存任何有关给定项的指数信息：所有内容都包含在数组中系数的位置中。

3. Numerical integration. Vectorisation
3. 数值积分。矢量化

The secular system (1) is close to its linear part
长期系统（1）接近其线性部分
α ˙ = − 1 Φ 1 α α ˙ = − 1 Φ 1 α alpha^(˙)=sqrt(-1)Phi_(1)alpha\dot{\alpha}=\sqrt{-1} \Phi_{1} \alpha
which is called the Laplace-Lagrange system. This system is easily integrated and gives some idea of the aspect of the solution. The solutions of (2) are linear combinations of periodic terms of frequencies ranging from about 50,000 years to several million years. The secular system (1) is very smooth, and we can use a very large stepsize for its numerical integration of about 500 years instead of about 0.5 days for the direct numerical integration of Newton’s law (and even less if we also include the Moon). The suppression of the short period terms allows to use a very large stepsize in the integration, but we must not forget that the new system (1) is only an approximation of second order of the secular system. The other drawback is that our secular system (1) is very large (153824 terms).
这被称为拉普拉斯-拉格朗日系统。该系统易于积分，并能在一定程度上反映解的性质。(2) 的解是周期项的线性组合，周期项的频率范围从约 5 万年到数百万年。长期系统 (1) 非常平滑，我们可以对其数值积分使用非常大的步长，约为 500 年，而不是像牛顿定律的直接数值积分那样使用约 0.5 天（如果包含月球，步长会更小）。抑制短周期项使得积分可以使用非常大的步长，但我们不能忘记，新系统 (1) 只是长期系统的二阶近似。另一个缺点是我们的长期系统 (1) 非常庞大（153824 个项）。
With such a huge second hand member, the best results will be obtained with high order numerical integration routines. We used a well experienced predictor-corrector Adams routine D02CBF from the NAG library with varying order and stepsize and with an output every 500 yr . The use of this routine on a CRAY-1 computer imposed to make some experiments for the precision of the integration over a long time. A typical try is made by integrating a well known circular function of frequency 25 s / yr 25 s / yr 25s//yr25 \mathrm{~s} / \mathrm{yr} and 50 s / yr 50 s / yr 50s//yr50 \mathrm{~s} / \mathrm{yr} (period of about 52,000 and 26 , 000 yr 26 , 000 yr 26,000yr26,000 \mathrm{yr} ) over
由于二手数据量如此巨大，使用高阶数值积分程序将获得最佳结果。我们使用了 NAG 库中一个经验丰富的预估-校正 Adams 程序 D02CBF，该程序具有可变的阶数和步长，每 500 年输出一次。在 CRAY-1 计算机上使用该程序进行了一些实验，以验证其长期积分的精度。一个典型的尝试是，对一个频率为 25 s / yr 25 s / yr 25s//yr25 \mathrm{~s} / \mathrm{yr} 和 50 s / yr 50 s / yr 50s//yr50 \mathrm{~s} / \mathrm{yr} （周期约为 52,000 和 26 , 000 yr 26 , 000 yr 26,000yr26,000 \mathrm{yr} ）的著名圆函数进行积分。
10 million yr. The evaluation of the error is then simply made by difference with the known value of the function. The relative accumulated error goes from 5 × 10 − 8 5 × 10 − 8 5xx10^(-8)5 \times 10^{-8} to 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} after 10 million years which means from 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} to 8 × 10 − 7 8 × 10 − 7 8xx10^(-7)8 \times 10^{-7} after 20 million years. This precision is good enough for our computation and will ensure that there are no problems coming from the numerical integration itself. We must mention that on a CRAY computer, computation are made with only 14 digits which is rather low compared to the 16 usual digits on other computers.
1000 万年。然后，误差的评估只需通过与已知函数值的差值即可。相对累积误差在 1000 万年后从 5 × 10 − 8 5 × 10 − 8 5xx10^(-8)5 \times 10^{-8} 变为 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} ，这意味着在 2000 万年后从 2 × 10 − 7 2 × 10 − 7 2xx10^(-7)2 \times 10^{-7} 变为 8 × 10 − 7 8 × 10 − 7 8xx10^(-7)8 \times 10^{-7} 。这个精度对于我们的计算来说已经足够好，并且能够确保数值积分本身不会出现问题。必须指出的是，在 CRAY 计算机上，计算仅使用 14 位数字，与其他计算机上通常的 16 位数字相比，这个数字相当低。
The main part of the computing time in the numerical integration is used in the evaluation of the 153824 terms of the second member of (1). Typically, the values of z i , ζ i z i , ζ i z_(i),zeta_(i)z_{i}, \zeta_{i} are given and a standard term of degree 5 is of the form
数值积分中计算时间的主要部分用于计算式（1）中第二个元素的153824个项。通常， z i , ζ i z i , ζ i z_(i),zeta_(i)z_{i}, \zeta_{i} 的值是给定的，并且5次标准项的形式为
C × X i 1 X i 2 X i 3 X i 4 X i 5 C × X i 1 X i 2 X i 3 X i 4 X i 5 C xxX_(i_(1))X_(i_(2))X_(i_(3))X_(i_(4))X_(i_(5))C \times X_{i_{1}} X_{i_{2}} X_{i_{3}} X_{i_{4}} X_{i_{5}}
where C C CC is a real coefficient, and X i k X i k X_(i_(k))X_{i_{k}} any of the complex variables z i , ζ i , z ¯ i , ζ ¯ i z i , ζ i , z ¯ i , ζ ¯ i z_(i),zeta_(i), bar(z)_(i), bar(zeta)_(i)z_{i}, \zeta_{i}, \bar{z}_{i}, \bar{\zeta}_{i}.
其中 C C CC 是实系数， X i k X i k X_(i_(k))X_{i_{k}} 是任意复变量 z i , ζ i , z ¯ i , ζ ¯ i z i , ζ i , z ¯ i , ζ ¯ i z_(i),zeta_(i), bar(z)_(i), bar(zeta)_(i)z_{i}, \zeta_{i}, \bar{z}_{i}, \bar{\zeta}_{i} 。
In a CRAY-1, computations are vectorised if the most inner loop is long, the optimum being obtained by a multiple of 64 which is the length of the accumulators. The difficulty is that the vectorisation is not possible with indirect addressing (the elements to be multiplied need to be physically in the accumulators before the multiplication). The only way we found to deal with this problem was to choose first the monomial type (corresponding to the first index in TT1, TT2 and TT3). The inner loop is then made over all the different combinations of planets involved in all the different equations. This mean an inner loop of 28 for the 2-planets monomials of TT2 and of 21 × 8 = 168 21 × 8 = 168 21 xx8=16821 \times 8=168 for the 3-planets monomials of TT3. The vectorisation obtained is then very efficient and reached the possible limit of the CRAY. The program was tested also on the NAS9080 computer of the CIRCE. The NAS9080 is a conventional computer with about the same speed as the scalar computation on the CRAY. We gained a factor of 3 with the CRAY due to the very rapid direct memory of the CRAY, and the vectorisation of the program added a factor of 5 , which means a total factor of 15 .
在 CRAY-1 中，如果最内层循环较长，则计算将进行矢量化，最佳值为累加器长度 64 的倍数。难点在于，矢量化无法通过间接寻址实现（待乘元素在乘法运算前必须物理存储在累加器中）。我们发现解决这个问题的唯一方法是首先选择单项式类型（对应于 TT1、TT2 和 TT3 中的第一个索引）。然后，对所有不同方程式中涉及的所有行星组合进行内循环。这意味着 TT2 的 2 行星单项式的内循环长度为 28，TT3 的 3 行星单项式的内循环长度为 21 × 8 = 168 21 × 8 = 168 21 xx8=16821 \times 8=168 。获得的矢量化非常高效，达到了 CRAY 的可能极限。该程序也在 CIRCE 的 NAS9080 计算机上进行了测试。 NAS9080 是一台传统计算机，其速度与 CRAY 上的标量计算速度大致相同。由于 CRAY 拥有非常快速的直接内存，因此 CRAY 的计算速度提高了 3 倍，而程序的矢量化又提高了 5 倍，总速度达到了 15 倍。
I should mention again that this vectorisation was only possible because of the very simple organization of the series. As no term was neglected, all the different series were alike and once a monomial was selected, we could loop on all the possible combinations of the planets involved. We obtain the somewhat surprising fact that keeping more terms leads to faster computations!
我应该再次强调，这种矢量化之所以能够实现，是因为该序列的组织非常简单。由于没有忽略任何项，所有不同的序列都相同，一旦选定一个单项式，我们就可以循环遍历所有可能涉及的行星组合。我们得到了一个令人惊讶的事实：保留更多项可以加快计算速度！
The total evaluation of the second member of (1) takes 0.06 second. The numerical integration of the secular system over 10 million years took about 30 minutes which means about 20,000 evaluations of the second member (about one every 500 yr ). The integration was made first 10 million years backwards, and then 20 million years forwards for a total computer time of 1.5 hour. Of course, there is nothing which prevents us from integrating
(1) 式中第二个元素的总计算耗时 0.06 秒。对 1000 万年的长期系统进行数值积分大约需要 30 分钟，这意味着对第二个元素进行了大约 20,000 次计算（大约每 500 年一次）。积分首先向后推 1000 万年，然后向前推 2000 万年，总共耗时 1.5 小时。当然，没有什么可以阻止我们进行积分
the system over a much longer time ( 1.5 hours is very reasonable), except that our system is an approximation of the real system, and we need first to analyse our results before thinking of extending the integration.
该系统需要更长的时间（ 1.5 小时非常合理），但我们的系统是真实系统的近似值，我们需要先分析我们的结果，然后才能考虑扩展集成。
The output consists of all the different variables z i , ζ i z i , ζ i z_(i),zeta_(i)z_{i}, \zeta_{i} every 500 years from − 10 , 000 , 000 − 10 , 000 , 000 -10,000,000-10,000,000 years to 20 , 000 , 000 20 , 000 , 000 20,000,00020,000,000 years.
输出包括从 − 10 , 000 , 000 − 10 , 000 , 000 -10,000,000-10,000,000 年到 20 , 000 , 000 20 , 000 , 000 20,000,00020,000,000 年每 500 年出现一次的所有不同变量 z i , ζ i z i , ζ i z_(i),zeta_(i)z_{i}, \zeta_{i} 。

4. Fourier analysis of the solution
4. 解的傅里叶分析

If we integrate the secular system (1) in an analytical manner as in Duriez (1979) or Laskar (1984), we obtain a quasi-periodic solution where the arguments are linear integer combinations of the basic frequencies of the system. As was outlined in Laskar (1984), in the Birkhoff normalization of the system (1), serious problems of convergence appear due to the presence of numerous small divisors for the 8 planet systems. This leads us to prefer a numerical integration of the system which avoids all these problems, and provides an accurate solution, although not very utilizable.
如果我们像 Duriez (1979) 或 Laskar (1984) 那样，以解析的方式对长期系统 (1) 进行积分，我们会得到一个准周期解，其中自变量是系统基频的线性整数组合。正如 Laskar (1984) 所述，在系统 (1) 的 Birkhoff 正则化中，由于 8 个行星系统存在大量小因子，因此会出现严重的收敛问题。因此，我们更倾向于对系统进行数值积分，这样可以避免所有这些问题，并提供精确的解，尽管实用性不强。
With the Fourier analysis of the solution, we try to obtain a more handy representation of the solution, similar to the result of an analytical solution.
通过对解进行傅里叶分析，我们试图获得更方便的解表示，类似于解析解的结果。
Different methods have been already used for a similar purpose by Applegate et al., (1986) and Carpino et al. (1986). The main problem is to get a very good resolution in frequency, which is not given by a usual Fourier transform. Our concern was also not to anticipate at all on the form of the solution, even if we could have an idea of the solution by the knowledge given by analytical integration. This excludes recognition of combination of frequencies and least-squares adjustments during the frequency analysis like in Carpino et al. (1986). All the solutions for the different variables and planets are analysed separately, and only the analysis of the results will show the existence of common frequencies in the solutions. The reason to do this is that in fact if we can forecast the possible results when the normalization process converges correctly, very little is known in case of bad convergence.
Applegate 等人（1986）和 Carpino 等人（1986）已经使用不同的方法实现类似的目的。主要问题是获得非常好的频率分辨率，而这通常的傅里叶变换无法实现。我们关注的重点是根本不预测解的形式，即使我们可以通过解析积分获得的知识对解有所了解。这排除了在频率分析过程中识别频率组合和最小二乘调整的可能性，就像 Carpino 等人（1986）的方法一样。所有不同变量和行星的解都被分别分析，只有对结果的分析才能揭示解中存在共同的频率。这样做的原因是，事实上，如果我们能够预测正则化过程正确收敛时的可能结果，那么在收敛性不佳的情况下，我们所知甚少。
The only assumption is that we search a solution on a quasiperiodic form, i.e. for each planet and variable, a solution in the form
唯一的假设是，我们寻找准周期形式的解，即对于每个行星和变量，形式为
α = ∑ i = 1 N A i exp ⁡ − 1 ( v i t + ϕ i ) α = ∑ i = 1 N   A i exp ⁡ − 1 v i t + ϕ i alpha=sum_(i=1)^(N)A_(i)exp sqrt(-1)(v_(i)t+phi_(i))\alpha=\sum_{i=1}^{N} A_{i} \exp \sqrt{-1}\left(v_{i} t+\phi_{i}\right)
where A i , v i , ϕ i A i , v i , ϕ i A_(i),v_(i),phi_(i)A_{i}, v_{i}, \phi_{i} are real numbers, and we will try to represent the solution with the minimum possible number of arguments. The principle of the process is simple. We make first a FFT (Fast Fourier Transform) of the solution. Then we search for the term of maximum amplitude of the FFT. Around this term, we search for the value of v v vv which gives the maximum of the power spectrum
其中 A i , v i , ϕ i A i , v i , ϕ i A_(i),v_(i),phi_(i)A_{i}, v_{i}, \phi_{i} 为实数，我们将尝试用尽可能少的参数来表示解。这个过程的原理很简单。我们首先对解进行快速傅里叶变换 (FFT)。然后，我们寻找 FFT 的最大振幅项。围绕该项，我们寻找 v v vv 的值，使得功率谱达到最大值。
ϕ ( v ) = ∫ α ( t ) exp ⁡ − 1 ( − v t ) d t ϕ ( v ) = ∫ α ( t ) exp ⁡ − 1 ( − v t ) d t phi(v)=int alpha(t)exp sqrt(-1)(-vt)dt\phi(v)=\int \alpha(t) \exp \sqrt{-1}(-v t) d t
this fine tuning of the frequency is made with a precision of 0.00002 arcsec / year 0.00002 arcsec / year 0.00002arcsec//year0.00002 \mathrm{arcsec} / \mathrm{year}. Once the frequency is found, the power spectrum gives the amplitude of the first term which is then subtracted from the solution. The process is then started again to search for the following term. As the two terms exp ⁡ − 1 ( v 1 t ) exp ⁡ − 1 v 1 t exp sqrt(-1)(v_(1)t)\exp \sqrt{-1}\left(v_{1} t\right) and exp ⁡ − 1 ( v 2 t ) exp ⁡ − 1 v 2 t exp sqrt(-1)(v_(2)t)\exp \sqrt{-1}\left(v_{2} t\right) are not orthogonal, on the second step, the amplitude of each term will not be given by the corresponding power spectrum like in regular Fourier analysis. We have then to orthog-
这种频率的微调精度为 0.00002 arcsec / year 0.00002 arcsec / year 0.00002arcsec//year0.00002 \mathrm{arcsec} / \mathrm{year} 。一旦找到频率，功率谱就会给出第一项的振幅，然后将其从解中减去。然后再次开始该过程以寻找下一项。由于两个项 exp ⁡ − 1 ( v 1 t ) exp ⁡ − 1 v 1 t exp sqrt(-1)(v_(1)t)\exp \sqrt{-1}\left(v_{1} t\right) 和 exp ⁡ − 1 ( v 2 t ) exp ⁡ − 1 v 2 t exp sqrt(-1)(v_(2)t)\exp \sqrt{-1}\left(v_{2} t\right) 不正交，因此在第二步中，每个项的振幅将不会像常规傅里叶分析那样由相应的功率谱给出。因此，我们必须进行正交-
onalize our basis of function by the Gauss algorithm, and then compute the amplitude of each term by computing the inner products
通过高斯算法对函数基进行函数化，然后通过计算内积来计算每个项的振幅
∫ α ( t ) f i ( t ) ― d t ∫ α ( t ) f i ( t ) ¯ d t int alpha(t) bar(f_(i)(t))dt\int \alpha(t) \overline{f_{i}(t)} d t
where the f i f i f_(i)f_{i} are the functions of the orthogonal basis made with the e i = exp ⁡ − 1 ( v i t ) e i = exp ⁡ − 1 v i t e_(i)=exp sqrt(-1)(v_(i)t)e_{i}=\exp \sqrt{-1}\left(v_{i} t\right) (the f i ( t ) f i ( t ) f_(i)(t)f_{i}(t) are linear combinations of the e i ( t ) ) e i ( t ) {:e_(i)(t))\left.e_{i}(t)\right). We have a decomposition of the solution on the orthogonal basis ( f i ( t ) ) i = 1 , N f i ( t ) i = 1 , N (f_(i)(t))_(i=1,N)\left(f_{i}(t)\right)_{i=1, N} which gives then the decomposition on the nonorthogonal basis ( exp ⁡ − 1 ( v i t ) ) i = 1 , N exp ⁡ − 1 v i t i = 1 , N (exp sqrt(-1)(v_(i)t))_(i=1,N)\left(\exp \sqrt{-1}\left(v_{i} t\right)\right)_{i=1, N}.
其中 f i f i f_(i)f_{i} 是用 e i = exp ⁡ − 1 ( v i t ) e i = exp ⁡ − 1 v i t e_(i)=exp sqrt(-1)(v_(i)t)e_{i}=\exp \sqrt{-1}\left(v_{i} t\right) 构成的正交基函数（ f i ( t ) f i ( t ) f_(i)(t)f_{i}(t) 是 e i ( t ) ) e i ( t ) {:e_(i)(t))\left.e_{i}(t)\right) 的线性组合）。我们对正交基 ( f i ( t ) ) i = 1 , N f i ( t ) i = 1 , N (f_(i)(t))_(i=1,N)\left(f_{i}(t)\right)_{i=1, N} 上的解进行分解，然后由此得出非正交基 ( exp ⁡ − 1 ( v i t ) ) i = 1 , N exp ⁡ − 1 v i t i = 1 , N (exp sqrt(-1)(v_(i)t))_(i=1,N)\left(\exp \sqrt{-1}\left(v_{i} t\right)\right)_{i=1, N} 上的分解。
The advantage of this method is that there is no suppositions on the analytical form of the solution and all the different solutions for the different planets are treated in a completely independent way. The information on the analytical form of the solution will be given by the existence of similar features in all the solutions, and by the recognition of all the basic frequencies from one solution to the other. No least square adjustments are made which could group together frequencies that are close and then prevent further analysis. Of course, the method has its drawbacks, one of which is that the time used for the orthogonalization of the basis of functions grows quadratically with the number of terms. For more than 40 terms, most of the time was used for this orthogonalization.
该方法的优点在于无需对解的解析形式进行任何假设，并且不同行星的所有不同解都以完全独立的方式处理。解的解析形式的信息将通过所有解中存在的相似特征以及从一个解到另一个解的所有基频的识别来提供。无需进行最小二乘调整，因为最小二乘调整可能会将接近的频率归为一类，从而妨碍进一步分析。当然，该方法也有其缺点，其中之一就是用于函数基正交化的时间会随着项数的平方增长。对于超过 40 个项的情况，大部分时间都用于这种正交化。
The solutions for the variables z i , ζ i z i , ζ i z_(i),zeta_(i)z_{i}, \zeta_{i} of the different planets are given in Tables 6-13. The presentation is different for the inner planets and for the outer planets. As it will be detailed further, the solutions for the inner planets converge poorly and the terms are not easily identified. They are given in Tables 6-9 with 80 terms ordered by decreasing amplitude. The solutions for the outer planets (Tables 10 − 13 10 − 13 10-1310-13 ) contain 51 terms ordered by decreasing amplitude. The amplitude of the terms decreases rapidly and several terms are easily identified. In this case the combination of fundamental frequencies corresponding to the term is given in column 1 . The residuals between the numerical integration and the series of the solutions are given in Table 5. The fundamental frequencies are easily identified and we can check their internal precision by comparison of the values obtained in the solutions for all the 8 planets. The values of these frequencies are gathered in Table 3.
表 6-13 给出了不同行星的变量 z i , ζ i z i , ζ i z_(i),zeta_(i)z_{i}, \zeta_{i} 的解。内行星和外行星的呈现方式不同。正如将进一步详细说明的那样，内行星的解收敛性较差，且项不易识别。表 6-9 给出了它们，其中 80 个项按振幅递减顺序排列。外行星的解（表 10 − 13 10 − 13 10-1310-13 ）包含 51 个项，按振幅递减顺序排列。项的振幅迅速下降，有几个项很容易识别。在这种情况下，与该项相对应的基频组合在第 1 列给出。表 5 给出了数值积分与解系列之间的残差。基频很容易识别，我们可以通过比较所有 8 颗行星的解中获得的值来检查它们的内部精度。这些频率的值汇总在表 3 中。

5. Accuracy of the solution NGT
5. NGT 解决方案的准确性

In every work on analytical theory with perturbation methods one of the most difficult aspects is the analysis of the accuracy of the final solution. Truncations are often made, usually without any real evaluation of the remainder. One basic reason for this is that the problem is too difficult and no definite results have been obtained on the convergence of the series in use. In our solution NGT, we have computed all the terms up to the second order with respect to the masses and degree 5 in eccentricity and inclinations. With these limitations, we obtain very accurate and controlled results ( P 1 P 1 P1P 1 ), but some estimation needs to be done on higher orders and higher degrees contributions. There exists no way to make comparisons on the whole solution, but the classical Le Verrier type semi-analytical theories give a very efficient comparison at the origin of time of the solution. The solution VSOP82 of Bretagnon (1982) is the most developed solution of this kind so far. It is a semi-analytical solution where the secular terms are expressed in polynomial form, and thus valid only over a few thousand years. The contribution of all
在每一项运用扰动方法的解析理论研究中，最困难的方面之一就是分析最终解的精度。通常会进行截断，但通常不会对余数进行任何实际计算。其中一个基本原因是问题过于困难，而且尚未获得所用级数收敛性的明确结果。在我们的解 NGT 中，我们计算了所有关于质量以及偏心率和倾角的五阶项，直至二阶。在这些限制条件下，我们获得了非常精确且可控的结果 ( P 1 P 1 P1P 1 )，但需要对更高阶和更高阶的贡献进行一些估算。目前无法对整个解进行比较，但经典的勒威耶型半解析理论可以在解的时间原点处提供非常有效的比较。Bretagnon (1982) 的解 VSOP82 是迄今为止此类解中最完善的。这是一个半解析解，其中长期项以多项式形式表示，因此仅在几千年内有效。所有