Efficient fourth order symplectic integrators for near-harmonic separable Hamiltonian systems

0.4 MB·English

by  Kristian Mads Egeris Nielsen

Efficient fourth order symplectic integrators
for near-harmonic separable Hamiltonian systems
KristianMadsEgerisNielsen∗
Abstract
Efficientfourthorder symplecticintegrators are proposedfor numericalintegration ofseparable Hamiltonian systems H(p,q) = T(p)+
V(q). Symmetricsplitting coefficients with five to nine stages are obtained by higher order decomposition of the simple harmonic oscillator.
TheperformanceofthemethodsisevaluatedforvariousHamiltoniansystems: Integrationerrorsarecomparedtothoseofacclaimedintegrators
composedbyS.Blanesetal.(2013),W.Kahanetal.(1999)andH.Yoshida(1990).Numericaltestsindicatethattheintegratorsobtainedinthis
5 paperperformsignificantlybetterthanpreviousintegratorsforcommonHamiltoniansystems.
1
0
Keywords: Hamiltoniansystem,symplecticintegration,simpleharmonicoscillator,optimizedfourthorder
2

b ntroduction undamental theorems
I. I II. F
e
F
Explicitsymplecticintegrators(SIs,singularSI)arenumer- Consider the general separable Hamiltonian H(q,p) =

9 ical integration schemes for separable Hamiltonian systems. V(q)+T(p)with p˙ =−∂H/∂qandq˙=∂H/∂pwhereqand p

SIsarewidelyusedinfieldsofcelestialmechanics,accelerator arecanonicalcoordinatesconventionallylabeledpositionand
]
h physics,moleculardynamics,and quantum chemistrydueto momentum respectively. Bydefining z = (q,p), the timeevo-
p their structure-preservingpropertiesand simpleimplementa- lutionofzcanbeexpressedby[2]:
- tion.
nt z(τ)=eτ(A+B)z(0) (1)
a
u Numerous SIs with increased accuracy and stability have where A and B are non-commutative operators associated
q been published since the first discovery of the fourth order with T and V, and τ is the time step. A set of real numbers
[ SI [1]. The perhaps most acclaimed splitting composition is exist(c ,c ,...,c )and (d ,d ,...,d ),referredtoassplittingco-
1 2 k 1 2 k
2 derivedbyH.YoshidawhoprovedthatSIsofarbitrarilyhigh efficients, that satisfy ∑ki=1ci = ∑ki=1di = 1 suchthat the time
even order exist and can be constructed using Lie algebra evolutionof zcan beapproximatedbyaproductofexponen-
v
5 [2]. Other more recent studies on the topic of high stability tialfunctions[2]:
SIs include force gradient schemes [3] and symplectic correc-
4
k
3 tors [4][5]. The integration orderis often associated with the eτ(A+B) =∏eciτAediτB+o τN+1 (2)
4 efficiency of numerical integration schemes. Nth order inte-
0 grators have their integration error reduced by a factor kN i=0 (cid:16) (cid:17)
1. when the time step is reducedby a factor k, making high or- where o(τN+1) denotes an error on the N+1 order term. It
0 derintegratorshighlyfavorableforsmalltimesteps. However, iscommonpracticetoexpandtheexponentialtermsbypower
5 higher order error terms cannot be neglected for large time series with respect to τ. Lie algebra can be applied to de-
1 steps. Symplecticcorrectorshave beenproposedforeliminat- termine the splitting coefficients which satisfy Eq. 2 exactly
v: inghighordererrorsofSIs,althoughthesemethodsareonly up to a given order. The time evolution of z(0) → z(τ) or
i competitivefor specificsystems. The force-gradientcomposi- (q0,p0)→(qk,pk)isexplicitlycomputableas[2]:
X
r tgieonnehraalssbeepeanrpabrolepoHseadmailstoanmiaonresyesffitecmiesntbsuptlirtteiqnugisrechseemvaelufoar- qi =qi−1+τci∂∂Tp (pi−1) (3)
a tionoftheforcegradient[3]. Theforce-gradientcompositions pi = pi−1−τdi∂∂Vq (qi)
benefit from all-positive splitting coefficients meaning that
for i = 1,2,...,k. Note that the intermediate steps (q,p ),j =
j j
onlyforwardstepsaretakeninphasespaceandaretherefore
1,2,...,k−1, referredto assubsteps, donotrepresentthe time
referred to as a forward SIs. Non-gradient splitting methods
evolution of (q ,p ). SIs are time reversibleif the symplectic
0 0
ofordershigherthantwoalwaysinvolvenegativecoefficients
coefficients are symmetric[2]. Time reversibilityimplies that
[5]. ToconstructmoreefficienthighorderSIsithasbeenpro-
z(0) → z(τ)can be exactlyrevertedby z(0) ← z(−τ). This
posedthat the absolute valueof negative coefficients and the
property is embraced throughout this paper. The symmetric
sumofpositivecoefficientsshouldbeminimized[5].
compositionintheformofEq. 3withd =0isreferredtoas
k
ABA[5]. Iftherolesof AandBinEq. 2areswitchedthetime
This paper presents a series of efficient fourth order SIs
evolutionisexplicitlycomputableas:
which are ideal for numerical integration of particle systems
with quadratic kinetic energy. Additionally, some of the pre- pi = pi−1−τdi∂∂Vq (qi−1) (4)
sented integrators exhibit higher order accuracy for systems qi =qi−1+τci∂∂Tp (pi)
thatresembletheharmonicoscillator. Thepresentedsplitting
coefficientsgenerallyhavesmallabsolutevaluesandarethere- fori=1,2,...,k. ThesymmetriccompositionintheformofEq.
forereferredtoasnear-forwardSIs. 4 with c = 0 is referred to as BAB. Table 1 summarizes the
k
∗[email protected]
1
symmetricstructuresoftheABA-andBAB-composition.Note by the entries Q in Table 2where h denotes the intermedi-
h,λ
thatthecanddcoefficientsareconsistentlyusedinconnection atestepnumber. ThisdecompositionfollowsdirectlyfromEq.
with the ∂T/∂p and ∂V/∂q term respectively. The computa- 6. A similar decomposition can be constructed for the ABA
tionalcostofSIsisapproximatelyproportionaltothenumber composition. Thepositioncontributionζλ foreachorderλis
ofderivativeevaluations(stages). Thenumber ofstagesisde- introducedas:
fined as s = k−1 for symmetric coefficients due to the first
same as last property, which implies that the first derivative ζλτλ ≡Qs,λτλ,λ=0,1,...,λmax (7)
evaluationisthe sameas the lastderivativeevaluationofthe wheresisthenumberofstages. Thenewpositioncanthenbe
previousiteration. calculatedasthesumofallpositioncontributions:
Table1:ABAandBABcompositionforsymmetriccoefficients qb =λ∑maxτλζλ+o τλmax+1 (8)
λ=0 (cid:16) (cid:17)
ABA BAB
Supposethataknownfunctionq(t)existssuchthatq(t0)=qa
d1 =dk−1 c1 =ck d1 =dk c1 =ck−1 and q(t0+τ) = qb+o(τλmax+1) then qb is expressed by the
d2 =dk−2 c2 =ck−1 d2 =dk−1 c2 =ck−2 powerseries:
. . . .
. . . .
. . . . q(t +τ)
0
ddkk−−21 ck−1 dk−1 cckk−−21 =qa+ q′(1t!0)τ1+ q′′2(!t0)τ2+...+ q(λλmmaxa)x(!t0)τλmax+o τλmax+1
dk =0 ck dk ck =0 =ζ0+ζ1τ1+ζ2τ2+...+ζλmaxτλmax+o τλmax+1 (cid:16) (cid:17)
(cid:16) (cid:17) (9)
IftheSIisof Nthorderthen
ethods
III. M
q(λ)(t )
Thissectionfirstdescribesadecompositionmethodwhich λ! 0 τλ =ζλτλ,λ=0,1,...,N (10)
followsdirectlyfromthecompositionofthesimpleharmonic for any function q(t) that qualifies as a time evolution in
oscillator. It is then shown that this decomposition can be H(p,q)=T(p)+V(q). Thesimpleharmonicoscillator(SHO):
used to satisfy general fourth order conditions and higher
p2 1
orderconditionsforthesimpleharmonicoscillator. H= + kq2 (11)
2m 2
Asingleiterationfrom(qa,pa)to(qb,pb)usingthes-stage isenforcedbyapplyingthecondition Q¯h,λ = ∂V/∂q(Qh,λ)in
BAB-composition(Eq. 4)canbeexpandedinto: Table2. TheSHOyieldsthewell-knownanalyticaltimeevolu-
tionofq:
p1 = pa−τd1∂∂Vq (qa)
q1=qa+τc1∂∂Tp (p1) q(t)=acos kt +bsin kt (12)
p = p −τd ∂V (q ) rm ! rm !
2 1 2∂q 1
... (5) Using Eq. 10 where ζ0,ζ1,...,ζ6 are evaluated in Table 2 the
following identities must be true for symmetric fourth order
qs =qs−1+τcs∂∂Tp (ps)=qb splittingcoefficients:
ps+1= ps−τds+1∂∂Vq (qs)= pb
h1 =ζ0 =qa
Suppose now that ∂T(p) ∂p = p m which applies to many 1 kb=ζ = pa (2c +c )
1! m 1 m 1 2
classicalsystemsincludingthesimpleharmonicoscillator. Eq.
(cid:14) (cid:14) −q1 ka=ζ =−kqa (d +d )(2c +c )
5thenyields: 2!m 2 m 1 2 1 2 (13)
3
−1 k 2b=ζ =−2c d kpa (c +c )
qq12==qq1a(+1τ+c1Cpm2a)+−τC22Dqa1m+−τ12∂∂DVq2(mq−a)1∂∂Vq (q1) 41!3!m(cid:16)k m2(cid:17)a=ζ4 =3 c1d2km2q12a (22mc21d11+2c22d1+c2d2)
q =q (1+C )−C q +τ2D m−1∂V (q ) (cid:16) (cid:17)
3 2 3 3 1 3 ∂q 2 Eliminationoftheconstantsaandbyieldstheidentities:
. (6)
.
. 1=2(d +d )(2c +c )
1 2 1 2
qs =qs−1(1+Cs)−Csqs−2+τ2Dsm−1∂∂Vq (qs−1) 1=12c1d2(c1+c2)/(2c1+c2) (14)
ps+1=mqs−τcqss−1 −τds+1∂∂Vq (qs) 1=24c1d2(2c1d1+2c2d1+c2d2)
ItisreadilyfoundthattheonlyrealsolutionsatisfyingEq. 14
where Cn = cn+1/cn and Dn = −cndn. The first same as last
istheoneoriginallyfoundby[1]:
propertyallowssubstitutionof pa by ps+1 yieldinganexpres-
sionthatisindependentof p. Eq. 6yieldssignificantlylarger d =d = 1 = 1 −d = 1 −d
truncation errors than Eq. 5 and is therefore not suitable for c1=c4= 2(21−21/=3) 1(12−c2) 2 3 (15)
practicalimplementation. However,Eq. 6allowsseparationof 1 3 2−21/3 2 2
τλtermsforanyorderλ=0,1,...,λmaxwhereλmaxisthemax- A closerexamination of Table 2reveals that m, k and qa only
imumnumberoferrortermswhichappearbydecomposition. appear as scalars of ζ0,ζ1,...,ζλmax the τ0, and that τ2 and τ4
Forinstanceconsiderthethreestagedecompositionexpressed separationscanbeexpressedasinTable3.
2
Table2:BABdecompositionforsymmetriccoefficientswhereC¯h=Ch+1,Q¯h,λ=m−1∂V/∂q(Qh,λ)andQh,λdenotesthetablevalueatindices(h,λ)
hτλ τ0 τ1 τ2 τ3 τ4 τ5 τ6
0 qa
1 Q0,0 c1pam−1 D1Q¯0,0
2 Q C¯ −C Q Q C¯ Q C¯ +D Q¯ D Q¯ D Q¯
1,0 2 2 0,0 1,1 2 1,2 2 2 1,0 2 1,1 2 1,2
3 Q C¯ −C Q Q C¯ −C Q Q C¯ −C Q +D Q¯ Q C¯ +D Q¯ Q C¯ +D Q¯ D Q¯ D Q¯
2,0 3 3 1,0 2,1 3 3 1,1 2,2 3 3 1,2 3 2,0 2,3 3 3 2,1 2,4 3 3 2,2 3 2,3 3 2,4
Table3:Decomposition of even orders expressed in terms of the splitting objectives κ and κ are always 0 if ∑k c = ∑k d = 1.
0 1 i=1 i i=1 i
coefficients.mandqaaresetto1 Additionally these identities eliminate one c and d variable
from the minimization procedure. The minimization is per-
hτλ τ0 τ2 τ4 formed using an adaptive simplex method [7] minimizing
0 1 only one objective κmax = max(κ1,κ2,...,κλ). Solutions are
1 1 provided in the following section with minimization criteria
1 1 − ∑ ci ∑ dj κmax < 1e−128 calculated using software-implemented ar-
i=1 j=1
bitrary precision arithmetic (2048bit, 256bit exponent). Four
2 i 1 1
2 1 − ∑ ci ∑ dj − ∑ Qm,2dm+1 ∑ ck unique sets of splitting coefficients are presented in the fol-
i=1 j=1 m=1 k=m+1 lowing section. The non-unique solutions are picked among
. . . .
.. .. .. .. thousands of solutions based on the sum of the absolute val-
h i h−1 h−1 uesofthesplittingcoefficientsandhigherorderSHOerrors.
h 1 − ∑ ci ∑ dj − ∑ Qm,2dm+1 ∑ ck
i=1 j=1 m=1 k=m+1
Symmetriccoefficientsomitoddordererrorterms[2],im-
plyingthatthefollowingidentitiesmustbetruetosatisfythe esults
IV. R
fourthconditionsfortheSHO:
ζ = q′′(0) =−1 =− ∑h c ∑i d A. Obtainedsplitting coefficients
2 2! 2! i j
i=1 j=1
(16)
ζ4= q′′4′′(!0) = 41! =−mh∑−=11Qm,2dm+1k=∑mh+1ck in TSapbllietti4nigncaorerfafiycifeonrtmsawtitfhorficvoentvoenniiennetsitmagpelsemareenptaretisoenntfeodr
variousprogramminglanguages.
It has previously been shown that Eq. 16 is in fact general
second and fourth order criteria [6]. This implies that split-
TwodifferentkindsofBABsolutionsareevaluated:
tingcoefficients satisfyingEq. 10to atleastfourthorderalso
satisfyEq. 2toatleastfourthorder. Compositionswithmore
• BAB: Solution to Eq. 10 using BAB decomposition(see
stagesandhigherorderSHOcriteriacanbeexploitedforsig-
Table2).
nificantly reducing higher order error terms. In this sense,
the higher order SHO criteria are used only as a constraint • BAB’:SolutiontoEq. 10usingABAdecompositionwith
for obtaining efficient splitting coefficients for more general coefficientsc0 =cs+1=0inEq. 3.
Hamiltoniansystems.
The difference between the BAB and BAB’ solutions is re-
Analyticalevaluationofsplittingcoefficientsformorethan
vealedbyEq. 5wheretheerrorof ps+1 remainsunaddressed
three stages quickly becomes hopeless. Instead a numerical
for the BAB solutionas opposedto BAB’ solutionwhere the
routine has been developed to satisfy the multi-objective op-
error of qs+1 remains unaddressed. Thus the BAB and BAB’
timization problem with variables c and d and minimization
solutions possess different constraints which may influence
objectiveκ:
theirperformanceineitherdirection.
κλ ≡ q(λ)(ti)−λ!ζλ(c,d) =0,λ=1,2,...,λH (17)
BAB and BAB’ solutions should be implemented accord-
(cid:12) (cid:12)
where λH(cid:12)(cid:12)is referred to as the(cid:12)(cid:12)SHO order constraint and κλ ingtoEq. 4. UniquesolutionssatisfyboththeABAandBAB
forλ=1,2,...,λ aretheminimizationobjectives. Numerical decompositions and can therefore be implemented as either,
H
experiments suggest that for λ > 5 every additional stage althoughtheperformanceofthepresentedsolutionsdepends
H
yields two increments of λmax and a single increment of λH ontheappliedscheme.
forwhichatleastonerealsolutionexists.
The naming convention of the obtained methods is
TheBABdecompositionisnumericallyevaluatedwiththe [scheme]s[stages]o[order]H, for instance BAB’s9o7H for a BAB-
choice of constants in Eq. 12: a = −b = 1 and m = k = 1 optimizedschemewith nine stagesand seventh harmonicor-
corresponding to the initial conditions qa = −pa = 1. The der.
3
B. Numericaltests
−2 Ruth [ABA]
s5odr4 [ABA]
Common methods used to benchmark SIs for integration
−4 ABA104 [ABA]
ofautonomousHamiltoniansystemsinclude[4][5][8]:
ABA864 [ABA]
r
ro −6 ABA1064 [ABA]
r
• ThemaximumHamiltonianerror: max(|∆H(t)/H0|) n e Yosh s7o6 A [ABA]
a −8
• ThemeanHamiltonianerror: mean(|∆H(t)/H |) ni
0 o
• Theaccumulatederrorofintegrals milt−10
a
H
m −12
where H istheinitialHamiltonian. Eachofthesebenchmark u
0 m
methodsareappliedthroughoutthissection. xi−14
a
m
ABAs5o6H A [ABA]
The following 10 symmetric SIs are benchmarked for three −16 BABs7o7H [BAB]
autonomous Hamiltonian systems: The simple harmonic os-
BAB’s8o7H [BAB]
cillator,theHénonHeilessystemandtheKeplerproblem: −18 BAB’s9o7H [BAB]

−3 −2.5 −2 −1.5 −1 −0.5
• Ruth[ABA]3-stage,4thorder[1] time per evaluation
• s5odr4[ABA]5-stage,4thorder[9]
• ABA104[ABA]7-stage,10-4generalizedorder[5]
• ABA864[ABA]7-stage,8-6-4generalizedorder[5] −2
• ABA1064[ABA]8-stage,10-6-4generalizedorder[5]
• Yoshs7o6A[ABA]7-stage,6thorder[2] −4
• ABAs5o6H[ABA]5-stage,4thorder
• BABs7o7H[BAB]7-stage,4thorder −6
• BAB’s8o7H[BAB]7-stage,4thorder or
r
• BAB’s9o7H[BAB]9-stage,4thorder n er −8
a
ni
o−10
where [ABA] and [BAB] denote suggested ABA- and BAB- milt
implementation respectively. It should be noted that the
a
generalizedorderschemesaredesignedforperturbedHamil- H−12
n
tonian systems H = HA+ǫHB with ǫ ≪ 1. It should also a
e−14
be noted Yosh s7o6 A yields the smallest errors of the three m
stage-optimized sixthorder integrators presentedby [2], and
−16
thatthisintegratorishereevaluatedtoreferencehigherorder
accuracyandstability.
−18
Numerical tests are performed using double precision
−3 −2.5 −2 −1.5 −1 −0.5
floating point format with compensated summation. Com-
time per evaluation
pensated summation drastically reduces truncation errors
for long-term symplectic integration (see [9] and references
Figure1:Integrator benchmark using the SHO in the time interval t ∈
therein). AllnumericalintegrationsareperformedusingZym-
[0,500]. Top: Maximum Hamiltonian error. Bottom: Average
plectic (v.1.01.00) from Zymplectic.com, which is a numerical relativeHamiltonianerror.Thedashedlinesshowbenchmarkpro-
frameworkdesignedfornumericalintegrationandbenchmark files for integrators of corresponding color using extended preci-
ofseparableHamiltoniansystems. sion. Thehorizontalaxescorrespondtoτ/s. Allaxesareinbase
10logarithmicunits.Legendappliestobothgraphs
Thesimpleharmonicoscillator
Consider the Hamiltonian of the one-dimensional SHO (Eq. The ABA- and BAB-decompositionfrom the previoussec-
11) with k = m = 1 and initial conditions (p,q) = (0,1). tioncanbeutilizedtodetermineerrorcontributionsofdiffer-
Figure 1 shows a benchmark profile of the 10 integrators for ent orders. More precisely, κλ can be evaluated for all λ to
theSHO.Everypointoneachlinecorrespondstoacomplete resolvethemagnitudeofhigherordererrorsfortheSHO.This
integration in a given time interval. The slopes of the lines canbeusedtocharacterizehigherordererrortermswhichof-
read the integrator order. It is noticable that the integrators tendivergerapidlyforhighorderstage-optimizedintegrators.
ABAs5o6H A, BABs7o7H, BAB’s8o7H and BAB’s9o7H be- While these error terms strictly apply to the SHO, they still
haveassixthorderintegrators. covergeneralequationspresentintheexpansionofEq. 2.
4
Q : Ruth
ABA Ruth [ABA]
−2
QBAB: Ruth s5odr4 [ABA]
Q : ABAs5o6H A ABA104 [ABA]
ABA
Q : ABAs5o6H A −4 ABA864 [ABA]
BAB r
106QABA: BAB’s8o7H erro AYBosAh1 0s76o46 [ AAB [AA]BA]
Q : BAB’s8o7H n −6
BAB a
Q : ABA1064 ni
ABA o
104QBAB: ABA1064 milt −8
Q : Yosh s7o6 A a
ABA H
Q : Yosh s7o6 A m −10
BAB
u
m
102 xi
a−12
λ m
κ ABAs5o6H A [ABA]
BABs7o7H [BAB]
−14 BAB’s8o7H [BAB]
0
10 BAB’s9o7H [BAB]

−3.5 −3 −2.5 −2 −1.5 −1 −0.5
time per evaluation
−2
10

−2
−4
10 −4
4 6 8 10 12 14 16 18
λ
r
o
r −6
r
e
Figure2:SHOdecompositionofselectedintegratorsforboththeABA-(cir- n
clemarker)andBAB-composition(solidline). Linesandmarkers a
ni −8
ofthesamecolorintersectexactlyatevenorders o
milt
Figure2showspositionerrortermsofselectedintegratorsfor a−10
H
the SHO. Notice that the obtained BAB methods have more n
a
error terms when implemented as ABA methods. Further- e
m−12
more the seventh harmonic order of BAB’s8o7H eliminates
theseventhordererrortermofthepositioncoordinateforthe
BAB implementation as opposed to the sixthharmonic order −14
integrator ABAs5o6H A. The largest integer value of λ with
no integration error reads the SHO integration order. Stage-
−16
optimizedintegratorsofhigherordersgenerallyhaverapidly
−3.5 −3 −2.5 −2 −1.5 −1 −0.5
diverginghigher ordererror contributions. This implies that
time per evaluation
theseintegratorsarenotsuitableforintegratingtheSHOwith
largetimesteps.
Figure3:IntegratorbenchmarkusingtheHénonHeilessysteminthetime
interval t ∈ [0,500]. Top: Maximum Hamiltonian error. Bot-
TheHénonHeilessystem tom:AveragerelativeHamiltonianerror. Thehorizontalaxescor-
TheHénon-HeilessystemisoneofthemoststudiedHamilto- respondtoτ/s.Allaxesareinbase10logarithmicunits.Legend
niansystemsduetoitsapplicationsincelestialmechanics,its appliestobothgraphs
non-integrable properties,chaoticdynamics and fractal struc-
tures. Considerthe Hamiltonian of the Hénon-Heiles system
TheKeplerproblem
withinitialconditions[qx,qy,px,py]=[0.3,0.0,0.0,0.4]:
Symplectic integration has been widely used for high accu-
H= p2x+2p2y + qx2+2qy2 +qx2qy− 31qy3 (18) racy simulations of the Solar System (see [8] and references
therein). IthasbeenproposedthatMercuryduetoitsfastor-
where in this case H = 1/8, inside the chaotic domain. Fig- bital periodand higheccentricityislimitingthe performance
ure3shows abenchmark profileof the 10 integratorsforthe of SIs when integrating the Solar System [8]. Consider the
HénonHeilessystem.TheHénonHeilessystemcanbeconsid- HamiltonianoftheplanarSun-Mercurysystem:
eredasaperturbedSHOforsmallvaluesofqxandqy. Forthis
reasonitappears that especiallyBAB’s8o7Hand BAB’s9o7H
exhibitasixthorderdescentoferrorforlargeτ. H= kpk2 −GmSmM ≈ kpk2 −GmSmM (19)
2M kqk 2mM kqk
5
where M isthereducedmass, G isthe gravitational constant,
Ruth [ABA]
and m and m are the masses of the Sun and Mercury re-
S M s5odr4 [ABA]
spectively. The planar Sun-Mercurysystemis hereevaluated −4.5
ABA104 [ABA]
usingdata fromnssdc.gsfc.nasa.gov wherethe aphelionischo-
ABA864 [ABA]
senastheinitialSun-Mercurydistanceyieldingtheminimum or −5 ABA1064 [ABA]
orbitalvelocity. rr
e Yosh s7o6 A [ABA]
n
a−5.5
Figure 5 shows a benchmark profile of the 10 integrators ni
o
for the Sun-Mercury system. It is noticable that ABA104 milt −6
shows the best performance in regardto the maximum error
a
and BAB’s9o7Hshows the bestperformanceinregardto the H
m
averageerror. −6.5
u
m
xi
The accumulated error in a two body system (Sun-Mercury) ma −7
ABAs5o6H A [ABA]
manifestsasasecularadvanceoftheorbitalellipsesimilarto
BABs7o7H [BAB]
the precession of perihelion induced by general relativity or
−7.5
BAB’s8o7H [BAB]
the gravitational pull of other planets. The mathematical na-
BAB’s9o7H [BAB]
ture of the secular SI error is described by [10]. The secular
advance of the orbital ellipse can be evaluated through the −8
0.4 0.6 0.8 1 1.2 1.4
rotationoftheLaplace-Runge-LenzvectorA:
time per evaluation

A=p×L−rˆ (20)
−5.5
whererˆ =r/r, L= r×p,pisthemomentumvectorandris
−6
thepositionvectorwithlengthr. r
o
r
r
e
n −6.5
a
ABAs5o6H A [ABA] ni
o
104 BABs7o7H [BAB] milt −7
BAB’s8o7H [BAB] a
H
ry) BAB’s9o7H [BAB] an −7.5
u 3 e
nt10 m
e −8
c
n
a
uli102 −8.5
J
c/
e
s −9
c
ar101 0.4 0.6 0.8 1 1.2 1.4
dt| ( Rs5uothd r[4A B[AAB]A] time per evaluation
θ/
d ABA104 [ABA] Figure5:Integrator benchmark using the planar Sun-Mercury system in
|100 ABA864 [ABA] thetimeintervalt ∈ [0,10yr]. Top: MaximumHamiltonian er-
ror. Bottom: AveragerelativeHamiltonianerror. Thehorizontal
ABA1064 [ABA]
axescorrespondtoτ/sinunitsofhours. Allaxesareinbase10
Yosh s7o6 A [ABA] logarithmicunits.Legendappliestobothgraphs
10−1
2
10 Theorbitalprecessionadvanceslinearlywithtimet[10]. The
orbital periods/(τ/s) rotation θ is introduced as the angle of A implying that θ(t)
is a linear function and is chosen such that θ(0) = 0. Fig-
Figure4:Advanceofperihelion(absoluterotationvelocity)asafunctionof ure 5 shows the perihelion precession dθ/dt which has been
thenumberofevaluationsperorbit. Linearfit95%confidencein-
obtained numerically through linear regression of θ(t) over
tervalsoftheperihelionadvancearesignificantlysmallerthanthe
dotsize. TheintegratorsABAs5o6H,ABA104andBAB’s8o7H numerousorbitalperiods. ComparisonofFigure4and5sug-
achieve approximately the same performance, obfuscating their gests that the accumulated error is in better agreement with
plots theaverageerrorthanthemaximumerror.
6
0.2 0.4 0.6 0.8 1 1.2 1.4 0.05 0.1 0.15 0.2 0.25 0.3

(A) (C)
1 0.4
0.3
0.8
0.2
0.1
0.6
0
p p
0.4 −0.1
−0.2
0.2
−0.3
0 −0.4
−0.5
−0.2
(B) (D)
1 0.4
0.3
0.8
0.2
0.1
0.6
0
p p
0.4 −0.1
−0.2
0.2
−0.3
0 −0.4
−0.5
−0.2
−0.2 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1
q q
Figure6:PhasespaceillustrationofthesubstepsperformedovertwoiterationsbyYoshs7o6A(A)(C)andBABs7o7H(B)(D).Theintegratorsareappliedto
theSHO(A)(B)withinitialconditions(q,p)=(0,1)andtimestepτ=π/4,andtheHénon-Heilessystem(C)(D)inthey-plane(qx= px =0)
withinitialconditions (q,p) =(qy,py) = (0.4,0.4)andtimestepτ =2. BluecirclesandXsindicatetheinitialandfinalpositionsrespectively
foreachiteration. Bluedotsindicatesubstepswithpositivec andd coefficientforcalculatingthenextsubstep. Reddotsandgreendotsindicate
i i
substepswithnegative c andd coefficientsrespectivelyforcalculatingthenextsubstep. Agreendotwithinareddotindicatesthatbothc and
i i i
d arenegative. ThecontourcolorsdisplaytheHamiltonianwithvaluesgivenbythecolorbars. (C)showsanunstableintegrationwhereHisnot
i
conserved.(A)(B)and(D)showstableintegrationswhichpreserveHindefinitely
iscussion
V. D evitablycausehigherrorsultimatelywaivingsymplecticprop-
ertiesofSIs. Higherorderintegratorsaregenerallymorecon-
Athoroughbenchmarkanalysisof10differentintegrators strained and usually involve several negative splitting coeffi-
has been performed for the SHO, the Hénon Heiles system cientswithlargeabsolutevalues. Theperformancebreakpoint
and the Keplerproblem. All numerical tests suggestthat the of higher order integrators can be visualized by considering
obtainedintegratorsinspiteofanincreasednumberofstages the executed substeps. Figure 6 illustrates the substeps per-
greatlyandconsistentlyimprovetheintegrationaccuracyand formedbytwo sevenstageSIs (Yoshs7o6Aand BABs7o7H)
stability. Infact, thenewschemesoutperformstheacclaimed for the SHO and Hénon-Heiles system. Figure6 reveals that
stage-optimized integrator Ruth by several orders of magni- Yosh s7o6 A (A) (C) performs large forward and backward
tude. The higher order behavior for SHO-perturbed systems stepscomparedtoBABs7o7H(B)(D).Thestepsaretoolarge
is an additional perk of the obtained integrators. This prop- in (C) for the Hamiltonian to be preserved. BABs7o7H and
ertymay be particularlyusefulin fieldsofmoleculardynam- othersplittingcoefficientsobtainedinthispapergenerallyper-
ics where symplectic correctors may not be suitable. Opti- formsmallnear-forwardsteps.Notethatsomeofthemethods
mized fourth order schemes appear to be useful for achiev- presentedinTable4onlyhavesmallcoefficientsforeitherdor
ing high numerical stability allowing large time steps with c. Forthisreasontheintegratorperformancemaychangedras-
minimal integration error. High stability is particularly im- ticallyiftherolesofTandVareswitched. Theperformanceof
portant for systems where the complexity of H varies with theobtainedintegratorsgenerallydeclinesforsystemswhere
time: CloseencountersinthegravitationalN-bodysystemin- T(p)deviatessignificantlyfromasquaredkineticenergy.
7
onclusion
VI. C
It has been found that harmonic order conditions higher
Efficient fourth order splitting coefficients have been ob- thansevenoftenintroducemorenegativesplittingcoefficients.
tained by satisfying higher order conditions for the simple A largernumber of stages presentsmoredegreesof freedom
harmonicoscillator. Numericaltestsindicatethattheobtained resulting in more inappropriate solutions. Additional con-
integratorsoutperformthethreestagefourthorderintegrator straints should be pursued to obtain efficient splitting coef-
by R. Ruth as well as the more optimized integrators by W. ficientswithmorethanninestages.
Kahan and S. Blanes. The obtained integrators also compare The presented decomposition can be used to satisfy gen-
favorably with higher order integrators by H. Yoshida for eral fourth order criteria and higher order conditions specif-
moderatetimesteps. Theintegratorsareparticularlysuitable ically for the simple harmonic oscillator. This work suggests
for particle systemsthat requirehighnumerical stability and thatgeneralpurposehighaccuracy,highstabilityexplicitsym-
accuracy at large time steps. The best performance is gen- plectic integrators can be pursued by approaching multiple
erally achieved by the methods BAB’s9o7H, BAB’s8o7H and higherorderconditions. Thishasbeenachievedherebysatis-
ABAs5o6HA. fyingindividualhighorderentries.
eferences
R [6] Optimized fifth ordersymplecticintegrators for orbital
problems. Tselios, Kostas and Simos, T.E. 1, 2013, Re-
[1] Fourth-ordersymplecticintegration.Forest,Etienneand vistaMexicanadeAstronomíayAstrofísica,Vol.49,pp.
Ruth,RonaldD.1,1990,PhysicsD,Vol.43,pp.105-117. 11-24.
[2] Construction of higher order symplectic integrators.
[7] ImplementingtheNelder-Meadsimplexalgorithmwith
Yoshida, Haruo. 5-7, 1990, Physics Letters A, Vol. 150,
adaptive parameters. Gao, Fuchang and Han, Lixing.
pp.262-268.
1, 2012, Computational Optimization and Applications,
[3] Forward Symplectic Integrators for Solving Gravita- Vol.51,pp.259-277.
tional Few-Body Problems. Chin, Siu A. and Chen, C.
R.3-4,2005,CelestialMechanicsandDynamicalAstron- [8] Highprecisionsymplecticintegrators forthe Solar Sys-
omy,Vol.91,pp.301-322. tem. Farrés,Ariadna, etal.2, 2013, CelestialMechanics
andDynamicalAstronomy,Vol.116,pp.141-174.
[4] Wisdom, J., Holman, M. and Touma, J. Symplectic
Correctors. [book auth.] Jerrold E. Marsden, George
[9] Composition constants for raising the orders of uncon-
W. Patrick and William F. Shadwick. Integration Algo-
ventional schemes for ordinary differential equations.
rithms and Classical Mechanics. s.l. : Providence, RI,
Kahan,WilliamandLi,Ren-Cang.219,1997,Mathemat-
1996,Vol.10,pp.217-244.
icsofComputation,Vol.66,pp.1089-1099.
[5] Newfamiliesofsymplecticsplittingmethodsfornumer-
icalintegrationindynamicalastronomy.Blanes,Sergio, [10] Physics of symplectic integrators: Perihelion advances
et al. June 2013, Applied Numerical Mathematics, Vol. and symplectic corrector algorithms. Chin, Siu A. 3,
68,pp.58-72. 2007,PhysicsreviewE,Vol.75.
8
Table4:Obtainedsplittingcoefficientswith77digits.Themethodsareorganizedaccordingtotheirstagecount.
UniqueABAs5o6H,A,BandC
d(1) = 0.1558593591762168313166117535752091422239663993391011462498104831549442591694
d(2) =-0.0070254990919573173514483364758218294773716640092220571342056284758867609611
c(1) =-0.6859195549562166768601873150414759494319985863677163820719179393682014399373
c(2) = 0.9966295909529363159571451429325843698583459772292551181721475637244006507927
d(1) = 0.4020196038964999834667409950496227775945673320979099323902806525851620445492
d(2) = 0.5329396856308538150258772262086702929451721575835842834460326556965220312130
c(1) = 0.9110842375676615218574607388486783304139753525628699898390474132061253024968
c(2) = 0.1740059542332660799009374186088931171982348451547482386207462271424421679090
d(1) = 0.1868565631155112597511173758337610451623768791420295598869906080256347098408
d(2) = 0.5520581660514781484261043096825685955052553493857487316732455515112095793516
c(1) = 0.5642486163110637621453746447826190031465518453448216443978248524452914255263
c(2) =-0.2393627021773294286793711975145735718917010075899623225091609656425715483488
d(3)=0.5-d(1)-d(2);d(4)=d(3);d(5)=d(2);d(6)=d(1);
c(3)=1-2*(c(1)+c(2));c(4)=c(2);c(5)=c(1);
UniqueBABs6o7H(moreexist),BABs6o5HandBAB’s6o5H
d(1) = 0.0832701092493097690276300822599156817795619881080575430174826369500044839553
d(2) = 0.3997273690963360211284395920007795550575060531634793748020207288976344439468
d(3) =-0.0541842778124726964199287659702152862181671805554302053695494244408226729818
c(1) = 0.2475471587650765967910125296669232190787926795528258860075742877866898482465
c(2) = 0.5446579217808193419580029125986805136192611468678745304306198457355253495088
d(1) = 0.0658831533161155021794371297629949214211270641450367882108652454225238292357
d(2) =-0.6711629060948253965117521242801468651670183829736696004743743104248379033034
d(3) = 0.9736703100725350498414312651550857191131218932308788320806762063004096320895
c(1) = 0.2265023974336291596186923088995152371194987043433278784212774210853229477086
c(2) =-0.0047799986678794678665602622568725658855054645768977416258774175978957628102
d(1) = 0.0650508268637574949487516678539036744380576078003520231297474895927182047842
d(2) =-0.3948051939117155639582651907195511796839512131373933326629623631591978696178
d(3) = 0.6918498547904058960782554213200966000604457266088855079657967408955821538731
c(1) = 0.2328962665845291347812910553597276545034489573682700034501734659308763580659
c(2) =-0.0111617638003721094728940473306267483522869816097800738075975236192041340802
d(4)=1-2*(d(1)+d(2)+d(3));d(5)=d(3);d(6)=d(2);d(7)=d(1);
c(3)=0.5-c(1)-c(2);c(4)=c(3);c(5)=c(2);c(6)=c(1);
BABs7o7HandBAB’s7o6H
d(1) = 0.0638745574250616045658401356462756092272737349204789877616691621039130680037
d(2) =-0.0650239777505938311516598494765811300128929849501107531440553736739769298894
d(3) = 0.2509446105745547370613575645855473357282136355718617210088090794709222342775
c(1) = 0.2752781729059777393394978710448690782125215018949186085075325605348526197756
c(2) =-0.0843138705589167473554015820986490036832890668438279781819362930106920807542
c(3) = 0.1674497222006475614401177016323447087805836086414469568091358611098423440220
d(1) = 0.0522155297747848201407012160969040693245471580104797248381281194964273517726
d(2) =-0.0824972558529561412131911937717420514162728339681056503508469680313691406287
d(3) = 0.3285541797987193353601113204079269672646845923663727576986276602617960257026
c(1) = 0.2487563308365098625528031803769571289196558939258433219240690943143969198069
c(2) =-0.0651011247076581799932061212576878177123945470202647856511921757791017775052
c(3) = 0.2480624780675545152650672751613106579864581926645260137078906816928505888862
d(4)=0.5-d(1)-d(2)-d(3);d(5)=d(4);d(6)=d(3);d(7)=d(2);d(8)=d(1);
c(4)=1-2*(c(1)+c(2)+c(3));c(5)=c(3);c(6)=c(2);c(7)=c(1);
BAB’s8o7H
d(1) = 0.0538184115480034769403763798524605188562842390760879592632218376015166638395
d(2) = 0.1648743326910472361014809085317059425299121141031052090901977952513984878990
d(3) = 0.3895399407808198068744134256203146340834631254960864069726823667050364522355
d(4) =-0.2288957415563594299572505173565338312542463595622272333825061768110435645417
c(1) = 0.1486140577445185629163082471176700173109512976367237631150576219945233462284
c(2) = 0.1071986675806227950500566279939336794589433458464489776124879870581484936262
c(3) =-0.0149646736494517061945681450558142918831874360003431672632178480031079700216
d(5)=1-2*(d(1)+d(2)+d(3)+d(4));d(6)=d(4);d(7)=d(3);d(8)=d(2);d(9)=d(1);
c(4)=0.5-c(1)-c(2)-c(3);c(5)=c(4);c(6)=c(3);c(7)=c(2);c(8)=c(1);
BAB’s9o7H
d(1) = 0.0464929004396589154281717058427105561306160230440930588914036807441235817244
d(2) = 0.1549010127028879927850680477816652638346460615901974901213193690401204696252
d(3) = 0.3197054828735917137611074311771339117602994884245091220333340037841616085048
d(4) =-0.1929200088157132136865513532391282410293753210475133631464188500663304857888
c(1) = 0.1289555065927298176557065467802633438775379080212831185779306825670371511433
c(2) = 0.1090764298548827040268039227200943338187149719339317536310302288046641781422
c(3) =-0.0138860356804715144111581981849964201100030653749527555344377031679795959892
c(4) = 0.1837549745641803566768357217228586277331494085368674804908537743649129597425
d(5)=0.5-d(1)-d(2)-d(3)-d(4);d(6)=d(5);d(7)=d(4);d(8)=d(3);d(9)=d(2);d(10)=d(1);
c(5)=1-2*(c(1)+c(2)+c(3)+c(4));c(6)=c(4);c(7)=c(3);c(8)=c(2);c(9)=c(1);
9