作者君在作品相關中其實已經解釋過這個問題。

　　不過仍然有人質疑——“你說得太含糊了”，“火星軌道的變化比你想象要大得多！”

　　那好吧，既然作者君的簡單解釋不夠有力，那咱們就看看嚴肅的東西，反正這本書寫到現在，嚷嚷著本書BUG一大堆，用初高中物理在書中挑刺的人也不少。

　　以下是文章內容:

　　Long-term and stability of orbits in our system

　　Abstract

　　We present the results of very long-term numerical of orbital over 109 -yr time-spans including all nine A quick of our numerical data shows that the ， at least in our simple dynamical model， seems to be quite stable even over this very long time-span. A closer look at the lowest-frequency using a low-pass filter shows us the potentially diffusive character of terrestrial ， especially that of Mercury. The of the of Mercury in our is qualitatively to the results from Jacques Laskar's theory (e.g. emax～ 0.35 over ～± 4 Gyr). However， there are no apparent increases of or in any orbital elements of the ， which may be revealed by still longer-term numerical We have also performed a couple of trial including of the outer five over the of ± 5 × 1010 yr. The result indicates that the three major resonances in the Neptune–Pluto system have been maintained over the 1011-yr time-span.

　　1

　　 of the problem

　　The of the stability of our system has been debated over several hundred years， since the era of Newton. The problem has attracted many famous mathematicians over the years and has a central role in the development of non-linear dynamics and chaos theory. However，

we do not yet have a definite answer to the of whether our system is stable or not. This is partly a result of the fact that the of the term ‘stability’ is vague when it is used in to the problem of in the system. Actually it is not easy to give a clear， rigorous and physically meaningful of the stability of our system.　　Among many of stability， here we adopt the Hill ( 1993): actually this is not a of stability， but of instability. We define a system as unstable when a close somewhere in the system， starting from a certain initial (Chambers， Wetherill & Boss 1996; Ito & Tanikawa 1999). A system is defined as experiencing a close when two bodies approach one another within an area of the larger Hill radius. Otherwise the system is defined as being stable. Henceforward we state that our system is dynamically stable if no close happens during the age of our system， about ±5 Gyr. Incidentally， this may be by one in which an of any orbital crossing between either of a pair of takes This is because we know from experience that an orbital crossing is very likely to lead to a close in and systems (Yoshinaga， Kokubo & Makino 1999). Of course this statement cannot be simply applied to systems with stable orbital resonances such as the Neptune–Pluto system.

　　 studies and aims of this research

　　In to the vagueness of the concept of stability， the in our system show a character typical of dynamical chaos (Sussman & Wisdom 1988， 1992). The cause of this chaotic is now partly understood as being a result of resonance (Murray & Holman 1999; Lecar， Franklin & Holman 2001). However， it would require integrating over an ensemble of systems including all nine for a covering several 10 Gyr to thoroughly understand the long-term of orbits， since chaotic dynamical systems are characterized by their strong dependence on initial

　　From that point of view， many of the long-term numerical included only the outer five (Sussman & Wisdom 1988; Kinoshita & Nakai 1996). This is because the orbital of the outer are so much longer than those of the inner four that it is much easier to follow the system for a given At present， the longest numerical published in journals are those of Duncan & Lissauer (1998). Although their main target was the effect of post-main-sequence mass loss on the stability of orbits， they performed many covering up to ～1011 yr of the orbital of the four jovian The initial orbital elements and masses of are the same as those of our system in Duncan & Lissauer's paper， but they decrease the mass of the Sun gradually in their numerical experiments. This is because they consider the effect of post-main-sequence mass loss in the paper. Consequently， they found that the crossing time-scale of orbits， which can 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 jovian remain stable over 1010 yr， or perhaps longer. Duncan & Lissauer also performed four experiments on the orbital of seven (Venus to Neptune)， which cover a span of ～109 yr. Their experiments on the seven are not yet comprehensive， but it seems that the terrestrial also remain stable during the ， maintaining almost

　　On the other hand， in his semi-analytical theory (Laskar 1988)， Laskar finds that large and can appear in the and of the terrestrial ， especially of Mercury and Mars on a time-scale of several 109 yr (Laskar 1996). The results of Laskar's theory should be confirmed and investigated by fully numerical

　　In this paper we present preliminary results of six long-term numerical on all nine orbits， covering a span of several 109 yr， and of two other covering a span of ± 5 × 1010 yr. The total time for all is more than 5 yr， using several dedicated PCs and One of the fundamental of our long-term is that system seems to be stable in terms of the Hill stability above， at least over a time-span of ± 4 Gyr. Actually， in our numerical the system was far more stable than what is defined by the Hill stability : not only did no close happen during the ， but also all the orbital elements have been confined in a narrow both in time and frequency domain， though are stochastic. Since the purpose of this paper is to exhibit and overview the results of our long-term numerical ， we show typical example figures as evidence of the very long-term stability of system For readers who have more specific and deeper interests in our numerical results， we have prepared a webpage ( )， where we show raw orbital elements， their low-pass filtered results， of elements and momentum deficit， and results of our simple time–frequency analysis on all of our

　　In 2 we briefly our dynamical model， numerical method and initial used in our 3 is devoted to a of the quick results of the numerical Very long-term stability of system is apparent both in and orbital elements. A rough of numerical errors is also given. 4 goes on to a of the longest-term of orbits using a low-pass filter and includes a of momentum deficit. In 5， we present a set of numerical for the outer five that spans ± 5 × 1010 yr. In 6 we also discuss the long-term stability of the and its possible cause.

　　2 of the numerical

　　(本部分涉及比較複雜的積分計算，作者君就不貼上來了，貼上來了起點也不一定能成功顯示。)

　　2.3 Numerical method

　　We utilize a -order Wisdom–Holman symplectic map as our main method (Wisdom & Holman 1991; Kinoshita， Yoshida & Nakai 1991) with a special start-up procedure to reduce the error of angle variables，‘warm start’(Saha & Tremaine 1992， 1994).

　　The stepsize for the numerical is 8 d throughout all of the nine (N±1，2，3)， which is about 1/11 of the orbital of the innermost (Mercury). As for the of stepsize， we partly follow the numerical of all nine in Sussman & Wisdom (1988， 7.2 d) and Saha & Tremaine (1994， 225/32 d). We rounded the decimal part of the their stepsizes to 8 to make the stepsize a multiple of 2 in order to reduce the of round-off error in the processes. In to this， Wisdom & Holman (1991) performed numerical of the outer five orbits using the symplectic map with a stepsize of 400 d， 1/10.83 of the orbital of Jupiter. Their result seems to be enough， which partly justifies our method of determining the stepsize. However， since the of Jupiter (～0.05) is much smaller than that of Mercury (～0.2)， we need some care when we compare these simply in terms of stepsizes.

　　In the of the outer five (F±)， we fixed the stepsize at 400 d.

　　We adopt Gauss' f and g in the symplectic map together with the third-order Halley method (Danby 1992) as a solver for Kepler The number of maximum we set in Halley's method is 15， but they never reached the maximum in any of our

　　The interval of the data output is 200 000 d (～547 yr) for the of all nine (N±1，2，3)， and about 8000 000 d (～21 903 yr) for the of the outer five (F±).

　　Although no output filtering was done when the numerical were in process， we applied a low-pass filter to the raw orbital data after we had completed all the See 4.1 for more detail.

　　2.4 Error

　　2.4.1 errors in total energy and momentum

　　 to one of the basic properties of symplectic integrators， which conserve the physically conservative quantities well (total orbital energy and momentum)， our long-term numerical seem to have been performed with very small errors. The averaged errors of total energy (～10?9) and of total momentum (～10?11) have remained nearly constant throughout the (Fig. 1). The special startup procedure， warm start， would have reduced the averaged error in total energy by about one order of magnitude or more.

　　 numerical error of the total momentum δA/A0 and the total energy δE/E0 in our numerical ± 1，2，3， where δE and δA are the absolute change of the total energy and total momentum， respectively， andE0andA0are their initial values. The horizontal unit is Gyr.

　　Note that different operating systems， different mathematical libraries， and different hardware architectures result in different numerical errors， through the in round-off error handling and numerical algorithms. In the upper panel of Fig. 1， we can this in the numerical error in the total momentum， which should be rigorously preserved up to machine-ε

　　2.4.2 Error in longitudes

　　Since the symplectic maps preserve total energy and total momentum of N-body dynamical systems inherently well， the degree of their may not be a good measure of the of numerical ， especially as a measure of the error of ， i.e. the error in longitudes. To estimate the numerical error in the longitudes， we performed the following procedures. We compared the result of our main long-term with some test ， which span much shorter but with much higher than the main For this purpose， we performed a much more with a stepsize of 0.125 d (1/64 of the main ) spanning 3 × 105 yr， starting with the same initial as in the N?1 We consider that this test provides us with a ‘pseudo-true’ of orbital Next， we compare the test with the main ， N?1. For the of 3 × 105 yr， we see a difference in mean anomalies of the Earth between the two of ～0.52°(in the case of the N?1 ). This difference can be to the value ～8700°， about 25 of Earth after 5 Gyr， since the error of longitudes increases linearly with time in the symplectic map. ， the longitude error of Pluto can 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. at the raw data

　　In this we briefly review the long-term stability of orbital through some snapshots of raw numerical data. The orbital of indicates long-term stability in all of our numerical : no orbital crossings nor close between any pair of took

　　3.1 General of the stability of orbits

　　First， we briefly look at the general character of the long-term stability of orbits. Our interest here focuses on the inner four terrestrial for which the orbital time-scales are much shorter than those of the outer five As we can see clearly from the orbital shown in Figs 2 and 3， orbital of the terrestrial differ little between the initial and final part of each numerical ， which spans several Gyr. The solid lines denoting the present orbits of the lie almost within the swarm of dots even in the final part of (b) and (d). This indicates that throughout the entire the almost of orbital remain nearly the same as they are at present.

　　Vertical view of the four inner orbits (from the z -axis ) at the initial and final parts of the ±1. The axes units are au. The xy - is set to the invariant of system total 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 denote the present orbits of the four terrestrial (taken from DE245).

　　The of and orbital for the inner four in the initial and final part of the N+1 is shown in Fig. 4. As expected， the character of the of orbital elements does not differ significantly between the initial and final part of each ， at least for Venus， Earth and Mars. The elements of Mercury， especially its ， seem to change to a significant extent. This is partly because the orbital time-scale of the is the shortest of all the ， which leads to a more rapid orbital than other ; the innermost may be nearest to instability. This result appears to be in some agreement with Laskar's (1994， 1996) that large and appear in the and of Mercury on a time-scale of several 109 yr. However， the effect of the possible instability of the orbit of Mercury may not fatally affect the global stability of the whole system owing to the small mass of Mercury. We will briefly the long-term orbital of Mercury later in 4 using low-pass filtered orbital elements.

　　The orbital of the outer five seems rigorously stable and quite over this time-span (see also 5).

　　3.2 Time–frequency maps

　　Although the exhibits very long-term stability defined as the non-existence of close events， the chaotic nature of dynamics can change the and amplitude of orbital gradually over such long time-spans. Even such slight of orbital in the frequency domain， in the case of Earth， can potentially have a significant effect on its surface climate system through (cf. Berger 1988).

　　To give an overview of the long-term change in in orbital ， we performed many fast Fourier (FFTs) along the time axis， and superposed the resulting to draw two- time–frequency maps. The specific approach to drawing these time–frequency maps in this paper is very simple – much simpler than the wavelet analysis or Laskar's (1990， 1993) frequency analysis.

　　Divide the low-pass filtered orbital data into many fragments of the same length. The length of each data segment should be a multiple of 2 in order to apply the FFT.

　　Each fragment of the data has a large part: for example， when the ith data begins from t=ti and ends at t=ti+T， the next data segment ranges from ti+δT≤ti+δT+T， where δT?T. We continue this until we reach a certain number N by which tn+T reaches the total length.

　　We apply an FFT to each of the data fragments， and obtain n frequency diagrams.

　　In each frequency diagram obtained above， the strength of can be by a grey-scale (or colour) chart.

　　We perform the ， and connect all the grey-scale (or colour) charts into one graph for each The horizontal axis of these new graphs should be the time， i.e. the starting times of each fragment of data (ti， where i= 1，…， n). The vertical axis represents the (or frequency) of the of orbital elements.

　　We have adopted an FFT because of its overwhelming speed， since the amount of numerical data to be into frequency components is terribly huge (several tens of Gbytes).

　　A typical example of the time–frequency map created by the above procedures is shown in a grey-scale diagram as Fig. 5， which shows the of in the and of Earth in N+2 In Fig. 5， the dark area shows that at the time indicated by the value on the abscissa， the indicated by the ordinate is stronger than in the lighter area around it. We can from this map that the of the and of Earth only changes slightly over the entire covered by the N+2 This nearly trend is qualitatively the same in other and for other ， although typical frequencies differ by and element by element.

　　4.2 Long-term exchange of orbital energy and momentum

　　We very long- and exchange of orbital energy and momentum using filtered elements L， G， H. G and H are equivalent to the orbital momentum and its vertical component per unit mass. L is to the orbital energy E per unit mass as E=?μ2/2L2. If the system is completely linear， the orbital energy and the momentum in each frequency bin must be constant. Non-linearity in the system can cause an exchange of energy and momentum in the frequency domain. The amplitude of the lowest-frequency should increase if the system is unstable and breaks down gradually. However， such a symptom of instability is not prominent in our long-term

　　In Fig. 7， the total orbital energy and momentum of the four inner and all nine are shown for N+2. The upper three panels show the long- of total energy (denoted asE- E0)， total momentum ( G- G0)， and the vertical component ( H- H0) of the inner four from the low-pass filtered elements.E0， G0， H0 denote the initial values of each quantity. The absolute difference from the initial values is plotted in the panels. The lower three panels in each figure showE-E0，G-G0 andH-H0 of the total of nine The shown in the lower panels is virtually entirely a result of the massive jovian

　　Comparing the of energy and momentum of the inner four and all nine ， it is apparent that the amplitudes of those of the inner are much smaller than those of all nine : the amplitudes of the outer five are much larger than those of the inner This does not mean that the inner terrestrial subsystem is more stable than the outer one: this is simply a result of the smallness of the masses of the four terrestrial compared with those of the outer jovian Another thing we notice is that the inner subsystem may unstable more rapidly than the outer one because of its shorter orbital time-scales. This can be seen in the panels denoted asinner 4 in Fig. 7 where the longer- and are more apparent than in the panels denoted astotal 9. Actually， the in theinner 4 panels are to a large extent as a result of the orbital of the Mercury. However， we cannot neglect the from other terrestrial ， as we will see in subsequent

　　4.4 Long-term coupling of several neighbouring pairs

　　Let us see some individual of orbital energy and momentum expressed by the low-pass filtered elements. Figs 10 and 11 show long-term of the orbital energy of each and the momentum in N+1 and N?2 We notice that some form apparent pairs in terms of orbital energy and momentum exchange. In ， Venus and Earth make a typical pair. In the figures， they show negative in exchange of energy and positive in exchange of momentum. The negative in exchange of orbital energy means that the two form a closed dynamical system in terms of the orbital energy. The positive in exchange of momentum means that the two are simultaneously under certain long-term Candidates for perturbers are Jupiter and Saturn. Also in Fig. 11， we can see that Mars shows a positive in the momentum to the Venus–Earth system. Mercury exhibits certain negative in the momentum versus the Venus–Earth system， which seems to be a caused by the of momentum in the terrestrial subsystem.

　　It is not clear at the moment why the Venus–Earth pair exhibits a negative in energy exchange and a positive in momentum exchange. We may possibly this through observing the general fact that there are no terms in semimajor axes up to -order theories (cf. Brouwer & Clemence 1961; & 1998). This means that the orbital energy (which is directly to the semimajor axis a) might be much less affected by perturbing than is the momentum exchange (which to e). Hence， the of Venus and Earth can be disturbed easily by Jupiter and Saturn， which results in a positive in the momentum exchange. On the other hand， the semimajor axes of Venus and Earth are less likely to be disturbed by the jovian Thus the energy exchange may be limited only within the Venus–Earth pair， which results in a negative in the exchange of orbital energy in the pair.

　　As for the outer jovian subsystem， Jupiter–Saturn and Uranus–Neptune seem to make dynamical pairs. However， the strength of their coupling is not as strong compared with that of the Venus–Earth pair.

　　5 ± 5 × 1010-yr of outer orbits

　　Since the jovian masses are much larger than the terrestrial masses， we treat the jovian system as an independent system in terms of the study of its dynamical stability. Hence， we added a couple of trial that span ± 5 × 1010 yr， including only the outer five (the four jovian plus Pluto). The results exhibit the rigorous stability of the outer system over this long time-span. Orbital (Fig. 12)， and of and (Fig. 13) show this very long-term stability of the outer five in both the time and the frequency domains. Although we do not show maps here， the typical frequency of the orbital of Pluto and the other outer is almost constant during these very long-term ， which is demonstrated in the time–frequency maps on our webpage.

　　In these two ， the numerical error in the total energy was ～10?6 and that of the total momentum was ～10?10.

　　5.1 Resonances in the Neptune–Pluto system

　　Kinoshita & Nakai (1996) integrated the outer five orbits over ± 5.5 × 109 yr . They found that four major resonances between Neptune and Pluto are maintained during the whole ， and that the resonances may be the main causes of the stability of the orbit of Pluto. The major four resonances found in research are as follows. In the following ，λ denotes the mean longitude，Ω is the longitude of the ascending node and ? is the longitude of Subs P and N denote Pluto and Neptune.

　　Mean resonance between Neptune and Pluto (3:2). The critical argument θ1= 3 λP? 2 λN??P librates around 180° with an amplitude of about 80° and a of about 2 × 104 yr.

　　The argument of of Pluto ωP=θ2=?P?ΩP librates around 90° with a of about 3.8 × 106 yr. The dominant of the and of Pluto are synchronized with the of its argument of This is anticipated in the theory constructed by Kozai (1962).

　　The longitude of the node of Pluto referred to the longitude of the node of Neptune，θ3=ΩP?ΩN， and the of this is equal to the of θ2 When θ3 zero， i.e. the longitudes of ascending nodes of Neptune and Pluto ， the of Pluto maximum， the minimum and the argument of 90°. When θ3 180°， the of Pluto minimum， the maximum and the argument of 90° again. Williams & Benson (1971) anticipated this type of resonance， later confirmed by ， Nobili & Carpino (1989).

　　An argument θ4=?P??N+ 3 (ΩP?ΩN) librates around 180° with a long ，～ 5.7 × 108 yr.

　　In our numerical ， the resonances (i)–(iii) are well maintained， and of the critical arguments θ1，θ2，θ3 remain during the whole (Figs 14–16 ). However， the fourth resonance (iv) appears to be different: the critical argument θ4 alternates and over a 1010-yr time-scale (Fig. 17). This is an interesting fact that Kinoshita & Nakai's (1995， 1996) shorter were not able to disclose.

　　6

　　What kind of dynamical mechanism maintains this long-term stability of the system? We can immediately think of two major features that may be responsible for the long-term stability. First， there seem to be no significant lower-order resonances (mean and ) between any pair among the nine Jupiter and Saturn are close to a 5:2 mean resonance (the famous ‘great inequality’)， but not just in the resonance zone. Higher-order resonances may cause the chaotic nature of the dynamical ， but they are not so strong as to destroy the stable within the lifetime of the real system. The feature， which we think is more important for the long-term stability of our system， is the difference in dynamical distance between terrestrial and jovian subsystems (Ito & Tanikawa 1999， 2001). When we measure by the mutual Hill radii (R_)， among terrestrial are greater than 26RH， whereas those among jovian are less than 14RH. This difference is directly to the difference between dynamical features of terrestrial and jovian Terrestrial have smaller masses， shorter orbital and wider dynamical They are strongly perturbed by jovian that have larger masses， longer orbital and narrower dynamical Jovian are not perturbed by any other massive bodies.

　　The present terrestrial system is still being disturbed by the massive jovian However， the wide and mutual among the terrestrial renders the disturbance ineffective; the degree of disturbance by jovian is O(eJ)(order of magnitude of the of Jupiter)， since the disturbance caused by jovian is a forced having an amplitude of O(eJ). Heightening of ， for example O(eJ)～0.05， is far from sufficient to provoke instability in the terrestrial having such a wide as 26RH. Thus we assume that the present wide dynamical among terrestrial (> 26RH) is probably one of the most significant for maintaining the stability of the system over a 109-yr time-span. Our detailed analysis of the between dynamical distance between and the instability time-scale of system is now on-going.

　　Although our numerical span the lifetime of the system， the number of is far from sufficient to fill the initial phase space. It is necessary to perform more and more numerical to confirm and examine in detail the long-term stability of our dynamics.

　　——以上文段引自 Ito， T.& Tanikawa， K. Long-term and stability of orbits in our System. Mon. Not. R. Astron. Soc. 336， 483–500 (2002)

　　這只是作者君參考的一篇文章，關於太陽系的穩定性。

　　還有其他論文，不過也都是英文的，相關課題的中文文獻很少，那些論文下載一篇要九美元(《Nature》真是暴利)，作者君寫這篇文章的時候已經回家，不在檢測中心，所以沒有數據庫的使用權，下不起，就不貼上來了。