作者君在作品相关中其实已经解释过这个问题。
不过仍然有人质疑。
那么作者君在此列出相关参考文献中的一篇开源论文。
以下是文章内容：
Long-term iions and stability of pary orbits in our Sor system
Abstract
We present the results of very long-term numerical iions of pary orbital motions over 109 -yr time-spans including all nine ps. A quispe of our numerical data shows that the pary motion， at least in our simple dynamical model， seems to be quite stable evehis very long time-span. A closer look at the lowest-frequency osciltions using a low-pass filter shows us the potentially diffusive character of terrestrial pary motion， especially that of Mercury. The behaviour of the etricity of Mercury in our iions is qualitatively simir to the results from Jacques Laskar's secur perturbation theory (e.g. emax～ 0.35 over ～± 4 Gyr). However， there are no apparent secur increases of etricity or ination in any orbital elements of the ps， which may be revealed by still loerm numerical iions. We have also performed a couple of trial iions including motions of the outer five ps over the duration of ± 5 × 1010 yr. The result indicates that the three major resonances in the une–Pluto system have been maintained over the 1011-yr time-span.
1 Introdu
1.1Definition of the problem
The question of the stability of our Sor system has beeed over several hundred years， sihe era of on. The problem has attracted many famous mathematis over the years and has pyed a tral role in the development of non-linear dynamid chaos theory. However， we do not yet have a definite ao the question of whether our Sor system is stable or not. This is partly a result of the fact that the definition of the term ‘stability’ is vague when it is used iion to the problem of pary motion in the Sor system. Actually it is not easy to give a clear， rigorous and physically meaningful definition of the stability of our Sor system.
Among many definitions of stability， here t the Hill definition (Gdman 1993): actually this is not a definition of stability， but of instability. We define a system as being unstable when a close enter occurs somewhere in the system， starting from a certain initial figuration (Chambers， Wetherill & Boss 1996; Ito & Tanikawa 1999). A system is defined as experieng a close enter when two bodies approae another within an area of the rger Hill radius. Otherwise the system is defined as being stable. Henceforward we state that our pary system is dynamically stable if no close enter happens during the age of our Sor system， about ±5 Gyr. Ially， this definition may be repced by one in whi occurrence of any orbital crossiweeher of a pair of pakes pce. This is because we know from experiehat an orbital crossing is very likely to lead to a close enter iary and protopary systems (Yoshinaga， Kokubo & Makino 1999). Of course this statement ot be simply applied to systems with stable orbital resonances such as the une–Pluto system.
1.2Previous studies and aims of this research
In addition to the vagueness of the cept of stability， the ps in our Sor system show a character typical of dynamical chaos (Sussman & Wisdom 1988， 1992). The cause of this chaotic behaviour is now partly uood as being a result of resonance overpping (Murray & Holman 1999; Lecar， Franklin & Holman 2001). However， it would require iing over an ensemble of pary systems including all nine ps for a period c several 10 Gyr to thhly uand the long-term evolution of pary orbits， since chaotiamical systems are characterized by their strong dependen initial ditions.
From that point of view， many of the previous long-term numerical iions included only the outer five ps (Sussman & Wisdom 1988; Kinoshita & Nakai 1996). This is because the orbital periods of the outer ps are so much lohan those of the inner four phat it is much easier to follow the system fiven iion period. At present， the lo numerical iions published in journals are those of Dun & Lissauer (1998). Although their main target was the effect of post-main-sequenass loss oability of pary orbits， they performed many iions c up to ～1011 yr of the orbital motions of the four jovias. The initial orbital elements and masses of ps are the same as those of our Sor system in Dun & Lissauer's paper， but they decrease the mass of the Sun gradually in their numerical experiments. This is because they sider the effect of post-main-sequenass loss in the paper. sequently， they found that the crossing time-scale of pary orbits， which be a typical indicator of the instability time-scale， is quite sensitive to the rate of mass decrease of the Sun. When the mass of the Sun is close to its present value， the jovias remain stable over 1010 yr， or perhaps longer. Dun & Lissauer also performed four simir experiments on the orbital motion of seves (Venus to une)， which cover a span of ～109 yr. Their experiments on the seves are not yet prehensive， but it seems that the terrestrial ps also remain stable during the iion period， maintaining almur osciltions.
Oher hand， in his accurate semi-analytical secur perturbation theory (Laskar 1988)， Laskar finds that rge and irregur variations appear in the etricities and inations of the terrestrial ps， especially of Mercury and Mars on a time-scale of several 109 yr (Laskar 1996). The results of Laskar's secur perturbation theory should be firmed and iigated by fully numerical iions.
In this paper we present preliminary results of six long-term numerical iions on all nine pary orbits， c a span of several 109 yr， and of two other iions c a span of ± 5 × 1010 yr. The total epsed time for all iions is more than 5 yr， using several dedicated Pd workstations. One of the fual clusions of our long-term iions is that Sor system pary motioo be stable in terms of the Hill stability mentioned above， at least over a time-span of ± 4 Gyr. Actually， in our numerical iions the system was far more stable than what is defined by the Hill stability criterion: not only did no close enter happen during the iion period， but also all the pary orbital elements have been fined in a narrion both in time and frequenain， though pary motions are stochastic. Sihe purpose of this paper is to exhibit and overview the results of our long-term numerical iions， we show typical example figures as evidence of the very long-term stability of Sor system pary motion. For readers who have more specifid deeper is in our numerical results， repared a webpage (access )， where we show raw orbital elements， their low-pass filtered results， variation of Deunay elements and angur momentum deficit， as of our simple time–frequenalysis on all of our iions.
Iion 2 we briefly expin our dynamical model， numerical method and initial ditions used in our iions. Se 3 is devoted to a description of the quick results of the numerical iions. Very long-term stability of Sor system pary motion is apparent both iary positions and orbital elements. A rough estimation of numerical errors is also giveion 4 goes on to a discussion of the loerm variation of pary orbits using a low-pass filter and includes a discussion of angur momentum deficit. Iion 5， we present a set of numerical iions for the outer five phat spans ± 5 × 1010 yr. Iion 6 we also discuss the long-term stability of the pary motion and its possible cause.
2 Description of the numerical iions
（本部分涉及比较复杂的积分计算，作者君就不贴上来了，贴上来了起点也不一定能成功显示。）
2.3 Numerical method
We utilize a sed-order Wisdom–Holman symplectic map as our main iiohod (Wisdom & Holman 1991; Kinoshita， Yoshida & Nakai 1991) with a special start-up procedure to reduce the truncation error of angle variables，‘warm start’(Saha & Tremaine 1992， 1994).
The stepsize for the numerical iions is 8 d throughout all iions of the nine ps (N±1，2，3)， which is about 1/11 of the orbital period of the innermost p (Mercury). As for the determination of stepsize， we partly follow the previous numerical iion of all nine ps in Sussman & Wisdom (1988， 7.2 d) and Saha & Tremaine (1994， 225/32 d). We rouhe decimal part of the their stepsizes to 8 to make the stepsize a multiple of 2 in order to reduce the accumution of round-off error in the putation processes. Iion to this， Wisdom & Holman (1991) performed numerical iions of the outer five pary orbits using the symplectic map with a stepsize of 400 d， 1/10.83 of the orbital period of Jupiter. Their result seems to be accurate enough， which partly justifies our method of determining the stepsize. However， sihe etricity of Jupiter (～0.05) is much smaller than that of Mercury (～0.2)， we need some care when we pare these iions simply in terms of stepsizes.
Iegration of the outer five ps (F±)， we fixed the stepsize at 400 d.
t Gauss' f and g funs in the symplectic map together with the third-order Halley method (Danby 1992) as a solver for Kepler equations. The number of maximum iteratio in Halley's method is 15， but they never reached the maximum in any of our iions.
The interval of the data output is 200 000 d (～547 yr) for the calcutions of all nine ps (N±1，2，3)， and about 8000 000 d (～21 903 yr) for the iion of the outer five ps (F±).
Although no output filtering was done when the numerical iions were in process， lied a low-pass filter to the raw orbital data after we had pleted all the calcutions. See Se 4.1 for more detail.
2.4 Error estimation
2.4.1 Retive errors in total energy and angur momentum
Acc to one of the basic properties of symplectitegrators， which serve the physically servative quantities well (total orbital energy and angur momentum)， our long-term numerical iioo have been performed with very small errors. The averaged retive errors of total energy (～10?9) and of total angur momentum (～10?11) have remained nearly stant throughout the iion period (Fig. 1). The special startup procedure， warm start， would have reduced the averaged retive error in total energy by about one order of magnitude or more.
Retive numerical error of the total angur momentum δA/A0 and the total energy δE/E0 in our numerical iionsN± 1，2，3， where δE and δA are the absolute ge of the total energy and total angur momentum， respectively， andE0andA0are their initial values. The horizontal unit is Gyr.
hat different operating systems， different mathematical libraries， and different hardware architectures result in different numerical errors， through the variations in round-off error handling and numerical algorithms. In the upper panel of Fig. 1， we reize this situation in the seumerical error ial angur momentum， which should be rigorously preserved up to mae-ε precision.
2.4.2 Error iary longitudes
Sihe symplectic maps preserve total energy and total angur momentum of N-body dynamical systems ily well， the degree of their preservation may not be a good measure of the accuracy of numerical iions， especially as a measure of the positional error of ps， i.e. the error iary longitudes. To estimate the numerical error in the pary longitudes， we performed the following procedures. We pared the result of our main long-term iions with some test iions， which span much shorter periods but with much higher accuracy than the main iions. For this purpose， we performed a much more accurate iion with a stepsize of 0.125 d (1/64 of the main iions) spanning 3 × 105 yr， starting with the same initial ditions as in the N?1 iion. We sider that this test iion provides us with a ‘pseudo-true’ solution of pary orbital evolutio， we pare the test iion with the main iion， N?1. For the period of 3 × 105 yr， we see a differen mean anomalies of the Earth betweewo iions of ～0.52°(in the case of the N?1 iion). This difference be extrapoted to the value ～8700°， about 25 rotations of Earth after 5 Gyr， sihe error of longitudes increases linearly with time in the symplectic map. Simirly， the longitude error of Pluto be estimated as ～12°. This value for Pluto is much better than the result in Kinoshita & Nakai (1996) where the difference is estimated as ～60°.
3 Numerical results – I. G the raw data
In this se we briefly review the long-term stability of pary orbital motion through some snapshots of raw numerical data. The orbital motion of ps indicates long-term stability in all of our numerical iions: no orbital crossings nor close enters between any pair of pook pce.
3.1 General description of the stability of pary orbits
First， we briefly look at the general character of the long-term stability of pary orbits. Our i here focuses particurly on the inner four terrestrial ps for which the orbital time-scales are much shorter than those of the outer five ps. As we see clearly from the pnar orbital figurations shown in Figs 2 and 3， orbital positions of the terrestrial ps differ little between the initial and final part of eaumerical iion， which spans several Gyr. The solid lines denoting the present orbits of the ps lie almost within the swarm of dots even in the final part of iions (b) and (d). This indicates that throughout the eegration period the almur variations of pary orbital motion remain nearly the same as they are at present.
Vertical view of the four inner pary orbits (from the z -axis dire) at the initial and final parts of the iionsN±1. The axes units are au. The xy -pne is set to the invariant pne of Sor system total angur momentum.(a) The initial part ofN+1 ( t = 0 to 0.0547 × 10 9 yr).(b) The final part ofN+1 ( t = 4.9339 × 10 8 to 4.9886 × 10 9 yr).(c) The initial part of N?1 (t= 0 to ?0.0547 × 109 yr).(d) The final part ofN?1 ( t =?3.9180 × 10 9 to ?3.9727 × 10 9 yr). In each panel， a total of 23 684 points are plotted with an interval of about 2190 yr over 5.47 × 107 yr . Solid lines in each panel dehe present orbits of the four terrestrial ps (taken from DE245).