作者君在作品相关中其实已经解释过这个问题。
不过仍然有人质疑“你说得太含糊了”，“火星轨道的变化比你想象要大得多！”
那好吧，既然作者君的简单解释不够有力，那咱们就看看严肃的东西，反正这本书写到现在，嚷嚷着本书bug一大堆，用初高中物理在书中挑刺的人也不少。
以下是文章内容：
long-term integrations and stabilityplaary orbitsour solar system
abstract
we present the resultsvery long-term numerical integrationsplaary orbital motions over 109 -yr time-spans including all nin inspectionour numerical data shows that the plaary motion，leastour simple dynamical model， seemsbe quite stable even over this very lon lookthe lowest-frequency oscillations using a low-pass filter showsthe potentially diffusive characterterrestrial plaary motion， especially that o behaviourthe eccentricitymercuryour integrationsqualitatively similarthe results from jacques laskar's secular perturbation theory (e.g. emax 0.35 over ± 4 gyr). however， there areapparent secular increaseseccentricityinclinationany orbital elementsthe plas， which mayrevealedstill longer-term numerical i have also performed a coupletrial integrations including motionsthe outer five plas over the duration± 5 x 1010 yr. the result indicates that the three major resonancesthe neptunepluto system have been maintained over the 1011-yr time-span.
1 introduction
1.1definitionthe problem
the questionthe stabilityour solar system has been debated over several hundred years， since the era o problem has attracted many famous mathematicians over the years and has played a central rolethe developmentnon-linear dynamics and chao，do not yet have a definite answerthe questionwhether our solar systemstable opartly a resultthe fact that the definitionthe term ‘stability’vague whenis usedrelationthe problemplaary motionthe solais not easygive a clear， rigorous and physically meaningful definitionthe stabilityour solar system.
among many definitionsstability， hereadopt the hill definition (gladman 1993): actually thisnot a definitionstability， but define a systembeing unstable when a close encounter occurs somewherethe system， starting from a certain initial configuration (chambers， wetherill & boss 1996; ito & tanikawa 1999). a systemdefinedexperiencing a close encounter when two bodies approach one another withinareathe larger hil the systemdefinedbeinstate that our plaary systemdynamically stableno close encounter happens during the ageour solar system， about ±， this definition mayreplacedonewhichoccurrenceany orbital crossing between eithera pairplas takebecauseknow from experience thatorbital crossingvery likelyleada close encounterplaary and protoplaary systems (yoshinaga， kokubo & makino 1999).course this statement cannotsimply appliedsystems with stable orbital resonances suchthe neptunepluto system.
1.2previous studies and aimsthis research
in additionthe vaguenessthe conceptstability， the plasour solar system show a character typicaldynamical chaos (sussman & wisdom 1988， 1992). the causethis chaotic behaviournow partly understoodbeing a resultresonance overlapping (murray & holman 1999; lecar， franklin & holman 2001). however，would require integrating overensembleplaary systems including all nine plas for a period covering severalgyrthoroughly understand the long-term evolutionplaary orbits， since chaotic dynamical systems are characterizedtheir strong dependenceinitial conditions.
from that pointview， manythe previous long-term numerical integrations included only the outer five plas (sussman & wisdom 1988; kinoshita & nakai 1996). thisbecause the orbital periodsthe outer plas aremuch longer than thosethe inner four plas thatis much easierfollow the system for a given integratio present， the longest numerical integrations publishedjournals are thoseduncan & lissauer (1998). although their main target was the effectpost-main-sequence solar mass lossthe stabilityplaary orbits， they performed many integrations coveringto 1011of the orbital motionsthe four jovia initial orbital elements and massesplas are the samethoseour solar systemduncan & lissauer's paper， but they decrease the massthe sun graduallytheir numerical because they consider the effectpost-main-sequence solar mass lossth， they found that the crossing time-scaleplaary orbits， which cana typical indicatorthe instability time-scale，quite sensitivethe ratemass decreaseth the massthe suncloseits present value， the jovian plas remain stable over 1010 yr，perhap & lissauer also performed four similar experimentsthe orbital motionseven plas (venusneptune)， which cover a span109 yr. their experimentsthe seven plas are not yet prehensive， butseems that the terrestrial plas also remain stable during the integration period， maintaining almost regular oscillations.
on the other hand，his accurate semi-analytical secular perturbation theory (laskar 1988)， laskar finds that large and irregular variations can appearthe eccentricities and inclinationsthe terrestrial plas， especiallymercury and marsa time-scaleseveral 109(laskar 1996). the resultslaskar's secular perturbation theory shouldconfirmed and investigatedfully numerical integrations.
in this paperpresent preliminary resultssix long-term numerical integrationsall nine plaary orbits， covering a spanseveral 109 yr， andtwo other integrations covering a span± 5 x 1010 yr. the total elapsed time for all integrationsmore than 5 yr， using several dedicated pcs and wthe fundamental conclusionsour long-term integrationsthat solar system plaary motion seemsbe stabletermsthe hill stability mentioned above，least over a time-span± ，our numerical integrations the system was far more stable than whatdefinedthe hill stability criterion: not only didclose encounter happen during the integration period， but also all the plaary orbital elements have been confineda narrow region bothtime and frequency domain， though plaary motionsthe purposethis paperto exhibit and overview the resultsour long-term numerical integrations，show typical example figuresevidencethe very long-term stabilitysolar system plaar readers who have more specific and deeper interestsour numerical results，have prepared a webpage (access )， whereshow raw orbital elements， their low-pass filtered results， variationdelaunay elements and angular momentum deficit， and resultsour simple timefrequency analysisallour integrations.
in section 2briefly explain our dynamical model， numerical method and initial conditions usedour i 3devoteda descriptionthe quick resultsthe numerical i long-term stabilitysolar system plaary motionapparent bothplaary positions and orbita estimationnumerical errorsals 4 goesto a discussionthe longest-term variationplaary orbits using a low-pass filter and includes a discussionangular mometu section 5，present a setnumerical integrations for the outer five plas that spans ± 5 x 1010 yr.section 6also discuss the long-term stabilitythe plaary motion and its possible cause.
2 descriptionthe numerical integrations
（本部分涉及比较复杂的积分计算，作者君就不贴上来了，贴上来了起点也不一定能成功显示。）
2.3 numerical method
we utilize a second-order wisdomholman symplectic mapour main integration method (wisdom & holman 1991; kinoshita， yoshida & nakai 1991) with a special start-up procedurereduce the truncation errorangle variables，‘warm start’(saha & tremaine 1992， 1994).
the stepsize for the numerical integrations8 d throughout all integrationsthe nine plas (n±1，2，3)， whichabout 1/11the orbital periodthe innermost pla (mercury).for the determinationstepsize，partly follow the previous numerical integrationall nine plassussman & wisdom (1988， 7.2 d) and saha & tremaine (1994， 225/32 d).rounded the decimal partthe their stepsizes8make the stepsize a multiple2orderreduce the accumulationround-off errorthe putatio relationthis， wisdom & holman (1991) performed numerical integrationsthe outer five plaary orbits using the symplectic map with a stepsize400 d， 1/10.83the orbital period o result seemsbe accurate enough， which partly justifies our methoddetermining th， since the eccentricityjupiter (0.05)much smaller than thatmercury (0.2)，need some care whenpare these integrations simplytermsstepsizes.
in the integrationthe outer five plas (f±)，fixed the stepsize400 d.
we adopt gauss' f and g functionsthe symplectic map together with the third-order halley method (danby 1992)a solver for keple numbermaximum iterationssethalley's method15， but they never reached the maximumanyour integrations.
the intervalthe data output200 000 d (547 yr) for the calculationsall nine plas (n±1，2，3)， and about 8000 000 d (21 903 yr) for the integrationthe outer five plas (f±).
althoughoutput filtering was done when the numerical integrations wereprocess，applied a low-pass filterthe raw orbital data afterhad pleted all the c section 4.1 for more detail.
2.4 error estimation
2.4.1 relative errorstotal energy and angular momentum
accordingohe basic propertiessymplectic integrators， which conserve the physically conservative quantities well (total orbital energy and angular momentum)， our long-term numerical integrations seemhave been performed with very smal averaged relative errorstotal energy (109) andtotal angular momentum (1011) have remained nearly constant throughout the integration period (fig. 1). the special startup procedure， warm start， would have reduced the averaged relative errortotal energyabout one ordermagnitudemore.
relative numerical errorthe total angular momentum δa/a0 and the total energy δe/e0our numerical integrationsn± 1，2，3， whereandare the absolute changethe total energy and total angular momentum， respectively， ande0anda0are their initia horizontal unitgyr.
note that different operating systems， different mathematical libraries， and different hardware architectures resultdifferent numerical errors， through the variationsround-off error handling and numerica the upper panel o，can recognize this situationthe secular numerical errorthe total angular momentum， which shouldrigorously preservedto machine-e precision.
2.4.2 errorplaary longitudes
since the symplectic maps preserve total energy and total angular momentumn-body dynamical systems inherently well， the degreetheir preservation may nota good measurethe accuracynumerical integrations， especiallya measurethe positional errorplas， i.e. the errorplaar estimate the numerical errorthe plaary longitudes，performed the followin pared the resultour main long-term integrations with some test integrations， which span much shorter periods but with much higher accuracy than the main i this purpose，performed a much more accurate integration with a stepsize0.125 d (1/64the main integrations) spanning 3 x 105 yr， starting with the same initial conditionsin the consider that this test integration provideswith a ‘pseudo-true’ solutionplaary orbita，pare the test integration with the main integration， n1. for the period3 x 105 yr，see a differencemean anomaliesthe earth between the two integrations0.52°(in the casetheintegration). this difference canextrapolatedthe value 8700°， aboutrotationsearth after 5 gyr， since the errorlongitudes increases linearly with timethe symplecti， the longitude errorpluto canestimated12°. this value for plutomuch better than the resultkinoshita & nakai (1996) where the differenceestimated60°.
3 numerical results i. glancethe raw data
in this sectionbriefly review the long-term stabilityplaary orbital motion through some snapshotsraw numerica orbital motionplas indicates long-term stabilityallour numerical integrations:orbital crossings nor close encounters between any pairplas took place.
3.1 general descriptionthe stabilityplaary orbits
first，briefly lookthe general characterthe long-term stabilityplaar interest here focuses particularlythe inner four terrestrial plas for which the orbital time-scales are much shorter than thosethe outer fivcan see clearly from the planar orbital configurations shownfigs 2 and 3， orbital positionsthe terrestrial plas differ little between the initial and final parteach numerical integration， which spans severa solid lines denoting the present orbitsthe plas lie almost within the swarmdots eventhe final partintegrations (b) and (d). this indicates that throughout the entire integration period the almost regular variationsplaary orbital motion remain nearly the samethey arepresent.
vertical viewthe four inner plaary orbits (from the z -axis direction)the initial and final partsthe integrationsn±1. the axes units are au. the-planesetthe invariant planesolar system total angular momentum.(a) the initial part ofn+1 ( t = 00.0547 x9 yr).(b) the final part ofn+1 ( t = 4.9339 x84.9886 x9 yr).(c) the initial partn1 (t= 00.0547 x 109 yr).(d) the final part ofn1 ( t =3.9180 x93.9727 x9 yr).each panel， a total23 684 points are plotted withintervalabout 2190over 5.47 x 107. solid lineseach panel denote the present orbitsthe four terrestrial plas (taken from de245).
the variationeccentricities and orbital inclinations for the inner four plasthe initial and final partthe integration n+1shown i expected， the characterthe variationplaary orbital elements does not differ significantly between the initial and final parteach integration，least for venus， earthelementsmercury， especially its eccentricity， seemchangea significanpartly because the orbital time-scalethe plathe shortestall the plas， which leadsa more rapid orbital evolution than other plas; the innermost pla maynearest result appearsbesome agreement with laskar's (1994， 1996) expectaions that large and irregular variations appearthe eccentricities and inclinationsmercurya time-scaleseveral 109 yr. however， the effectthe possible instabilitythe orbitmercury may not fatally affect the global stabilitythe whole plaary system owingthe small mass o will mention briefly the long-term orbital evolutionmercury latersection 4 using low-pass filtered orbital elements.
the orbital motionthe outer five plas seems rigorously stable and quite regular over this time-span (see also section 5).
3.2 timefrequency maps
although the plaary motion exhibits very long-term stability definedthe non-existenceclose encounter events， the chaotic natureplaary dynamics can change the oscillatory period and amplitudeplaary orbital motion gradually over such lon such slight fluctuationsorbital variationthe frequency domain， particularlythe caseearth， can potentially have a significant effectits surface climate system through solar insolation variation (cf. berger 1988).
to giveoverviewthe long-term changeperiodicityplaary orbital motion，performed many fast fourier transformations (ffts) along the time axis， and superposed the resulting periodgramsdraw two-dimensional timefrequenc specific approachdrawing these timefrequency mapsthis papervery simple much simpler than the wavelet analysislaskar's (1990， 1993) frequency analysis.
divide the low-pass filtered orbital data into many fragmentsthe sam lengtheach data segment shoulda multiple2orderapply the fft.
each fragmentthe data has a large overlapping part: for example， when the ith data begins from t=ti and endst=ti+t， the next data segment ranges from ti+δt≤ti+δt+t， where δt?t.continue this division untilreach a certain number nwhich tn+t reaches the total integration length.
we applyffteachthe data fragments， and obtain n frequency diagrams.
in each frequency diagram obtained above， the strengthperiodicity canreplaceda grey-scale (or colour) chart.
we perform the replacement， and connect all the grey-scale (or colour) charts into one graph for each horizontal axisthese new graphs shouldthe time， i.e. the starting timeseach fragmentdata (ti， where i= 1，…， n). the vertical axis represents the period (or frequency)the oscillationorbital elements.
we have adoptedfft becauseits overwhelming speed， since the amountnumerical databe deposed into frequency ponentsterribly huge (several tensgbytes).
a typical examplethe timefrequency map createdthe above proceduresshowna grey-scale diagram a， which shows the variationperiodicitythe eccentricity and inclinationearthn+2 fig. 5， the dark area shows thatthe time indicatedthe valuethe abscissa， the periodicity indicatedthe ordinatestronger thanthe lighter area around it.can recognize from this map that the periodicitythe eccentricity and inclinationearth only changes slightly over the entire period coveredthe n+2 nearly regular trendqualitatively the sameother integrations and for other plas， although typical frequencies differ plapla and elementelement.
4.2 long-term exchangeorbital energy and angular momentum
we calculate very long-periodic variation and exchangeplaary orbital energy and angular momentum using filtered delaunay elements l， g， h. g and h are equivalentthe plaary orbital angular momentum and its vertical ponent per uni relatedthe plaary orbital energy e per unit masse=μ2/2l2.the systempletely linear， the orbital energy and the angular momentumeach frequency bin must bthe plaary system can causeexchangeenergy and angular momentumthe frequenc amplitudethe lowest-frequency oscillation should increasethe systemunstable and breaks dow， such a symptominstabilitynot prominentour long-term integrations.
i， the total orbital energy and angular momentumthe four inner plas and all nine plas are shown for integration n+2. the upper three panels show the long-periodic variationtotal energy (denoted ase- e0)， total angular momentum ( g- g0)， and the vertical ponent ( h- h0)the inner four plas calculated from the low-pass filtered delaunay ， g0，denote the initial valueseac absolute difference from the initial valuesplottedth lower three panelseach figure showe-e0，g-g0 andh-h0the totalnin fluctuation shownthe lower panelsvirtually entirely a resultthe massive jovian plas.
paring the variationsenergy and angular momentumthe inner four plas and all nine plas，is apparent that the amplitudesthosethe inner plas are much smaller than thoseall nine plas: the amplitudesthe outer five plas are much larger than thosethe inne does not mean that the inner terrestrial plaary subsystemmore stable than the outer one: thissimply a resultthe relative smallnessthe massesthe four terrestrial plas pared with thosethe outer jovia thingnoticethat the inner plaary subsystem may bee unstable more rapidly than the outer one becauseits shorter orbital canseenthe panels denoted asinner 4 i the longer-periodic and irregular oscillations are more apparent thanthe panels denoted astotal 9. actually， the fluctuationstheinner 4 panels area large extenta resultthe orbital variationth，cannot neglect the contribution from other terrestrial plas，we will seesubsequent sections.
4.4 long-term couplingseveral neighbouring pla pairs
letsee some individual variationsplaary orbital energy and angular momentum expressedthe low-pass filtered delaunaandshow long-term evolutionthe orbital energyeach pla and the angular momentumn+1 andi notice that some plas form apparent pairstermsorbital energy and angular momentu particular， venus and earth make a typica the figures， they show negative correlationsexchangeenergy and positive correlationsexchangeangula negative correlationexchangeorbital energy means that the two plas form a closed dynamical systemtermsthe orbita positive correlationexchangeangular momentum means that the two plas are simultaneously under certain long-termfor perturbers are jupiteri，can see that mars shows a positive correlationthe angular momentum variationthe venuseart exhibits certain negative correlationsthe angular momentum versus the venusearth system， which seemsbe a reaction causedthe conservationangular momentumthe terrestrial plaary subsystem.
itnot clearthe moment why the venusearth pair exhibits a negative correlationenergy exchange and a positive correlationangular momentu may possibly explain this through observing the general fact that there aresecular termsplaary semimajor axesto second-order perturbation theories (cf. brouwer & clemence 1961; boccaletti & pucacco 1998). this means that the plaary orbital energy (whichdirectly relatedthe semimajor axis a) mightmuch less affectedperturbing plas thanthe angular momentum exchange (which relatese). hence， the eccentricitiesvenus and earth can bedisturbed easilyjupiter and saturn， which resultsa positive correlationthe angular momentu the other hand， the semimajor axesvenus and earth are less likelybe disturbedthe jovia the energy exchange maylimited only within the venusearth pair， which resultsa negative correlationthe exchangeorbital energythe pair.
as for the outer jovian plaary subsystem， jupitersaturn and uranusneptune seemmake dynamica， the strengththeir couplingnotstrong pared with thatthe venusearth pair.
5 ± 5 x 1010-yr integrationsouter plaary orbits
since the jovian plaary masses are much larger than the terrestrial plaary masses，treat the jovian plaary systeman independent plaary systemtermsthe studyits dynamica，added a coupletrial integrations that span ± 5 x 1010 yr， including only the outer five plas (the four jovian plas plus pluto). the results exhibit the rigorous stabilitythe outer plaary system over this lon configurations (fig. 12)， and variationeccentricities and inclinations (fig. 13) show this very long-term stabilitythe outer five plasboth the time and the frequencdo not show maps here， the typical frequencythe orbital oscillationpluto and the other outer plasalmost constant during these very long-term integration periods， whichdemonstratedthe timefrequency mapsour webpage.
in these two integrations， the relative numerical errorthe total energy was 106 and thatthe total angular momentum was 1010.
5.1 resonancesthe neptunepluto system
kinoshita & nakai (1996) integrated the outer five plaary orbits over ± 5.5 x 109. they found that four major resonances between neptune and pluto are maintained during the whole integration period， and that the resonances maythe main causesthe stabilitythe orbit o major four resonances foundprevious research are a the following description，λ denotes the mean longitude，Ωthe longitudethe ascending node and the longitude o p and n denote pluto and neptune.
mean motion resonance between neptune and pluto (3:2). the critical argument θ1= 32 λnp librates around 180° withamplitudeabout 80° and a libration periodabout 2 x 104 yr.
the argumentperihelionpluto wp=θ2=pΩp librates around 90° with a periodabout 3.8 x 106 yr. the dominant periodic variationsthe eccentricity and inclinationpluto are synchronized with the librationits argument oanticipatedthe secular perturbation theory constructedkozai (1962).
the longitudethe nodepluto referredthe longitudethe nodeneptune，θ3=ΩpΩn， circulates and the periodthis circulationequalthe periodθbees zero， i.e. the longitudesascending nodesneptune and pluto overlap， the inclinationpluto bees maximum， the eccentricity bees minimum and the argumentperihelion bees 90°. whenbees 180°， the inclinationpluto bees minimum， the eccentricity bees maximum and the argumentperihelion bees 90° & benson (1971) anticipated this typeresonance， later confirmedmilani， nobili & carpino (1989).
an argument θ4=pn+ 3 (ΩpΩn) librates around 180° with a long period， 5.7 x 108 yr.
in our numerical integrations， the resonances (i)(iii) are well maintained， and variationthe critical arguments θ1，θ2，θ3 remain similar during the whole integration period (figs 1416 ). however， the fourth resonance (iv) appearsbe different: the critical argumentalternates libration and circulation over a 1010-yr time-scale (fig. 17). thisan interesting fact that kinoshita & nakai's (1995， 1996) shorter integrations were not abledisclose.
6 discussion
what kinddynamical mechanism maintains this long-term stabilitythe plaary system?can immediately thinktwo major features that mayresponsible for the long-ter， there seembesignificant lower-order resonances (mean motion and secular) between any pair among the nin and saturn are closea 5:2 mean motion resonance (the famous ‘great inequality’)， but not justthe resonanc resonances may cause the chaotic naturethe plaary dynamical motion， but they are notstrongto destroy the stable plaary motion within the lifetimethe real sola second feature， whichthinkmore important for the long-term stabilityour plaary system，the differencedynamical distance between terrestrial and jovian plaary subsystems (ito & tanikawa 1999， 2001). whenmeasure plaary separationsthe mutual hill radii (r_)， separations among terrestrial plas are greater than 26rh， whereas those among jovian plas are less tha differencedirectly relatedthe difference between dynamical featuresterrestrial and jovia plas have smaller masses， shorter orbital periods and wider dynamica are strongly perturbedjovian plas that have larger masses， longer orbital periods and narrower dynamica plas are not perturbedany other massive bodies.
the present terrestrial plaary systemstill being disturbedthe massive jovia， the wide separation and mutual interaction among the terrestrial plas renders the disturbance ineffective; the degreedisturbancejovian plaso(ej)(ordermagnitudethe eccentricityjupiter)， since the disturbance causedjovian plasa forced oscillation havingamplitudeo(ej). heighteningeccentricity， for example o(ej)0.05，far from sufficientprovoke instabilitythe terrestrial plas having such a wide separation aassume that the present wide dynamical separation among terrestrial plas (> 26rh)probably ohe most significant conditions for maintaining the stabilitythe plaary system over a 109-y detailed analysisthe relationship between dynamical distance between plas and the instability time-scalesolar system plaary motionnow on-going.
although our numerical integrations span the lifetimethe solar system， the numberintegrationsfar from sufficientfill the initial phasnecessaryperform more and more numerical integrationsconfirm and examinedetail the long-term stabilityour plaary dynamics.
以上文段引自 ito， t.& tanikawa， k. long-term integrations and stabilityplaary orbitsour sola， 483500 (2002)
这只是作者君参考的一篇文章，关于太阳系的稳定性。
还有其他论文，不过也都是英文的，相关课题的中文文献很少，那些论文下载一篇要九美元（《nature》真是暴利），作者君写这篇文章的时候已经回家，不在检测中心，所以没有数据库的使用权，下不起，就不贴上来了。