On variational and symplectic time integrators for Hamiltonian systems

On

variational

and

symplectic

time

integrators

for

Hamiltonian

systems

E. Gagarina

a

,

∗

,

V.R. Ambati

a

,

S. Nurijanyan

a

,

J.J.W. van der Vegt

a

,

O. Bokhove

b

,

∗

a_{MathematicsofComputationalScienceGroup,Dept.ofAppliedMathematics,UniversityofTwente,P.O.Box217,7500AE,Enschede,} The Netherlands

b_{SchoolofMathematics,UniversityofLeeds,LS2 9JT,Leeds,UnitedKingdom}

a

r

t

i

c

l

e

i

n

f

o

a

b

s

t

r

a

c

t

Articlehistory:

Received18August2015

Receivedinrevisedform22October2015 Accepted22November2015

Availableonline2December2015

Keywords:

Nonlinearwaterwaves FiniteelementGalerkinmethod (Non-)autonomous variationalformulation Symplectictimeintegration

Various systems in nature have a Hamiltonian structure and therefore accurate time integrators for those systemsare of great practical use. In thispaper, afinite element methodwillbeexploredtoderivesymplectictimesteppingschemesfor(non-)autonomous systemsinasystematicway.The techniqueused isavariationaldiscontinuousGalerkin finite element method in time. This approach provides a unified framework to derive known and new symplectic time integrators. An extended analysis for the new time integratorswillbeprovided.The analysisshowsthatanovelthirdordertimeintegrator presented in this paper has excellent dispersion properties. These new time stepping schemesarenecessarytogetaccurateandstablesimulationsof(forced)waterwavesand othernon-autonomousvariationalsystems,whichweillustrateinournumericalresults. ©2015TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCC

BYlicense(http://creativecommons.org/licenses/by/4.0/).

1. Introduction

The dynamicsof various physical phenomena, such asthe movement of pendulums, planets, or waterwaves can be described inavariationalframework. Thedevelopment ofvariationalprinciplesforclassical mechanicstracesbackto Eu-ler, Lagrange, and Hamilton;an overview ofthis historycan be found in [1,19]. Thisapproach allows to expressall the dynamics of a systemin asingle functional– the Lagrangian – whichis an action integral. Hamiltonianmechanics isa reformulationofLagrangianmechanicswhichprovidesaconvenientframeworktostudythesymmetrypropertiesofa sys-tem. ThisisexpressedbyNoether’stheoremwhichestablishesthedirectconnectionbetweenthesymmetrypropertiesof Hamiltoniansystemsandconservationlaws.Whenoneapproximatesthesystemnumerically,itisadvantageoustopreserve theHamiltonianstructurealsoatthediscretelevel.GiventhatHamiltoniansystemsareabundantinnature,theirnumerical approximationisthereforeatopicofsignificantrelevance.

ThetopicoftimeintegrationofHamiltoniansystemshasseenarigorousdevelopmentduringthepastdecades.Thishas resultedinsymplectictimeintegratorsthathavebeendesignedtopreservetheHamiltonianstructureforanapproximated system.MoreinformationonsymplecticintegrationofHamiltoniansystemscanbefoundin

[12]

and

[16]

.Accordingtothe

*

Correspondingauthors.

E-mailaddresses:[email protected](E. Gagarina),[email protected](V.R. Ambati),[email protected](S. Nurijanyan),

[email protected](J.J.W. van der Vegt),[email protected](O. Bokhove).

http://dx.doi.org/10.1016/j.jcp.2015.11.049

0021-9991/©2015TheAuthors.PublishedbyElsevierInc.ThisisanopenaccessarticleundertheCCBYlicense (http://creativecommons.org/licenses/by/4.0/).

reviewin

[20]

,variationaltimeintegratorscanbeproventobesymplecticbyconstructionundercertainconditions.Inthe variationalapproach,onecanderivewell-knownschemesornewschemes,asreportedin

[10,17,22,26]

.

One of the mostactively developed numerical methods is the finite element method withtime discontinuous basis functions. Such a discontinuous Galerkin methodin time was studied in [5–7]. In particular, the variationalGalerkin fi-nite element methodintime was implementedin [11,20,22]. Integratorsfor variationalprinciplesincludingforcing were developedamongothersin

[20]

.

Thestudyofthevariationaltimeintegratorsthatarediscussedinthispaperwasmotivatedbyourresearchinto numer-icaldiscretizationsoftheequationsdescribingnonlinearwaterwaves,whichhaveaHamiltonianstructure.TheHamiltonian structure of nonlinear gravity free surface water waveswas shownby Zakharov [29].In [9] (and[8]), we constructed a numericalwaterwavetankbasedonadiscreteMiles’variationalprinciplewithforcing.Thecomputationalmethodhadto fulfillanumberofrequirementsinordertoperformthesimulationsofthelaboratoryexperimentspresentedin

[14]

.First, theexperimentswithwhichwecomparednumericalresultsrunoveralongperiodoftime,coveringalsoanextensive spa-tialdomain,whichmeantthatconservationofmassandboundedfluctuationsofenergyatthediscretelevelareimportant. Second,inoneofthetestcasesin

[9]

(also

[8]

)afocussingwavewascreatedwithawave heightfivetimesexceedingthe ambient,incomingwaveheight.Therefore numericalstabilityhadtobeensured. Third,apiston wavemaker wasusedin theexperiments,whichgivesexternalforcingtothesystemandresultsinanon-autonomousHamiltoniansystem.Itis non-trivialtomeetall theserequirements,inparticularincombinationwitha spatialdiscretizationofthe nonlinearpotential flowwaterwaveequations.Wethereforeinvestigatevariationaltimeintegratorsindetailinthatpaper.

In the construction of variational time integrators we choose to work with a discontinuous Galerkin finite element method in time, because it allows additional freedom to develop new symplectic integrators, also for non-autonomous Hamiltoniansystems.Inthismethodthetimedomainissplitintoelements,andineachelementthevariablesare approxi-matedwithpolynomialexpansions.Thediscretevariationalprincipleisthenobtainedbysubstitutionoftheapproximations in the continuum variational principle, provided one calculates the contributions of the delta functionsthat arise. After evaluation ofthe variational principle by calculating its variations, one obtains the discrete system of equations forthe variationaltimeintegrationmethod.

Acrucialpartofthediscontinuousfiniteelementmethodisthenumericalfluxarisingfromthebefore-mentioneddelta functions.Ithastobeintroduced toestablisha connectionbetweenthetimeelements sinceineachelementonlyalocal polynomialapproximationisused.Westudya speciallychosen numericalflux,somewhat inspiredbytheoneexploredin

[24] and

[25]

,whichinturnwas basedonthetheoryofnon-conservative productsdevelopedin

[4]

and

[28]

.Thefluxis alsoimplementedinthenumericaldiscretizationforthenonlinearpotentialflowwaterwaveequations,discussedin

[8,9]

andgivesexcellentresults.

Combiningthediscretevariationalprinciplewiththederived numericalflux,weobtainaunifiedapproachtoconstruct timeintegratorsfor(non-)autonomous Hamiltoniansystems.Wehavederivedbothwell-knownandnovelsymplectic time steppingschemesoffirst,secondandthirdorderaccuracy.Anextensiveanalysisforthediscretizationsisprovided,including alinearstability analysisandaninvestigation ofthesymplectic natureoftheschemes.The novelthirdorder schemewe have developedis shownto have improveddispersion properties. Furthermore,within this approachwe derive andtest time stepping schemes for non-autonomous Hamiltoniansystems, such asforced anddamped oscillators. We study the approximationofforcinganddampingtermsconsideringthediscreteversionsofcertainintegrals.

In

[10]

,a unifieddiscontinuous Galerkin frameworkfortime integration isgivenusing weighted residualsfor general nonlinearordinary differentialequations.New examplesaregivenoftwo(implicit)symplecticRunge–Kuttamethods, sat-isfyingthe symplectic condition by construction.Our workis complementary aswe derivediscrete variational principles by staying withinthe variationalframework, whichvariationsubsequently yieldsthe nonlinearalgebraic setofequations to be solved. While our schemes are always variational by construction, the symplectic andstability conditions need to be verified. These conditionscan be met depending on the choice of quadraturewe useand thechoice ofthe free pa-rameter(s) in thejump conditions. This complementsZhao and Wei’s work [10], also because we developed our stable, variationalandsymplecticintegratorsforapplicationstowaterwaveproblems,sincethesecanthenbeexpressedasspace andtime discretevariationalprinciples

[9,27]

. Extensionsofthe researchpresented herearefound inBokhoveand Kalo-girou[3].

Theoutlineofthispaperisasfollows.InSection2weintroducethedynamicsofaHamiltoniansystemwithonedegree offreedom. Next,Section 3introduces thevariational discontinuousGalerkin method.Herewe obtain secondorder time stepping schemes, which can be formulated in the form of well-known symplectic schemes. In Section 4 a novel time integratorofthirdorderaccuracyintimeisderivedandanalyzed.Weanalyzeitslinearstabilityandsymplecticconditions. InSection 5wepresentnumericalresultsfornon-autonomoussystemsby extendingthewell-knownsymplectic schemes wederivedinthiscomplementaryvariationalframework.Ournovelthirdordertimeintegratorisusedtosimulatenonlinear waterwavesinSection6.ConclusionsaredrawninSec.7.Moredetailsontheanalysisofthetimeintegrationmethodsare presentedintwoappendices.

where p

(

t

)

isthemomentum,q

(

t

)

thecoordinate,andH

(

p

,

q

)

theenergyorHamiltonian.Thevariationalprincipleforthe systemis

δ

L

:= δ

p

L

+ δ

q

L

=

0,where

δ

p

L

and

δ

q

L

arethefunctionalderivativeswithrespectto p and q. Thefunctional

derivativesaredefinedas

δ

L

=

lim →0 1





L



p

+



δ

p

,

q

+



δ

q



−

L



p

,

q



.

(2)

ConsiderthevariationoftheLagrangianfunctional

(1)

δ

L

(

p

,

q

)

=

T



0



_dq dt

δ

p

+

p d

δ

q dt

−

∂

H

∂

p

δ

p

−

∂

H

∂

q

δ

q



dt

=

0

.

(3)

Afterintegrationby partsof

(3)

,usingtheendpointconditions

δ

q

(

0

)

= δ

q

(

T

)

=

0,andthearbitrarinessofthevariations, thedynamicsofaHamiltoniansystememergesasfollows:

δ

p

:

dq dt

−

∂

H

∂

p

=

0 and

δ

q

:

dp dt

+

∂

H

∂

q

=

0

,

(4)

withinitialconditions p

(

t

=

0

)

=

p0 andq

(

t

=

0

)

=

q0.

3. VariationaldiscontinuousGalerkintimediscretizationofaHamiltoniansystem

3.1. Discretefunctional

The accurate numerical approximation ofequations (4)is a well-developed subjectof research. Geometric integrators target thepreservationofsystempropertiesatthe discretelevelandsymplecticintegratorsensurethat theapproximated systemisalsoHamiltonian

[12]

.WestarttoapproximatetheLagrangian

(1)

ratherthanapproximatingequations

(4)

.

Toapproximatethefunctional

(1)

forasingle-degree-of-freedomsystemintime,wefirstdividethetimedomain

[

0

,

T

]

intoN finitetimeintervalsIn

= (

tn

,

tn+1

)

,n

=

0

,

. . . ,

N

−

1,withlength



t

=

tn+1

−

tn.Eachtimeinterval In hasaconstant

lengthandisrelatedtoareferencedomain

ˆ

I

= (−

1

,

1

)

throughthemappingFn,definedas

Fn

: ˆ

I

→

In

:

τ

→

t

=

1 2



tn

(

1

−

τ

)

+

tn+1

(

1

+

τ

)



.

(5)

Second,functions

(

p

,

q

)

areapproximatedelement-wisewithapolynomialexpansion.Ineveryelementbasisfunctionsare continuous,butcouldbediscontinuousacrosselementboundaries.Thefiniteelementspaceisdefinedas

Vτ

:= {

vτ

|

vτ

∈

L2

(

[

0

,

T

]),

vτ



_I

n

∈

Ps

(

In

),

∀

In

,

n

=

0

,

· · · ,

N

−

1

},

(6)

where Ps is the space ofpolynomials of degree s. Eachvariable in finite element In is now approximated as pτ

∈

V τ ,

qτ

∈

V τ andrepresentedinatimeslab

(

tn

,

tn+1

)

as

pτ

=

pi

ϕ

i

,

qτ

=

qi

ψ

i

,

(7)

withbasis functions

ϕ

i_,

_ψ

i_{,which donot necessarilyhaveto coincide,andexpansion coefficients p}i_,qi_{.We willin}

par-ticularconsiderLagrangeandBernsteinpolynomials.Also,theEinsteinsummationconventionwillbeused,whichimplies summationoverrepeatingindices.

Inordertoaccountforthediscontinuities attheelementboundaryweintroducetherightandlefttracesofavariable

d

(

t

)

attimetnas

dn,+

:=

lim

↓0d

(

tn

+



)

and d

n,−

_:=

_lim

↓0d

(

tn

−



),

(8)

respectively.Attimeleveltn thefollowingjump

[[·]]

andaverage

{

{·}

}

βα operatorscanthenbedefined

[[

d

]]



_

tn

=

dn,−

−

dn,+ and

{{

d

}}

β_α



_

tn

=

Fig. 1. Piecewise linear approximation in time.

where

α

,

β

are arbitraryrealnumbers with

α

+ β =

1 and

α

,

β

≥

0.The discrete functional fortheHamiltonian system

(1)isnow obtainedafterintroducingpolynomial approximationsof p andq intothefunctional andsplitting ofthetime integralintoasummationovertimeintervals

L

τ

₍

_pτ

_,

_qτ

₎

_:=

N−1



n=0 tn+1



tn



pτdq τ dt

−

H

(

p τ

_,

_qτ

₎



_dt

₋

N−1



n=−1

[[

qτ

]]{{

pτ

}}

β_α





tn+1

.

(10)

Thelasttermin

(10)

comesasanumericalflux.1 _{Thefluxisrequiredtoestablishaconnectionbetweenthetimeelement}

boundaries,asweusediscontinuousbasisfunctions.

There is still some freedom in the definition of the discrete Lagrangian (10). In the first place we need to specify the polynomial order ofthe basis functions. As in generalit is not possible ordesirable to compute the integral ofthe Hamiltoniananalytically,wealsohavetochooseaquadraturerule

[20]

.In

[13]

itisproventhattheorderofthequadrature rulewillgiveanupperboundfortheattainableorderofthetimeintegrationscheme.Finally,thechoiceoftheweight

α

in thefluxtermalsoleavessomefreedomtoderiveandoptimizedifferenttimeintegrationmethods.

The derivation ofa time integration scheme isnow performedas follows.We choose a polynomial approximation of orders.Toobtainaschemeoforders

+

1 weshouldtakeaquadraturerulewithanorderofaccuracythatisatleasts

+

1

[13].Aftertheschemeisderived,alinearstabilityanalysiswillbeperformedtofindappropriateweights

α

inthefluxterm. Forournovelschemefurtheranalysisisdonetoshowthattheschemeissymplectic.

3.2. Secondordervariationaltimediscretization:Störmer–Verlet

We have rederivedand reinterpreted the symplectic Euler scheme asa genuine discontinuous Galerkin schemewith piecewiseconstantpolynomialsforp andq,ineachtime slab,aswellasthe(modified)midpointschemeusinga discon-tinuousGalerkinschemeinwhichp andq areexpandedusingquadraticpolynomialsandbothevaluatedatthemid-points

pn+1/2 _andqn+1/2 _{oftheintegrationintervalonly(see,}

_[8]

_).

Nextweshowhowthiscanbedoneforaschemeofsecondorderaccuracybecauseweuseanextensionoftheresulting Störmer–Verlet schemefor the non-autonomoussystems considered in §5. To obtain a second order time integrator we needtochooselinearbasisfunctionsinthepolynomialexpansions.Variousoptionsfortheapproximationofthevariables andintegralsarepossible.Since wewantto derivethe knownStörmer–Verletschemeinournewapproach,wetake the trapezoidalruleforp andthemidpointrulefor q,yielding

tn+1



tn H

(

pτ

,

qτ

)

dt

=



t 2



H

(

pn,+

,

qn+1/2

)

+

H

(

pn+1,−

,

qn+1/2

)



,

(11)

wheretheexpansioncoefficientqn+1/2_{isthevalueinthemiddleofthetimeinterval.Weapproximatethevariablesineach} timeelementInas pτ

=

t n+1

₋

_t



t p n,+

₊

t

−

tn



t p n+1,−

_,

_(12a) qτ

=

2

(

t

−

t n

₎



t q n+1/2

₊

tn

+

tn+1

−

2t



t q n,+

_,

_(12b)

see

Fig. 1

.Introducing

(11)

–(12)inthediscretefunctional

(10)

,weobtain

Fig. 2. Piecewise linear approximation in time.

L

τ

₍

_pτ

_,

_qτ

_,

_t

₎

₌

N

_

−1 n=0

(

pn,+

+

pn+1,−

)(

qn+1/2

−

qn,+

)

−



t 2



H

(

pn,+

,

qn+1/2

)

+

H

(

pn+1,−

,

qn+1/2

)



−

N−1



n=−1

(

2qn+1/2

−

qn,+

−

qn+1,+

)(

α

pn+1,−

+ β

pn+1,+

).

(13)

Applying thevariational principle

δ

L

τ

=

0, usingthearbitrariness ofthe variations andend-pointconditions

δ(

2q−1/2

₋

q−1,+

)

:= δ

q0,−

=

0

,

δ

p0,−

=

0,

δ

qN,+

=

0

,

δ

pN,+

=

0,weobtainthefollowingvariationaltimeintegrator

δ

qn,+

: (β −

1

)

pn,+

+ (

α

−

1

)

pn+1,−

+ β

pn+1,+

+

α

pn,−

=

0

,

(14a)

δ

pn,+

:

qn+1/2

=

qn,+

+ β(

2qn−1/2

−

qn−1,+

−

qn,+

)

+



t 2

∂

H

(

pn,+

_,

_qn+1/2

₎

∂

pn,+

,

(14b)

δ

qn+1/2

:

2

β

pn+1,+

+ (

2

α

−

1

)

pn+1,−

=

pn,+

−



t 2



_∂

_H

₍

_pn,+

_,

_qn+1/2

₎

∂

qn+1/2

+

∂

H

(

pn+1,−

_,

_qn+1/2

₎

∂

qn+1/2



,

(14c)

δ

pn+1,−

:

qn+1/2

−

α

(

2qn+1/2

−

qn,+

−

qn+1,+

)

=

qn,+

+



t 2

∂

H

(

pn+1,−

,

qn+1/2

)

∂

pn+1,−

.

(14d)

In Appendix Bweshow thatthisschemeisstableforweights

α

∈ [

0

.

5

,

1

]

andcoincides underthischoice withthe well-knownsecondorderStörmer–Verlettimeintegrator.Itisinterestingtonotethatwestartwithtwodiscontinuousvariables, butattheendthereisnojumpinp atthetimelevelstn,n

=

0

,

. . . ,

N anymore.Moreover,thepresenceofajumpintheq

variabledoesnotresultinaschemedifferentfromtheclassicalone(see

Appendix B

).

An adjoint Störmer–Verlet schemecan be obtainedvia a swapped linearapproximation of the variables in each ele-ment In qτ

=

t n+1

₋

_t



t q n,+

₊

t

−

tn



t q n+1,−

_,

_(15a) pτ

=

2

(

t

−

t n

₎



t p n+1/2

₊

tn

+

tn+1

−

2t



t p n,+

_,

_(15b)

see

Fig. 2

,andatrapezoidalruleforq andamidpointruleforp

tn+1



tn H

(

pτ

,

qτ

)

dt

=



t 2



H

(

pn+1/2

,

qn,+

)

+

H

(

pn+1/2

,

qn+1,−

)



.

(16)

Analogously totheprevious case,itcanbeshownthat thisweightedschemeisstablefor

α

∈ [

0

,

0

.

5

]

andareformulation allowstorecoverawellknownformulationofthescheme,asgivenin

[12]

.

4. Newvariationaltimeintegrators

Anaturalquestioniswhetheritispossibletofindnewtimeintegratorsfollowingthevariationalapproachdescribedin theprevious section.Inthissectionwe willconstructanewtime integrationschemewiththirdorderaccuracyusingthe variationalapproach.Wewillinvestigateconvergence,symplecticity,linearstabilityanddispersionproperties.

4.1. Thirdorderscheme

Fig. 3. Piecewise quadratic approximation in time. pτ

=

(

t n

₊

_tn+1

₋

_2t

₎₍

_tn+1

₋

_t

₎



t2 p n,+

₊

4

(

t

−

tn

)(

tn+1

−

t

)



t2 p n+1/2

₋

(

tn

+

tn+1

−

2t

)(

t

−

tn

)



t2 p n+1,−

_,

_(17a) qτ

=

(

t n

₊

_tn+1

₋

_2t

₎₍

_tn+1

₋

_t

₎



t2 q n,+

₊

4

(

t

−

tn

)(

tn+1

−

t

)



t2 q n+1/2

₋

(

tn

+

tn+1

−

2t

)(

t

−

tn

)



t2 q n+1,−

_,

_(17b)

ascanbeseenin

Fig. 3

.Therearethreeindependentexpansioncoefficientspresentinthepolynomialrepresentation:pn,+

attheleftboundaryoftheinterval, pn+1/2 atthemiddleoftheintervaland pn+1,− attherightboundaryoftheinterval andsimilarlyfor q.WeapproximatetheintegralovertheHamiltonianwithSimpson’sthirdorderquadratureruleas

tn+1



tn H

(

pτ

,

qτ

)

dt

∼

₌



t 6



H

(

pn,+

,

qn,+

)

+

4H

(

pn+1/2

,

qn+1/2

)

+

H

(

pn+1,−

,

qn+1,−

)



.

(18)

Thediscretefunctional

(10)

thenbecomes

L

τ

₍

_pτ

_,

_qτ

₎

₌

N

_

−1 n=0 tn+1



tn pτdq τ dt dt

−



t 6



H

(

pn,+

,

qn,+

)

+

4H

(

pn+1/2

,

qn+1/2

)

+

H

(

pn+1,−

,

qn+1,−

)



−

N

_

−1 n=−1

(

qn+1,−

−

qn+1,+

)(

α

pn+1,−

+ β

pn+1,+

).

(19)

This choice of quadrature rule and basis functions corresponds to the choices made in [22] for the scheme called “ P 2N3Q 4Lob” in their paper (keeping this acronym). Nevertheless, the resulting scheme is different due to the choice ofthenumericalflux.

Applyingthevariationalprinciple

δL

τ

₌

_{0,usingthearbitrarinessofthevariations andend-pointconditions}

_δ

_q0,−

₌

_0,

δ

p0,−

₌

_0,

_δ

_qN,+

₌

₀

_,

_δ

_pN,+

₌

_{0,wegetanewthirdorderscheme}

qn+1/2

= (

1

−

3 2

β)

q n,+

₊

3 2

β

q n,−

₊



t 4



_∂

_H

₍

_pn,+

_,

_qn,+

₎

∂

pn,+

+

∂

H

(

pn+1/2

_,

_qn+1/2

₎

∂

pn+1/2



,

(20a) pn+1/2

= (

1

−

3 2

α

)

p n,+

₊

3 2

α

p n,−

₋



t 4



_∂

_H

₍

_pn,+

_,

_qn,+

₎

∂

qn,+

+

∂

H

(

pn+1/2

_,

_qn+1/2

₎

∂

qn+1/2



,

(20b) qn+1,−

=

qn,+

+ 

t

∂

H

(

p n+1/2

_,

_qn+1/2

₎

∂

pn+1/2

,

(20c) pn+1,−

=

pn,+

− 

t

∂

H

(

p n+1/2

_,

_qn+1/2

₎

∂

qn+1/2

,

(20d)

α

qn+1,+

= −

1 6q n,+

₊

2 3q n+1/2

_{+ (}

_α

₋

1 2

)

q n+1,−

₊



t 6

∂

H

(

pn+1,−

,

qn+1,−

)

∂

pn+1,−

,

(20e)

β

pn+1,+

= −

1 6p n,+

₊

2 3p n+1/2

_{+ (β −}

1 2

)

p n+1,−

₋



t 6

∂

H

(

pn+1,−

_,

_qn+1,−

₎

∂

qn+1,−

.

(20f)

Equations

(20a)

and

(20b)

havetobesolvedtogetherasanimplicitsystem.Theotherfourequationsareexplicit.

4.1.1. Linearstability

ToperformthelinearstabilityanalysisweintroducetheHamiltonianoftheharmonicoscillatorH

(

p

,

q

)

=

1₂

ω

2_q2

₊

1 2p

2 intoscheme

(20)

.WiththeprescribedHamiltonianwecanformulatescheme

(20)

inaform

Fig. 4. Modulusofeigenvaluesλiofsystem(20)foralinearharmonicoscillator.Theshadowedareacorrespondstotheregionwherethemodulusofall theeigenvaluesisequaltoone,i.e.,wheretheschemeisstable.

Fig. 5. Modulus of the eigenvalues for the linear system(22)withα=0.5 – solid line versusλ=1 – dotted line.

Toensure thatthenumericalsolutionisbounded, itisrequiredfortheeigenvalues

λ

i ofthematrix A to belessthanor

equal toone inmodulus

[12]

.Thecharacteristicpolynomial ofthematrix A isaquartic functionin

λ

.Since,theexplicit formoftheeigenvaluesiscomplicatedwe thereforecompute thestabilitycriterianumericallyinMatlab.We takearange of weights

α

, a range of frequencies multiplied by time step

ω

t, substitute them into the expressions for the exact solutions ofa quartic function, andshow thedomains wherethe maximum in modulus eigenvalueis smaller than one

maxi=1,...,4

|λ

i

|

≤

1.Onecan seein Fig. 4that theonly stability regionlocated near

ω

t

=

0 is theone corresponding to

α

=

0

.

5.

For

α

=

0

.

5 andtheHamiltonianoftheharmonicoscillatorthetransformationmatrixbecomesequalto

A

=

⎛

⎜

⎜

⎜

⎝

−

3 5g

−

12g 5t 1

−

g

−

ω

2



t 1−g 4 12g 5tω2

−

3₅g



t 1−g 4 1

−

g 1

−

g

−

ω

2

_

_t1−g 4 1 15g

(

t2

ω

2

−

9

)

4g 15t

(

t2

ω

2

−

9

)



t1−₄g 1

−

g

−

4g 15tω2

(

t2

ω

2

−

9

)

₁₅1g

(

t2

ω

2

−

9

)

⎞

⎟

⎟

⎟

⎠

,

(22)

where g

=

_5_t2_ωt22ω₊2₁₆.The eigenvalues are computednumerically, see

Fig. 5

,which isa two-dimensional slice ofthe plot in

Fig. 4

withafixedweight

α

=

0

.

5.Withinthe(grey)stability domain,alltheeigenvaluesareequaltoone inmodulus, whichleadstoanabsenceofdissipationinthenumericalscheme.Thenumericalcalculationoftheeigenvaluesprovidesthe linearstability condition

|

ω



t

|

≤

1

.

757,whichisslightlymorerestrictive thanthecriteriafortheStörmer–Verletscheme. ComparedwiththeschemeP 2N3Q 4Lob in

[22]

,weseethatbumpsarepresentinthelinearstabilitydomainforbothcases, butthestability restrictionfor P 2N3Q 4Lob islessrestrictiveandequals

|

ω



t

|

≤

2

√

2.Astheauthorsmentionin

[22]

,the scheme P 2N3Q 4Lob isimplicitandall theequationshaveto besolved simultaneously, whichmayposeamore difficult numerical challenge. The third orderscheme proposed heredoes not possessthisshortcoming, as willbe demonstrated further.

Sincethevalue

α

=

0

.

5 givesthebeststabilityresults,thetimeintegrationscheme

(20)

becomes qn+1/2

=

3 4q n,−

₊

1 4q n,+

₊



t 4



_∂

_H

₍

_pn,+

_,

_qn,+

₎

∂

pn,+

+

∂

H

(

pn+1/2

_,

_qn+1/2

₎

∂

pn+1/2



,

(23a) pn+1/2

=

3 4p n,−

₊

1 4p n,+

₋



t 4



_∂

_H

₍

_pn,+

_,

_qn,+

₎

∂

qn,+

+

∂

H

(

pn+1/2

,

qn+1/2

)

∂

qn+1/2



,

(23b)

Table 1

ConvergenceinL∞-,L2_{-normsoftheerrorinq foraharmonicoscillator,andL}∞_{-errorintheenergyofthesystem.}

Initialconditionsareq0,−_=q0,+_{= −0.001, p}0,−_=p0,+_=0,ω2₌₀

.1,thebiggesttimestep t=1,thefinaltime

T=40.Theprecisionsettosolvetheimplicitsystemis=10−12_.

Step q q Energy

L∞-error order L2_{-error order}

H-error order

t 4.1204E−6 – 6.1191E−6 – 1.3489E−9 –

t/2 5.1937E−7 2.9879 7.0457E−7 3.1185 1.6565E−10 3.0256

t/4 6.5058E−8 2.9970 8.6159E−8 3.0317 2.0613E−11 3.0065

t/8 8.1366E−9 2.9992 1.0710E−8 3.0081 2.5735E−12 3.0017

t/16 1.0172E−9 2.9998 1.3369E−9 3.0020 3.2171E−13 2.9999

t/32 1.2716E−10 3.0000 1.6705E−10 3.0005 4.0211E−14 3.0001

t/64 1.5895E−11 3.0000 2.0879E−11 3.0001 5.0263E−15 3.0000

qn+1,−

=

qn,+

+ 

t

∂

H

(

p n+1/2

_,

_qn+1/2

₎

∂

pn+1/2

,

(23c) pn+1,−

=

pn,+

− 

t

∂

H

(

p n+1/2

_,

_qn+1/2

₎

∂

qn+1/2

,

(23d) qn+1,+

=

4 3q n+1/2

₋

1 3q n,+

₊



t 3

∂

H

(

pn+1,−

,

qn+1,−

)

∂

pn+1,−

,

(23e) pn+1,+

=

4 3p n+1/2

₋

1 3p n,+

₋



t 3

∂

H

(

pn+1,−

,

qn+1,−

)

∂

qn+1,−

.

(23f)

The accuracy of the scheme (23) is verified for a linear systemwith Hamiltonian H

=

1₂p2

+

1₂

ω

2_q2_{. The results in}

Table 1showthatthescheme

(23)

isthirdorderaccurate.WedefinetheL2_errorasL2

err

=



N n=1



tn tn−1



qτ

−

qex

(

t

)



2 dt for thepolynomial approximationqτ definedin

(17b)

andtheexact solutionqex oftheautonomoussystemderived in(A.8).

Thefinalresultsoftheschemearethevaluesqn+1,+_{, p}n+1,+_{;thereforetheL}∞_–errorand

_

_{H –erroraredefinedas} L∞_err

=

max 1≤n≤N

|

q n,+

₋

_q ex

(

tn

)|,

(24a)



Herr

=

max 1≤n≤N



H

(

pn,+

,

qn,+

)



−

min 1≤n≤N



H

(

pn,+

,

qn,+

)



.

(24b) 4.1.2. Dispersionrelation

Thedispersionanalysisisperformedfollowing

[15]

.Inordertocomputethedispersionpropertiesofthenewvariational timeintegrationscheme

(23)

,wetaketheFouriermodes(i.e.,wenowtake

λ

=

e−iωˆt _suchthat

_|λ|

₌

₁₎

⎛

⎜

⎜

⎝

pn,− qn,− pn,+ qn,+

⎞

⎟

⎟

⎠ =

⎛

⎜

⎜

⎝

ae−iωˆnt be−iωˆnt ce−iωˆnt de−iωˆnt

⎞

⎟

⎟

⎠

(25)

andsubstitutetheseintothesystem

(23)

,whichgivesalinearsystem

e−iωˆt

(

a

,

b

,

c

,

d

)

T

=

A

(

a

,

b

,

c

,

d

)

T

.

(26)

Ingeneral, we are notinterested inthe resultsforthe variables pn+1,−,qn+1,−,asthose areintermediate data.What mattersisarelationforthevariables pn+1,+_,qn+1,+_{.Inthecontinuumcase,asshownin}

_{Appendix A}

_{,theexactdispersion}

relationisc

= ±

id

ω

,cf.

(A.7)

.Nevertheless,astherearefourequations,allofthemhavetobeincludedintheanalysis.We manipulatesystemofequations

(26)

symbolicallywithMatlab tofindtherelations betweena, b,c andd;and, toobtain thenumericalfrequency

ω

ˆ

.

Forthedispersionanalysisweintroducethevariables

x

=



1

−



t 6

_ω

6 36

(

4

− 

t2

_ω

2

₎

2

,

y

=



t3

_ω

3 6

(

4

− 

t2

_ω

2

₎

.

(27)

Notice,thatx

→

1

,

y

→

0,when



t

→

0.Thefirstfamilyofsolutionsof

(26)

is a

= −

ib

ω

,

c

= −

id

ω

,

b

= −



x

+

iy



d

,

(28a)

ˆ

ω

=

1



tarccos

−

4

(

4

− 

t2

_ω

2

₎

_x

_{+ }

_t2

_ω

2

₍

t2_ω2 6

−

3

)



t2

_ω

2

₊

₁₆

.

(28b)

Fig. 6. Adifferencebetweennumericalandexactfrequenciesω− ˆω.ThereddashlinepresentsthedispersionrelationfortheStörmer–Verletschemeand the

blue

solidlinethedispersionrelation(29b)forthethirdorderscheme(23).

Wesee,thatthesecondrelationin

(28)

istheexactone.Thenumericalfrequencyis,however,wrongascanbeverified bytakingthelimit



t

→

0 in

(28b)

.

Thesecondfamilyofsolutionsis a

= −

ib

ω

,

c

= −

id

ω

,

b

=



x

−

iy



d

,

(29a)

ˆ

ω

=

1



tarccos 4

(

4

− 

t2

ω

2

₎

_x

_{+ }

_t2

_ω

2

₍

t2_ω2 6

−

3

)



t2

_ω

2

₊

₁₆

.

(29b)

The relationc

= −

id

ω

istheexactdispersion relation,andthenumericalfrequencyconvergestothe exactone

ω

ˆ

→

ω

if



t

→

0 in

(29b)

,see

Fig. 6

.

Thethirdfamilyofsolutionsis a

=

ib

ω

,

c

=

id

ω

,

b

=



x

+

iy



d

,

(30a)

ˆ

ω

= −

1



tarccos 4

(

4

− 

t2

_ω

2

₎

_x

_{+ }

_t2

_ω

2

₍

t2_ω2 6

−

3

)



t2

_ω

2

₊

₁₆

.

(30b)

Inthiscasethenumericalfrequencyconvergestothenegativevalue,

ω

ˆ

→ −

ω

when



t

→

0. Thefourthfamilyofsolutionsis

a

=

ib

ω

,

c

=

id

ω

,

b

=



−

x

+

iy



d

,

(31a)

ˆ

ω

= −

1



tarccos

−

4

(

4

− 

t2

ω

2

₎

_x

_{+ }

_t2

_ω

2

₍

t2_ω2 6

−

3

)



t2

_ω

2

₊

₁₆

.

(31b)

Herethenumericalfrequencydoesnotconvergetothecorrectvaluewhen



t

→

0.

Sincetherearetwoparasiticmodesweneedtobecarefulthattheydonotpollutethesolution.Fortheparasiticmodes therelationb

→ −

d holds.Oneimportantstepistochoosea

=

c andb

=

d,i.e.take p0,−

₌

_p0,+_andq0,−

₌

_q0,+_asinitial

conditions.Inthiscasetheparasiticmodesareeliminated.Forthenonlinearwaterwavecaseweverifiedtheconsequence ofthisinitialfilteringnumericallyin§6,whichappearstobesufficienttoensurestabilityinthenonlinearcase.Thisabsence ofaninitialjumpleadstosatisfactoryresults.

4.1.3. Symplecticity

Ifaone-degree-of-freedomHamiltoniansystemisapproximatedwithtwovariables pn andqn,thesymplecticity condi-tionisformulatedas

MTJ−1M

=

J−1

,

(32)

wheretheJacobianofthetransformationisdefinedas

M

=

∂(

q

n+1

_,

_pn+1

₎

∂(

qn

_,

_pn

₎

.

(33)

Thestructurematrix J ,equalto

J

=



0 I

−

I 0



,

(34)

Fig. 7. Evolutionofvariables(pn,+_,qn,+_{)inphasespaceoveronetimeperiod.System(23)withω=0.3162,t=0.2484 issolvedforfourinitialconditions.}

Theerrorinsolvingthelinearsystemintheimplicitstepof(23)issetto=10−15_.

where I istheidentity matrixwiththedimensionequal tothedegreesoffreedomofthe Hamiltoniansystem,relatesto theHamiltoniansystemas



q p



t

=

J

∇

H

=



0 1

−

1 0

 

Hq Hp



.

(35)

Here thesubscript

(

·)

t refers totakinga time derivative and Hq

≡ ∂

H

/∂

q, etc.Fordetails we referto

[12]

. Sincewe are

workingwith an extended system, i.e.a one-degree-of-freedom Hamiltoniansystem is approximatedwith fourvariables

pn,−

,

pn,+andqn,−

,

qn,+,thematricesinthesymplecticitycondition

(32)

arealsoextended.Wetake

M

=

∂(

q

n+1,−

_,

_qn+1,+

_,

_pn+1,−

_,

_pn+1,+

₎

∂(

qn,−

_,

_qn,+

_,

_pn,−

_,

_pn,+

₎

,

(36)

andthestructurematrix J ,definedasin

(34)

fora two-dimensionalsystem.TheHamiltoniansystemcanthenbewritten as

⎛

⎜

⎜

⎝

q− q+ p− p+

⎞

⎟

⎟

⎠

t

= (

J

)∇

H

=

⎛

⎜

⎜

⎝

0 0 0 1 0 0 1 0 0

−

1 0 0

−

1 0 0 0

⎞

⎟

⎟

⎠

⎛

⎜

⎜

⎝

Hq− Hq+ Hp− H_p+

⎞

⎟

⎟

⎠ .

(37)

WeverifiedinMatlabusingsymbolicmanipulationthatthecondition

(32)

isvalidforthenewsystem

(23)

.Thedetailsof thematrixM arenotgivenexplicitlyheresincetheexpressionsaretoolargetopresent.

For a Hamiltoniansystem with one degree of freedom the symplecticity condition means preservation of an area in phasespace

[12]

.Aphase–spaceplaneisshownin

Fig. 7

.A pointontheplaneisgivenas

(

pn,+

,

qn,+

)

,asourprimeinterest liesinthesevariables.We definefour initialconditionsinthe plane,that constitutea quadrilateralshape.Thenumerical approximationofaharmonicoscillatorsystemdevelopsaccordingtothescheme

(23)

.Thesystemcyclesovertimethrough thesamestates,whicharedefinedbytheinitialconditions.The areaoftherectangularshapeoscillatesaroundtheinitial state,see

Fig. 8

.

From a geometrical point ofview, the symplecticity condition (32) means in the generalcase withd

>

1 degrees of freedomthatthesumoftheorientedareasoftheprojectionsofashapein2d spaceontotheplanes

(

pi

_,

_qi

₎

_isconserved [12].Therefore

Figs. 7 and 8

donotprovesymplecticity,butillustratethepropertyoftheplus-variables pn,+

_,

_qn,+_tocycle

throughthesamestates.

Theconclusions onthenovelthirdorderscheme

(23)

areasfollows.Theschemeissymplectic,butthelinearstability conditionismorerestrictive thantheconditionfortheStörmer–Verletscheme.Theschemeisthird orderaccurate, which isverifiedbothforlinearandnonlinearproblems(see§6andalsoBokhoveandKalogirou

[3]

).

5. Applicationstonon-autonomoussystems

Inthefollowingsection,wedeveloptimeintegrationofnon-autonomousHamiltoniansystems.

5.1. Dampedoscillator

Consideradampedone-degree-of-freedomHamiltoniansystemwiththefollowingnon-autonomousvariationalprinciple, (see

[23]

and

[27]

forrelatedexamples):

Fig. 8. Oscillationintheareaofaquadrilateralshapeinphase–spaceshowninFig. 7fortheharmonicoscillatorcomputedwiththethirdorderaccurate symplectictimeintegrator(23).System(23)issolvedover10timeperiodswithω=0.3162,t=0.2484.Theprecisioninsolvingtheimplicitstepin(23)

issetto=10−15_. 0

= δ ˜

L

(

p

,

q

)

= δ

T



0



pdq dt

−

H

(

p

,

q

)



eγtdt

=

T



0



_dq dt

−

∂

H

∂

p



δ(

peγt

)

−



_d

₍

_peγt

₎

dt

+

∂

H

∂

qe γt



_δ

_qdt

₊

_peγt

_δ

_q

_|

T 0

,

(38)

with

γ

>

0 the coefficient of the linear momentum damping and end-point conditions

δ

q

(

0

)

= δ

q

(

T

)

=

0. Variation of

(38)withrespecttothevariables p andq yieldsHamilton’sequationsforthedampedoscillatorwhentheHamiltonianis

H

(

p

,

q

)

=

1₂p2

+

1₂q2 (andunitfrequency).Wefind

δ(

peγt

)

:

dq dt

=

∂

H

∂

p

=

p

,

(39a)

δ

q

:

dp dt

+

γ

p

= −

∂

H

∂

q

= −

q

.

(39b)

Withthecoordinatetransformationq

=

Q exp

(

−

γ

t

/

2

)

andp

=

P exp

(

−

γ

t

/

2

)

,thevariationalprinciplethenbecomes

0

= δ

T



0



PdQ dt

− ˜

H

(

P

,

Q

)



dt

,

(40)

withmodifiedHamiltonian

˜

H

(

P

,

Q

)

=

1 2P 2

₊

1 2Q 2

₊

1 2

γ

P Q

.

(41)

Hence, thesystem

(40)

isreformulatedasanautonomoussystem,viz.themodifiedHamiltonian

(41)

istime independent inthenewvariables P and Q .WecannowmonitorthismodifiedHamiltonianatthediscretelevel.

5.1.1. Discretization

Westartwithafirstorderapproximationofthefunctional

(38)

intime.Weconsiderthediscretevariationalprinciple

˜

L

τ

₍

_pτ

_,

_qτ

₎

_:=

N

_

−1 n=0 tn+1



tn



pτdq τ dt

−

H

(

p τ

_,

_qτ

₎



_eγt_dt

₋

N

_

−1 n=−1

[[

qτ

]]{{

pτeγt

}}

βα





tn+1

,

(42)

where the variables

(

p

,

q

)

are approximated as piecewise-constant per time interval

(

tn

,

tn+1

)

. As in [8], to obtain the well-knownformulationofthesymplecticEulerscheme,werenametheapproximationsasqτ

=

qn,+

,

pτ

=

pn+1,−andtake

α

=

1

,

β

=

0.Theexponentialtime-dependenttermresponsibleforthedampingisnewandwillbeapproximatedusingthe valueattherightboundaryofthetimeelementeγtn+1

. Substitutingtheseapproximationsinto

(42)

,weobtain

˜

L

τ

₍

_pτ

_,

_qτ

₎

_:=

N

_

−1 n=0

−

t H

(

pn+1,−

,

qn,+

)

eγtn+1

−

N

_

−1 n=−1

(

qn,+

−

qn+1,+

)(

pn+1,−eγtn+1

).

(43)

Fig. 9. ModifiedtotalenergyHexforadampedoscillatorversustime.Thediscretemodifiedenergyforapproximation(44a)isplottedasasolidline,the discretemodifiedenergyforapproximation(48)isplottedwithcrossesandthediscretemodifiedenergyforapproximation(47)isplottedversustime withcircles.

We take variations

δ ˜

L

τ

=

0 of (43) with respect to the variables pn+1,−

,

qn,+, use the arbitrariness of variations and end-pointconditions

δ

q0,−

=

0,

δ

p0,−

=

0,

δ

qN,+

=

0,toobtain

pn+1,−

=

pn,−

− 

t

∂

H

(

p n+1,−

_,

_qn,+

₎

∂

qn,+

−



1

−

eγ(tn−tn+1)



pn,−

,

(44a) qn+1,+

=

qn,+

+ 

t

∂

H

(

p n+1,−

_,

_qn,+

₎

∂

pn+1,−

,

(44b)

afterdivisionbytheexponentialtermeγtn+1.Weseethatthefirststepisimplicitandthesecondstepisexplicit. Theexponentialtermin

(44a)

needssomeclarification.Using

1

−

eγ(tn−tn+1)

=

1

−

e−γt

≈

γ



t

,

(45)

showsthatthediscreteequation

(44a)

isanapproximationtothedampedequation

˙

p

+

∂

H

∂

q

+

γ

p

=

0

.

(46)

Thisensuresthatthemodifiedtotalenergy

(41)

willbepreservedaswillbedemonstratednext.

5.1.2. Numericalresults

Consider theHamiltonian H

(

p

,

q

)

=

1₂p2

+

1₂q2.Wecompare threecasesforthedampedvariational principle,viz. ap-proximation

(44a)

andthemorestraightforward(forwardandbackwardEuler)approximationsforthedampingterm:

pn+1

=

pn

− 

t

∂

H

(

p n+1

_,

_qn

₎

∂

qn

− 

t

γ

p n

_,

₍₄₇₎ or pn+1

=

pn

− 

t

∂

H

(

p n+1

_,

_qn

₎

∂

qn

− 

t

γ

p n+1

_,

₍₄₈₎

whilethesecondequation

(44b)

inthesystemremains unchanged.Weverifytheabsenceofdriftinadiscreteversion of themodifiedenergy

(41)

ateachtimeleveltn,n

=

0

,

. . . ,

N,

Hex

(

pn,−

,

qn,+

,

tn

)

=

1 2



(

pn,−

)

2

+ (

qn,+

)

2

+

γ

pn,−qn,+



eγtn (49)

forthedifferentapproximations.Theresultsofthediscretemodifiedenergyarepresentedin

Fig. 9

forthedifferent approx-imations versus time. Wesee, that themodified energyfor approximation

(44a)

oscillatesaround the initial state,while the other two approximations reveal a drift in energy. The latter was expected for the forward Eulerapproximation of theforcing termin

(47)

.Nevertheless,theresultisunexpected forthe implicitscheme

(48)

,wheretheapproximation of thedampingtermisdone analogouslytothesymplectic Eulerscheme.Thereforeourvariationalapproachallowsalsothe derivationoftimeintegratorsthatresultinaboundedfluctuationofspecialintegrals.

Weconsiderthesimplecaseofaforcedoscillator.TheHamiltonianisH

(

p

,

q

)

=

₂1p2

+

1₂q2withforcingfunction f

(

p

,

q

,

t

)

=

−

a

γ

tq,wherea

γ

istheamplitude.Thecorrespondingsystemreads

dq

dt

=

p and

dp

dt

= −

q

+

a

γ

t

.

(51)

Thissystemisnon-autonomous,butwecanintroduceacoordinatetransformationtomakeitautonomous,asfollows

q

=

Q

+

a

γ

t

,

p

=

P

.

(52)

TheLagrangian

(50)

foraforcedharmonicoscillatoristhenreformulatedas

L

(

p

,

q

)

=

T



0



PdQ dt

−



1 2P 2

₊

1 2Q 2

₋

_a

_γ

_P



_dt

₊

1 6a 2

_γ

2_T3

_.

₍₅₃₎

AnimportantquantitytomonitoristheHamiltonianinthenewcoordinates

ˆ

H

=

1 2P 2

₊

1 2Q 2

₋

_a

_γ

_P

_,

₍₅₄₎

whichwillbeusedtocomputetheorderofaccuracyoftheforcedscheme.Theexactsolutiontotheproblem

(51)

is qex

=

q

(

0

)

cos

(

t

)

+



p

(

0

)

−

a

γ



sin

(

t

)

+

a

γ

t

.

(55)

Next,wewillconsiderthesymplecticStörmer–Verletschemetocomputetheevolutionofthisforcedsystem.

5.2.1. ForcedsymplecticStörmer–Verletscheme

Weconsiderthesecondorderapproximationwithaddedforcing,andinvestigatetheaccuracyfordifferentdiscretizations oftheforcingterm.WeconstructaStörmer–Verlettimesteppingmethodadjointtotheonediscussedindetailin§3.2.We take appropriate polynomial approximations (15), the mixedquadrature rule (16)and

α

=

0. Then we obtaina discrete approximationoftheLagrangian

(50)

formulatedasfollows

L

τ

₍

_pτ

_,

_qτ

₎

₌

N

_

−1 n=0

pn+1/2

(

qn+1,−

−

qn,+

)

−



t 2



H

(

pn+1/2

,

qn,+

)

+

H

(

pn+1/2

,

qn+1,−

)



−



t 2



f

(

pn+1/2

,

qn,+

,

t∗

)

+

f

(

pn+1/2

,

qn+1,−

,

t∗∗

)



−

N

_

−1 n=−1

(

qn+1,−

−

qn+1,+

)

pn+1,+

.

(56)

Computingthevariations,weobtainthefollowingsystemofequations

δ

pn,+

:

qn+1,+

=

qn+1,−

,

(57a)

δ

qn,+

:

pn+1/2

=

pn,+

−



t 2

∂

H

(

pn+1/2

_,

_qn,+

₎

∂

qn,+

−



t 2

∂

f

(

pn+1/2

_,

_qn,+

_,

_t∗

₎

∂

qn,+

,

(57b)

δ

pn+1/2

:

qn+1,−

=

qn,+

+



t 2



_∂

_H

₍

_pn+1/2

_,

_qn,+

₎

∂

pn+1/2

+

∂

H

(

pn+1/2

_,

_qn+1,−

₎

∂

pn+1/2



,

+



t 2



_∂

_f

₍

_pn+1/2

_,

_qn,+

_,

_t∗

₎

∂

pn+1/2

+

∂

f

(

pn+1/2

_,

_qn+1,−

_,

_t∗∗

₎

∂

pn+1/2



(57c)

δ

qn+1,−

:

pn+1,+

=

pn+1/2

−



t 2

∂

H

(

pn+1/2

_,

_qn+1,−

₎

∂

qn+1,−

−



t 2

∂

f

(

pn+1/2

_,

_qn+1,−

_,

_t∗∗

₎

∂

qn+1,−

,

(57d)