Time-symmetry, symplecticity and stability of Euler-Maclaurin and Lanczos-Dyche integration

# Time-symmetry, symplecticity and stability of Euler-Maclaurin and Lanczos-Dyche integration

Charalampos M. Markakis Michael F. O’Boyle Derek Glennon Khoa Tran Pablo Brubeck Roland Haas Hsi-Yu Schive Kōji Uryū DAMTP, Centre for Mathematical Sciences, University of Cambridge, Cambridge, UK School of Mathematical Sciences, Queen Mary University of London, London, UK Department of Physics, University of Illinois at Urbana-Champaign, Urbana, IL, USA NCSA, University of Illinois at Urbana-Champaign, Urbana, IL, USA Department of Mathematics, University of Illinois at Urbana-Champaign, Urbana, IL, USA Department of Mathematics, University of California San Diego, San Diego, CA, USA Argonne National Laboratory, Lemont, IL, USA Physics Department & Institute of Astrophysics, National Taiwan University, Taiwan Department of Physics, University of the Ryukyus, Senbaru, Nishihara, Okinawa, Japan
###### Abstract

Numerical evolution of time-dependent differential equations via explicit Runge–Kutta or Taylor methods typically fails to preserve symmetries of a system. It is known that there exists no numerical integration method that in general preserves both the energy and the symplectic structure of a Hamiltonian system. One is thus normally forced to make a choice. Nevertheless, a symmetric integration formula, obtained by Lanczos-Dyche via two-point Taylor expansion (or Hermite interpolation), is shown here to preserve both energy as well as symplectic structure for linear systems. This formula shares similarities with the Euler-Maclaurin formula, but is superconvergent rather than asymptotically convergent. For partial differential equations, the resulting evolution methods are unconditionally stable, i.e, not subject to a Courant-Friedrichs-Lewy limit. Although generally implicit, these methods become explicit for linear systems.

###### keywords:
Hamiltonian systems, Symplectic integration, Two-point Taylor series, Hermite interpolation, Padé approximant, Multi-derivative methods, Euler-Maclaurin formula, Trapezoidal rule, Hermite rule, Lotkin rule, Lanczos-Dyche formula
###### Pacs:
02.60.Lj, 45.20.Jj, 03.65.-w, 04.25.Dm
journal: Journal of Computational Physics

## 1 Introduction

In this paper we consider the numerical solution of initial value problems of the form

 dudt=f(t,u),u(t0)=u0 (1)

with where is an open set of . Eq. (1) may represent an ordinary differential equation (ODE) for the scalar function , or a system of ODEs, or (in a method of lines approach Alcubierre (2012)) a partial differential equation (PDE) for the vector with components constructed from discretizing a scalar field on the spatial grid-points .

To evolve this differential equation numerically, one may use the fundamental theorem of calculus to write it as an integral equation

 u(tν+1)=u(tν)+∫tν+1tνf(t,u(t))dt (2)

(where counts the sequential time steps of numerical integration) and compute the time integral via a suitable quadrature or other numerical integration rule Miyatake and Butcher (2016). Explicit integration schemes, such as Runge–Kutta methods, or Taylor series of the integrand expanded about , yield numerical integration formulas that violate symmetry under time reversal . Symmetries of the continuum system, such as conservation of energy or symplectic structure, are thus generally violated upon discretization.

It is known that there exists no numerical integration method that in general preserves both the energy and the symplectic structure of a Hamiltonian system Butcher et al. (2016); Chartier et al. (2006); Zhong and Marsden (1988). One is thus normally forced to make a choice between conserving one or the other. However, for quadratic Hamiltonians, we show that it is in fact possible to preserve both energy as well as symplectic structure (up to numerical roundoff error).

Our focus in this paper is to integrate via a single-step rule, originally due to Lanczos Lanczos (1956) and Dyche Dyche (1956), which is based on a two-point Taylor expansion about and , and is symmetric under time reversal . We will refer to this integration rule as the Lanczos-Dyche (LD) formula, and will study the degree to which it preserves symmetries of the continuum system. When compared with Runge–Kutta methods of the same order, these symmetric methods are shown to be efficient, and invariants of the motion are well–preserved over long time intervals, exactly for linear systems and approximately for non-linear systems.

Moreover, when using symmetric methods in a method of lines framework to evolve partial differential equations, the resulting scheme is shown to be unconditionally stable, that is, not subject to a Courant-Friedrichs-Lewy (CFL) condition, thus allowing arbitrarily large time steps in numerical evolution. In terms of accuracy, time steps cannot be arbitrarily large, but the high order of the methods, combined with the absence of a CFL limit, does allow significantly larger time steps compared to explicit methods. Moreover, while symmetric formulas are implicit, for problems that allow inverting and storing the evolution operator in memory, there is no computational overhead compared to explicit methods, as the symmetric formulas can become effectively explicit.

## 2 Numerical integration via Hermite interpolation

Let be a function. We wish to construct approximations to the definite integral

 I(f)=∫t2t1f(t)dt, (3)

for a small interval , given the function values and derivatives at the endpoints of this interval.

### 2.1 One-point Taylor expansion

A one-point Taylor expansion of the integrand about has the form

 f(t)=Pn(t)+Rn(t) (4)

where

 Pn(t)=n∑l=0(t−t1)ll!f(l)1 (5)

is the Taylor polynomial of order , with , and

 Rn(t)=(t−t1)n+1(n+1)!f(n+1)(τ) (6)

is the remainder, with . Integration of Eq. (4) yields the familiar formula

 I(f)=n∑l=1Δtll!f(l−1)1+Rn(f) (7)

with remainder

 Rn(f)=Δtn+1(n+1)!f(n)(τ) (8)

Eq. (1) can be used to express all derivatives in terms of Butcher (2000). This procedure can generally be automated using algebraic manipulation software. Brown (1977); Jeltsch and Kratz (1978). In particular, if Eq. (1) is a linear ODE, then formula (7) is equivalent to a Runge–Kutta method of order .

The above formula has the advantage of being explicit, that is, it only requires initial conditions at time to be specified. It is thus suitable for initial value problems such as Eq. (1). However, this property means that the formula is not symmetric under time reversal and thus fails to preserve symmetries of the differential equation (1). In what follows, we use two-point Hermite interpolation or, equivalently, a two-point Taylor expansion, to derive formulas symmetric under time reversal.

### 2.2 Two-point Taylor expansion

#### Trapezoidal rule (2nd order)

Given the function values and at the endpoints, one may approximate in the interval with a linear Lagrange interpolating polynomial,

 P(t)=t2−tΔtf1+t−t1Δtf2, (9)

uniquely determined by the two conditions and . Substituting this expression into Eq. (3) and integrating yields the trapezoidal rule:

 I(f)=Δt2(f1+f2)+O(Δt3) (10)

with the error term to be computed below.

#### Hermite’s rule (4th order)

If, in addition to the function values, the derivatives and are known at the end-points, one may approximate with a cubic Hermite interpolating polynomial , uniquely determined by the four conditions , , and . Using the method of undetermined coefficients Markakis and Barack (2014); Nozawa et al. (1998), one finds:

 P(t) = (t2−t)2(2t+t2−3t1)Δt3f1−(2t+t1−3t2)(t−t1)2Δt3f2 (11) + (t2−t)2(t−t1)Δt2f′1−(t2−t)(t−t1)2Δt2f′2.

Substituting this expression into Eq. (3) and integrating yields the Hermite rule:

 I(f)=Δt2(f1+f2)+Δt212(f′1−f′2)+O(Δt5) (12)

A detailed comparison between Hermite’s rule and Simpson’s rule may be found in Lampret (2004).

#### Lotkin’s rule (6th order)

If, in addition to the function values and first derivatives, the second derivatives and at the end-points are known, one may approximate with a quintic Hermite interpolating polynomial, uniquely determined by the six conditions , , , , and . The method of undetermined coefficients yields:

 P(t) = [6t2+3t2t+10t21+t22−5t1(3t+t2)](t2−t)3Δt5f1 (13) + [6t2+3t1t+t21+10t22−5(3t+t1)t2](t−t1)3Δt5f2 + (3t−4t1+t2)(t−t1)(t2−t)3Δt4f′1 + + (t2−t)2(t−t1)32Δt3f′′2+(t2−t)3(t−t1)22Δt3f′′1.

Substituting this expression into Eq. (3) and integrating yields a order generalization to the trapezoidal rule, which we will refer to as Lotkin’s Lotkin (1952) rule:

 I(f)=Δt2(f1+f2)+Δt210(f′1−f′2)+Δt3120(f′′1+f′′2)+O(Δt7) (14)

Unlike formulas such as Eq. (7) stemming from one-point Taylor series, one notices that Eqs. (12) and (14) have different coefficients, which depend on the order of the series truncation. This is a characteristic of two-point Taylor series, as shown below.

#### The Lanczos-Dyche formula (2nth order)

The above procedure can be continued up to order if the function values and derivatives , at the end-points are known for . The finite form of the two-point Taylor series is Hummel and Seebeck (1949); Dyche (1956)

 f(t)=Pn(t)+Rn(t) (15)

where

 Pn(t)=n−1∑l=0Δt−(2l+1)[(t−t1)(t−t2)]l[(t−t2)al+(t−t1)bl] (16)

is the two-point Hermite interpolation polynomial,

 al=[dldtlf(t)(t−t2)l]t=t1,bl=[dldtlf(t)(t−t1)l]t=t2 (17)

are constant coefficients, and

 Rn(t)=[(t−t1)(t−t2)]n(2n)!f(2n)(τ),τ∈[t1,t2] (18)

is the remainder Estes (1966). Substituting Eq. (15) into Eq. (3) and integrating yields the order generalization of the trapezoidal rule, which we will refer to as the Lanczos-Dyche111We derived this formula independently, but discovered Lanczos’ book and Dyche’s thesis when this work was near completion, and thus named it after these authors. Despite its remarkable accuracy and symmetry properties, little attention has been paid to the LD formula in the literature, and we found no references other than the ones above. Lanczos (1956); Dyche (1956) (cf. also Squire (1961, 1962); Gould and Squire (1963); Lambert and Mitchell (1963)) formula:

 I(f)=n∑l=1ClnΔtll![f(l−1)1+(−1)l−1f(l−1)2]+Rn(f) (19)

with coefficients

 Cln:=n!(2n−l)!(2n)!(n−l)!, (20)

and remainder

 Rn(f)=(−1)nn!2(2n+1)!(2n)!Δt2n+1f(2n)(τ). (21)

Eqs. (10), (12) and (14) are special cases of the Lanczos-Dyche (LD) formula (19), which forms the basis of this paper.

Higher-order rules easily follow from the above formula. For , the LD formula (19) gives the order rule

 I(f) = Δt2(f1+f2)+3Δt228(f′1−f′2)+Δt384(f′′1+f′′2) (22) + Δt41680(f′′′1−f′′′2)+O(Δt9).

For , the same formula gives the order rule

 I(f) = Δt2(f1+f2)+Δt29(f′1−f′2)+Δt372(f′′1+f′′2) (23) + Δt41008(f′′′1−f′′′2)+Δt530240(f(4)1+f(4)2)+O(Δt11).

and so forth.

As mentioned earlier, the expansion coefficients of the LD formula (19) depend on the order of the series truncation. This is in contrast to familiar formulas based on one-point Taylor series [cf. Eq. (7)] or the Euler-Maclaurin formula [cf. Eq. (24) below] which have fixed coefficients. This dependence on the order of truncation is a characteristic of two-point Taylor series, and is partly responsible for their rapid convergence.

In fact, one of the most interesting features of the LD formula (19) is its remainder term, given by Eq. (21). Inspection of this remainder term, as well as numerical applications outlined in the next sections, reveal that the LD formula is superconvergent: when the function derivatives of order are known, the formula is accurate to order (rather than order , like the usual Taylor formula (7)).

### 2.3 Comparison to the Euler-Maclaurin formula

It is worth comparing the Lanczos-Dyche formula (19) to the more widely known (albeit often asymptotically convergent) Euler-Maclaurin (EM) formula:

 I(f)=n∑l=1BlΔtll![f(l−1)1+(−1)l−1f(l−1)2]+Rn(f) (24)

where, for even , the remainder term is Epperson (2013)

 Rn(f)=−Bn+2Δtn+3(n+2)!f(n+2)(τ),τ∈(t1,t2). (25)

and denotes the Bernoulli numbers

 B1=12,B2=16,B3=0,B4=−130,B5=0,B6=142,... (26)

For and , the EM formula (24) coincides with the trapezoidal rule (10) and Hermite rule (12) respectively. However, the EM formula differs from the LD formula (19) at higher order. For , Eq. (24) yields the order rule:

 I(f)=Δt2(f1+f2)+Δt212(f′1−f′2)−Δt4720(f′′′1−f′′′2)+O(Δt7) (27)

which differs from the Lanczos-Dyche order rule (22), despite the fact that both formulas use up to third-order derivatives of .

We note that, according to Eq. (26), only odd-order derivatives appear in the EM formula (24), and these derivatives appear with opposite sign at the end-points. As a result, the EM formula can be used to construct a composite (multi-step) rule by summing up the contributions from a sequence of points (). For the composite EM formula, all derivative contributions cancel out, except for the values at the first and last point ( and ). The EM formula is thus often used to convert a sum to an integral Tao (2014, 2010). An important application is Riemann zeta-function regularization Hawking (1977); Elizalde and Romeo (1989); Fermi and Pizzocchero (2011), which is used in conformal field theory, renormalization and in fixing the critical spacetime dimension of string theory. An example of zeta-function regularization is the calculation of the vacuum expectation value of the energy of a particle field in quantum field theory. More generally, the zeta-function approach can be used to regularize the full energy-momentum tensor in curved spacetime Parker and Toms (2009).

On the other hand, the LD formula (19) contains both odd- and even-order derivatives, and the latter appear with the same sign at the end-points (e.g. and ). Thus, if one attempts to construct a composite (multi-step) formula, the even-order derivatives do not cancel out. For this reason, the LD formula is not readily222unless the function derivatives can be easily expressed in terms of the function values . suited for converting sums to integrals. Nevertheless, the LD formula is ideally suited for the purpose of numerical integration, i.e. conversion of integrals to sums. In this paper, restrict attention to single-step, rather than composite, formulas.

For the purpose of integration, the difference in accuracy between the LD formula (19) and the EM formula (24) is dramatic. Inspection of the respective remainder terms (21) and (25) reveals that, when the function derivatives of all orders up to are used, the EM formula is (formally) accurate to order whereas the LD formula is accurate to order . Further inspection of the remainder term (25) reveals that, because the Bernoulli numbers decay slowly for increasing , the EM formula is often only asymptotically convergent, and in practice may fail to converge for certain functions or large time-steps. On the other hand, the coefficients in the remainder term (21) decay faster than exponentially for increasing , and the LD formula is superconvergent.

A comparison of the error terms of the EM and LD formulas when computing the integral of the exponential function, shown in Fig. 1, confirms that the latter decreases faster than exponentially with the order of derivatives used.

## 3 Ordinary differential equations

### 3.1 Absolute Stability

To study stability of a time integration scheme, we consider a linear test equation of the form (1) with , , that is,

 dudt=λu,u(0)=1 (28)

and consider the relation between and arising from Eqs. (2), (3) and the relevant integration formula. For any one-step method, this relation has the form

 u(tν+1)=ζ(λΔt)u(tν). (29)

The function is the increment function. For the exact solution of Eq. (28), we have , and thus the exact increment function is

 ζ(μ)=eμ. (30)

A numerical method is considered absolutely stable if the increment function for that method satisfies the condition .

For a Runge–Kutta method of order , the time-stepping has the form of Eq. (29) with an increment function given by

 ζ(μ)=n∑l=01l!μl, (31)

which coincides with the Taylor expansion of order about of the exact increment function (30). These methods exhibit stability within the regions plotted in Fig. (a)a.

On the other hand, substituting the -order derivatives of Eq. (28) into the Lanczos-Dyche formula (19), yields an increment function

 ζ(μ)=∑nl=0Clnl!μl∑nl=0Clnl!(−μ)l, (32)

with coefficients given by Eq. (20). The above expression coincides with the Padé approximant of order to the exact function (30) expanded about . It has already been proven in Ehle (1973) that Padé methods with stabiltity function given by (32) are stable on the entire left half plane for all orders . Thus, the methods resulting from the Lanczos-Dyche formula (19) are A-stable (but not L-stable, as Eq. (32), does not satisfy the condition as ).

We note that implicit multistep methods can only be A-stable if their order is at most 2. This result, known as the second Dahlquist barrier Dahlquist (1956, 1961, 1963); Butcher (2000); Wanner (2006), restricts the usefulness of linear multistep methods for stiff equations. The optimal A-stable method is the trapezoidal rule, which is a special case of the Lanczos-Dyche formula for . However, the use of higher order derivatives allows one to exceed this barrier Brown (1975, 1977); Jeltsch (1976, 1978); Jeltsch and Kratz (1978); Jeltsch (1979); Gekeler (2006), and methods based on the Lanczos-Dyche formula are A-stable to all orders. The property of A-stability will prove extremely useful when applying this method to partial differential equations, as it will be seen that systems discretized via the method of lines using the Lanczos-Dyche formula are unconditionally stable, that is, these methods are not subject to a Courant-Friedrichs-Lewy condition on the time-step.

We notice that the Lanczos-Dyche formula (19) reduces to Padé methods for linear differential equations. However, Padé methods rely on a rational approximation of the exponential function , and are applicable only to linear systems. In contrast, methods resulting from the Lanczos-Dyche formula are generally applicable to both linear and non-linear systems (albeit they are implicit in the latter case). Moreover, even for linear systems, if the system is large, Padé methods provide an approximation to the large matrix . However, the LD formula (19) can be applied to each sub-equation separately, and then one can proceed to invert a linear sub-system of much smaller size, as will be seen below in Sec. 4.1 and in more detail in a companion paper.

### 3.2 Hamiltonian systems

Consider the Hamilton equations

 dqdt = ∂H∂p (33a) dpdt = −∂H∂q (33b)

where the Hamiltonian is a smooth function on an open set of .

To linear order in the Taylor expansion (7), the time evolution of this system could be accomplished via Euler integration

 qν+1 = qν+Δt∂H∂pν (34a) pν+1 = pν−Δt∂H∂qν, (34b)

where counts the time steps. For a typical Hamiltonian of the form

 H=12p2+V(q) (35)

the above discrete transformation is not canonical, as the Jacobian:

 J=∂(qν+1,pν+1)∂(qν,pν) (36)

does not preserve its unitary value, . Instead, the Euler integration (34) introduces an error of order in the Jacobian (36) of the transformation.

One can include higher order terms in the Taylor expansion so the error accumulates more slowly, but symplectic structure would still not be preserved. For example, a second order Runge-Kutta integration yields the integration rule

 qν+1 = qν+Δt∂H∂pν+Δt22ddt(∂H∂pν) (37a) pν+1 = pν−Δt∂H∂qν−Δt22ddt(∂H∂qν), (37b)

Since the partial derivatives of with respect to and are functions of those variables, the chain rule:

 dHdt=∂H∂t+˙p∂H∂p+˙q∂H∂q (38)

and the equations of motion (33a) may be used to remove any total time derivatives that appear at second order. A scheme of this sort would introduce an error in the Jacobian (36) of order .

Symplectic structure may be preserved if the dependence of the Hamiltonian on the future position and momentum at time is taken into account. For example, one may replace the Euler method by the momentum-Verlet method Apostolatos and Ioannou (2004)

 qν+1 = qν+Δt∂H∂pν (39a) pν+1 = pν−Δt∂H∂qν+1. (39b)

or the position-Verlet method

 qν+1 = qν+Δt∂H∂pν+1 (40a) pν+1 = pν−Δt∂H∂qν (40b)

In one of the two equations, the partial derivative of the Hamiltonian is evaluated at as opposed to . For the Hamiltonian (35), Substituting Eq. (39) or (40) Eq. (36) yields a Jacobian identically. Thus, both Verlet methods preserve symplectic structure. In accordance to Liouville’s theorem, phase-space volume is conserved by these methods, as demonstrated in Fig. 3. Although these methods yield closed trajectories in phase space, they need not conserve energy (cf. Fig. 3.). In addition, the resulting scheme will often be a system of nonlinear equations that must be solved numerically. We refer to the method with advanced momentum given by Eq. (40) as VP (momentum-Verlet) and the method with advanced position given by Eq. (39) as VQ (position-Verlet) Gans and Shalloway (2000); Rowley and Marsden (2002); Marsden and West (2001); Kraus (2017).

Using the second order Lanczos-Dyche (LD2) formula (or trapezoidal rule) to expand the solution to the Hamiltonian system (33a) yields the scheme

 qν+1 = qν+Δt2(∂H∂pν+∂H∂pν+1) (41a) pν+1 = pν−Δt2(∂H∂qν+∂H∂qν+1). (41b)

Note that this formula may be obtained by averaging the Verlet schemes (40) and (39). Intuitively, by inspecting the oval-shaped trajectories for VQ and VP in Fig. 3, one may expect the symmetric formula to (41) to yield circular trajectories and thus preserve symplectic structure. For linear systems, such as the harmomic osillator, this is indeed the case. Furthermore, the formula (41) is time-symmetric: the scheme does not change under the reversal . Since energy conservation is intimately related to time symmetry, by means of Noether’s theorem, we expect this scheme to preserve energy. We theoretically and numerically explore these conjectures in the following sections.

If the fourth order LD formula (LD4) is used, one obtains

 qν+1 = qν+Δt2(∂H∂pν+∂H∂pν+1)+Δt212ddt(∂H∂pν−∂H∂pν+1) (42a) pν+1 = pν−Δt2(∂H∂qν+∂H∂qν+1)−Δt212ddt(∂H∂qν−∂H∂qν+1). (42b)

We again expect this formula to preserve energy and symplecticity for the same reasons discussed above.

### 3.3 Linear systems: harmonic oscillator

As a first example of a linear Hamiltonian system, we consider the simple harmonic oscillator. The Hamiltonian for this problem (with an appropriate choice of units) is

 H=p22+q22 (43)

and the Hamilton equations are

 ˙q = p (44a) ˙p = −q (44b)

Since the system is linear, the integration schemes outlined in Section (3.2) can easily be written explicitly. We explicitly construct the first two LD schemes given by Eqs. (41) and (41). LD2 reads

 qν+1=qν+Δt2(pν+pν+1) (45a) pν+1=pν−Δt2(qν+qν+1) (45b)

 qν+1=qν+Δt2(pν+pν+1)−Δt212(qν−qν+1) (46a) pν+1=pν−Δt2(qν+qν+1)−Δt212(pν−pν+1) (46b)

Since the equations are linear, the two schemes may be explicitly solved for the advanced momentum and position . LD2 results in

 qν+1=qν−12Δt21+14Δt2qν+Δt1+14Δt2pν (47a) pν+1=pν−Δt1+14Δt2qν−12Δt2 1+14Δt2pν (47b)

while LD4 results in

 qν+1=qν−12Δt21+112Δt2(1+112Δt2)qν+Δt(1−112Δt2)1+112Δt2(1+112Δt2)pν (48a) pν+1=pν−Δt(1−112Δt2)1+112Δt2(1+112Δt2)qν−12Δt21+112Δt2(1+112Δt2)pν (48b)

In both cases, polynomial expressions in were written in Horner form, to reduce the round-off error from floating point arithmetic. Moreover, separating out the incremental additive change in and explicitly, also reduces round-off error from accumulated summation, and facilitates compensated summation if desired.

The Jacobians of these transformations may be computed in the usual way, and it is found that for both methods. Moreover, it can be shown that

 p2ν+12+q2ν+12=p2ν2+q2ν2 (49)

in both cases. Thus, the LD methods are exactly symplectic and energy conserving for this linear system. Symplecticity implies the numeric phase space trajectories should form closed curves. Energy conserving implies the Hamiltonian should have a constant value on the numeric trajectories; curves of constant energy are circles in this problem. We demonstrate this in Fig. (3). The Euler (RK1) and RK2 methods are neither symplectic nor energy conserving, so their phase portraits exhibit unphysical growth. The growth in RK2 is slower than Euler because the error term is of higher order. The two Verlet methods are both symplectic, so their numeric trajectories are closed curves, but they are not energy conserving, so they are not circles as expected for this problem. The two LD methods, second and fourth order, are both symplectic and energy conserving, so their numeric trajectories are circles in phase space as expected.

The LD methods were shown to be exactly energy conserving, so we numerically confirm this. The LD2 and LD4 methods were compared to RK2 and RK4 (the fourth order Taylor expansion). A step size of was used to evolve the system for 5000 periods. At each time step, the relative error between the numeric energy and exact energy was computed; this quantity is defined by

 δEE≡|E(tν)−E(0)|E(0) (50)

where

 E(tν)≡p2ν2+q2ν2. (51)

The RK methods were found to exhibit polynomial growth over long time scales while the LD methods conserved energy to machine precision. We note that it was necessary to use the Horner forms shown in equations (47) and (48) with compensated summation to achieve the accuracy shown in Fig. (4).

### 3.4 Nonlinear systems: anharmonic oscillator

Although the method was proved to be symplectic for linear systems, such as the harmonic oscillator, it is worth noting that this property does not generalize exactly to non-linear systems: the method is only approximately symplectic for non-linear systems. To see this, consider a general time-independent system:

 ˙p=−V′(q),    ˙q=p (52)

This system is integrated to find

 pν+1=pν−∫tν+1tνV′(q)dt (53)
 qν+1=qν+∫tν+1tνp dt. (54)

Let us apply the second order LD method (or trapezoidal integral rule) to approximate the above integrals by

 pν+1=pν−Δt2(V′(qν)+V′(qν+1)) (55)
 qν+1=qν+Δt2(pν+pν+1) (56)

This system can be rearranged to give equations for the advanced position and momentum :

 qν+1+Δt24V′(qν+1)=qν−Δt24V′(qν)+Δt pν (57)
 pν+1=pν−Δt2(V′(qν)+V′(qν+1)) (58)

The Jacobian for this method can be calculated as

 J=1+Δt24V′′(qν)1+Δt24V′′(qν+1) (59)

We see that the method is not symplectic unless , which is only true if is either a constant, linear, or quadratic function of position.

Nevertheless, we emphasize that the method is still superior to explicit methods. Consider the pendulum equation, with potential , as an example of an anharmonic oscillator. If we assume is small, then Eq. (59) may be expanded as

 J≃1−Δt24(cosqν+1−cosqν)=1+O(Δt3) (60)

Although the order of deviation from symplecticity is formally the same as the order of the method, the method is more accurate for this problem than explicit methods, such as the ordinary Runge-Kutta methods. This is because of two reasons: (i) For the LD methods, the error in the Jacobian is oscillatory and thus has an upper bound, in contrast to the RK methods, where the error grows in time without bound. (ii) As seen by inspecting the remainder terms of each formula, the error for the LD methods is significantly lower compared to RK methods of the same order. The difference increases by several orders of magnitude at high order, as demonstrated in the previous section.

We demonstrate this explicitly by comparing the second and fourth order LD methods to second and fourth order RK methods. To carry out the LD integration, equation (57) and the corresponding fourth order equation needed to be solved numerically. A simple Newton-Raphson algorithm was adequate. For the second order method we iterated the equality

 qi+1=qi−qi+Δt24sinqi−(qν+Δt pν+Δt24sinqν)1+Δt24cosqi (61)

with initialized to . We found that the algorithm converged rapidly, within three iterations, regardless of ’s magnitude. A similar procedure was followed for the fourth order LD method.

We plot the deviation from the true energy as a function of time from a numeric simulation in Fig. (5). We note that the deviation from the true energy is a periodic function of time for both LD methods and remains bounded, as opposed to explicit methods for which the deviation grows polynomially in time. This demonstrates that our LD methods are still superior to explicit methods for the pendulum equation, despite their non-exact symplecticity.

This can be understood by inspecting Eq. (60), which suggests that the deviation from symplecticity is formally of order , but oscillates sinusoidally with and . Each of these contributions averages to zero over one period, so the second order LD method can be considered “symplectic on average”. A similar result holds for the fourth order LD method, since the generalized force is sinusoidal.

### 3.5 Nonconservative Systems: Damped Oscillator

Having demonstrated the method’s utility for both linear and nonlinear conservative systems, we now turn to nonconservative systems. Although there is no longer a manifestly conserved quantity, such systems can still be Hamiltonian and thus preserve symplectic structure. We expect our method to more useful in these problems than an explicit method, such as Runge-Kutta, for this reason.

Consider the damped harmonic oscillator with unit mass and spring constant as a specific example:

 ¨q+γ˙q+q=0. (62)

The damping force cannot be expressed as the derivative of a potential, so this system does not conserve mechanical energy. However, one can still find a Lagrangian that gives rise to this equation Havas (1957); Leach (1977); Lemos (1979); Leubner (1981); Apostolatos and Ioannou (2004); McDonald (2015); Choudhuri et al. (2008):

 L=e−γt(˙q22−q22). (63)

One can then compute a Hamiltonian,

 H=p22e−γt+q22eγt, (64)

and canonical equations of motion,

 ˙p=−eγtq,   ˙q=e−γtp. (65)

Note that is the canonical momentum, distinct from the kinematic momentum in this problem. Although energy is not conserved, one can show that the following quantity is a constant of motion Qureshi et al. (2010):

 C=p22e−γt+γ2pq+q22eγt. (66)

Since this system is Hamiltonian, its symplectic structure should be preserved, and the constant should be conserved through the evolution. This provides a criterion to test numeric schemes. In Figure (6) we plot the relative error in as a function of time resulting from second and fourth order LD and RK methods. Like the pendulum equation, the LD scheme for this problem is not exactly conservative or symplectic, so exhibits oscillations. As mentioned earlier, the oscillations are bounded, and increasing the method’s order reduces the amplitude of these oscillations. Contrast this with the RK methods, in which grows polynomially in time.

The performance of the LD methods may be compared to that of other variational integrators, such as the ‘slimplectic’ integrator Tsang et al. (2015). For dissipative linear systems, such as the damped harmonic oscillator, numerical accuracy is essentially the same. For linear conservative systems, the LD methods conserve both energy and symplectic structure (i.e. the Jacobian 36 of the canonical transformation associated with motion) exactly, while most symplectic integrators typically conserve the Jacobian exactly but the energy is conserved approximately (in the sense that its growth is bound) Marsden and West (2001). On the other hand, for non-linear systems, LD methods conserve energy and the Jacobian approximately, but symplectic methods typically conserve the Jacobian exactly. A discussion of the connection between symmetric and symplectic integration may be found in Hairer et al. (2002).

## 4 Systems of differential equations

### 4.1 Ordinary differential equations: coupled harmonic oscillators

The properties demonstrated above generalize to systems of ordinary of partial differential equations. For example, let us consider a 1-dimensional system of coupled harmonic oscillators. The Hamiltonian can be written in matrix form as

 H(q,p)=12(M−1)abpapb+12Kabqaqb (67)

where and are constant matrices that depend on the mass and Hooke constant of each oscillator, and the indices are summed over the different oscillators. The Hamilton equations of motion,

 ˙qa = (M−1)abpb (68a) ˙pa = −Kabqb, (68b)

can be integrated in time to any order using the LD formula (19), and the equations of motion can be used to eliminate all time derivatives. This procedure, after solving the resulting linear system, yields an explicit scheme

 qν+1 = Aqν+Bpν (69a) pν+1 = Cqν+Dpν, (69b)

where the matrices can be computed using the methods outlined in Sec. 3.3. This example serves as a prelude to the symplectic integration of hyperbolic PDEs such as the wave equation. For example, upon discretization via the method of lines (cf. Sec. 4.2), the Klein-Gordon equation takes the form of a system of coupled ODEs identical to Eq. (68).

### 4.2 Partial differential equations: method of lines

As mentioned above, a common approach to solving systems of evolution equations is the method of lines. For concreteness, let us consider a first-order in time scalar differential equation in 1+1 dimensions

 ∂∂tu(t,x)=L(u(t,x)) (70)

where is a (not necessarily linear) spatial differential operator. If one uses finite-difference or pseudospectral methods for the spatial derivatives Bengt Fornberg (1996); Markakis and Barack (2014), one may rewrite the original partial differential equation as a coupled system of ordinary differential equations of the form

 du(t)dt=L(u(t)) (71)
 dui(t)dt=N∑j=0Lij(uj(t)) (72)

where is a vector constructed from the values of the scalar field in the spatial grid points and is a matrix (band-diagonal in the case of finite differences, or full in the case of pseudo-spectral methods) that couples the different grid points.

Eq. (72) can be time-stepped using the fundamental theorem of calculus, in accordance with Eq. (2). Thus, the methods applied in the previous sections to solve ordinary differential equations, can be generalized to partial differential equations. A detailed derivation and demonstration of the symmetry, symplecticity and stability properties afforded by the LD methods in such PDEs will be given in a companion paper. It will be shown that the method is particularly powerful for linear or quasi-linear PDEs, as it is unconditionally stable, and thus overcomes the Courant-Friedrichs-Lewy limit, and moreover preserves properties of the continuum system, such as energy or probability in the Schrödinger equation and other parabolic or hyperbolic partial differential equations.

## 5 Conclusions

To our knowledge, although Padé approximation methods have been known for linear PDEs, their derivation from the LD formula has not been brought forward before. More importantly, since Padé methods amount to a Padé expansion of the exponential function, their applicability is restricted to linear systems. In contrast, the LD formula only assumes that the problem has continuous derivatives up to some order, so it is equally valid for non-linear systems. Moreover, even in the context of linear systems, the LD methods are more general than Padé methods, as the former can be straightforwardly applied to each sub-equation separately, requiring inversion of matrices of smaller size. This is not straightforward with Padé methods. This property is particularly useful for integrating periodic Hamiltonian systems or wave-like equations (cf. Sec. 4.1), as will be demonstrated in a companion paper.

In the context of ODEs, the LD methods are examples of “multi-derivative methods”, as the schemes depend on derivatives of the function in question. The remarkable stability properties of such methods have been studied by Brown et al. Brown (1975, 1977); Jeltsch and Kratz (1978); Jeltsch (1979). But in addition, their time-symmetric nature gives them great promise for symplectic and energy-conserving integration. Our hope is to draw attention to the properties of these schemes by demonstrating their utility for integrating Hamiltonian systems and many other problems of physical interest. We have shown that the LD methods are exactly symplectic and energy conserving for linear problems and outperform explicit schemes for nonlinear and dissipative problems. Their remarkable accuracy for quadratic Hamiltonians makes them very useful for systems of linear or quasi-linear ODEs.

Moreover, when linear PDEs are discretized using the method of lines, a system of coupled linear ODEs is obtained. Many such problems are endowed with symplectic structure analogous to that outlined in Sec. 4.1. This suggests that an integration scheme with the same advantages as for ODEs (stability, symplecticity, and energy conservation) is often possible. In a companion paper, we will use the LD formula to obtain such schemes for prototypical hyperbolic or parabolic PDEs, such as the advection, wave, and Schrödinger equations. While LD methods certainly posses utility for ODEs, we will demonstrate that it is in the numerical integration of PDEs that these schemes are most powerful.

Finally, we emphasize that, although much of the discussion has focused on linear problems, the LD formula is equally valid nonlinear systems. It will be interesting to explore the use of such integration schemes in nonlinear partial differential equations.

## Acknowledgements

We thank Theocharis Apostolatos, Petros J. Ioannou, Leo Stein and George Pappas for valuable discussions. C.M. was supported by the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 753115. H.S. was supported in part by the Jade Mountain Young Scholar Award No 107V0201.

## References

• Alcubierre (2012) M. Alcubierre, Introduction to 3+1 Numerical Relativity, Oxford University Press, Oxford, 2012.
• Miyatake and Butcher (2016) Y. Miyatake, J. Butcher, A Characterization of Energy-Preserving Methods and the Construction of Parallel Integrators for Hamiltonian Systems, SIAM Journal on Numerical Analysis 54 (2016) 1993–2013.
• Butcher et al. (2016) J. C. Butcher, A. T. Hill, T. J. T. Norton, Symmetric general linear methods, BIT Numerical Mathematics 56 (2016) 1189–1212.
• Chartier et al. (2006) P. Chartier, E. Faou, A. Murua, An Algebraic Approach to Invariant Preserving Integators: The Case of Quadratic and Hamiltonian Invariants, Numerische Mathematik 103 (2006) 575–590.
• Zhong and Marsden (1988) G. Zhong, J. E. Marsden, Lie-Poisson Hamilton-Jacobi theory and Lie-Poisson integrators, Physics Letters A 133 (1988) 134–139.
• Lanczos (1956) C. Lanczos, Applied Analysis, Englewood Cliffs, N.J., Prentice-Hall, 1956.
• Dyche (1956) C. W. Dyche, On multiple point Taylor series expansions, MSc Thesis, Oregon State University, 1956.
• Butcher (2000) J. C. Butcher, Numerical methods for ordinary differential equations in the 20th century, Journal of Computational and Applied Mathematics 125 (2000) 1–29.
• Brown (1977) R. L. Brown, Some Characteristics of Implicit Multistep Multi-Derivative Integration Formulas, SIAM Journal on Numerical Analysis 14 (1977) 982–993.
• Jeltsch and Kratz (1978) R. Jeltsch, L. Kratz, On the stability properties of Brown’s multistep multiderivative methods, Numerische Mathematik 30 (1978) 25–38.
• Markakis and Barack (2014) C. Markakis, L. Barack, High-order difference and pseudospectral methods for discontinuous problems, arXiv:1406.4865 [gr-qc] (2014).
• Nozawa et al. (1998) T. Nozawa, N. Stergioulas, E. Gourgoulhon, Y. Eriguchi, Construction of Highly Accurate Models of Rotating Neutron Stars - Comparison of Three Different Numerical Schemes, Astronomy and Astrophysics Supplement Series 132 (1998) 431–454.
• Lampret (2004) V. Lampret, An Invitation to Hermite’s Integration and Summation: A Comparison between Hermite’s and Simpson’s Rules, Society for Industrial and Applied Mathematics 46 (2004) 311–328.
• Lotkin (1952) M. Lotkin, A New Integrating Procedure of High Accuracy, Journal of Mathematics and Physics 31 (1952) 29–34.
• Hummel and Seebeck (1949) P. M. Hummel, C. L. Seebeck, A Generalization of Taylor’s Expansion, The American Mathematical Monthly 56 (1949) 243–247.
• Estes (1966) R. H. L. Estes, Two-Point Taylor Series Expansions, Technical Report, 1966.
• Squire (1961) W. Squire, Some Applications of Quadrature by Differentiation, Journal of the Society for Industrial and Applied Mathematics 9 (1961) 94–108.
• Squire (1962) W. Squire, A Comment on the Approximate Calculation of Eigenvalues for Certain Second-Order Linear Differential Equations, The Aeronautical Journal 66 (1962) 722–722.
• Gould and Squire (1963) H. W. Gould, W. Squire, Maclaurin’s Second Formula and Its Generalization, The American Mathematical Monthly 70 (1963) 44–52.
• Lambert and Mitchell (1963) J. D. Lambert, A. R. Mitchell, The use of higher derivatives in quadrature formulae, The Computer Journal 5 (1963) 322–327.
• Epperson (2013) J. F. Epperson, An Introduction to Numerical Methods and Analysis, 2nd edition ed., Wiley, 2013.
• Tao (2014) T. Tao, Compactness and Contradiction, AMS, University of California, Los Angeles, CA, 2014.
• Tao (2010) T. Tao, The Euler-Maclaurin formula, Bernoulli numbers, the zeta function, and real-variable analytic continuation, https://terrytao.wordpress.com/2010/04/10/, 2010.
• Hawking (1977) S. W. Hawking, Zeta function regularization of path integrals in curved spacetime, Communications in Mathematical Physics 55 (1977) 133–148.
• Elizalde and Romeo (1989) E. Elizalde, A. Romeo, Expressions for the zeta-function regularized Casimir energy, Journal of Mathematical Physics 30 (1989) 1133–1139.
• Fermi and Pizzocchero (2011) D. Fermi, L. Pizzocchero, Local zeta regularization and the Casimir effect, Progress of Theoretical Physics 126 (2011) 419–434.
• Parker and Toms (2009) L. Parker, D. Toms, Quantum Field Theory in Curved Spacetime by Leonard Parker, 2009. doi:10.1017/CBO9780511813924.
• Ehle (1973) B. Ehle, A-Stable Methods and Padé Approximations to the Exponential, SIAM Journal on Mathematical Analysis 4 (1973) 671–680.
• Dahlquist (1956) G. Dahlquist, Convergence and stability in the numerical integration of ordinary differential equations, MATHEMATICA SCANDINAVICA 4 (1956) 33–53.
• Dahlquist (1961) G. Dahlquist, Stability and Error Bounds in the Numerical Integration of Ordinary Differential Equations, Zeitschrift Angewandte Mathematik und Mechanik 41 (1961) 267–268.
• Dahlquist (1963) G. G. Dahlquist, A special stability problem for linear multistep methods, BIT Numerical Mathematics 3 (1963) 27–43.
• Wanner (2006) G. Wanner, Dahlquist’s classical papers on stability theory 46 (2006) 671–683.
• Brown (1975) R. L. Brown, Jr., Multi-Derivative Numerical Methods for the Solution of Stiff Ordinary Differential Equations., PhD Thesis, University of Illinois at Urbana-Champaign, Champaign, IL, USA, 1975.
• Jeltsch (1976) R. Jeltsch, Note on A-stability of multistep multiderivative methods, BIT Numerical Mathematics 16 (1976) 74–78.
• Jeltsch (1978) R. Jeltsch, On the stability regions of multistep multiderivative methods, in: R. Bulirsch, R. D. Grigorieff, J. Schröder (Eds.), Numerical Treatment of Differential Equations, Lecture Notes in Mathematics, Springer Berlin Heidelberg, 1978, pp. 63–80.
• Jeltsch (1979) R. Jeltsch, A0-Stability and stiff stability of Brown’s multistep multiderivative methods, Numerische Mathematik 32 (1979) 167–181.
• Gekeler (2006) E. Gekeler, Discretization Methods for Stable Initial Value Problems, Springer, 2006.
• Apostolatos and Ioannou (2004) T. Apostolatos, P. Ioannou, Elements of Theoretical Mechanics, Leader Books, 2004.
• Gans and Shalloway (2000) J. Gans, D. Shalloway, Shadow mass and the relationship between velocity and momentum in symplectic numerical integration, Physical Review E 61 (2000) 4587–4592.
• Rowley and Marsden (2002) C. W. Rowley, J. E. Marsden, Variational integrators for degenerate Lagrangians, with application to point vortices, in: Proceedings of the 41st IEEE Conference on Decision and Control, 2002., volume 2, 2002, pp. 1521–1527 vol.2. doi:10.1109/CDC.2002.1184735.
• Marsden and West (2001) J. E. Marsden, M. West, Discrete mechanics and variational integrators, Acta Numerica 10 (2001) 357–514.
• Kraus (2017) M. Kraus, Projected Variational Integrators for Degenerate Lagrangian Systems, arXiv:1708.07356 [math-ph, physics:physics] (2017).
• Havas (1957) P. Havas, The range of application of the lagrange formalism — I, Il Nuovo Cimento 5 (1957) 363–388.
• Leach (1977) P. G. L. Leach, On a direct method for the determination of an exact invariant for the time-dependent harmonic oscillator, The ANZIAM Journal 20 (1977) 97–105.
• Lemos (1979) N. A. Lemos, Canonical approach to the damped harmonic oscillator, American Journal of Physics 47 (1979) 857–858.
• Leubner (1981) C. Leubner, Inequivalent lagrangians from constants of the motion, Physics Letters A 86 (1981) 68–70.
• McDonald (2015) K. McDonald, A Damped Oscillator as a Hamiltonian System, 2015.
• Choudhuri et al. (2008) A. Choudhuri, S. Ghosh, B. Talukdar, Symmetries and conservation laws of the damped harmonic oscillator, Pramana 70 (2008) 657.
• Qureshi et al. (2010) M. I. Qureshi, M. Rafat, S. I. Azad, The exact equation of motion of a simple pendulum of arbitrary amplitude: A hypergeometric approach, European Journal of Physics 31 (2010) 1485.
• Tsang et al. (2015) D. Tsang, C. R. Galley, L. C. Stein, A. Turner, Slimplectic Integrators: Variational Integrators for General Nonconservative Systems, The Astrophysical Journal Letters 809 (2015) L9.
• Hairer et al. (2002) E. Hairer, G. Wanner, C. Lubich, Symplectic Integration of Hamiltonian Systems, Springer Series in Computational Mathematics, Springer, Berlin, Heidelberg, 2002.
• Bengt Fornberg (1996) Bengt Fornberg, A Practical Guide to Pseudospectral Methods, Cambridge Monographs on Applied and Computational Mathematics, 1996.