
The numerical solution of an ordinary differential equation can be interpreted as the exact solution of a nearby modified equation. Investigating the behaviour of numerical solutions by analysing the modified equation is known as backward error analysis. If the original and modified equation share structural properties, then the exact and approximate solution share geometric features such as the existence of conserved quantities. Conjugate symplectic methods preserve a modified symplectic form and a modified Hamiltonian when applied to a Hamiltonian system. We show how a blended version of variational and symplectic techniques can be used to compute modified symplectic and Hamiltonian structures. In contrast to other approaches, our backward error analysis method does not rely on an ansatz but computes the structures systematically, provided that a variational formulation of the method is known. The technique is illustrated on the example of symmetric linear multistep methods with matrix coefficients.
Citation: Robert I McLachlan, Christian Offen. Backward error analysis for conjugate symplectic methods[J]. Journal of Geometric Mechanics, 2023, 15(1): 98-115. doi: 10.3934/jgm.2023005
[1] | Valentin Duruisseaux, Melvin Leok . Time-adaptive Lagrangian variational integrators for accelerated optimization. Journal of Geometric Mechanics, 2023, 15(1): 224-255. doi: 10.3934/jgm.2023010 |
[2] | Elias Maciel, Inocencio Ortiz, Christian E. Schaerer . A Herglotz-based integrator for nonholonomic mechanical systems. Journal of Geometric Mechanics, 2023, 15(1): 287-318. doi: 10.3934/jgm.2023012 |
[3] | Álvaro Rodríguez Abella, Melvin Leok . Discrete Dirac reduction of implicit Lagrangian systems with abelian symmetry groups. Journal of Geometric Mechanics, 2023, 15(1): 319-356. doi: 10.3934/jgm.2023013 |
[4] | Maulik Bhatt, Amit K. Sanyal, Srikant Sukumar . Asymptotically stable optimal multi-rate rigid body attitude estimation based on lagrange-d'alembert principle. Journal of Geometric Mechanics, 2023, 15(1): 73-97. doi: 10.3934/jgm.2023004 |
[5] | Jordi Gaset, Arnau Mas . A variational derivation of the field equations of an action-dependent Einstein-Hilbert Lagrangian. Journal of Geometric Mechanics, 2023, 15(1): 357-374. doi: 10.3934/jgm.2023014 |
[6] | Marcin Zając . The dressing field method in gauge theories - geometric approach. Journal of Geometric Mechanics, 2023, 15(1): 128-146. doi: 10.3934/jgm.2023007 |
[7] | Jacob R. Goodman . Local minimizers for variational obstacle avoidance on Riemannian manifolds. Journal of Geometric Mechanics, 2023, 15(1): 59-72. doi: 10.3934/jgm.2023003 |
[8] | Francesco Bonechi, Jian Qiu, Marco Tarlini . Generalised Kähler structure on $ \mathbb{C}P^2 $ and elliptic functions. Journal of Geometric Mechanics, 2023, 15(1): 188-223. doi: 10.3934/jgm.2023009 |
[9] | William Clark, Anthony Bloch . Existence of invariant volumes in nonholonomic systems subject to nonlinear constraints. Journal of Geometric Mechanics, 2023, 15(1): 256-286. doi: 10.3934/jgm.2023011 |
[10] | Xavier Rivas, Daniel Torres . Lagrangian–Hamiltonian formalism for cocontact systems. Journal of Geometric Mechanics, 2023, 15(1): 1-26. doi: 10.3934/jgm.2023001 |
The numerical solution of an ordinary differential equation can be interpreted as the exact solution of a nearby modified equation. Investigating the behaviour of numerical solutions by analysing the modified equation is known as backward error analysis. If the original and modified equation share structural properties, then the exact and approximate solution share geometric features such as the existence of conserved quantities. Conjugate symplectic methods preserve a modified symplectic form and a modified Hamiltonian when applied to a Hamiltonian system. We show how a blended version of variational and symplectic techniques can be used to compute modified symplectic and Hamiltonian structures. In contrast to other approaches, our backward error analysis method does not rely on an ansatz but computes the structures systematically, provided that a variational formulation of the method is known. The technique is illustrated on the example of symmetric linear multistep methods with matrix coefficients.
While the forward error of a numerical method compares the exact solution of an ODE with the numerical solution after one time-step h, to obtain qualitative statements about the long-term behaviour of numerical solutions to ODEs, it is helpful to consider a modified ODE whose exact solution closely approximates the numerical flow map at grid points. A modified equation can be obtained as an expansion of the numerical solution as a power series in the step-size h. Though the series does not converge in general, optimal truncation techniques have been established such that the flow of the modified system and the numerical method coincide at grid points up to exponentially small error terms. Computing and analysing the structural properties of modified equations is known as backward error analysis (BEA) (see, for instance, [3,§Ⅸ] or [6,§5]). Next to the analysis of long-term behaviour of numerical schemes, backward error analysis has been used to improve the initialisation of multi-step methods [2] as well as to improve physics informed machine learning techniques [13,12,10].
If a Hamiltonian ODE is discretised by a symplectic integrator, then any truncation of the modified equation is itself a Hamiltonian system with respect to the original symplectic structure and a modified Hamiltonian. These are also called Shadow Hamiltonians. The existence of a modified Hamiltonian or a modified Lagrangian is a key ingredient to obtain statements about long-term behaviour of symplectic method, such as oscillatory energy errors over exponentially long time intervals. Moreover, just as in the exact system, symplectic symmetries of the modified system yield conserved quantities for the modified dynamics by Noether's theorem. This explains why symplectic integrators behave well on completely integrable systems. A detailed discussion can be found in [3]. Vermeeren observed that backward error analysis for variational integrators can be done entirely on the Lagrangian side [15].
In contrast to symplectic methods, conjugate symplectic methods preserve a modified symplectic structure rather than the original symplectic structure. Conjugate symplectic methods share the excellent long-term behaviour of symplectic methods. Moreover, Noether's theorem applies such that symmetries of the modified system yield modified conserved quantities of the modified dynamics. When modified structures are explicitly known, explicit expressions of modified conserved quantities can be derived. This motivates the development of techniques to compute modified symplectic structures and Hamiltonians.
While traditional methods for the computation of modified Hamiltonians use an Ansatz (i.e. an educated guess) of the Hamiltonian as a power series and match terms, working with an Ansatz is challenging when a modified Hamiltonian and a modified symplectic structure need to be computed simultaneously: the components of matrices representing symplectic structures are in this context not constant but depend on the state space variables, fulfil a symmetry condition, and satisfy the Jacobi identity, which is a partial differential equation [3,§Ⅶ.2]. This makes finding a suitable Ansatz difficult.
A typical strategy [3,8] to obtain a structure preserving numerical method is to approximate the variational principle
δS=0,S(y)=∫tNt0L(y(t),˙y(t))dt,y(t0)=y0,y(tN)=yN | (1.1) |
which governs the exact Euler–Lagrange equations
ddt∂L∂˙y−∂L∂y=0 |
by a discrete variational principle
∇SΔ({yi}i)=0,SΔ({yi}i)=N−1∑i=0hLΔ(yi,yi+1). |
Since y0 and yN are fixed in the variations considered in (1.1), the gradient above is taken with respect to all inner grid points y1,…,yN−1. We obtain the discrete Euler–Lagrange equations
D1LΔ(yi,yi+1)+D2LΔ(yi−1,yi)=0,i=0,…,N−1 |
which yield approximations yi≈y(t0+ih) to an exact solution y. The term recursion is called a variational method. Indeed, the class of variational methods is equivalent to the class of symplectic integrators [3,8].
While backward error analysis for discrete Lagrangians LΔ(yi,yi+1) are established, discrete Lagrangians depending on several grid-points LΔ(yi,yi+1,…,yi+s) corresponding to multistep methods require different approaches because they are not symplectic but only preserve a modified symplectic structure. In other words, these methods are conjugate to symplectic methods. However, the modified symplectic structures or conjugacies, respectively, are given by a formal power series that might not be convergent. Although rigorous optimal truncation results are not available, we will demonstrate in numerical examples that truncations can be useful objects in the analysis of the numerical methods.
In the following, we will introduce blended backward error analysis to systematically compute modified Hamiltonian and modified symplectic structures. We will prove the following theorem which applies, for instance, to series expansions LΔ of consistent discrete Lagrangians LΔ(y(t),y(t+h),…,y(t+sh)).
Theorem 1.1. Consider a power series LΔ(y[∞]) in a formal variable h. The series depends on the jet y[∞]=(y,˙y,¨y,…) of a variable y such that any truncation only depends on a finite jet of y. Assume further that the truncation to zeroth order constitutes a regular Lagrangian L(y,˙y), i.e. (∂2L∂˙yi∂˙yj)i,j is invertible.
● There exists a 2nd order modified equation given as a formal power series in h such that for any N∈N a solution of the modified equation truncated to order O(hN) solves the Euler–Lagrange equations to L[N]Δ up to an error of order O(hN+1), where L[N]Δ is the truncation of LΔ to order O(hN).
● If we denote z=(y,˙y), there exists a symplectic structure matrix Jmod and a Hamiltonian Hmod given as formal power series in h such that solutions to Hamilton's equations
˙z=J[N]mod(z)−1∇H[N]mod(z) |
fulfil the modified equation. Here, J[N]mod and H[N]mod are truncations to order O(hN) of Jmod and Hmod, respectively such that L[N]Δ is regular, i.e. (∂2L[N]Δ∂y(M) i∂y(M) j)i,j is invertible, where y(M) is the highest derivative of L[N]Δ.
The technique will be illustrated for linear multistep methods with matrix valued coefficients (2). These occur, for instance, when in a system of coupled ODEs components of the differential equation are discretised separately with traditional linear multistep method, when multistep methods are stabilised [4,5,14], or when discretisation schemes of PDEs are analysed [9].
Moreover, we will analyse under which conditions modified Lagrangians exist: if the original equation is the Euler–Lagrange equation to a variational principle of the form (1.1) for a Lagrangian L(y,˙y), it is natural to ask, whether there exists a modified Lagrangian Lmod(y,˙y) such that the modified equation is governed by the Euler–Lagrange equations to Lmod. In contrast to classical variational integrators, for which Lmod(y,˙y) is known to exist and can be computed [15], for conjugate symplectic schemes Lmod may only exists in modified variables Lmod(˜y,˙˜y). We will prove that Lmod exists in the original variables (y,˙y) if Jmod has the form
Jmod=(∗∗∗0). |
In particular, this shows that for classical consistent, symmetric, linear multistep methods with matrix coefficients with a central force evaluation applied to second order equations Lmod exists in the original variables (y,˙y) if all coefficients are multiples of the identity matrix, i.e. we have a multistep method with scalar coefficients. However, in the general case of matrix coefficients Lmod only exists in modified variables (˜y,˙˜y).
The article is structured as follows: Section 2 illustrates the ideas of blended backward error analysis on linear multi-step methods. Section 3 shows how to compute the modified data Jmod and Hmod introduced in 1.1. The technique is then applied to linear multi-step methods with matrix valued coefficients in 4 and results are illustrated by numerical experiments. Additionally, for comparison of blended backward error analysis with classical backward error analysis, ?? contains an application of blended backward error analysis to a mechanical ordinary differential equation discretised with the Störmer-Verlet scheme. A formal proof of 1.1 is provided in 5. Finally, 6 discusses the existence of modified Lagrangians as formal power series and future research directions are suggested in 7.
To illustrate the idea of blended backward error analysis, we compute modified symplectic structures and Hamiltonians of linear multistep methods.
Consistent, symmetric linear multistep methods with a single force evaluation applied to the second order ODE
¨y=∇U(y(t)) | (2.1) |
take the form
s2∑j=1Aj(y(t−jh)−2y(t)+y(t+jh))=h2∇U(y(t)) | (2.2) |
with
s2∑j=1j2Aj=I. | (2.3) |
Relation (2.3) is coming from the consistency requirement [4]. These are s-step methods, where s is even. Here we allow matrix valued coefficients Aj [4,5,9,14], h is the step-size and I denotes the identity matrix. If the coefficients Aj are scalars, then the schemes constitute classical consistent, symmetric linear multistep methods. A series expansion of (2.2) in h is equivalent to a power series expression of the form
¨y=∇U(y(t))+N∑i=1hi˜gi(y(t),…,y(ai))+O(hN+1) | (2.4) |
for some h-independent expressions ˜gi in the ai-jet of y at t. Substituting 3rd and higher derivatives on the right hand side of (2.4) with derivatives of (2.4) itself iteratively yields an equation of the form
¨y=∇U(y(t))+N∑i=1higi(y(t),˙y(t))+O(hN+1) | (2.5) |
which is called the modified equation of method (2.2) applied to (2.1). We refer to [3] for optimal truncation techniques and a discussion of spurious solutions not covered by the considered modified system for the case of linear multistep methods with scalar coefficients. In the following, we will focus on the question which structural properties the modified equation (2.5) shares with the original ODE (2.1).
Variational Structure
The ODE (2.1) has first order variational structure as it is the Euler–Lagrange equation
ddt∂L∂˙y−∂L∂y=0 |
to the variational principle
δS=0forS(y)=∫L(y(t),˙y(t))dt |
with
L(y,˙y)=12‖˙y‖2+U(y). |
Moreover, there exists a variational principle for (2.2):
Lemma 2.1. Let the matrices Aj∈Rn×n be symmetric. For T>0 let T either be the circle T=R/TZ or the real line T=R. For y defined on T the variational principle
δSΔ=0forSΔ(y)=∫TLΔ(y(t),y(t+h),…,y(t+s2h))dt | (2.6) |
with
LΔ(y(t),y(t+h),…,y(t+s2h))=12h2(s2∑j=1⟨Aj(y(t+jh)−y(t)),y(t+jh)−y(t)⟩)+U(y(t)). |
implies the functional equation (2.2). Moreover, if T=[a,b] is an interval, then (2.6) implies (2.2) on the interval ∘T=[a+s2h,b−s2h]. Here we assume that the function space for y and the potential U are such that U∘y and ∇U∘y constitute square integrable functions.
Proof. Let Δτy denote the forward difference, i.e. (Δτy)(t)=y(t+τ)−y(t). Then ⟨y,Δτz⟩L2(T,Rn)=⟨Δ∗τy,z⟩L2(T,Rn) holds with Δ∗τ=Δ−τ on T if T∈{R,R/TZ} or on the interval [a+τ,b−τ] if T=[a,b]. The expression Δ∗τΔτy corresponds to the central difference
(Δ∗τΔτy)(t)=−y(t+τ)+2y(t)−y(t−τ). |
Let δ denote the variational derivative in the direction of a variation δy, i.e. δSΔ(y)=limϵ→01ϵ(SΔ(y+ϵδy)−SΔ(y)). Let δy∈C∞c(T,Rn), if T∈{R,R/TZ} and let δy∈C∞c(∘T,Rn), if T=[a,b]. We have
δSΔ(y)=12h2s2∑j=1δ⟨AjΔjhy,Δjhy⟩L2(T,Rn)+δ∫TU(y(t))dt=1h2s2∑j=1⟨AjΔjhy,Δjhδy⟩L2(T,Rn)+⟨∇U(y),δy⟩L2(T,Rn)=1h2s2∑j=1⟨AjΔ∗jhΔjhy+∇U(y),δy⟩L2(T,Rn). |
Now (2.2) follows from the fundamental lemma of the calculus of variations on T or ∘T, respectively.
To analyse structure preserving properties of the method (2.2), it might seem natural to seek a modified Lagrangian Lmod(y,˙y) given as a formal power series in the step-size h such that the modified variational principle
δSmod=0forSmod(y)=∫Lmod(y(t),˙y(t))dt |
covers smooth solutions of (2.6) up to any order in the step-size h. However, we show that although a 1st order Lagrangian Lmod covering the modified equations always exists as a power series, it only exist in modified variables (˜y,˙˜y) in the most general case. Even for simple methods, the existence of an expression in closed form for the change of coordinates from (y,˙y) to (˜y,˙˜y) can not be expected. This makes it difficult to compute Lmod using an ansatz.
Hamiltonian structure
Another approach is to work on the Hamiltonian side. The ODE (2.1) has the form of a Hamiltonian system
˙z(t)=J−1∇H(z(t)),H(z)=12‖˙y‖2−U(y),z=(y˙y). | (2.7) |
Here
J=(0−II0) | (2.8) |
is the standard symplectic structure.
In this paper we use a blended approach of the variational and Hamiltonian viewpoint to systematically compute a modified Hamiltonian system
˙z(t)=J−1mod(z(t))∇Hmod(z(t)) | (2.9) |
consisting of a modified Hamiltonian Hmod and a modified symplectic structure Jmod given as formal power series in h such that for a truncation to arbitrary order N the ODE (2.9) covers the modified equation (2.5) up to higher order terms. Here Jmod is a skew symmetric matrix which satisfies a Jacobi identity. The following theorem can be considered as an instance of 1.1 and will be proved in 5.
Theorem 2.2. Let N∈N denote the considered order of the series expansion of the matrix multistep method (2.2), where the coefficients are symmetric matrices. There exists a modified symplectic structure J[N]mod, which is O(h) close to J and a modified Hamiltonian H[N]mod, which is O(h) close to H, such that
˙z(t)=(J[N]mod(z))−1∇H[N]mod(z(t)) | (2.10) |
with coordinate z=(y,˙y) is equivalent to the modified equation (2.5) up to terms of order O(hN+1).
Here O(h)-closeness of J[N]mod and J means that the zeroth coefficient of the polynomial J[N]mod in the formal variable h is given by J. O(h)-closeness of H[N]mod and H has an analogous meaning.
We observe conditions under which a modified first order Lagrangian Lmod(y,˙y) exists in the original variable y by analysing the modified symplectic structure Jmod.
Theorem 2.3. If the matrices Aj in (2.2) are scalar multiples of the identity matrix, then there exists a modified Lagrangian L[N]mod depending on (y,˙y) that is O(h)-close to L(y,˙y) such that the modified variational principle
δS[N]mod=0forS[N]mod(y)=∫L[N]mod(y(t),˙y(t))dt |
yields the modified equation (2.5) up to terms of order O(hN+1).
Again, O(h)-closeness is to be interpreted in a formal sense, analogously to its meaning in 2.2.
Theorem 2.3 applies to traditional multistep methods with scalar coefficients. However, we will see that for general linear multistep methods with matrix-valued coefficients the existence of a first order modified Lagrangian in the original variable y cannot be expected. A proof of 2.3 is postponed to 6.
In this section we introduce a method to compute the modified data J[N]mod and H[N]mod of 2.2 such that (2.10) governs (2.5). We will then verify the validity of the construction method and prove 1.1 and 2.2 in 5.
Let L[N]Δ denote the series expansion of
LΔ(y(t),y(t+h),…,y(t+s2h)) |
to order N in the step-size h. The expression L[N]Δ depends on the M-jet of y at t for some M∈N. Notice that the order M variational principle
δS[N]Δ(y)=0withS[N]Δ(y)=∫L[N]Δ(y(t),˙y(t),…,y(M)(t))dt | (3.1) |
recovers (2.4) up to higher order terms. We first compute a high-dimensional Hamiltonian system defined on the 2M−1-jet space of y corresponding to the order M variational principle (3.1). The Hamiltonian principle is then reduced to a Hamiltonian system defined on the 1-jet space of y. It has the form (2.9) and covers (2.5) up to higher order terms.
To construct the high-dimensional Hamiltonian system, we use Ostrogradsky's Hamiltonian description of high-order Lagrangian systems [16]. For this we define variables
q=(y,˙y,…,y(M−1)) |
and for i=1,…,M
pji=M−i∑k=0(−1)kdkdtk∂L[N]Δ∂(yj)(k+i). |
Here the index j enumerates the components of y and ddt denotes the total derivative operator on the jet-space of y, which acts like
ddtρ(y,˙y,…,y(M))=M∑i=0⟨∇y(i)ρ,y(i+1)⟩. |
on a scalar valued function ρ defined on the M-jet space of y. The high dimensional Hamiltonian system consists of the Hamiltonian
H[N]=M∑i=1⟨pi,˙qi⟩−L[N]Δ, | (3.2) |
where all expressions are expressed in y[2M−1]=(y,˙y,…,y(2M−1)), and the symplectic structure matrix J[N]mod(y[2M−1]). The skew-symmetric matrix J[N]mod(y[2M−1]) is the representing matrix of the differential 2-form
Ω[N]=M∑i=1n∑j=1dpji∧dqij, | (3.3) |
where pji and qij are interpreted as functions in the variable y[2M−1] of the 2M−1-jet space, i.e. J[N]mod is the anti-symmetrised tensor product* ∧ of the gradients ∇y[2M−1] pij and ∇y[2M−1] qij summed over all indices.
*This corresponds to the command TensorWedge in Wolfram Mathematica.
To compute the modified Hamiltonian system on the 1-jet space with variable y[1]=(y,˙y), the variables y(2),…,y(2M−1) in the expression (3.2) for the Hamiltonian H[N](y[2M−1]) are repeatedly replaced by (2.5) until higher derivatives only occur in O(hN+1) terms. This yields H[N]mod(y,˙y). Similarly, we can consider pji and qij as functions of (y,˙y) truncating O(hN+1) terms. The matrix J[N]mod is then given as the representing matrix of the 2-form Ω[N] pulled to the 1-jet space of the variable y, i.e. interpreting y[1]=(y,˙y) as the only independent variables. Equivalently, J[N]mod is the anti-symmetrised tensor product ∧ of the gradients ∇y[1]pij and ∇y[1]qij summed over all indices. (As J[N]mod is constructed from a closed differential 2-form, it is automatically skew-symmetric and satisfies the Jacobi identity.) This completes the construction of the modified data.
The system (2.9) recovers (2.5) up to higher order terms as we will prove after a computational example.
As introduced in 2, consider the multistep method
A2y(t−2h)+A1y(t−h)−2(A1+A2)y(t)+A1y(t+h)+A2y(t+2h)=h2∇U(y(t)) | (4.1) |
with matrix coefficients in dimension n=2. By the consistency requirement, A2=1/4(I−A1). We obtain
J[4]mod=J+h2(J211−J221J2210)+h4(J411−J421J421J422) |
with
J211=(0−b1b10) |
where
b1=14(˙y2(a12(U(2,1)−U(0,3))+(a22−a11)U(1,2))+14(˙y1(−a12(U(2,1)−U(0,3))+(a22−a11)U(2,1))). |
Here and in the following U(k,l)=∂k+lU(y)∂ky1∂ly2.
As the expressions of higher order terms become quite complicated, we refer the reader to the Mathematica Notebooks of our accompanying source code [11]. However, as this will be relevant in the discussion later, we are reporting J422 for the special case that A1 is a diagonal matrix: we have
J422=(0−b2b20) |
with
b2=18(a211−a222+2(a22−a11))(U(2,1)˙y1+U(1,2)˙y2)ifA1=(a1100a22). |
The modified Hamiltonian H[4]mod for the general case is given as
H[4]mod(y,˙y)=H(y,˙y)+h2H2(y,˙y)+h4H4(y,˙y) |
with
H2(y,˙y)=124((˙y1)2(2(4−3a11)U(2,0)−6a12U(1,1))−2˙y2˙y1(3a12U(0,2)+(3a11+3a22−8)U(1,1)+3a12U(2,0))−6a22U(0,2)(˙y2)2−6a12U(1,1)(˙y2)2+(3a22−4)(U(0,1))2+3a11(U(1,0))2+6a12U(0,1)U(1,0)+8U(0,2)(˙y2)2−4(U(1,0))2). |
For further terms, we refer the reader to the Mathematica Notebooks of our accompanying source code. Hamilton's equations
˙z=(J[4]mod(z))−1∇H[4]mod(z) |
for the modified Hamiltonian system are equivalent to the modified equation (2.5) truncating terms of order O(h6).
Figure 1 shows a numerical experiment with a rotational invariant potential. The start values for the multistep formula were obtained using the fourth order modified equation. Trajectories computed with the multistep scheme look very regular. The quantities H[0]mod=H, H[2]mod, and H[4]mod evaluated along a trajectory show oscillatory energy error behaviour. Experiments with different values for the step-size h confirm the preservation of H[k]mod, k=0,2,4 up to truncation error. Initialising the multistep scheme with the fourth order modified equation, the effects of spurious solutions was minimised. However, as h is decreased, spurious solutions cause wriggles in the energy error H[k]mod−H(zinit), k=0,2,4 which eventually prevent further energy error decay.
If A2=αI with α≥14 and A1=I−4A2, then (4.1) corresponds to a classical stable† multistep scheme: the generating polynomial ρ to (4.1) is given as
†A multistep scheme for second order equations is stable if all roots of its generating polynomial lie in the closed unit disk and those on the circle are at most double zeros [3,XV.1.2].
ρ(ξ)=α+(1−4α)ξ−2(1−3α)ξ2+(1−4α)ξ3+αξ4. |
The polynomial ρ has a double root at 1 as well as the roots
ξ3=1−12α+√1−4α2α,ξ4=1−12α−√1−4α2α. |
Since α≥14, the roots ξ3 and ξ4 are complex conjugate to each other and lie on the unit circle such that the scheme is stable.
Moreover, LΔ and LΔ are rotationally invariant because A1 and A2 commute with rotation matrices. An application of Noether's theorem to L[N]Δ yields the following modified angular momentum:
I[N]=M∑m=1m−1∑k=0(−1)k⟨∇ykL[N]Δ,dRy(m−1−k)⟩,dR=(0−110) |
The integer M is the order of the highest derivative of y in L[N]Δ. For the truncation order N=4 we have M=5. Using (2.5) repeatedly, derivatives of y of order greater than two are replaced by terms in y, ˙y and O(h6)-terms. After truncation at order 4 in h we obtain the modified angular momentum I[4]mod(y,˙y). Figure 2 shows the conservation error of I=I[0]mod, I[2]mod, I[4]mod along a trajectory. We see oscillatory error behaviour. The oscillations shrink as the step-size h decreases until spurious solutions prevent further decrease.
For source code of our experiments refer to [11].
Let us prepare the proof of 1.1. We will use the language of Differential Geometry.
To motivate the following proposition and to fix notation recall the following fact: if (Z′,ω′,H′) is a Hamiltonian system with Hamiltonian vector field X′H′, ZΨ↪Z′ an embedding such that Ψ(Z) is a symplectic submanifold that is invariant under motions, i.e. X′H′(Ψ(z))∈TΨ(z)Ψ(Z) for all z∈Z, then the pull-back system (Z,ω,H) with ω=Ψ∗ω′, H=H′∘Ψ is a Hamiltonian system with Hamiltonian vector field XH. Here Ψ∗ω′ denotes the pull-back of ω′ along Ψ. The Hamiltonian vector fields XH and X′H′ relate by
Ψ∗XH=X′H′∘Ψ, |
i.e. for a motion γ with ˙γ=XH∘γ the curve γ′=Ψ∘γ is a motion of (Z′,ω′,H′), i.e. ˙γ′=X′H′∘γ′. Here Ψ∗ denotes the push-forward map, i.e. (Ψ∗XH)(z)=dΨ|z(XH(z)) for z∈Z.
In the following, we will adapt this statement to a setting, where the definition of ω′, H′, and Ψ′ contain a formal variable h and where Ψ(Z) is left invariant only up to higher order terms in h.
Definition 5.1 (Formal Hamiltonian system). Let N∈N be a truncation index, Z be a smooth manifold, h a formal variable, and H=∑Nk=0hkHk a formal polynomial whose coefficients are smooth maps Z→R. Further, let ω=∑Nk=0hkωk be a formal polynomial whose coefficients are 2-forms on Z. The collection (Z,ω,H) is called a formal Hamiltonian system with truncation index N if the following conditions are satisfied.
● The formal symplectic form ω is closed, i.e. dω=∑Nk=0hkdωk is a formal polynomial whose coefficients are 3-forms that are zero.
● The 2-form ω0 is non-degenerate.
The relation −dH=ω(X,⋅) defines a formal power series X=∑∞k=0hkXk, where Xk are vector fields on Z. The truncation XH:=∑Nk=0hkXk is called (formal) Hamiltonian vector field. The formal differential equation ˙γ=XH∘γ is called Hamilton's equation.
Proposition 5.2. Let N∈N be a truncation index, Z′ and Z be manifolds, and let (Z′,ω′,H′) be a formal Hamiltonian system with truncation index N. Consider a formal polynomial Ψ=∑Nk=0Ψkhk such that Ψk:Z→Z′ are smooth and Ψ0:Z→Z′ is an embedding. To z∈Z define formal tangent spaces
TΨ(z)Ψ(Z):={N∑k=0hkdΨk|z(v)|v∈TzZ}⊂N⨁k=0hkTz′kZ′, |
where z′:=∑Nk=0hkz′k:=Ψ(z).
● Assume that the Hamiltonian vector field is tangential to Ψ(Z), i.e. trk(X′H′(Ψ(z)))∈TΨ(z)Ψ(Z) for all z∈Z. Here trk denotes the truncation of O(hN+1)-terms.
● Assume that Ψ(Z) is a symplectic submanifold of Z′, i.e. for all z∈Z
TΨ(z)Ψ(Z)∩TΨ(z)Ψ(Z)⊥ω′=0. |
Here Ψ(Z)⊥ω′ denotes the symplectic complement
Tz′Ψ(Z)⊥ω′={v∈N⨁k=0hkTz′kZ′|trk(ω′(v,w))=0∀w∈Tz′Ψ(Z)}, |
where z′:=∑Nk=0hkz′k:=Ψ(z).
Then the pull-back system (Z,ω,H) with ω=Ψ∗ω′, H=H′∘Ψ constitutes a formal Hamiltonian system with truncation index N and
trk(Ψ∗XH)=X′H′∘Ψ. |
Proof. Let ω=∑Nk=0hkωk Since Ψ0 is an embedding, it is an immersion and we conclude that ω0 is non-degenerate as ω′0 is non-degenerate. Since pull-back and the differential d commute, ω is closed.
In the following calculations we identify two formal polynomials P1 and P2 with coefficients of the same type (real numbers, n-forms, smooth functions, …) if and only if P1−P2∈O(hN+1). For all z∈Z we have
−ω′Ψ(z)(dΨ|z(XH(z)),dΨ|z(⋅))=−(Ψ∗ω′)z(XH(z),⋅)=−ωz(XH(z),⋅)=dH|z=d(H′∘Ψ)|z=dH′|Ψ(z)∘dΨ|z=−ω′Ψ(z)(X′H′(Ψ(z)),dΨ|z(⋅)) |
⟹ω′Ψ(z)(dΨ|z(XH(z))−X′H′(Ψ(z)),dΨ|z(⋅))=0⟹dΨ|z(XH(z))−X′H′(Ψ(z))∈TΨ(z)Ψ(Z)⊥ω′. |
As the inclusion dΨ|z(XH(z))−X′H′(Ψ(z))∈TΨ(z)Ψ(Z) holds as well and Ψ(Z) is a symplectic submanifold, it follows that
dΨ|z(XH(z))=X′H′(Ψ(z)). |
We can now proceed to the proof of 1.1.
Proof. Step 1. Construction and validity of the modified equation. The truncated power series L[N]Δ is a formal polynomial of the form
L[N]Δ(y[M])=L(y,˙y)+N∑k=1hkLkΔ(y[M]), |
where h is the formal variable. Since the Lagrangian L is regular, the Euler–Lagrange equations to L[N]Δ
M∑j=0(−1)jddt(∇y(j)L[N]Δ(y[M]))=0 |
are equivalent to an ordinary differential equation of the form
¨y=g0(y,˙y)+N∑k=1hk˜gk(y[2M]) | (5.1) |
when truncating terms of order O(h[N+1]). The ordinary differential equation is of order 2M because L[N]Δ(y[M]) is regular as well. Iterative replacements of derivatives of order j≥2 by derivatives of (5.1) yield the modified equation
¨y=g0(y,˙y)+N∑k=1hkgk(y,˙y)+O(h[N+1]). | (5.2) |
Solutions to
¨y=g0(y,˙y)+N∑k=1hkgk(y[M]) | (5.3) |
fulfil (5.1) up to O(h[N+1]) terms by construction. This proves the first part of 1.1.
Step 2. Existence of Hamiltonian structure on a higher jet-space. Let Y denote the domain of y (smooth manifold), Z=Jet1(Y) the 1-jet space, and Z′=Jet2M−1(Y).
The iterative substitution procedure gives rise to a formal polynomial Ψ of maps Z→Z′ defined by y[1]↦y[2M−1], where y(j) for j≥2 is expressed as a formal polynomial of functions depending on y[1]=(y,˙y). The expression for y(j) is obtained by deriving (5.1) j−2 times, iteratively replacing derivatives of order greater than 2 by derivatives of (5.1) followed by a truncation of O(h[N+1]) terms.
In the remainder of this part of the proof we suppress the fact that we operate on formal polynomials as the following steps can also be done when h is substituted with a sufficiently small real number h>0.
In the following, we will show that there exists a Hamiltonian structure (Z′,ω′,H′) on Z′ whose motions correspond to the motions induced by the order M-Lagrangian L[N]Δ(y[M]). To compute the Hamiltonian structure (Z′,ω′,H′), we first construct a Hamiltonian system on T∗JetM−1(Y) which we will then pull back to Z′.
As the Lagrangian L[N]Δ(y[M]) is regular, by Ostrogradsky's principle for high order Lagrangians [16], there exists a transformation χ:Z′→T∗JetM−1(Y) between the jet variable y[2M−1]∈Z′ and variables (q,p)∈T∗JetM−1(Y) with
q=(y,˙y,…,y(M−1)) |
and
pi=M−i∑k=0(−1)kdkdtk∇y(k+i)L[N]Δ(y[M]),i=1,…,M |
such that with
Ω=M∑i=1n∑j=1dpji∧dqij,andH[N](q,p)=M∑k=1⟨pk,qk⟩−L[N]Δ |
the motions of (5.1) are exactly mapped to the motions of the Hamiltonian vector field XH[N] on T∗JetM−1(Y) with
−dH[N]=Ω(XH[N],⋅) | (5.4) |
by the transformation χ−1. Let ω′=χ∗Ω and H′=H[N]∘χ. Now (Z′,ω′,H′) is a Hamiltonian system on Z′ whose motions are exactly the motions induced by the Lagrangian structure. This completes step 2 of the proof.
Remark 5.3. In the setting of formal polynomials, Ω is a two form whose coefficients are formal polynomials or, alternatively, Ω is a formal polynomial whose coefficients are 2-forms. In this case, (5.4) defines a formal series XH[N] in h from which we truncate O(hN+1) terms. This is also done when defining X′H′ through −dH′=ω′(X′H′,⋅). The differential equation defined by X′H′ is then equivalent to the differential equation induced by the Lagrangian structure up to O(hN+1) terms.
Step 3. Pull-back of (Z′,ω′,H′) to Z. As the Lagrangian L is regular, the 2-form ω′0=∑nj=1dyj∧dpj, where p is obtained by the Legendre transform for L, is a symplectic form. The pull-back form ω=Ψ∗ω′ is a formal polynomial in h, whose coefficients are 2-forms. The zeroth coefficient ω0 coincides with ω′0. ω is, therefore, non-degenerate and Ψ(Z) is a symplectic submanifold of Z′ in the sense of 5.2. Moreover, by construction of H′ and Ψ the condition trk(X′H′(Ψ(z)))∈TΨ(z)Ψ(Z) for all z∈Z is fulfilled. Therefore, 5.2 completes the theorem's proof.
Proof of 2.2. The construction method of the modified data J[N]mod and H[N]mod coincides with the construction method verified in the proof of 1.1.
Proposition 6.1. All Hamiltonian systems
˙z(t)=ˉJ−1(z(t))∇ˉH(z(t)),z=(y˙y) |
with a fixed symplectic structure represented by ˉJ can be formulated as variational problems of the form
δˉS=0forˉS(y)=∫ˉL(y(t),˙y(t))dt | (6.1) |
if and only if ˉJ is of the form
ˉJ=(∗∗∗0n), | (6.2) |
where 0n denotes an n×n-dimensional zero matrix with n the dimension of y, (i.e. if the distribution spanned by the vector fields ∂∂˙y1,…,∂∂˙yn is Lagrangian).
Proof. Denote the domain of the variable y by Y, the 1-jet space over Y by Jet1(Y), and the symplectic form represented by the matrix ˉJ by ˉω. Consider a Hamiltonian ˉH:Jet1(Y)→R. If the distribution D spanned by ∂∂˙y1,…,∂∂˙yn is Lagrangian w.r.t. ˉω, then there exists a primitive ˉλ of ˉω with kernel D [7,Cor.15.7]. The 1-form ˉλ is of the form
ˉλ=n∑j=1ℓ(y,˙y)dyj. |
By Hamilton's principle, a curve γ:[tstart,tend]→Y is a Hamiltonian motion if the action functional
ˉS(γ)=∫γ(ˉλ−ˉHdt) |
is stationary w.r.t. variations of γ through smooth curves fixing the endpoints. Thanks to the special structure of ˉλ (absence of d˙yj-terms), the pullback of ˉλ−ˉHdt along a curve γ has the form
ˉL(y(t),˙y(t))dt |
with y(t) describing γ(t) in the coordinate y. Therefore, in coordinates, the variational principle has the form (6.1).
On the other hand, a variational principle of the form (6.1) with regular Lagrangian ˉL (i.e. invertible (∂2ˉL∂˙yi∂˙yj)i,j) can be formulated as a Hamiltonian system with Hamiltonian
ˉH(y,˙y)=⟨˙y,p(y,˙y)⟩−L(y,˙y),and withp(y,˙y)=∂ˉL∂˙y(y,˙y). |
The symplectic structure is given as
ˉω=dˉλ,withˉλ=n∑j=1pj(y,˙y)dyj. |
As the distribution D spanned by ∂∂˙y1,…,∂∂˙yn is in the kernel of a primitive of ˉω, the distribution is Lagrangian.
Remark 6.2. The strength of 6.1 lies in the assertion that ˉL is a first-order Lagrangian in the original variable y, i.e. it depends on (y,˙y) only. If ˉJ is not of the required form, then, by Darboux's theorem, we can perform a change of variables on Jet1(Y) such that ˉω is the standard symplectic form ∑dˉpi∧dˉqi and ˉL is a 1st-order Lagrangian in q, i.e. depends on (q,˙q) but not on higher derivatives in q. However, as q and p each depend on (y,˙y), the variables have lost their dynamical meaning. This is because the required change of variables on Jet1(Y) is not fibred, i.e. the jet-space structure is not preserved. Expressed in the original variable y, the Lagrangian ˉL then depends on (y,˙y,¨y), i.e. describes a higher order variational structure.
Remark 6.3. The computational example presented in 4 provides an example for which the modified symplectic structure is not of the form that is required in 6.1, unless A1 is of the form A1=αI, i.e. the method coincides with a multistep method with scalar coefficients. Indeed, if the considered method is a classical multistep method, i.e. all coefficients are scalar, then Lmod(y,˙y) exists by 2.3.
We now proceed to the proof of 2.3. We exploit that linear multistep methods can be interpreted as 1-step methods on the original phase space [3,§XV.2]. Here and below we refer to the theory of linear multistep methods for 2nd order ODEs.
Proof of 2.3. As proved in [1,§5], for the underlying 1-step method ϕ of a symmetric linear multistep method there exists a local diffeomorphism ψ such that ˜ϕ=ψ∘ϕ∘ψ−1 is symplectic with respect to the original symplectic structure ω=∑jdpj∧dqj. The conjugacy ψ is given as a P-series applied to the original Hamiltonian vector field X0 given by the right hand side of
˙q=p˙p=∇U(q). |
The P-series ψ is a formal power series that is in general not convergent. Conjugacy, pull-back and push-forward operations are to be interpreted in a formal sense. The map ϕ is symplectic with respect to the modified symplectic structure ωmod=ψ∗ω and is the time-h-flow of a vector field X, for which the flow equations correspond to a first order formulation of the modified equation (2.5). The ω-symplectic map ˜ϕ is the time-h-flow of the ψ-related vector field ˜X=ψ∗X∘ψ−1. By standard results on backward error analysis for symplectic integrators [3,§Ⅸ], ˜X is a Hamiltonian vector field w.r.t. the standard symplectic structure ω for a Hamiltonian H up to any order in the step-size h. Pulling back the Hamiltonian structure (ω,H) using ψ, we obtain a modified Hamiltonian system (ωmod,Hmod)=(ψ∗ω,H∘ψ) such that Hamilton's equations are equivalent to the modified equation (2.5).
Since ψ is a P-series in X0, the distribution D spanned by the vertical vector fields ∂∂p1,…,∂∂pn is Lagrangian w.r.t. ωmod=ψ∗ω. By 6.1, for any order N in the step-size h there exists a first order Lagrangian Lmod(y,˙y) in the original variable such that the variational principle
δ(∫Lmod(y(t),˙y(t))dt)=0 |
recovers the modified equation up to higher order terms in h.
Remark 6.4. The proof of 2.3 also shows that Lmod in 2.3 has the structure of an S-series applied to a P-series (see [1]). The modified Lagrangian Lmod can, thus, be computed with an ansatz as well. The modified data Hmod and Jmod can then be computed from Lmod by a Legendre transformation.
Motivated by optimal truncation results for modified equations, it would be interesting to analyse the convergence properties of modified symplectic structures, modified Hamiltonians, and modified Lagrangians. Moreover, in view of 6.4, a systematic description of the combinatorial structure of the modified quantities appears feasible.
We would like to thank Peter Hydon and Chris Budd for discussions on multisymplecticity and variational principles at the Isaac Newton Institute of the University of Cambridge in 2019. This research was supported by the Marsden Fund of the Royal Society Te Apˉarangi and the School of Fundamental Sciences, Massey University, Manawatˉu, New Zealand. The authors would like to thank the Isaac Newton Institute for Mathematical Sciences for support and hospitality during the programme Geometry, compatibility and structure preservation in computational differential equations (2019) when work on this paper was undertaken. RM would like to thank the Simons Foundation for their support. This work was supported by: EPSRC grant number EP/R014604/1. The authors acknowledge support from the European Union Horizon 2020 research and innovation programmes under the Marie Skodowska-Curie grant agreement No. 691070.
[1] |
P. Chartier, E. Faou, A. Murua. An algebraic approach to invariant preserving integators: The case of quadratic and Hamiltonian invariants. Numer. Math, 103 (2006), 575–590. https://doi.org/10.1007/s00211-006-0003-8 doi: 10.1007/s00211-006-0003-8
![]() |
[2] | C. L. Ellison, J. W. Burby, J. M. Finn, H. Qin, W. M. Tang, Initializing and stabilizing variational multistep algorithms for modeling dynamical systems, 2014, arXiv preprint arXiv. |
[3] | E. Hairer, C. Lubich, G. Wanner. Geometric Numerical Integration: Structure-Preserving Algorithms for Ordinary Differential Equations. Springer Series in Com- putational Mathematics, 2013 https://doi.org/10.1007/3-540-30666-8 |
[4] |
L. Jódar, J. L. Morera, E. Navarro, On convergent linear multistep matrix methods, Int. J. Comput. Math, 40 (1991), 211–219. https://doi.org/10.1080/00207169108804014 doi: 10.1080/00207169108804014
![]() |
[5] |
J. D. Lambert, S. Sigurdsson, Multistep methods with variable matrix coefficients, Siam. J. Numer. Anal, 9 (1972), 715–733. https://doi.org/10.1137/0709060 doi: 10.1137/0709060
![]() |
[6] |
B. Leimkuhler, S. Reich, Simulating Hamiltonian Dynamics, Cambridge University Press, 2 (2005). https://doi.org/10.1017/CBO9780511614118 doi: 10.1017/CBO9780511614118
![]() |
[7] |
P. Libermann, C.M. Marle, Symplectic manifolds and Poisson manifolds, In Symplectic Geometry and Analytical Mechanics., (1987) 89–184. https://doi.org/10.1007/978-94-009-3807-6-3 doi: 10.1007/978-94-009-3807-6-3
![]() |
[8] |
J.E. Marsden, M. West. Discrete mechanics and variational integrators. Acta Numerica, 10 (2001), 357–514. https://doi.org/10.1017/S096249290100006X doi: 10.1017/S096249290100006X
![]() |
[9] |
R. I. McLachlan, C. Offen, Backward error analysis for variational discretisa-tions of PDEs, J. Geom. Mech, 14 (2022), 447. https://doi.org/10.3934/jgm.2022014 doi: 10.3934/jgm.2022014
![]() |
[10] |
S. Ober-Blöbaum, C. Offen, Variational learning of Euler-Lagrange dynam- ics from data. J. Comput. Appl. Math, (2022), 114780. https://doi.org/10.1016/j.cam.2022.114780 doi: 10.1016/j.cam.2022.114780
![]() |
[11] | C. Offen, Release of GitHub repository Christian-Offen/BEAConjugateSymplectic, Available from: https://github.com/Christian-Offen/BEAConjugateSymplectic. |
[12] | C. Offen, S. Ober-Blöbaum, Symplectic integration of learned Hamiltonian systems. Poster presented at the workshop Theory of deep learning (9 August 2021 to 13 August 2021), Isaac Newton Institute, University of Cambridge, Cambridge, UK, 2021. Available from: https://github.com/Christian-Offen/symplectic-shadow-integration/raw/master/PosterWorkshopCambridge.pdf. |
[13] |
C. Offen, S. Ober-Blöbaum, Symplectic integration of learned Hamiltonian systems. Chaos: An Interdisciplinary Journal of Nonlinear Science, 32 (2022). https://doi.org/10.1063/5.0065913 doi: 10.1063/5.0065913
![]() |
[14] | S. Sigurdsson, Linear multistep methods with variable matrix coefficients, In Lecture Notes in Mathematics, (1971), 327-331. Springer Berlin Heidelberg. https://doi.org/10.1007/BFb0069468 |
[15] |
M. Vermeeren, Modified equations for variational integrators, Numer. Math, 137 (2017), 1001–1037. https://doi.org/10.1007/s00211-017-0896-4 doi: 10.1007/s00211-017-0896-4
![]() |
[16] | E. T. Whittaker, S. W. McCrae, Hamiltonian systems and their integral invariants, Cambridge Mathematical Library, Cambridge University Press, https://doi.org/10.1017/CBO9780511608797.012 |
1. | Sina Ober-Blöbaum, Christian Offen, Variational learning of Euler–Lagrange dynamics from data, 2023, 421, 03770427, 114780, 10.1016/j.cam.2022.114780 |