High order explicit Runge-Kutta pairs for
ephemerides of the Solar System and the Moon
PHILIPW.SHARP y
DepartmentofMathematics,UniversityofAuckland,PrivateBag92019,
Auckland,NEWZEALAND.
Abstract. NumericallyintegratedephemeridesoftheSolarSystemandtheMoonrequirevery
accurate integrations of systemsof secondorder ordinarydierential equations. Wepresent a
newfamilyof8-9explicitRunge-Kuttapairsandassesstheperformanceoftwonew8-9pairson
theequationsusedtocreatetheephemerisDE102. Partofthisworkistheintroductionofthese
equationsasatestproblemforintegratorsofinitialvalueordinarydierentialequations.
1. Introduction
An ephemeris of the planetsand the Moon consists of tabularinformation from
which accurate positions and velocities of the celestial bodies can be calculated
for any value of astronomical time on a prescribed range. Modern ephemerides
typically contain accurate values of the position and velocity at equally spaced
astronomicaltimes,andthecoeÆcientsofChebyshevpolynomialsforinterpolation
betweenthevalues.
The information in an ephemeris can beobtained by numericallyintegrating a
system of ordinary dierential equations that model all signicant gravitational
attractionsbetweenthebodies. Totakefull advantageoftheaccuracyof modern
astronomicalobservationsandtodistinguishbetweencompetinganalyticaltheories
forthemotionofthebodies,theglobalerrorintheintegrationsmustbeverysmall.
Anothercharacteristicoftheintegrationsisthattheyoftenspanalargeintervalof
astronomicaltime,necessitating manyintegrationsteps.
The accumulated round-o error in an integration will in general grow as an
integrationproceeds. Iftheintegrationisdonein doubleprecisionarithmetic,the
accumulated round-o error may be far larger than the required accuracy. This
diÆcultycan beovercomebyusing80-bitarithmeticor evenquadrupleprecision.
Theordinarydierentialequationsforephemeridesarenon-stiandhenceexplicit
Runge-Kutta (ERK) pairs are suitable methods for performing the integrations.
Pairsconsistofformulaeoforderspandq,whereq<pandistypicallyp 1. The
computational eort requiredto advance astep with a pair can bemeasured by
thenumberofderivativeevaluations,knownasstages,performedonthestep. For
conciseness,werefertoapairofsstagesasans-stageq ppair.
y
ThisworkwassupportedbytheUniversityofAucklandResearchCommittee.
Ofthe many ERKpairsavailable,the13-stage7-8pairofPrinceand Dormand
[6]hasproventobeaseÆcientasanyotheronmanyproblemswhenusingdouble
precisionarithmetic, exceptpossiblyforlowaccuracyrequirements.
Inparticular,thepairisnoticeablymoreeÆcientthan8-9pairs. Weinvestigate
whetherthisresultholdsfornumericallyintegratedephemerides. Insectiontwo,we
summarisethederivationoftwofamiliesof8-9pairs,one ofwhichisanewfamily,
andpresentanearoptimalpairfromeachfamily. Insectionthree,wecomparethe
performanceofthetwonewpairsandthe7-8pairof PrinceandDormandonthe
modelequationsofDE102. Weendin sectionfour withadiscussionofourwork.
2. Order ninepairs
2.1. Denitions
Considertheinitialvalueproblem
y 0
=f(x;y); y(x
0 )=y
0
; (1)
where 0
d=dx,f :RR n
!R n
andthesolutiony(x)issuÆcientlydierentiable.
The8-9 ERKpairsweinvestigate haves-stagesand generatean ordernine ap-
proximationy
i
andanordereightapproximationyb
i toy(x
i
),i=1;2;:::;according
to
y
i
= y
i 1 +h
s
X
j=1 b
j f
j
; (2)
b y
i
= y
i 1 +h
s
X
j=1 b
b
j f
j
; (3)
whereh=x
i x
i 1 and
f
j
=f(x
i 1 +hc
j
;y
i 1 +h
j 1
X
k =1 a
jk f
k
); j=1;:::;s; (c
1
=0):
Werefertoc
j
;j =1;:::;s,astheabscissae,a
ij
;j =1;:::;i 1;i=2;:::;s,asthe
interiorweights,b
j ,
b
b
j
, j=1;:::;s,astheexteriorweights,and totheabscissae,
theinteriorweightsand theexteriorweightscollectively as thecoeÆcientsof the
pair. Toensure theone stepnature of the pairs, we restrictthe abscissaeto the
interval[0;1].
WhenthecoeÆcientsofthepairarechosensothat y
i andyb
i
areordernine and
eight respectively, somecoeÆcients are available asfree parameters, leadingto a
familyofpairsandnotauniquepair. Individualpairsfromthisfamilyareobtained
by assigningvalues to thefree parameters. Since weare interestedin doingvery
accurateintegrations,wehavechosenthevaluessothat thelocal errorintroduced
Thiserrorforthestepfrom x
i 1 to x
i
canbewrittenas(for,example[3],page
151)
h 10
X
t2T
10 e
10 (t)D
10
(t)+O(h 11
);
where T
10
is theset ofrootedtrees oforder ten,e
10
(t)isthe principalerrorcoef-
cient fortree t,and D
10
(t) is the elementary dierentialfor t. The elementary
dierentialisformedfromthepartialderivativesoff withrespecttox andy and
evaluatedat(x
i 1
;y
i 1 ).
TheprincipalerrorcoeÆcientfortreetcanbewritten as
e
10 (t)=
(t)
10!
((t) s
X
k =1 b
k
k
(t) 1);
where (t) and (t) are positive integers, and
k
(t) is a function of the interior
weightsandabscissae.
Numerical experiments have shown (see, for example [6]) that decreasing the
size of theprincipal errorcoeÆcientswill in generalimprovetheeÆciency of the
method. Hence, we choose the free parameters so that the error coeÆcients are
closetotheirminimumvalue.
WeusetwomeasuresofthesizeoftheprincipalerrorcoeÆcients
E 2
10
=
"
X
t2T
10 e
2
10 (t)
#
1=2
; E 1
10
= max
t2T10 fje
10
(t)jg: (4)
2.2. Sixteen stages
Verner [9] deriveda familyof 16-stage8-9 pairs with c
2 , c
5 , c
9 , c
10 , c
11 , c
13 , c
14
and a
11;6
asfreeparameters(Tosimplify what follows,wehaveinterchangedthe
coeÆcientsforthefourteenandsixteenthstages,thiscanbedonewithoutchanging
thepropertiesofthepairs.) Theordernineformulainthepairsusestherstfteen
stages and the order eight formulauses all sixteen stages. ThecoeÆcients b
j
; b
b
j ,
j=2;:::;7,b
16 ,
b
b
14 and
b
b
15
areidenticallyzero.
VernerpresentedthecoeÆcientsofapairfromthisfamilywhichhadc
2
=1=12,
c
5
= (2+2 p
6)=15, c
9
= 1=2, c
10
= 1=3, c
11
= 1=4, c
13
= 5=6, c
14
= 1=6 and
a
11;6
=0. ThepairhasE 2
10
=6:110 5
andE 1
10
=3:110 5
,andhasbeenused
whencomparingthenumericalperformanceof8-9pairswithpairsofotherorders.
However,the pairwas intended asanillustration ofthederivationand notasan
optimalornearoptimalpair.
Toassessinaproblem-independentwayifthe8-9familyofVernercontainsmore
eÆcientpairs,andifso,howmuchmoreeÆcient,weperformedaminimisationof
E 2
10
overthefreeparameters,subjectto theconstraintthat thecoeÆcientsofthe
selectingapairfromafamilyandisintendedtopreventtheselectionofapairwith
poorround-oerrorproperties. AlthoughnoonevalueofM isused, itisoften20
or30andwechose30.
Weperformedtheminimisation usingan interactivegridsearchand obtaineda
minimumvalueforE 2
10
of7:510 7
whenworkingwithagridsizeof0:001. The
pairweobtainedhadc
2
=0:020,c
5
=0:311,c
9
=0:312,c
10
=0:105,c
11
=0:587,
c
13
=0:879,c
14
=0:916 and a
11;10
= 0:150 (as amatter of preference wehave
useda
11;10
inplaceofa
11;6
asafreeparameter). Thealgorithmsin[9]canbeused
tondtheremainingcoeÆcients. ThepairhasE 1
10
=2:810 7
.
AslightlysmallervalueofE 2
10
ispossibleifasmallergridsizeisused, butsince
thenumberofderivativeevaluationsvariesapproximatelyastheninthrootofE 2
10 ,
thegainin eÆciencyis small. AsignicantlysmallervalueofE 2
10
, approximately
twice assmall,ispossibleiftheabscissaearenotconstrainedtotheinterval[0;1],
butthis choicemeansthepairisnolongeraonestepmethod.
An estimate of the eÆciency of the new pair relativeto that of Verner can be
calculatedbyusing E 2
10
forthe twopairs. Todothis, weassume theglobalerror
for a xed t and stepsize is proportional to E 2
10
. The relative eÆciency is then
estimatedas
3:110 5
7:510 7
1=9
=1:63::::
Thissuggeststhenewpairwillbeapproximately63percentmoreeÆcientthanthe
pairofVerneratsmallstepsizes,raisingthepossibilityofitbeingcompetitivewith
pairsofotherorders.
2.3. Seventeen stages
TheworkofSharpandSmart[7]for4-5and5-6ERKpairsshowsagainineÆciency
ispossibleifanextrastageisused toform thepair. Theextrastage meansmore
freeparametersareavailable,permittingasmallervalueofE 2
10
,but thisis atthe
expenseofincreasingbyonethenumberoffunction evaluationsrequiredtotakea
step.
ToinvestigateifagainineÆciencywaspossiblefor8-9pairs,wederivedafamily
of 17-stage 8-9 pairs. The family has six more free parameters (three abscissae,
threeinteriorweights)thanthe16-stage8-9pairs.
Weimpose,asiscommonlydoneforhighorderERKpairs(see, forexample[9],
[11]),thefollowingconditionsonthecoeÆcientsofthepair
c k +1
i
k+1
= i 1
X
j=1 a
ij c
k
j
; k=0;:::;
i
1; i=1;:::;s; (5)
a
ij
= 0; if
i
>
j
+1; j=1;:::;i 1; i=1;:::;s: (6)
Theimposition of theseconditionsreduces thenumberof independentorder con-
Theconditionscanberepresentedbythestage-ordervector=[
1
;
2
;:::;
s 1 ]
T
.
The16-stagepairshave =[5;1;2;3;3;4;4;5;5;5;5;5;5;5;5]
T
;toobtain forthe
17-stagepairs, onepositiveintegerlessthanve mustbeinserted. Weexamined
threechoices andfoundthat insertinga4after thesecond 4togive
=[5;1;2;3;3;4;4;4;5;5;5;5;5;5;5;5]
T
ledtothelargestnumberoffreeparameters.
With specied,thederivationissimilartothatforthe16-stagepairs,themain
dierencebeingfewerconstraintsontheabscissaefortherstninestages. Wetook
c
2 ,c
5 ,c
6 ,c
7 ,c
9 ,c
10 ,c
11 ,c
12 ,c
14 ,c
15 ,a
8;7 ,a
11;10 ,a
12;10 ,a
12;11
asfreeparameters;
other choices are possible, but the number of free parametersremains the same.
Theabscissaec
3 ,c
4 ,c
8 ,c
16 andc
17
areconstrainedas
c
3
= 2
3 c
4
; c
4
= 3c
6 4c
5
4c
6 6c
5 c
6
;
c
8
=c
9 20c
6 c
7 15c
6 c
9 15c
7 c
9 +12c
2
9
5(3c 2
9 4c
6 c
9 +6c
6 c
7 4c
7 c
9 )
; c
16
=c
17
=1:
(7)
Theexpressionforc
13
isthesameasforc
12
inthe16-stagepairsexceptc
8 ,c
9 ,c
10
andc
11
arereplacedbyc
9 ,c
10 ,c
11 andc
12
respectively.
WeperformedaminimisationofE 2
10
forthenewfamilyusinganinteractivegrid
search and steepest descent (a grid search by itself wasimpracticable because of
thelarge numberof free parameters)and obtained apairwith E 2
10
= 1:010 7
andE 1
10
=3:610 8
. Thevalueof thefreeparameterstofourdecimalplacesare
c
2
=0:0757,c
5
=0:3617,c
6
=0:4139,c
7
=0:1074,c
9
=0:7607,c
10
=0:6068,c
11
=
0:1531,c
12
=0:8333,c
14
=0:9733,c
15
=0:9888,a
8;7
= 0:0001,a
11;10
= 0:0078,
a
12;10
=0:0067and a
12;11
= 0:0026. Equations(7) together with given above
andthealgorithmsin[9]canbeusedtondtheremainingcoeÆcients.
Inasimilarwaytothat forthetwo16-stagepairs, E 2
10
canbeused toestimate
therelativeeÆciencyofthenew16-stageand17-stagepairs. Weget
7:510 7
1:010 7
1=9
16
17
=1:18:::;
where the factor(16/17) is the ratioof the numberstages. Hence weexpect the
new17-stagepairtobeabout18percentmoreeÆcientthanthenew16-stagepair
forsmallstepsizes.
2.4. Generalised
Thefamiliesof8-9pairsdescribedintheprevioussub-sectionarereadilygeneralised
toinclude eitheroneortwoextrafreeparameters.
Onegeneralisationistoreplace b
b
j
,j=1;:::;sbytheconvexlinearcombination
b +(1 )
b
b . This substitution is equivalent to making one of thepreviously
identicallyzero b
bafreeparameter. Thelocalerrorestimateforthepairischanged,
but since b
j
, j =1;:::;s 1 remain thesame, the principal errorcoeÆcients of
theordernineformulaeandhenceE 2
10
(andE 1
10
)areunchanged.
Thesecond generalisation isbased onatransformationobtainedby Verner [10]
fortwofamilies of8-stage5-6ERKpairs. VernershowedthefamilyofPrinceand
Dormand[6]whichhasc
2 ,c
3 , c
5 ,c
6 ,b
8 and
b
b
7
asfreeparameterscanbeobtained
fromthefamilyofVerner[9]whichhasc
2 ,c
3 ,c
5 andc
6
asfreeparametersusinga
simpletransformationonthelasttworowsofinteriorweights.
Thistransformationgeneralises(Verner,privatecommunication)tootherfamilies
ofpairs, includingthe8-9pairsin thispaper. This meansb
16 and
b
b
15
,previously
zeroin the16-stage8-9pairs, andb
17 and
b
b
16
,previouslyzeroin the17-stage8-9
pairs,canbefreeparameters.
The introduction of these twofree parameters changes the local error estimate
and the principal error coeÆcientsof the order nine formula. However, asis the
casefor the 5-6pairs in [9], the change in E 2
10 and E
1
10
is smallfor near optimal
pairs.
3. DE102
Newhall, Standishand Williams [5] presented DE102, an ephemeris of the Solar
SystemandtheMoon,obtainedbyintegratingasystemof33secondorderordinary
dierentialequationsoftheform
y 00
=f(t;y;y 0
): (8)
The system (8) consists of equations of motion for the nine planets, the Moon
and three equations for the lunar physical librations. The motion of the Sun is
found from thedenition ofthe Solar System barycentre. Theequations contain
contributionsfrom point-massinteractions,gureeectsfor Earthand theMoon,
Earth tides and perturbations from the ve asteroids (1) Ceres, (2) Pallas, (4)
Vesta,(7)Irisand(324)Bamberga.
Thecalculationsrequired foroneevaluation of thesecond derivativefor(8) are
described in Figure 1. A fuller descriptionis given in [5] and byinference in the
programDE118i.ARCofMoshier,availableontheinternet.
ThemodelequationsusedinDE102canbegeneralisedinanumberofways. For
example,termsmodellingthedeformationoftheMoon'ssurfacebytheEarthand
perturbationsfromotherasteroidscanbeincluded. However,themodelequations
of DE102 have proven suÆciently accurate and renements to DE102 (see, for
example [8]) have been in the coordinate systems used, and in the observations
usedtodene theinitial conditionsandphysicalconstants.
1. Initialise
a)Calculatetheheliocentricpositionandvelocityfortheasteroidsandtransformtoap-
proximatebarycentricvalues. Thesevaluesarecorrectedoncethecorrectpositionofthe
Sunisknown.
b)Calculatethedistancebetweenthebodies. ThedistancesinvolvingtheSunorasteroids
areestimates only. These distancesarecorrected oncethe correctpositionofthe Sunis
known.
c)Usexed-pointiterationtondthecorrectpositionandvelocityoftheSunandaster-
oids,thencorrectthedistancesinvolvingtheSunorasteroids.
d)Calculatethecubeofthedistancesbetweenallbodies.
2. Point-massacceleration
a)CalculatetheNewtonianaccelerationofallbodies.
b)Calculatethepost-NewtonianaccelerationoftheplanetsandtheMoon.
3. Figure ofthe Moon
a)Formtherotationmatrixforthetransformationfromspacetobodycoordinates.
b) Calculatethe eectof the point-mass Earthonthe lunargureand add this to the
lunaracceleration.
c) Calculate the torque on the Moon due to the gravitational interaction between the
lunargureandtheexternalpoint-massEarth.
d)Theaccelerationfromb)inducesanaccelerationintheEarth-addthistotheEarth's
acceleration.
e)Calculatetheeectofthepoint-massSunonthelunargureandaddthistothelunar
acceleration.
f)CalculatethetorqueontheMoonduetothegravitationalinteractionbetweenthelunar
gureandtheexternalpoint-massSun.
g)Calculatetheaccelerationofthelibrationangles.
4. Figure ofthe Earth
a)Calculatetheeectof thepoint-massMoonontheEarth'sgureandaddthis tothe
Earth'sacceleration.
b)Theaccelerationfrom a)inducesanacceleration inthe Moon-addthis tothe lunar
acceleration.
c)Calculatethecontributiontothe acceleration ofthe Moonand theEarthdueto the
Earthtides.
d) Calculatethe eectof thepoint-mass Sun onthe Earth'sgure andadd this to the
Earth'sacceleration.
Theaccelerationsinthissectionareadjustedfortheprecessionandnutationoftheequinox
andobliquityoftheecliptic.
Figure1. Asummaryofthecalculationsrequiredforoneevaluationofthesecondderivativein
themathematicalmodelofDE102.
4. Numericalexperiments
Weconducted numericaltestsofthetwonew8-9pairsand the7-8pairofPrince
andDormandonthemodelequationsdescribedintheprevioussection. Theresults
areillustratedbelow. Thepairsare denotedbyPD78(Princeand Dormand7-8),
P16(new16-stage)andP17(new17-stage).
A computerwhich performed quadruple precision in hardware was unavailable
and henceweusedthethe multiprecisionFortran90packageMPFUN90 ofBailey
[1], withtheprecisionlevelset at35 digits,approximatelythat of quadruplepre-
cision. Themultiprecision version ofour programwas270times slowerthanour
double precisionversion which makes the use of MPFUN90 impractical for long
integrations.
−24 −22 −20 −18 −16 −14
2.6 2.7 2.8 2.9 3 3.1 3.2 3.3 3.4 3.5 3.6
Problem DE102: integration interval is [0,20]
Base 10 logarithm of the norm of the end−point global error
Base 10 logarithm of the number of derivative evaluations
−−− 13−stage ... 16−stage
− 17−stage
Figure2. Alog-loggraphofthe numberof derivativeevaluationsagainstthe normof theend-
pointglobalerrorfor DE102witha integrationintervalof20. Prince andDormand 7-8pair-
dashedline,new16-stagepair-dottedline,new17-stagepair-solidline.
ThecoeÆcients ofthe7-8 pairasspeciedin [6] areaccurate to approximately
18digits. Werecalculatedthe coeÆcientsin100digitarithmetic usingthevalues
ofthefreeparametersin [6]and usedthesecoeÆcients,roundedto35digits. The
global errorin anumerical solutionwas obtainedby calculatinga moreaccurate
solutionandtakingthedierencebetweenthetwosolutions.
The rst exampleis foran integration interval of 20 and local error tolerances
of 10 i
, i = 14;:::; 22. Figure 2contains the log-loggraph of the number of
derivative evaluations against the norm of the end-point global error (the data
pointshavebeenjoinedforclarity).
PairP17ismoreeÆcientthanP16suggestingtheeÆciencyisimprovedbyadding
astage.ThegainineÆciencyvariesfrom15to20percent,ingoodagreementwith
that predicted using E 2
10
. The pairs P16 and P17 are more eÆcient than PD8
forglobal errorssmallerthan(approximately) 10 16
,and 10 18:5
respectively. In
additionandascanbeexpectedfromtheorderofthepairs,theeÆciencyofthe8-9
pairsrelativetothe7-8 pairincreasesasthe globalerrordecreases. Forexample,
P17 is 16 percent moreeÆcient for a global error of 10 20
and 29 percent more
eÆcientforaglobalerrorof10 22
.
Thesecondexampleisforanintegrationintervalof50usingthesamelocalerror
tolerances as in the rst example. The results are given in Figure 3. P16 was
−24 3 −22 −20 −18 −16 −14 −12
3.1 3.2 3.3 3.4 3.5 3.6 3.7 3.8 3.9 4
Problem DE102: integration interval [0,50]
Base 10 logarithm of the norm of the end−point global error
Base 10 logarithm of the number of derivative evaluations
−−− 13−stage
− 17−stage
Figure3. Alog-loggraphofthe numberof derivativeevaluationsagainstthe normof theend-
pointglobalerrorfor DE102witha integrationintervalof50. Prince andDormand 7-8pair-
dashedline,new17-stagepair-solidline.
excluded because ourtest resultssuch as those in Figure 2 showed P16 wasless
eÆcientthanP17forthelocalerrortoleranceswewereusing.
TheeÆciencyofP17relativetoPD78asafunctionoftheglobalerrorissimilarto
Theglobalerrorsare largerthan in therstexample,aresultwhich isconsistent
withalargerintervalofintegration.
5. Discussion
The main aim of our work was to investigate if 8-9 explicit Runge-Kutta pairs
were more eÆcient than lower order pairs, principally 7-8 pairs, for numerically
integratedephemerides. Wederivedanewfamilyof8-9pairs,obtainednearoptimal
8-9 pairsfrom this familyand anexisting one, and comparedthe performance of
thesepairsandthe7-8pairofPrinceandDormandonthemodelequationsofthe
ephemerisDE102.
Our testing showed the8-9 pairs were usually moreeÆcientthan the 7-8 pair.
Thegainin eÆciency wasnotlarge,but giventheamountof CPU timerequired
toproduceephemerides, thegainissignicant. Ourtestingalso showedthat near
optimal 17-stage 8-9 pairs can be more eÆcient than near optimal 16-stage 8-9
pairs.
AspartofthisworkweintroducedthemodelequationsofDE102asatestproblem
for integrators of initial value ordinary dierential equations. This problem, in
addition to being a realistic one, has several interesting numerical aspects. For
example, the position and velocity of the Sun is found by solving a system of
nonlinear(algebraic)equations. Asin[5],weusedxedpointiteration;thequestion
arisesasto whetherthereisabetterwayto solvetheequations.
References
1. D.H.Bailey,AFortran-90basedmultiprecisionSystem,RNRTechnicalReportRNR-94-013,
NASScienticComputationBranch,NASAAmesResearchCenter,January,1995.
2. E.Fehlberg,Classicalfth-,sixth-,seventh-,andeighth-orderRunge-Kuttaformulaswithstep-
sizecontrol,NASATechnicalReportNASATRR-287(1968),82pages.
3. E.Hairer,S.P.Nrsett,G.Wanner,SolvingordinarydierentialequationsI:nonstiprob-
lems,Springer-Verlag,1987.
4. S.L.Moshier,Comparison of a 7000-year lunar ephemeriswith analyticaltheory, Astron.
Astrophys.262(1982),613-616.
5. X.X.Newhall,E.M.Standish,J.G.Williams,DE102: anumericallyintegratedephemerisof
theMoonandplanetsspanningforty-fourcenturies,Astron.Astrophys.125(1983),150-167.
6. P.J. Prince and J.R. Dormand, High-order embedded Runge-Kutta formulae, J. Comput.
Appl.Math.7(1981),67-76.
7. P.W.Sharp,E.Smart,ExplicitRunge-Kuttapairswithonemorederivativeevaluationthan
theminimum,SIAMJ.Sci.Comput.14(1993),338-348.
8. E.M.Standish,X.X.Newhall,J.G.Williams,W.F.Folkner,W.F,JPLPlanetaryandLunar
Ephemerides,DE403/LE403,JPLIOM314.10-127,1995.
9. J.H. Verner, Explicit Runge-Kutta methods with estimates of the local truncation error,
SIAMJ.Num.Anal.15(1978),772-790.
10. J.H.Verner,AcontrastofsomeRunge-Kuttaformulapairs,SIAMJ.Num.Anal.27(1990),
1332-1344.
11. J.H.Verner,HighorderexplicitRunge-Kuttapairswithlowstageorder,App.Num.Math.
22(1996),345-357.
