Frequency analysis and the representation of slowly diffusing planetary solutions

0.52 MB·

by  Yanning Fu

Astronomy&Astrophysicsmanuscriptno.naffdifN_c (cid:13)cESO2017
January26,2017
Frequency analysis and the representation of slowly diffusing
planetary solutions
Y.N.Fu1 andJ.Laskar2
1 PurpleMountainObservatory,ChineseAcademyofSciences,2WestBeijingRoad,Nanjing210008,China
e-mail:[email protected]
2 ASD,IMCCE-CNRSUMR8028,ObservatoiredeParis,UPMC,77Av.Denfert-Rochereau,75014Paris,France
e-mail:[email protected]
Received;accepted
7
1 ABSTRACT
0
2 Context. Over short time intervals planetary ephemerides have been traditionally represented in analytical form as finite sums of
periodictermsorsumsofPoissontermsthatareperiodictermswithpolynomialamplitudes.Nevertheless,thisrepresentationisnot
n
welladaptedfortheevolutionoftheplanetaryorbitsinthesolarsystemovermillionofyearsastheypresentdriftsintheirmain
a
J frequencies,duetothechaoticnatureoftheirdynamics.
Aims. The aim of the present paper is to develop a numerical algorithm for slowly diffusing solutions of a perturbed integrable
4
Hamiltoniansystemthatwillapplytotherepresentationofthechaoticplanetarymotionswithvaryingfrequencies.
2
Methods.Bysimpleanalyticalconsiderations,wefirstarguethatitispossibletorecoverexactlyasinglevaryingfrequency.Then,

afunctionbasisinvolvingtime-dependentfundamentalfrequenciesisformulatedinasemi-analyticalway.Finally,startingfroma
]
M numericalsolution,arecursivealgorithmisusedtonumericallydecomposethesolutiononthesignificantelementsofthefunction
basis.
I Results.Simpleexamplesshowthatthisalgorithmcanbeusedtogivecompactrepresentationsofdifferenttypesofslowlydiffusing
h. solutions.Asatestexample,weshowhowthisalgorithmcanbesuccessfullyappliedtoobtainaverycompactapproximationofthe
p La2004solutionoftheorbitalmotionoftheEarthover40Myr([-35Myr,5Myr]).Thisexamplehasbeenchosenasthissolutionis
- widelyusedforthereconstructionoftheclimatesofthepast.
o
r Key words. Celestial mechanics – Ephemerides – Chaos – Methods: numerical– planet and satellites : dynamical evolution and
t stability
s
a
[

1. Introduction order to compare the output of the numerical integrations with
1
v the quasiperiodic solutions of the perturbative methods, Laskar
Beforethecomputerages,thelongtermsolutionsfortheplan-
8 (1988) introduced the frequency analysis method that allowed
0etaryorbitswerederivedbyperturbationmethods,andobtained toobtaininaveryefficientwayapreciseapproximationofthe
in the form of sum of periodic terms. The first of such solu-
1 numericalsolutioninquasiperiodicform(e.g.Laskar(2005)).
7tions was obtained by Lagrange (1782) who only considered
Butalthoughthegoalwastosearchforthemostpreciselong
0the known planets of the time, that is the planets that are visi-
term solution for the planetary orbits, one of the outcomes of
.blebynakedeye(Mercury,Venus,Earth,Mars,JupiterandSat-
1 thesecomputationswastodemonstratethatthesolarsystemmo-
urn).ThiswaslateronimprovedbyLeVerrier(1840,1841)who
0 tionischaotic(Laskar1989,1990).Asaconsequence,thesolu-
took also Uranus into account. These long term solutions have
7 tions are not quasiperiodic, although they can be approximated
1revealedtobeoffundamentalimportancefortheunderstanding
byaquasiperiodicexpressionoveralimitedtimeofafewmil-
:of the past climate of the Earth, when it was understood that
v lionofyears(Laskar1988,1990;Laskaretal.2004).
the changes of the orbit of the Earth induce also some change
i In the present work, we derive a more adapted strategy for
Xin its obliquity and in the insolation at the Earth surface (Mi-
theslowlydiffusingtrajectoriesofadynamicalsystemthatwill
lankovitch1941)(foradetailedreview,seeLaskaretal.(2004)).
r
aWiththeadventofcomputers,twodifferentapproachesbecome beverywellsuitedtotheconstructionofcompactformsforthe
longtimebehaviouroftheplanetaryorbits.Weintroduceanal-
possible. Computer algebra allowed to extend the perturbation
gorithm,thatisderivedfromthefrequencyanalysis,butwherea
methods (e.g. Bretagnon 1974), and as the computer speed in-
slowvariationofthefrequenciesandamplitudesisadded.Here,
creased,directintegrationsofthefullSolarsystembecomepos-
by a slowly diffusing solution, we mean a solution that, while
sible(Quinnetal.1991;Sussman&Wisdom1992).Meanwhile,
experiencing significant frequency drifts in the whole consid-
Laskar(1985,1986,1988)developedamixedstrategy,withan
ered time interval, is nearly quasiperiodic in time subintervals.
analyticalaveragingoftheplanetaryequationsobtainedbyper-
Similar to frequency analysis, a key step of our algorithm is to
turbation methods using dedicated computer algebra. This ana-
constructfrequency-dependentfunctionbasisonwhichthecon-
lytical averaging was then followed by a numerical integration
sideredsolutionwillbedecomposed.
of the averaged system with high order multistep method. In
In section 2, after reminding some fundamental results on
Sendoffprintrequeststo:J.Laskar frequency analysis, we consider a single term model and show
Articlenumber,page1of14
A&Aproofs:manuscriptno.naffdifN_c
that it is possible to recover exactly its varying frequency as a recovering KAM solutions. Actually, most relevant works are
functionoftime.Thisisimportantwhenoneisinterestedinthe aboutfrequencydrifts.Inparticular,thevariationsofthefunda-
detailsofhowasolutiondiffusesinthefrequencyspace.Butfor mentalfrequenciesofthesecularsolarsystemover200Myrwas
constructingacompactsolutionrepresentationwithareasonably used to study the chaotic behavior of the solar system (Laskar
highprecision,smallflickerswithoutsignificantcumulativeef- 1990).
fectscanbesmoothedoutfromavaryingfrequency.Therefore, Letusfirstrecallbrieflythetheoremabouttheconvergence
it is preferable to use a model with limited number of parame- oftheNAFFalgorithm(forthecompleteversionwithproof,see
ters, e.g. low-order polynomials, to approximate the frequency. (Laskar2005)).Withasetof N appropriatevariables,thecom-
Todoso,wesampletheaveragefrequenciesoveraslidingtime plex function (1) describing a degree of freedom of an analytic
interval. For a slowly diffusing solution, it is assumed that all KAMsolutioncanbewrittenas
fundamental frequencies can be sampled in this way by using
theNAFFalgorithm(e.g.(Laskar2005)). (cid:88)
Insection3,ageneralalgorithmofrepresentingaslowlydif- z(t)=eßν1t+ akeß(cid:104)k,ν(cid:105)t (|ak|<1), (2)
fusingsolutionisdesigned,whereafrequency-dependentfunc- k∈ZN−(1,0,...,0)
tionbasisisconstructedbasedontheChebyshevapproximations
whereν isafundamentalfrequency.Wehavethen(e.g.Laskar
of fundamental frequencies. The basis functions with signifi- 1
1999)
cant but non-fundamental frequencies are generated according
Theorem1LetνT bethevalueofσ∈Rthatmaximizesthe
to the assumption that, at any instant, a slowly diffusing solu- 1
functionφ(σ)=|(cid:104)z(t),eßσt(cid:105)χ|,whereχ=χ(t/T)>0isaweight
tionisnearlyquasi-periodicwiththesmoothedfundamentalfre- T
quenciesatthatinstant,namelyallmainfrequenciesareintegral function,and(cid:104)·,·(cid:105)χT theinnerproductoftwocomplexfunctions
linearcombinationsofthefundamentalfrequencies.Thiseffec- oftdefinedas
tively avoids the difficulties in determining the frequencies of 1 (cid:90) T
long-periodtermsand/orgroupsofneighboring-periodterms. (cid:104)f(t),g(t)(cid:105)χ = f(t)g¯(t)χ(t/T)dt, (3)
T 2T
Twosimpleexamplesaregiveninsection4.Inthefirstex- −T
ample,aweaklydissipatedsystemisconsidered,andtherepre- wethenhavelim νT =ν .(cid:3)
sentedsolutiondiffusesbecauseofdissipation.Inthesecondex- T→∞ 1 1
Similartothefundamentalfrequencyν ,anyothermainfre-
1
ample,aHamiltoniansystemofdegree1.5isconsidered,anda
quencycanberecoveredbysearchingitsneighborhoodorRfor
solutionstartingfromanobviousresonanceoverlapzoneisrep-
thevalueofσthatmaximizesφ(σ)definedwiththeremaining
resented. To show the flexibility of our algorithm in represent-
z(t),thatis,theoriginalz(t)minustherecoveredterms.
ingsuchasolution,weexcludeallpossiblelibrationfrequencies A diffusion solution is then characterized by fundamental
fromthesetoffundamentalfrequenciesusedinconstructingthe
frequencyvariations(Laskaretal.1992;Laskar1993).Regard-
functionbasis.
ing to the determination of a varying frequency, let’s first con-
Applicationstotherepresentationofephemeridesoftheso-
siderthefollowingsimplestcase,
lar system bodies are provided in section 5. As examples, the
eccentricityandtheinclinationoftheEarth,asgivenbythelong- f(t)=a(t)eß(cid:82)0tν(τ)dτ, (4)
termnumericalsolutionLa2004(Laskaretal.2004),arerepre-
sented.Forsucharealisticsolution,itisnaturaltotakeintoac- whereν(t)anda(t)>0arerealintegrablefunctions.Writing
counttherestrictionsontherepresentationmodelfromprevious
raensaullytssi,se(.eg..gt.hLeaismkaprorettaanlt.l2ib0r0a4t)io,nsofrtehqautethnecireessufrlotimngnruempreersiecna-l φ(σ(t))=(cid:12)(cid:12)(cid:12)(cid:12)(cid:104)f(t),eß(cid:82)0tσ(τ)dτ(cid:105)χT(cid:12)(cid:12)(cid:12)(cid:12)=(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)21T (cid:90)−TTa(t)eßθ(t)χ(t/T)dt(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12), (5)
tationscanbeascloseaspossibletothephysicalmodel.
whereθ(t) = (cid:82)t(ν(τ)−σ(τ))dτbelongstoC0,thesetofcontin-
0
uousfunctiononR,wehavethefollowingtheorem
2. Frequencyanalysisandtime-dependent
Theorem2ForanygivenT >0andν(t)∈C0,thefunctional
frequencies φ(σ(t))hasoneandonlyonemaximuminC0,whichisattained
atσ(t)≡ν(t).(cid:3)
ForaKAMsolutionofadynamicalsystem(Kolmogorov1954;
Proof.Theright-handsideof(5)canbewrittenas
Arnol’d1963;Moser1962),withfundamentalfrequencyvector
dνo=m(νo1f,t.h..e,νsNy)st∈emR,Nt,hwehmeoretioNniosfthaenynugmivbeenrdoefgdreeegroefesfroefedfroeme-, (cid:12)(cid:12)(cid:12)(cid:12) 1 (cid:90) T((cid:112)a(t))((cid:112)a(t)eßθ(t))χ(t/T)dt(cid:12)(cid:12)(cid:12)(cid:12)=(cid:12)(cid:12)(cid:12)(cid:104)(cid:112)a(t), (cid:112)a(t)e−ßθ(t)(cid:105)χ(cid:12)(cid:12)(cid:12).
withassociatedvariablez(t),canbedescribedincomplexform (cid:12)(cid:12)2T −T (cid:12)(cid:12) (cid:12) T(cid:12)
as
By using the Cauchy-Bunyakovsky-Schwarz inequality1, it is
easytodeducefromthisformulathat
(cid:88)
z(t)= akeß(cid:104)k,ν(cid:105)t, (1) (cid:112) (cid:112) 1 (cid:90) T
k∈ZN φ(σ(t))≤(cid:107) a(t)(cid:107)×(cid:107) a(t)e−ßθ(t))(cid:107)= a(t)χ(t/T)dt
2T
−T
wherek = (k ,...,k )isthefrequencyindexvector,a complex
amplitude,an1d(cid:104)·,·(cid:105)Ndenotestheusualinnerproductokftworeal where (cid:107) · (cid:107) denotes the norm induced by the inne√r product
vectors. √(3). Moreover, the equality is true if and only if a(t) and
a(t)e−ßθ(t) are linearly dependent, i.e. θ(t) is a constant. This
The NAFF (Numerical Analysis of Fundamental Frequen-
cies) algorithm was designed to numerically recover signifi- 1 Givenanytwovectors f andginaninnerproductspace,theCauchy-
cant terms of (1) from a numerical sample set of z(t), e.g. Bunyakovsky-Schwarzinequalitywrites(cid:104)f,g(cid:105)≤(cid:107)f(cid:107)×(cid:107)g(cid:107),where(cid:104)·,·(cid:105)
an ephemeris obtained by numerical integration (Laskar 1988). denotesinnerproductand(cid:107)·(cid:107)theinducednorm.Also,theequalityis
However,theapplicationofNAFFisnotlimitedtonumerically trueifandonlyif f andgarelinearlydependent.
Articlenumber,page2of14
Fu&Laskar:Frequencyanalysis&solutionrepresentation
constantis0,sinceθ(0)=0bydefinition.Theabovearguments Table1.Summaryofprocedureparameters.
imply that φ(σ(t)) has a unique maximum attained at θ(t) ≡ 0.
Meaning Venue
Sinceν(t)iscontinuousandthesearchedσ(t)isalsocontinuous,
θ(t) ≡ 0 is equivalent to σ(t) ≡ ν(t). Therefore, it is concluded
foranygivenT >0thatφ(σ(t))hasoneandonlyonemaximum M Degreeofthepolynomial (7)
attainedatσ(t)≡ν(t).(cid:3) n approximatingthenthfundamentalfrequency
Theorem 2 implies that, even in the case of a varying fre-
Maximumindexnumber
quency,itisstillpossibletorecoverexactlythefrequency,asa Kn associatedwiththenthfundamentalfrequency (9)
functionoftime.Now,letusconsiderthefollowingmoregeneral
complexfunction Degreeofthepolynomial
L (10)
f(t)=a1(t)eß(cid:82)0tν1(τ)dτ+ (cid:88) ak(t)eß(θk+(cid:82)0t(cid:104)k,ν(τ)(cid:105)dτ), (6) k approximatingtheamplitudewithindexk
J Maximumnumberofrepresentationterms
k∈ZN−(1,0,...,0)
where a (a ) ∈ S with S a subspace of real function space,
1 k a a δ Absolutetruncationerror (15)
θ ∈ R, and {ν }N ∈ S with S a linear subspace of the real
k n n=1 ν ν
integrable function space. As an extension of (2), this function δ Relativetruncationerror (16)
r
inheritsanimportanttimevaryingcharacterof(2),i.e.allphase
incrementsaredescribedbyasinglefundamentalfrequencyvec-
torν(t).Wewillmaketheheuristicassumptionthat,undersuit-
able conditions, a similar result as Theorem 1 holds for more
generalexpressionwithvaryingfrequenciesasin(6). where c ∈ R is a Chebyshev coefficient, x = 2(t−τ1) − 1 ∈
Assumption:If(cid:107)ν1T −ν1(cid:107)issufficientlysmall,whereν1T ∈Sν [−1,1] am,nnormalized time, and T (x) the Chebys(hτΛe−vτ1)polyno-
and(cid:107)·(cid:107)isanappropriatelydefinedT-dependentnorm(e.g.the m
one induced by the inner product (3)), and σ(t) = νT(t) (lo- mial2 of degree m. We then construct numerically a frequency-
cally)maximizesthefunctionalφ(σ) = |(cid:104)f(t),eß(cid:82)0tσ(τ)dτ(cid:105)1χT|,then dtheepecnodnesindtefruendcstoioluntiboansiwsiBll(bseeedetchoemnpexotsesdu.bsection), on which
limT→∞ν1T(t)=ν1(t). The final step, i.e. decomposing zn(t) on B, is the same
as that of the NAFF algorithm (Laskar 1999), except for the
searchedfunctionbasesforsignificantterms.Thefunctionbases
3. Representationofslowlydiffusingsolutions
are {eßωkt,ωk ∈ R} for NAFF and B for the presently described
From now on, we will restrict ourselves to slowly diffusing so- procedure.
lutions of an ordinary differential equation system obtained by
slightly perturbing an integrable Hamiltonian system. We de-
note by (I,θ) = {I ,θ }N the action-angle variables of the in- 3.2. Representationmodel
n n n=1
tegrableHamiltoniansystem,whichwillbeusedtoexpressthe Byvariationofparameters,(1)becomes
ephemerisofadiffusingsolution{zn(t)≡ In(t)eßθn(t)}nN=1. zv(t)= (cid:88) ak(t)eß(cid:82)0t(cid:104)k,ν(τ)(cid:105)dτ. (8)
3.1. Representationprocedure k∈ZN
Now,let
We start from a sample set of {z (t)}N . Suppose that the sam-
ples are given, respectively, at gnridnp=1oints from t to t with ZˆN ={k :|k |≤ K ∈Z+}N (9)
0 1 n n n n=1
fixed time step h, which is much smaller than the minimum of
beatruncatedsetofthefrequencyindexvectork.And,foreach
the fundamental periods. In accordance with the condition of
slow diffusion, we assume that the solution is close to quasi- k∈ZˆN,let
periodic in below mentioned time subintervals (of [t0,t1]) with (cid:88)Lk
lengthmuchlargerthanthemaximumofthefundamentalperi- ak(t)≈ a(cid:102)l,kTl(x(t)) (10)
ods.Theseroughlystatedpreconditionsarerequiredbecausewe l=0
willuseNAFFalgorithmtoestimatethechangingfundamental
frequencies. be a Chebyshev expansion of the amplitude valid on [τ1,τΛ],
As the first step, we apply NAFF algorithm to obtain fun- where a(cid:102)l,k ∈ C. These lead us from (8) to the following repre-
sentationmodelofz (t)
damental frequency samples. For each given degree (n), this v
aevlgeonrliythsmpacweidlltibmeeaspupbliinetde,rvraelssp{e[cτtλiv−eld2y,,τtλo+znd2(]t})Λλ=s1amofp[let0s,to1v]e=r zv(t)≈ (cid:88) (cid:88)Lk a(cid:102)l,kTl(x(t))eß(cid:82)0t(cid:104)k,ν(τ)(cid:105)dτ. (11)
[τ1−d2,τΛ+d2].Foreachofthesesubintervals,thefirstrecovered k∈ZˆN l=0
frequencywillbetakenastheaveragedvalueofthefundamental
FromtheChebyshevapproximationsofthefundamentalfre-
frequencyν overthesamesubinterval.The N averagedfunda-
n quencies, it is easy to obtain by integration the Chebyshev ex-
mentalfrequenciesobtainedinthiswayarethentakenasthein-
pansion representing the phase increment associated with the
stantonesatτ ,resultinginthefundamentalfrequencysamples
{τ ,ν (τ )}Λ ,λn=1,...N. frequency (cid:104)k,ν(τ)(cid:105), from the value, as is preferred, at the mid-
λ Sencoλndλly=,1wefitforeachgivendegreethefrequencysamples dleofthetimeinterval[τ1,τΛ],
toaChebyshevexpansionvalidon[τ1,τΛ], (cid:90) t (cid:104)k,ν(τ)(cid:105)dτ=ϕ +ϕ (t) (k∈ZˆN) (12)
0k k
ν (t)=(cid:88)Mn c T (x) (n=1,...,N), (7) (τ1+τΛ)/2
n m,n m 2 TheexplicitexpressionsofT (x)uptom=15aregivenin(A.8).
m=0 m
Articlenumber,page3of14
A&Aproofs:manuscriptno.naffdifN_c
whereϕ ∈ R,andϕ (t)gathersallChebyshevpolynomialsof
0k k
degreelargerthan0.Withϕ andϕ (t),(11)canbewrittenas
0k k
(cid:88) (cid:88)Lk
zv(t)≈ al,kTl(x(t))eßϕk(t), (13)
k∈ZˆN l=0
whereal,k = a(cid:102)l,keßϕ0k ∈C issimplythecoordinateofzv(t)when
itisdecomposedonthefunctionbasis
B={bl,k}={Tl(x(t))eßϕk(t)}. (14)
For the obtained representation, it should be noted that,
thoughwedecomposeasolutiononthebasisB,andcorrespond-
ingly, express its representation by (13), ϕ (t) will not be used
k
to specify the representation, {ν (t)}N will be used instead. In
n n=1
other words, the representation will be specified by the coor-
dinates a , and the Chebyshev coefficients of {ν (t)}N . This
l,k n n=1
choice is made for the following two reasons. One is that the
representation could otherwise be unnecessarily cumbersome.
The other reason is that, from the Chebyshev coefficients of
Fig.1.Phasetrajectoryofthedissipatedsolutionspecifiedby(18),(19)
{ν (t)}N ,thoseofϕ (t)canbeeasilyobtainedaccordingto(12).
n n=1 k and(20). Alsoshownisthephaseorbitof theunperturbedpendulum
Here, one should also be reminded of the fact that the solution system,namely(18)withε=0,whichpassesthroughtheinitialphase
representation,becauseofitsdependenceontheChebyshevap- pointofthedissipatedsolution.Thedashedlinedepictsthespepratrix
proximations of νn(t) and ak(t), is valid on [τ1,τΛ] rather than oftheunperturbedpendulum.
[t ,t ].
0 1
Tocompletethedescriptionoftherepresentationmodel,we
point out that the integers {M ,K }N and {L } , introduced
n n n=1 k k∈ZˆN
respectively in (7), (9) and (10), are necessary parameters for
defining a particular representation procedure. Of course, the
number of these parameters can be reduced by requiring that
some or all of these parameters take the same sufficiently large
value, e.g. L = L for all k ∈ ZˆN, at the price of unnecessar-
k
ily increasing the basis dimension. Another important point to
mention is that, for terminating the procedure, one needs to set
beforehand the maximum number of representation terms (J)
and/or the required precision. The required precision are spec-
ified by using either the absolute truncation error δ or the rela-
tive truncation error δ . Correspondingly, the procedure will be
r
terminatedif
(cid:107)z−z (cid:107)<δ (15)
v
or
(cid:107)z−z (cid:107)/(cid:107)z(cid:107)<δ , (16)
v r Fig.2.Thechangingfundamentalfrequencyν(t)ofthedissipatedso-
wherethemodule(cid:107)·(cid:107)isinducedbytheinnerproduct(3)with lutionspecifiedby(18),(19)and(20).AlsoshownisaChebyshevap-
proximationofν(t).
prescribedχ≡1.Inthefollowing,theabove-mentionedparam-
eters,assummarizedinTable1,willbereferredtoasprocedure
parameters.
4.1. Anexampleofdissipatedsolution
Considerthefollowingweaklydissipatedsystem
4. Examples
Inthefollowingtwosubsections,wewillillustrateourrepresen- d2θ =sin(θ)−εdθ (18)
tation procedure with two slowly diffusing solutions, of which dt2 dt
onediffusesduetoadissipativeperturbation,andtheotherdue
where
toitschaoticnature.Inordertoillustratetheconvergenceprop-
ertyofthisprocedure,itisconvenienttowritetheresultingrep-
ε=10−5. (19)
resentationwithasingletermindex,i.e.,
(cid:88)J (cid:88)J If we nullify the dissipative term −εdθ in (18), the system is a
zv(t)= zj(t)≡ ajTl(j)(x(t))eßϕk(j)(t) (17) simplependulumwithHamiltonian Hdt(I) = I2 +cos(θ),where
j=1 j=1 I = θ˙. Take as an example the solutio0n z(t) ≡2I(t)eßθ(t) of (18)
where the terms are arranged in the same order as they are ob- withthefollowinginitialconditions
tainedintheprocedure,whichcorrespondsroughlytotheorder
ofdecreasing|a |. t =0,θ =0,I =1. (20)
j 0 0 0
Articlenumber,page4of14
Fu&Laskar:Frequencyanalysis&solutionrepresentation
Due to the presence of the dissipative perturbation, the phase Table2.ThecoefficientsoftheChebyshevexpansion(cid:80)9m=0cmTm(x(t))
orbit starting from (θ ,I ) decays gradually away from the un- approximating ν(t), the changing fundamental frequency of the dissi-
0 0
patedsolutionspecifiedby(18),(19)and(20).
perturbedorbitpassingthroughthesamephasepoint.Thisdis-
sipationeffectissignificantinthelongrun,asisshowninFig.1.
m c
Asasolutionofasystemof1degreeoffreedomwithadissi- m
pativeperturbation,thedecayingz(t)hasachangingfundamen- 0 +1.450265×100
talfrequencyν,whichinheritsthecharacteristicfrequencyofthe 1 −1.032502×10−1
unperturbedsystemandvariesasthesolutiondecays.
2 −7.924442×10−4
By numerical integration, an ephemeris of z(t) is obtained
at {t = ih : i = 0,...,65535,h = 0.15}. We then apply NAFF 3 −1.163806×10−4
i
to 129 evenly spaced time subintervals with length d = 2046h 4 −7.914948×10−6
andmidpoints{τ = d +496(λ−1)h}129,respectively,toobtain
λ 2 λ=1 5 −6.369786×10−7
asamplesetofthechangingfrequency,{τ ,ν(τ )}129.Thisfre-
quencysamplesetisshowninFig.2.Alsosλhownλisλa=1Chebyshev 6 −5.975314×10−8
expansionapproximatingν(t).Thisexpansionisofdegree9and 7 −3.609371×10−9
expressedas
8 +1.771340×10−9
(cid:88)9 9 +8.159405×10−9
ν(t)≈ c T (x), (21)
m m
m=0
wheretheChebyshevpolynomialsT (x)aregivenin(A.8),and
m
x(t)= τ21(2t9−−ττ11)−1isthetimeafterbeingnormalizedfrom[τ1,τ129] Table3.Theleading10termsof(25),arepresentationofthedissipated
tothe[−1,1]intervalonwhichtheChebyshevpolynomialsare solution z(t) specified by (18), (19) and (20), which is obtained with
defined.Thevaluesofthecoefficients{cm}9m=0arelistedinTable {M,K,L,δr}={9,10,9,10−5}
2. We write the Chebyshev expansion approximating the phase
increment, from the value at x = 0 (i.e. t = τ1+τ129), associated .
2
withνintwoparts j l(j) k(j) Arg(a ) |a |
j j
(cid:82)t ν(τ)dτ = τ129−τ1 (cid:82) x(t)(cid:80)9 c T (x)dx 1 0 1 −1.431544 1.378074489
(τ1+τ129)/2 2 0 m=0 m m (22) 2 0 2 0.278103 0.622837454
=ϕ T (x)+ϕ(x). 3 0 3 1.987348 0.159698128
0 0
4 1 1 1.748566 0.116187680
In this equation, the time scaling factor τ129−τ1 accounts for the
2 5 0 4 −2.587002 0.032622276
changeofintegrationvariablefromthephysicaltimetothenor-
malisedtime,ϕ0 =−τ1292−τ1ϕ(0)isarealconstant,and 6 0 -1 −1.709654 0.017747854
7 1 2 0.133848 0.028210547
τ −τ (cid:88)10
ϕ(x)= 129 1 C T (x) (23) 8 1 3 1.931736 0.027791194
m m
2
m=1 9 0 5 −0.878547 0.005904572
gathersallChebyshevpolynomialsofdegreelargerthan0with3 10 1 4 −2.628911 0.009823726
(cid:40)
C = 2c02−c2 (m=1), (24)
m cm−1−cm+1 (m=2,...,10).
2m
wherec =c =0. |cm/c0|<6×10−4 form>1.Buttheprecisionofthisrepresen-
10 11
tationisseveralorderslessprecisethanthepreviousone,asseen
Thiscompletesthedeterminationofϕ(x),whichwillbeused
inthefollowing.Withtheprocedureparameters{M,K,L,δ } = bycomparingthetwopanelsofFig.3.Thisisbecausethewhole
r
{9,10,9,10−5},weobtaina37-termrepresentation considered time interval is long, and so, the small discrepancy
inν(t)canresultinsignificantphaseerrors.Ontheotherhand,
(cid:88)37 2(t−τ ) however, the fact that there is no terms with either |k| > 8 or
z37(t)= ajTl(j)(x)eßϕk(j)(x) with x= τ −τ1 −1 (25) l > 4 in the 37-term representation indicates that the represen-
j=1 129 1 tationmodelpracticallyconvergeswithrespecttoK andL.Our
representation procedure also converges rapidly with respect to
ofthesolutionz(t),whereTl(j)(x)istheChebyshevpolynomial J.ThisisillustratedinFig.4,forwhichtheprocedureparameters
ofdegreel(j),andϕk(j)(x) = k(j)ϕ(x)thephaseincrementasso- excludingJarefixedas{M,K,L}={9,10,9}.
ciatedwiththefrequencyk(j)ν.Thedataneededforspecifying
thefirst10representationtermsaregiveninTable3.
Now, we discuss how the procedure parameters affect the 4.2. AnexampleofchaoticsolutionsofHamiltoniansystem
precisionoftheresultingrepresentation.Forthis,let’sfirstcon-
sider a 150-term representation obtained by resetting M = 1. ConsiderthefollowingHamiltoniansystem,
Thisresettingintroducesonlyasmalldiscrepancyinν(t),since
I2
3 forthegeneralcase,see(A.7). H(I,θ,t)= +εcos(θ)[1+cos(ν1t)+cos(ν2t)], (26)
2
Articlenumber,page5of14
A&Aproofs:manuscriptno.naffdifN_c
Table5.Theleading40termsofthe100-termrepresentationz (t),asexpressedin(31).
3
No. (cid:104)k,f(cid:105) Abs(a )×106 Arg(a )(degree)
k k
1 g 18984 −68.812
5
2 g 16088 95.535
2
3 g 13041 17.345
4
4 g 9042 −54.984
3
5 g 4314 −175.273
1
6 g −r 2583 −92.093
4 1
7 g −r 2415 15.391
3 1
8 g +r +r 2377 −136.939
3 1 2
9 g +s −s 1934 −128.504
4 3 4
10 g 1498 160.676
6
11 2g −g 1393 82.417
1 5
12 g +g −g 1372 178.968
3 4 6
13 g −r 1298 140.905
1 2
14 g −s +s 1282 −87.475
3 3 4
15 g −r 1156 −127.879
2 2
16 g −2r 1153 15.739
4 1
17 g +r 1085 49.430
1 2
18 g −g +g −r 1028 −40.633
2 3 6 2
19 g −r 946 149.088
4 2
20 g +r 942 123.828
2 2
21 −g +2g −g +s 916 82.353
1 4 5 3
22 −g +2g −2r 903 63.075
3 4 1
23 g +s −s 824 150.057
3 3 4
24 g +2r 816 145.441
3 1
25 g +r 806 165.483
4 1
26 g +s −s −r 756 113.428
4 3 4 1
27 g −s +s +r 712 −121.353
1 3 4 2
28 g −s +s 697 −14.943
4 3 4
29 2g −g +2r 605 −21.586
3 4 2
30 g 577 −146.073
7
31 g +r 573 52.634
4 2
32 g +r 504 −21.947
3 2
33 g +g −g +r 432 152.449
1 5 7 1
34 −g +g +g −r 383 −28.070
1 2 5 2
35 g +r 275 −4.253
3 1
36 g −r 144 89.319
3 2
37 g −r −r 121 −36.562
2 1 2
38 g −r +r 104 −105.175
2 1 2
39 −g +g +g 92 −173.117
1 4 5
40 g −r 87 99.123
1 1
where Bytheanalysisofthefrequencymap(e.g.Laskar1999),de-
ε=5×10−3,ν1 =1.5,ν2 = π. (27) fithneepdhaasseI0p→ointν3 withθ0 = 0andshowninFig.5,weknowthat
2
Asolutionofthissystemhasthreefundamentalfrequencies.The t =0,θ =0,I =1.535, (28)
0 0 0
firsttwoaresimplytheforcingfrequenciesν andν ,whichare
1 2
non-commensurable. The other, denoted as ν , can be approxi- liesinachaoticzoneformedbyresonanceoverlap.Andtheso-
3
matedinthesimilarwayasintheprevioussubsection. lutionz(t)startingfromthispointischaotic.
Articlenumber,page6of14
Fu&Laskar:Frequencyanalysis&solutionrepresentation
Table6.Theleading40termsofthe100-termrepresentationζ (t),asexpressedin(31).
3
No. (cid:104)k,f(cid:105) Abs(b )×106 Arg(b )(degree)
k k
1 s 13774 107.587
5
2 s 8666 −62.318
3
3 s 4647 96.756
4
4 s 4085 27.817
1
5 s 3312 80.364
2
6 g −g +s 2745 −167.132
3 4 4
7 s +2r 2041 −44.496
2 2
8 s +r 1543 125.410
2 2
9 g −g +s 1530 −137.669
3 4 3
10 s +s −s −r 1469 91.798
1 3 4 2
11 s −r 1450 23.071
2 2
12 s +r 1417 −116.942
1 2
13 s 1333 110.029
6
14 s 889 9.186
7
15 s +r 646 −50.277
2 1
16 s 641 26.053
8
17 s −r 613 −178.905
3 1
18 s +s −s 532 −177.568
1 3 4
19 s −2r 518 101.955
1 2
20 s −r 484 106.762
3 2
21 g −g +s +r 481 −21.125
3 4 3 1
22 g −g +s +r 445 −5.831
3 4 2 1
23 s −s +s 364 72.700
2 3 4
24 s −s +s −r 344 31.282
2 3 4 2
25 s −r 341 20.734
2 1
26 s −2r 320 −51.578
2 2
27 g −g +s 315 −118.127
3 4 1
28 g −s +s +s −s 304 172.667
3 1 6 7 8
29 s −r 293 149.314
1 1
30 −s +s −s +2s 293 −76.858
4 6 7 8
31 s +2r −r 292 134.167
1 1 2
32 −g +s +s +r 285 −170.886
3 5 8 1
33 −g +g +s 268 47.758
3 4 2
34 −g +g +s −r 258 −27.836
3 4 2 2
35 −g +g +s 244 −18.592
3 4 4
36 g +s −s +r 230 109.121
4 6 8 1
37 g −g +s −r 221 −91.732
3 4 4 1
38 s +s −s −2r 217 173.446
2 5 8 1
39 −g +g +s −r 198 24.257
3 4 2 1
40 s +2r 78 77.412
1 1
An ephemeris of z(t) is obtained by the symplectic inte- callydifferenterrorsourcesinthepresentwayofapproximating
grator SBABc4 (Laskar & Robutel 2001) at {t = ih : i = a changing frequency. One is related to the NAFF process that
i
0,...,4095,h = 0.186058}. We then apply NAFF to 129 evenly givesthesamplesofthefrequency,whiletheotherrelatedtothe
spaced time intervals with length d = 512h and midpoints fitting process that leads to a Chebyshev approximation of the
{τ = d +28(λ−1)h}129,respectively.Thisgivesthesampleset frequency.Accordingly,weconsidertheChebyshevapproxima-
λ 2 λ=1
{τ ,ν }129, partly shown in Fig.6 together with a Chebyshev tionassufficientlygoodifitdeviatesfromthefrequencysample
λ 3,λ λ=1
approximation. It should be noted that there are two intrinsi- setmuchlessthanthesampleuncertainties.AsshowninFig.6,
Articlenumber,page7of14
A&Aproofs:manuscriptno.naffdifN_c
Table 4. The coefficients of Chebyshev expansion (cid:80) c T approxi-
m m m
mating the major fundamental frequencies (g ,...,g ,s ,..,s ) (arcsec
1 8 1 8
yr−1)ofthesolarsystemoverthetimeintervalfrom-35Myrto5Myr
withoriginatJ2000.Alsolistedaretheusedconstantlibrationfrequen-
ciesr andr (arcsecyr−1).
1 2
m cm m cm
0 5.59436858 0 −5.61412432
1 −0.02916946 1 0.01897975
2 0.00259949 2 0.00616405
3 −0.00858191 3 −0.00967633
4 −0.00538351 4 −0.00048289
5 −0.00216691 5 −0.00405450
6 −0.00311006 6 −0.00319098
7 0.00137912 7 0.00116236
g1 8 −0.00225709 s1 8 −0.00048593
9 0.00318074 9 0.00207233
10 0.00431717 10 0.00194425
11 −0.00288184 11 −0.00103824
12 −0.00189708 12 −0.00117501
13 0.00127879 13 0.00006311
14 0.00065723 14 0.00030964
15 −0.00017362 15 −0.00037496
0 7.45660678 0 −7.07313208 Fig.3.Errorsoftwodifferentrepresentationsofthedissipatedsolution
1 −0.00276205 1 0.04218433 specified by (18), (19) and (20). In the case of the upper panel, our
2 0.00177713 2 −0.00001315
3 0.00056064 3 −0.00225169 representation procedure is terminated with a 37-term representation,
4 −0.00002096 4 0.00386256 becausethisrepresentationreachestheprecisionrequirementδr =10−5.
5 0.00090782 5 −0.00280878 Inthecaseofthelowerpanel,withsligtlylargererrorinthefrequency
6 0.00096229 6 0.00247914 approximation,therepresentationisabout4orderslessprecise,though
7 −0.00024908 7 0.00109449
g2 8 0.00042855 s2 8 −0.00064733 thenumberofrepresentationtermsisincreasedto150.
9 −0.00071228 9 0.00229722
10 −0.00091971 10 −0.00000725
11 0.00061668 11 −0.00193901
12 0.00047006 12 0.00088350
13 −0.00024111 13 0.00053257
14 −0.00015053 14 0.00024897
15 0.00007778 15 0.00024128
0 17.36445990 0 −18.84810087
1 0.01585020 1 −0.00165582
2 −0.00496995 2 0.00007105
3 0.00540088 3 −0.00140634
4 0.00496032 4 0.00050744
5 −0.00226103 5 0.00039039
6 −0.00004476 6 0.00135948
7 0.00006943 7 0.00065306
g3 8 0.00007190 s3 8 −0.00041702
9 −0.00037475 9 −0.00026792
10 −0.00059252 10 −0.00018687
11 0.00025574 11 0.00004492
12 0.00035119 12 0.00013288
13 −0.00020600 13 0.00001794
14 −0.00000042 14 −0.00005746
15 0.00011279
0 17.91086281 0 −17.75496646
1 0.02044603 1 0.00670633
2 −0.00851446 2 −0.00717162
3 0.00812180 3 0.00313019
4 0.00866961 4 0.00807398 Fig.4.Convergencepropertyoftheprocedureofrepresentingthedis-
5 −0.00482118 5 −0.00498282
6 −0.00137335 6 −0.00116638 sipatedsolutionspecifiedby(18),(19)and(20).Theresidualofrepre-
7 0.00011721 7 0.00058378 sentationisplottedagainstthenumberJofrepresentationterms,where
g4 8 0.00083465 s4 8 0.00021726 otherprocedureparametersarefixedas{M,K,L}={9,10,9}.
9 −0.00041179 9 0.00016650
10 −0.00109334 10 −0.00060643
11 0.00069484 11 0.00014173
12 0.00045529 12 0.00039520 thisrequirementcanbemetwiththeChebyshevpolynomialof
1134 −−00..0000002194014720 1134 −−00..0000002046650908 degreeM3 =9.
15 0.00008189 15 0.00012665 Testcalculationsshowthatourrepresentationprocedurecan-
g5 0 4.25745185 s5 0 0.00000015 not lead to an acceptable representation of the solution on the
0 28.24498422 0 −26.34785292
1 0.00010582 1 0.00000053 whole sampling time interval. There are two possible reasons
g6 2 0.00011695 s6 2 −0.00000514 forthis.
3 −0.00002698
4 −0.00001157 Themostintrinsicreasonwouldbethatoursolutionexperi-
g7 01 −30..0008070905021476 s7 0 −2.99252583 ences passages into or out of resonance zones. Such a passage
isassociatedwithashiftbetweencirculationandlibrationofthe
0 0.67302182 0 −0.69173649
g8 s8 1 0.00000001 correspondingresonanceangle,andso,withoccurrenceordis-
r1 0 0.251085 r2 0 0.117222 appearance of certain terms. If some of these terms are signifi-
cantenoughinthewholeconsideredtimeinterval(τ ,τ ),then
1 129
therewouldbenowaytogetanyacceptablenon-piecewiserep-
Articlenumber,page8of14
Fu&Laskar:Frequencyanalysis&solutionrepresentation
Fig.5.Fundamentalfrequencymapofthesystemspecifiedby(26)and Fig. 7. Error of a 100-term representation of the chaotic solution
(27). In the chaotic zone formed by resonance overlap lies the phase specified by (26), (27) and (28). This representation is obtained with
point (θ0,I0) = (0,1.535), which will be chosen as the initial phase {M1,M2,M3,K1,K2,K3,Lk,J}={0,0,9,10,10,10,9,100}.
pointoftheconsideredchaoticsolution.
Fig.8.Theconvergencepropertyoftheprocedureofrepresentingthe
chaoticsolutionspecifiedby(26),(27)and(28).Theresidualofrepre-
Fig.6.Thechangingfundamentalfrequencyν3(t)ofthechaoticsolu- sentationisplottedagainstthenumber J ofrepresentationterms.The
tionspecifiedby(26),(27)and(28).Thesamplesofthisfrequencyis otherprocedureparametersarefixedas{M ,M ,M ,K ,K ,K ,L } =
1 2 3 1 2 3 k
computedbyapplyingNAFFalgorithmoveraslidingtimeinterval,and {0,0,9,10,10,10,9}.
theerrorsareestimatedastheirrespectivedifferencesfromthefrequen-
ciesofthequasiperiodicapproximationofthesolutionoverthesame
slidingtimeinterval.AlsoshownisaChebyshevappoximationofν3(t). interval,inthesensethateveryimportantlibrationfrequencyis
notfarfromatleastoneelementofthefrequencyset.
Settingtheprocedureparametersas
resentation. Therefore, we will restrict ourselves to the shorter
timeinterval(τ1,τ50). {M1,M2,M3,K1,K2,K3,Lk,J}={0,0,9,10,10,10,9,100},
Anotherpossiblereasonisthatthereisoneormoresignifi-
weobtaina100-termrepresentationofoursolution.Fig.7shows
cantlibrationfrequencies,whicharenottakenintoconsideration
that the errors of this representation are typically of order less
whenwegeneratethefunctionbasisB.Whileitiseasytomake than10−4I .And,Fig.8illustrates,inthesamewayasinthepre-
a necessary extension of B in order to include known libration 0
vioussubsection,theconvergencepropertyoftherepresentation
frequencies(seesection5),itisnotthatstraightforwardtoiden-
procedure.
tify and sample these frequencies (Laskar 1990). To show the
flexibilityofouralgorithm,wewillnotsearchforanylibration
frequencyandmakethecorrespondingextensionofB.Theflex-
5. Applicationtoplanetaryephemerides
ibility comes from the fact that, if we choose reasonably large
values of our procedure parameters K ’s, then the resulting set Numerical integration is now an efficient way of constructing
n
of frequencies would be dense enough over a large frequency ephemeridesofthesolarsystembodieswithhighprecision.For
Articlenumber,page9of14
A&Aproofs:manuscriptno.naffdifN_c
Following Laskar et al. (2004), we compute the fundamen-
tal frequencies of the secular solar system by applying the
NAFF algorithm over time intervals of length 20 Myr for the
proper modes (z•,...,z•,ζ•,...,ζ•), and 50 Myr for the proper
1 4 1 4
modes (z•,...,z•,ζ•,...,ζ•), respectively4. From -35 Myr to 5
5 8 5 8
Myr with step 0.1 Myr, we generate the samples of the funda-
mental frequencies (g ,...,g ,s ,...,s ) corresponding to these
1 8 1 8
proper modes. The resulting nominal value of s is about
5
0.00000015arcsecyr−1.Theotherfrequencysamplesareshown
respectivelyinthepanelsofFig.9,wheretheerrorsareestimated
as the difference between the values of a frequency computed
fromtheassociatedephemerisanditsquasiperiodicapproxima-
tion. Also shown in this figure are the Chebyshev approxima-
tionsofthesefundamentalfrequencies.AlloftheseChebyshev
approximations, the coefficients of which are listed in Table 4,
are obtained respectively by truncating the ones of degree 15.
Fig. 10. Convergence property of the two representation procedures
Thetruncationcriterionisroughlythatthediscrepancyinafun-
leadingrespectivelytothe100-termrepresentationsofz (t)andζ (t).
3 3 damental frequency should not induce an error in phase larger
TheresidualofrepresentationisplottedagainstthenumberJofrepre-
than2π/104overseveraltensofmillionyears.
sentationterms.
There are two important libration frequencies, i.e. r =
1
0.251085 arcsecyr−1 of the resonance argument 2((cid:36)• −(cid:36)•)−
4 3
practicalapplications,however,itcanbeusefultorepresentan- (Ω•−Ω•)andr =0.117222arcsecyr−1of2((cid:36)•−(cid:36)•)−(Ω•−
4 3 2 1 5 1
alytically,andthusinacontinuousway,thesediscretesolutions. Ω•).Toincludethemasadditionalfundamentalfrequencies,we
2
Theserepresentationscanbedoneintheformofsegmented expressamainfrequencyas
Chebyshev expansions, like the ones representing the classical
planetary ephemerides as INPOP (e.g. Fienga et al. 2008), or (cid:88)8 (cid:88)2
other generally applied approximation models without physi- ω= (migi+nisi)+ kjrj, (30)
cal basis. A drawback of these representations is that they re- i=1 j=1
quire large amount of data. In order to get compact represen-
Thed’Alembertcharacteristic,i.e.(cid:80)8 (m+n)=1,willbeused
tations, Chapront (1995) uses a model of Poisson series, with i=1 i i
toexcludenon-physicalfrequencyindexvectors.
fixed main frequencies, that were obtained with the NAFF al-
To show to what degree our algorithm is practically use-
gorithm.Thoughthismodelalreadyinvolvessomelong-termor
long-period-termeffectsbyallowingPoissonterms,andasthus, ful,wediscussherethefollowingtwo100-termrepresentations
((cid:60) forshort),
workswellwithplanetaryephemeridesoffiveouterplanetsover 100
afewhundredyears,itdoesnottakeintoconsiderationthefre-
(cid:88)100 (cid:88)100
quencydrifts,andcannotbeusedovermillionofyears. z3(t)= ak(j)eßϕk(j)(t), ζ3(t)= bk(j)eßϕk(j)(t), (31)
The frequency drifts are important over a few tens of mil-
j=1 j=1
lionyears,asshownbyLaskar(1990).Thealgorithmdeveloped
in the present paper should thus be more appropriate in repre- where, with f = (g ,...,g ,s ,...,s ,r ,r ) and k =
1 8 1 8 1 2
senting ephemerides spanning this long time interval. It is thus (m ,...,m ,n ,...,n ,k ,k ), the phase increment ϕ (t) =
interesting to test whether the present algorithm can be used to (cid:82)t(cid:104)1k,f(τ)8(cid:105)dτ1is ass8oci1ate2d with the main frequency (cid:104)k,kf(cid:105). The
represent over such a time scale the long-term numerical solu- w0aytheresidualsdecreasewiththeincreasingnumberofrepre-
tionofmajorsolarsystembodies,(e.g.Laskaretal.2004).For sentationtermsinthetworepresentationproceduresof(cid:60) are
100
this, we apply our algorithm to the eccentricity and inclination showninFig.10.
variablesoftheEarth,i.e., Theleading40termsofthesetworepresentationsaregiven,
respectively, in Tables 5 and 6, where the terms are reordered
z =e exp(i(cid:36) ) and ζ =sin(i /2)exp(iΩ ). (29)
3 3 3 3 3 3 accordingtotheirrealamplitudesanddataareroundedtoacon-
venientnumberofdigits.Comparingtheseresultswiththeones
To be more precise, e and i are the eccentricity and inclina-
3 3 given in tables 4 and 5 of Laskar (1990), we find that they are
tion, respectively, of the instantaneous orbit of the Earth-Moon
coherentwitheachotherinthesensethatallthemainfrequen-
barycenter, and (cid:36) and Ω are, respectively, the longitudes of
3 3 cies explicitly identified previously can be found in the present
theperihelionandofthenodeofthesameorbitwithrespectto
tables.
thefixedJ2000.0equatorialreferencesystem.
Theon-lineelectronicfilesassociatedwiththepresentpaper,
The chaotic behaviour of the terrestrial orbit certainly limit
i.e.z3R100.datandzeta3R100.dat,arethefullversionofTable
thetimespanoverwhichthisorbitcanbepreciselydetermined,
5 and Table 6, respectively. They are plain text tables provid-
butz andζ fromLa2004arereliableandpreciseatleastover
3 3 ingalltermsof(cid:60) intheformof(31).Togetherwiththedata
thetimespan[−35Myr,+5Myr](Laskaretal.2011).Therefore, 100
presented in Table 4 for computing f, these two tables can be
werestrictourrepresentationstothistimespan.Testcalculations
used to computethe eccentricity and inclination variables from
showthatouralgorithmcanleadtocompactandpreciserepre-
(cid:60) . The errors of (cid:60) as solution representation are shown
sentationsforbothdegreesoffreedom.Withdifferentprocedure 100 100
inFig.11.Fromthisfigure,weexpectthatthepresentalgorithm
parameters, the resulting representations contain very different shouldbeefficientinproducingcompactandpreciserepresenta-
terms.Thisconfirmstheflexibilityofthepresentalgorithm.
tionsoflong-termephemeridesofmajorsolarsystembodies.
In order to give representations as physical as possible,
we will resort to the knowledge we have for the solution. 4 See(Laskar1990)forthedefinitionofthepropermodes.
Articlenumber,page10of14