Global dynamics of the chemostat with different removal rates and variable yields

1.   Introduction

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

which governs the exact Euler–Lagrange equations

by a discrete variational principle

Since $ y_0 $ and $ y_N $ are fixed in the variations considered in (1.1), the gradient above is taken with respect to all inner grid points $ y_1, \ldots, y_{N-1} $. We obtain the discrete Euler–Lagrange equations

which yield approximations $ y_i \approx y(t_0+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_\Delta(y_i, y_{i+1}) $ are established, discrete Lagrangians depending on several grid-points $ L_\Delta(y_{i}, y_{i+1}, \ldots, y_{i+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 $ \mathcal{L}_\Delta $ of consistent discrete Lagrangians $ L_\Delta(y(t), y(t+h), \ldots, y(t+sh)) $.

Theorem 1.1. Consider a power series $ \mathcal{L}_\Delta(y^{[\infty]}) $ in a formal variable $ h $. The series depends on the jet $ y^{[\infty]} = (y, \dot{y}, \ddot{y}, \ldots) $ 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, \dot{y}) $, i.e. $ \left(\frac{{\partial}^2 L}{{\partial} \dot{y}^i {\partial} \dot{y}^j}\right)_{i, j} $ is invertible.

● There exists a 2nd order modified equation given as a formal power series in $ h $ such that for any $ N \in {\mathbb{N}} $ a solution of the modified equation truncated to order $ \mathcal{O}(h^N) $ solves the Euler–Lagrange equations to $ \mathcal L^{[N]}_\Delta $ up to an error of order $ \mathcal{O}(h^{N+1}) $, where $ \mathcal L^{[N]}_\Delta $ is the truncation of $ \mathcal{L}_\Delta $ to order $ \mathcal{O}(h^N) $.

● If we denote $ z = (y, \dot{y}) $, there exists a symplectic structure matrix $ J_{\mathrm{mod}\, } $ and a Hamiltonian $ H_{\mathrm{mod}\, } $ given as formal power series in $ h $ such that solutions to Hamilton's equations

fulfil the modified equation. Here, $ J^{[N]}_{\mathrm{mod}\, } $ and $ H^{[N]}_{\mathrm{mod}\, } $ are truncations to order $ \mathcal{O}(h^N) $ of $ J_{\mathrm{mod}\, } $ and $ H_{\mathrm{mod}\, } $, respectively such that $ \mathcal{L}^{[N]}_\Delta $ is regular, i.e. $ \left(\frac{{\partial}^2 \mathcal{L}^{[N]}_\Delta}{{\partial} {y^{(M)\ \ }}^i {\partial} {y^{(M)\ \ }}^j}\right)_{i, j} $ is invertible, where $ y^{(M)} $ is the highest derivative of $ \mathcal{L}_\Delta^{[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, \dot{y}) $, it is natural to ask, whether there exists a modified Lagrangian $ L_{\mathrm{mod}\, }(y, \dot{y}) $ such that the modified equation is governed by the Euler–Lagrange equations to $ L_{\mathrm{mod}\, } $. In contrast to classical variational integrators, for which $ L_{\mathrm{mod}\, }(y, \dot{y}) $ is known to exist and can be computed ^[15], for conjugate symplectic schemes $ L_{\mathrm{mod}\, } $ may only exists in modified variables $ L_{\mathrm{mod}\, }(\tilde{y}, \dot{\tilde{y}}) $. We will prove that $ L_{\mathrm{mod}\, } $ exists in the original variables $ (y, \dot{y}) $ if $ J_{\mathrm{mod}\, } $ has the form

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 $ L_{\mathrm{mod}\, } $ exists in the original variables $ (y, \dot{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 $ L_{\mathrm{mod}\, } $ only exists in modified variables $ (\tilde{y}, \dot{\tilde{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 $ J_{\mathrm{mod}\, } $ and $ H_{\mathrm{mod}\, } $ 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.

2.   Application of blended backward error analysis to linear multistep methods with matrix coefficients

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

take the form

with

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 $ A_j $ ^[4,5,9,14], $ h $ is the step-size and $ I $ denotes the identity matrix. If the coefficients $ A_j $ 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

for some $ h $-independent expressions $ \tilde g_i $ in the $ a_i $-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

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

to the variational principle

with

Moreover, there exists a variational principle for (2.2):

Lemma 2.1. Let the matrices $ A_j \in {\mathbb{R}}^{n \times n} $ be symmetric. For $ T > 0 $ let $ \mathbb{T} $ either be the circle $ \mathbb{T} = {\mathbb{R}} / T {\mathbb{Z}} $ or the real line $ \mathbb{T} = {\mathbb{R}} $. For $ y $ defined on $ \mathbb{T} $ the variational principle

with

implies the functional equation (2.2). Moreover, if $ \mathbb{T} = [a, b] $ is an interval, then (2.6) implies (2.2) on the interval $\mathop {\mathbb{T}}\limits^\circ = [a+\frac{s}{2}h, b-\frac{s}{2}h] $. Here we assume that the function space for $ y $ and the potential $ U $ are such that $ U \circ y $ and $ \nabla U \circ y $ constitute square integrable functions.

Proof. Let $ \Delta_{\tau} y $ denote the forward difference, i.e. $ (\Delta_{\tau} y)(t) = y(t+\tau)-y(t) $. Then $ \langle y, \Delta_{\tau} z \rangle_{L^2(\mathbb{T}, {\mathbb{R}}^n)} = \langle \Delta_{\tau}^\ast y, z \rangle_{L^2(\mathbb{T}, {\mathbb{R}}^n)} $ holds with $ \Delta_{\tau}^\ast = \Delta_{-\tau} $ on $ \mathbb{T} $ if $ \mathbb{T} \in \{{\mathbb{R}}, {\mathbb{R}} / T{\mathbb{Z}}\} $ or on the interval $ [a+\tau, b-\tau] $ if $ \mathbb{T} = [a, b] $. The expression $ \Delta_{\tau}^\ast \Delta_{\tau} y $ corresponds to the central difference

Let $ \delta $ denote the variational derivative in the direction of a variation $ \delta y $, i.e. $ \delta S_\Delta (y) = \lim_{{\epsilon} \to 0} \frac{1}{{\epsilon}} (S_\Delta(y+{\epsilon} \delta y)-S_\Delta(y)) $. Let $ \delta y \in \mathcal{C}_c^\infty(\mathbb{T}, {\mathbb{R}}^n) $, if $ \mathbb{T} \in \{{\mathbb{R}}, {\mathbb{R}} / T{\mathbb{Z}}\} $ and let $ \delta y \in \mathcal{C}_c^\infty(\mathop {\mathbb{T}}\limits^\circ, {\mathbb{R}}^n) $, if $ \mathbb{T} = [a, b] $. We have

Now (2.2) follows from the fundamental lemma of the calculus of variations on $ \mathbb{T} $ or $ \mathop {\mathbb{T}}\limits^\circ $, respectively.

To analyse structure preserving properties of the method (2.2), it might seem natural to seek a modified Lagrangian $ L_{{\mathrm{mod}\, }}(y, \dot{y}) $ given as a formal power series in the step-size $ h $ such that the modified variational principle

covers smooth solutions of (2.6) up to any order in the step-size $ h $. However, we show that although a 1st order Lagrangian $ L_{{\mathrm{mod}\, }} $ covering the modified equations always exists as a power series, it only exist in modified variables $ (\tilde y, \dot{\tilde 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, \dot y) $ to $ (\tilde y, \dot{\tilde y}) $ can not be expected. This makes it difficult to compute $ L_{{\mathrm{mod}\, }} $ 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

Here

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

consisting of a modified Hamiltonian $ H_{\mathrm{mod}\, } $ and a modified symplectic structure $ J_{\mathrm{mod}\, } $ 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 $ J_{\mathrm{mod}\, } $ 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 \in {\mathbb{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]}_{\mathrm{mod}\, } $, which is $ \mathcal{O}(h) $ close to $ J $ and a modified Hamiltonian $ H^{[N]}_{\mathrm{mod}\, } $, which is $ \mathcal{O}(h) $ close to $ H $, such that

with coordinate $ z = (y, \dot y) $ is equivalent to the modified equation (2.5) up to terms of order $ \mathcal{O}(h^{N+1}) $.

Here $ O(h) $-closeness of $ J^{[N]}_{\mathrm{mod}\, } $ and $ J $ means that the zeroth coefficient of the polynomial $ J^{[N]}_{\mathrm{mod}\, } $ in the formal variable $ h $ is given by $ J $. $ O(h) $-closeness of $ H^{[N]}_{\mathrm{mod}\, } $ and $ H $ has an analogous meaning.

We observe conditions under which a modified first order Lagrangian $ L_{{\mathrm{mod}\, }}(y, \dot{y}) $ exists in the original variable $ y $ by analysing the modified symplectic structure $ J_{\mathrm{mod}\, } $.

Theorem 2.3. If the matrices $ A_j $ in (2.2) are scalar multiples of the identity matrix, then there exists a modified Lagrangian $ L^{[N]}_{\mathrm{mod}\, } $ depending on $ (y, \dot{y}) $ that is $ \mathcal{O}(h) $-close to $ L(y, \dot{y}) $ such that the modified variational principle

yields the modified equation (2.5) up to terms of order $ \mathcal{O}(h^{N+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.

3.   Computation of modified hamiltonian structure

In this section we introduce a method to compute the modified data $ J^{[N]}_{\mathrm{mod}\, } $ and $ H^{[N]}_{\mathrm{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 $ \mathcal{L}^{[N]}_\Delta $ denote the series expansion of

to order $ N $ in the step-size $ h $. The expression $ \mathcal{L}^{[N]}_\Delta $ depends on the $ M $-jet of $ y $ at $ t $ for some $ M \in {\mathbb{N}} $. Notice that the order $ M $ variational principle

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

and for $ i = 1, \ldots, M $

Here the index $ j $ enumerates the components of $ y $ and $ \frac{{\mathrm{d}}}{{\mathrm{d}} t} $ denotes the total derivative operator on the jet-space of $ y $, which acts like

on a scalar valued function $ \rho $ defined on the $ M $-jet space of $ y $. The high dimensional Hamiltonian system consists of the Hamiltonian

where all expressions are expressed in $ y^{[2M-1]} = (y, \dot y, \ldots, y^{(2M-1)}) $, and the symplectic structure matrix $ \mathcal{J}^{[N]}_{\mathrm{mod}\, }(y^{[2M-1]}) $. The skew-symmetric matrix $ \mathcal{J}^{[N]}_{\mathrm{mod}\, }(y^{[2M-1]}) $ is the representing matrix of the differential 2-form

where $ p_i^j $ and $ q^i_j $ are interpreted as functions in the variable $ y^{[2M-1]} $ of the $ 2M-1 $-jet space, i.e. $ \mathcal{J}^{[N]}_{\mathrm{mod}\, } $ is the anti-symmetrised tensor product^* $ \wedge $ of the gradients $ \nabla_{y^{[2M-1]}}\ \ p_j^i $ and $ \nabla_{y^{[2M-1]}} \ \ q^i_j $ summed over all indices.

^*This corresponds to the command ${\texttt{TensorWedge}}$ in Wolfram Mathematica.

To compute the modified Hamiltonian system on the 1-jet space with variable $ y^{[1]} = (y, \dot{y}) $, the variables $ y^{(2)}, \ldots, y^{(2M-1)} $ in the expression (3.2) for the Hamiltonian $ \mathcal H^{[N]}(y^{[2M-1]}) $ are repeatedly replaced by (2.5) until higher derivatives only occur in $ \mathcal{O}(h^{N+1}) $ terms. This yields $ H^{[N]}_{\mathrm{mod}\, }(y, \dot{y}) $. Similarly, we can consider $ p_i^j $ and $ q^i_j $ as functions of $ (y, \dot{y}) $ truncating $ \mathcal{O}(h^{N+1}) $ terms. The matrix $ J^{[N]}_{\mathrm{mod}\, } $ is then given as the representing matrix of the 2-form $ \Omega^{[N]} $ pulled to the 1-jet space of the variable $ y $, i.e. interpreting $ y^{[1]} = (y, \dot{y}) $ as the only independent variables. Equivalently, $ J^{[N]}_{\mathrm{mod}\, } $ is the anti-symmetrised tensor product $ \wedge $ of the gradients $ \nabla_{y^{[1]}} p_j^i $ and $ \nabla_{y^{[1]}} q^i_j $ summed over all indices. (As $ J^{[N]}_{{\mathrm{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.

4.   Computational example

As introduced in 2, consider the multistep method

with matrix coefficients in dimension $ n = 2 $. By the consistency requirement, $ A_2 = 1/4(I-A_1) $. We obtain

with

where

Here and in the following $ U^{(k, l)} = \frac{{\partial}^{k+l}U(y)}{{\partial}^k y_1 {\partial}^l y_2} $.

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 $ J^4_{22} $ for the special case that $ A_1 $ is a diagonal matrix: we have

with

The modified Hamiltonian $ H^{[4]}_{\mathrm{mod}\, } $ for the general case is given as

with

For further terms, we refer the reader to the Mathematica Notebooks of our accompanying source code. Hamilton's equations

for the modified Hamiltonian system are equivalent to the modified equation (2.5) truncating terms of order $ \mathcal{O}(h^6) $.

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]}_{\mathrm{mod}\, } = H $, $ H^{[2]}_{\mathrm{mod}\, } $, and $ H^{[4]}_{\mathrm{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]}_{\mathrm{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]}_{\mathrm{mod}\, }-H(z_{\mathrm{init}}) $, $ k = 0, 2, 4 $ which eventually prevent further energy error decay.

If $ A_2 = \alpha I $ with $ \alpha \ge \frac{1}{4} $ and $ A_1 = I-4A_2 $, then (4.1) corresponds to a classical stable^† multistep scheme: the generating polynomial $ \rho $ 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].

The polynomial $ \rho $ has a double root at 1 as well as the roots

Since $ \alpha \ge \frac{1}{4} $, the roots $ \xi_3 $ and $ \xi_4 $ are complex conjugate to each other and lie on the unit circle such that the scheme is stable.

Moreover, $ L_\Delta $ and $ \mathcal{L}_\Delta $ are rotationally invariant because $ A_1 $ and $ A_2 $ commute with rotation matrices. An application of Noether's theorem to $ \mathcal{L}^{[N]}_\Delta $ yields the following modified angular momentum:

The integer $ M $ is the order of the highest derivative of $ y $ in $ \mathcal{L}^{[N]}_\Delta $. 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 $, $ \dot{y} $ and $ \mathcal{O}(h^6) $-terms. After truncation at order 4 in $ h $ we obtain the modified angular momentum $ I^{[4]}_{\mathrm{mod}\, }(y, \dot{y}) $. Figure 2 shows the conservation error of $ I = I^{[0]}_{\mathrm{mod}\, } $, $ I^{[2]}_{\mathrm{mod}\, } $, $ I^{[4]}_{\mathrm{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].

5.   Validity of the construction method

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', \omega', H') $ is a Hamiltonian system with Hamiltonian vector field $ X'_{H'} $, $ Z \mathop {\hookrightarrow}\limits^\Psi Z' $ an embedding such that $ \Psi(Z) $ is a symplectic submanifold that is invariant under motions, i.e. $ X'_{H'}(\Psi(z)) \in T_{\Psi(z)} \Psi(Z) $ for all $ z \in Z $, then the pull-back system $ (Z, \omega, H) $ with $ \omega = \Psi^\ast \omega' $, $ H = H' \circ \Psi $ is a Hamiltonian system with Hamiltonian vector field $ X_H $. Here $ \Psi^\ast \omega' $ denotes the pull-back of $ \omega' $ along $ \Psi $. The Hamiltonian vector fields $ X_H $ and $ X'_{H'} $ relate by

i.e. for a motion $ \gamma $ with $ \dot{\gamma} = X_H \circ \gamma $ the curve $ \gamma' = \Psi \circ \gamma $ is a motion of $ (Z', \omega', H') $, i.e. $ \dot{\gamma}' = X'_{H'} \circ \gamma' $. Here $ \Psi_\ast $ denotes the push-forward map, i.e. $ (\Psi_\ast X_H)(z) = {\mathrm{d}} \Psi|_z(X_H(z)) $ for $ z \in Z $.

In the following, we will adapt this statement to a setting, where the definition of $ \omega' $, $ H' $, and $ \Psi' $ contain a formal variable $ h $ and where $ \Psi(Z) $ is left invariant only up to higher order terms in $ h $.

Definition 5.1 (Formal Hamiltonian system). Let $ N \in {\mathbb{N}} $ be a truncation index, $ Z $ be a smooth manifold, $ h $ a formal variable, and $ H = \sum_{k = 0}^N h^k H_k $ a formal polynomial whose coefficients are smooth maps $ Z \to {\mathbb{R}} $. Further, let $ \omega = \sum_{k = 0}^{ N} h^k \omega_k $ be a formal polynomial whose coefficients are 2-forms on $ Z $. The collection $ (Z, \omega, H) $ is called a formal Hamiltonian system with truncation index $ N $ if the following conditions are satisfied.

● The formal symplectic form $ \omega $ is closed, i.e. $ {\mathrm{d}} \omega = \sum_{k = 0}^{ N} h^k {\mathrm{d}} \omega_k $ is a formal polynomial whose coefficients are 3-forms that are zero.

● The 2-form $ \omega_0 $ is non-degenerate.

The relation $ -{\mathrm{d}} H = \omega(X, \cdot) $ defines a formal power series $ X = \sum_{k = 0}^\infty h^k X_k $, where $ X_k $ are vector fields on $ Z $. The truncation $ X_H : = \sum_{k = 0}^N h^k X_k $ is called (formal) Hamiltonian vector field. The formal differential equation $ \dot{\gamma} = X_H \circ \gamma $ is called Hamilton's equation.

Proposition 5.2. Let $ N \in {\mathbb{N}} $ be a truncation index, $ Z' $ and $ Z $ be manifolds, and let $ (Z', \omega', H') $ be a formal Hamiltonian system with truncation index $ N $. Consider a formal polynomial $ \Psi = \sum_{k = 0}^N \Psi_k h^k $ such that $ \Psi_k \colon Z \to Z' $ are smooth and $ \Psi_0 \colon Z \to Z' $ is an embedding. To $ z \in Z $ define formal tangent spaces

where $ z' : = \sum_{k = 0}^N h^k z_k' : = \Psi(z) $.

● Assume that the Hamiltonian vector field is tangential to $ \Psi(Z) $, i.e. $ {\mathrm{trk}} (X'_{H'}(\Psi(z))) \in T_{\Psi(z)}\Psi(Z) $ for all $ z \in Z $. Here $ {\mathrm{trk}} $ denotes the truncation of $ \mathcal{O}(h^{N+1}) $-terms.

● Assume that $ \Psi(Z) $ is a symplectic submanifold of $ Z' $, i.e. for all $ z \in Z $

Here $ \Psi(Z)^{\bot_{\omega'}} $ denotes the symplectic complement

where $ z' : = \sum_{k = 0}^N h^k z_k' : = \Psi(z) $.

Then the pull-back system $ (Z, \omega, H) $ with $ \omega = \Psi^\ast \omega' $, $ H = H' \circ \Psi $ constitutes a formal Hamiltonian system with truncation index $ N $ and

Proof. Let $ \omega = \sum_{k = 0}^N h^k\omega_k $ Since $ \Psi_0 $ is an embedding, it is an immersion and we conclude that $ \omega_0 $ is non-degenerate as $ \omega_0' $ is non-degenerate. Since pull-back and the differential $ {\mathrm{d}} $ commute, $ \omega $ is closed.

In the following calculations we identify two formal polynomials $ P_1 $ and $ P_2 $ with coefficients of the same type (real numbers, $ n $-forms, smooth functions, $ \ldots $) if and only if $ P_1-P_2 \in \mathcal{O}(h^{N+1}) $. For all $ z \in Z $ we have

As the inclusion $ {\mathrm{d}} \Psi|_z(X_H(z))-X'_{H'}(\Psi(z)) \in T_{\Psi(z)}\Psi(Z) $ holds as well and $ \Psi(Z) $ is a symplectic submanifold, it follows that

We can now proceed to the proof of 1.1.

Proof. Step 1. Construction and validity of the modified equation. The truncated power series $ \mathcal{L}_\Delta^{[N]} $ is a formal polynomial of the form

where $ h $ is the formal variable. Since the Lagrangian $ L $ is regular, the Euler–Lagrange equations to $ \mathcal{L}_\Delta^{[N]} $

are equivalent to an ordinary differential equation of the form

when truncating terms of order $ \mathcal{O}(h^{[N+1]}) $. The ordinary differential equation is of order $ 2M $ because $ \mathcal{L}_\Delta^{[N]}(y^{[M]}) $ is regular as well. Iterative replacements of derivatives of order $ j\ge 2 $ by derivatives of (5.1) yield the modified equation

Solutions to

fulfil (5.1) up to $ \mathcal{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 = {\mathrm{Jet}}^1(Y) $ the 1-jet space, and $ Z' = {\mathrm{Jet}}^{2M-1}(Y) $.

The iterative substitution procedure gives rise to a formal polynomial $ \Psi $ of maps $ Z \to Z' $ defined by $ y^{[1]} \mapsto y^{[2M-1]} $, where $ y^{(j)} $ for $ j\ge 2 $ is expressed as a formal polynomial of functions depending on $ y^{[1]} = (y, \dot{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 $ \mathcal{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', \omega', H') $ on $ Z' $ whose motions correspond to the motions induced by the order $ M $-Lagrangian $ \mathcal{L}_\Delta^{[N]}(y^{[M]}) $. To compute the Hamiltonian structure $ (Z', \omega', H') $, we first construct a Hamiltonian system on $ T^\ast {\mathrm{Jet}}^{M-1}(Y) $ which we will then pull back to $ Z' $.

As the Lagrangian $ \mathcal{L}_\Delta^{[N]}(y^{[M]}) $ is regular, by Ostrogradsky's principle for high order Lagrangians ^[16], there exists a transformation $ \chi \colon Z' \to T^\ast {\mathrm{Jet}}^{M-1}(Y) $ between the jet variable $ y^{[2M-1]}\in Z' $ and variables $ (q, p) \in T^\ast {\mathrm{Jet}}^{M-1}(Y) $ with

and

such that with

the motions of (5.1) are exactly mapped to the motions of the Hamiltonian vector field $ X_{\mathcal{H}^{[N]}} $ on $ T^\ast {\mathrm{Jet}}^{M-1}(Y) $ with

by the transformation $ \chi^{-1} $. Let $ \omega' = \chi^\ast \Omega $ and $ H' = \mathcal{H}^{[N]} \circ \chi $. Now $ (Z', \omega', 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, $ \Omega $ is a two form whose coefficients are formal polynomials or, alternatively, $ \Omega $ is a formal polynomial whose coefficients are 2-forms. In this case, (5.4) defines a formal series $ X_{\mathcal{H}^{[N]}} $ in $ h $ from which we truncate $ \mathcal{O}(h^{N+1}) $ terms. This is also done when defining $ {X'}_{H'} $ through $ -{\mathrm{d}} H' = \omega'({X'}_{H'}, \cdot) $. The differential equation defined by $ {X'}_{H'} $ is then equivalent to the differential equation induced by the Lagrangian structure up to $ \mathcal{O}(h^{N+1}) $ terms.

Step 3. Pull-back of $ (Z', \omega', H') $ to $ Z $. As the Lagrangian $ L $ is regular, the 2-form $ \omega_0' = \sum_{j = 1}^n{\mathrm{d}} y_j \wedge {\mathrm{d}} \mathfrak{p}^j $, where $ \mathfrak{p} $ is obtained by the Legendre transform for $ L $, is a symplectic form. The pull-back form $ \omega = \Psi^\ast\omega' $ is a formal polynomial in $ h $, whose coefficients are 2-forms. The zeroth coefficient $ \omega_0 $ coincides with $ \omega_0' $. $ \omega $ is, therefore, non-degenerate and $ \Psi(Z) $ is a symplectic submanifold of $ Z' $ in the sense of 5.2. Moreover, by construction of $ H' $ and $ \Psi $ the condition $ {\mathrm{trk}} (X'_{H'}(\Psi(z))) \in T_{\Psi(z)}\Psi(Z) $ for all $ z \in Z $ is fulfilled. Therefore, 5.2 completes the theorem's proof.

Proof of 2.2. The construction method of the modified data $ J^{[N]}_{\mathrm{mod}\, } $ and $ H^{[N]}_{\mathrm{mod}\, } $ coincides with the construction method verified in the proof of 1.1.

6.   On the existence of modified lagrangians

Proposition 6.1. All Hamiltonian systems

with a fixed symplectic structure represented by $ \bar{J} $ can be formulated as variational problems of the form

if and only if $ \bar{J} $ is of the form

where $ 0_n $ denotes an $ n\times n $-dimensional zero matrix with $ n $ the dimension of $ y $, (i.e. if the distribution spanned by the vector fields $ \frac{{\partial}}{{\partial} \dot{y}_1}, \ldots, \frac{{\partial}}{{\partial} \dot{y}_n} $ is Lagrangian).

Proof. Denote the domain of the variable $ y $ by $ Y $, the 1-jet space over $ Y $ by $ \mathrm{Jet}^1(Y) $, and the symplectic form represented by the matrix $ \bar{J} $ by $ \bar{\omega} $. Consider a Hamiltonian $ \bar{H}\colon \mathrm{Jet}^1(Y) \to {\mathbb{R}} $. If the distribution $ \mathcal{D} $ spanned by $ \frac{{\partial}}{{\partial} \dot{y}_1}, \ldots, \frac{{\partial}}{{\partial} \dot{y}_n} $ is Lagrangian w.r.t. $ \bar{\omega} $, then there exists a primitive $ \bar{\lambda} $ of $ \bar{\omega} $ with kernel $ \mathcal{D} $ [7,Cor.15.7]. The 1-form $ \bar{\lambda} $ is of the form

By Hamilton's principle, a curve $ \gamma\colon [t_{\mathrm{start}}, t_{\mathrm{end}}] \to Y $ is a Hamiltonian motion if the action functional

is stationary w.r.t. variations of $ \gamma $ through smooth curves fixing the endpoints. Thanks to the special structure of $ \bar{\lambda} $ (absence of $ {\mathrm{d}} \dot{y}_j $-terms), the pullback of $ \bar{\lambda} - \bar{H} {\mathrm{d}} t $ along a curve $ \gamma $ has the form

with $ y(t) $ describing $ \gamma(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 $ \bar{L} $ (i.e. invertible $ \left(\frac{{\partial}^2\bar{L}}{{\partial} \dot{y}_i {\partial} \dot{y}_j }\right)_{i, j} $) can be formulated as a Hamiltonian system with Hamiltonian

The symplectic structure is given as

As the distribution $ \mathcal{D} $ spanned by $ \frac{{\partial}}{{\partial} \dot{y}_1}, \ldots, \frac{{\partial}}{{\partial} \dot{y}_n} $ is in the kernel of a primitive of $ \bar{\omega} $, the distribution is Lagrangian.

Remark 6.2. The strength of 6.1 lies in the assertion that $ \bar{L} $ is a first-order Lagrangian in the original variable $ y $, i.e. it depends on $ (y, \dot{y}) $ only. If $ \bar{J} $ is not of the required form, then, by Darboux's theorem, we can perform a change of variables on $ \mathrm{Jet}^1(Y) $ such that $ \bar{\omega} $ is the standard symplectic form $ \sum {\mathrm{d}} \bar{p}_i \wedge {\mathrm{d}} \bar{q}_i $ and $ \bar{L} $ is a 1st-order Lagrangian in $ q $, i.e. depends on $ (q, \dot{q}) $ but not on higher derivatives in $ q $. However, as $ q $ and $ p $ each depend on $ (y, \dot{y}) $, the variables have lost their dynamical meaning. This is because the required change of variables on $ \mathrm{Jet}^1(Y) $ is not fibred, i.e. the jet-space structure is not preserved. Expressed in the original variable $ y $, the Lagrangian $ \bar{L} $ then depends on $ (y, \dot{y}, \ddot{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 $ A_1 $ is of the form $ A_1 = \alpha 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 $ L_{\mathrm{mod}\, }(y, \dot{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 $ \phi $ of a symmetric linear multistep method there exists a local diffeomorphism $ \psi $ such that $ \tilde \phi = \psi \circ \phi \circ \psi^{-1} $ is symplectic with respect to the original symplectic structure $ \omega = \sum_j {\mathrm{d}} p_j \wedge {\mathrm{d}} q_j $. The conjugacy $ \psi $ is given as a $ P $-series applied to the original Hamiltonian vector field $ X_0 $ given by the right hand side of

The $ P $-series $ \psi $ 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 $ \phi $ is symplectic with respect to the modified symplectic structure $ \omega_{\mathrm{mod}\, } = \psi^\ast \omega $ 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 $ \omega $-symplectic map $ \tilde{\phi} $ is the time-$ h $-flow of the $ \psi $-related vector field $ \tilde{X} = \psi_\ast X \circ \psi^{-1} $. By standard results on backward error analysis for symplectic integrators [3,§Ⅸ], $ \tilde X $ is a Hamiltonian vector field w.r.t. the standard symplectic structure $ \omega $ for a Hamiltonian $ H $ up to any order in the step-size $ h $. Pulling back the Hamiltonian structure $ (\omega, H) $ using $ \psi $, we obtain a modified Hamiltonian system $ (\omega_{\mathrm{mod}\, }, H_{\mathrm{mod}\, }) = (\psi^\ast \omega, H \circ \psi) $ such that Hamilton's equations are equivalent to the modified equation (2.5).

Since $ \psi $ is a $ P $-series in $ X_0 $, the distribution $ \mathcal D $ spanned by the vertical vector fields $ \frac{{\partial}}{{\partial} p_1}, \ldots, \frac{{\partial}}{{\partial} p_n} $ is Lagrangian w.r.t. $ \omega_{\mathrm{mod}\, } = \psi^\ast \omega $. By 6.1, for any order $ N $ in the step-size $ h $ there exists a first order Lagrangian $ L_{\mathrm{mod}\, }(y, \dot{y}) $ in the original variable such that the variational principle

recovers the modified equation up to higher order terms in $ h $.

Remark 6.4. The proof of 2.3 also shows that $ L_{\mathrm{mod}\, } $ in 2.3 has the structure of an $ S $-series applied to a $ P $-series (see ^[1]). The modified Lagrangian $ L_{\mathrm{mod}\, } $ can, thus, be computed with an ansatz as well. The modified data $ H_{\mathrm{mod}\, } $ and $ J_{\mathrm{mod}\, } $ can then be computed from $ L_{\mathrm{mod}\, } $ by a Legendre transformation.

7.   Future work

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.

Acknowledgements

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${\bar{{\rm{a}}}}$rangi and the School of Fundamental Sciences, Massey University, Manawat${\bar{{\rm{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.