Modified equations for variational integrators applied to Lagrangians linear in velocities

Variational integrators applied to degenerate Lagrangians that are linear in the velocities are two-step methods. The system of modified equations for a two-step method consists of the principal modified equation and one additional equation describing parasitic oscillations. We observe that a Lagrangian for the principal modified equation can be constructed using the same technique as in the case of non-degenerate Lagrangians. Furthermore, we construct the full system of modified equations by doubling the dimension of the discrete system in such a way that the principal modified equation of the extended system coincides with the full system of modified equations of the original system. We show that the extended discrete system is Lagrangian, which leads to a construction of a Lagrangian for the full system of modified equations.


Introduction
An important technique to study the long-time behavior of numerical integrators is backward error analysis. This consists in finding a modified equation, a perturbation of the original differential equation whose solutions exactly interpolate the numerical solutions. When a modified equation has been found, one can study the behavior of the numerical solutions by comparing two differential equations, rather than comparing a differential equation with a difference equation. Several long-time (near) conservation laws for symplectic integrators can be proved this way. For a detailed introduction to modified equations we refer to [7,Chapter IX].
In [15] we considered modified equations for variational integrators in the case of non-degenerate Lagrangians. We gave a construction for a modified Lagrangian, which produces the modified equation as its Euler-Lagrange equation up to a truncation error of arbitrarily high order. Although the construction was new, the claim that modified equations for variational integrators are Lagrangian was not. This follows by Legendre transformation from the well-known fact that modified equations for symplectic integrators are Hamiltonian. The construction of a modified Lagrangian was combined in [2,3] with the idea of modifying integrators [1] to construct variational integrators of improved convergence order.
In this work we extend our previous construction to the case of degenerate Lagrangians that are linear in velocities. In this context the Legendre transformation is not invertible, so the fact that the modified equation is Lagrangian cannot be inferred in the same way from the theory of symplectic integrators. We consider Lagrangians L : where α : R N → R N , H : R N → R, and the brackets , denote the standard scalar product. Variational integrators for such Lagrangians were studied for example in [12] and [14]. An important role will be played by the matrices A(q) = α (q) = ∂α i (q) ∂q j i,j=1,...,N and A skew (q) = A(q) T − A(q).
We assume that A skew (q) is invertible, then the Euler-Lagrange equation for L is given byq where q is considered to be a column vector and H (q) is the row vector of partial derivatives of H with respect to q 1 , . . . , q N . In contrast to the case of non-degenerate Lagrangians, this is a first order ODE. A well-known example where a Lagrangian of the form (1) arises is the dynamics of point vortices in the plane. We will discuss this example in detail in Section 6.2. Another reason to study this class of Lagrangians is that its extension to PDEs covers several important equations. For example, the nonlinear Schrödinger equation is the Euler-Lagrange equation of a Lagrangian whose kinetic term is linear in the time-derivatives (see e.g. [13,Section 2.1.]). Perhaps the most general application of Lagrangians that are linear in velocities is the variational formulation in phase space of mechanics, where L : T T * Q ∼ = R 4N → R is given by L(p, q,ṗ,q) = p ,q − H(p, q).
Its Euler-Lagrange equations are Hamilton's canonical equationṡ Note that even though A = 0 0 I 0 is singular in this case, the assumption that A skew is invertible still holds. Like many concepts in classical mechanics, the variational principle in phase space dates back to the 19th century [11,Chapter XXIX]. A modern treatment can be found for example in [5,, and an application to geometric integration in [8].
The construction of modified Lagrangians for variational integrators, which we introduced in [15], carries over to the case of degenerate Lagrangians that are linear in velocities. However, there is a catch. The original differential equation is of first order for the Lagrangians considered here, but the difference equation produced by a variational integrator is of second order. Hence, in this context, variational integrators are two-step methods and parasitic solutions can occur.
In Section 2 we present two variational integrators which will be the protagonists of all examples discussed in this work. In Section 3 the essentials of the theory of modified equations for multi-step methods are reviewed. In Section 4 we summarize the construction of modified Lagrangians from [15] and in Section 5 we will present a method to extend it to the full system of modified equations. In Section 6 we look at some example systems.

A note on notation
As mentioned above, we use the convention that a derivative with respect to a column vector yields a row vector. In particular, this means that the derivative of the scalar product of two column vectors is calculated as Later on we will be taking higher derivatives of vectors with respect to other vectors, resulting in a zoo of tensors. We want to avoid heavy notations using indices, like If the tensor involved is symmetric, we will use the notations instead. If the tensor is of first or second order, we will often write these expressions as matrix multiplication, Ax and x T By.
We will also use the inner product notation A T , x as an alternative to Ax. This allows us to emphasize one particular pairing in a product of more than two tensors. Using these notations interchangeably allows us to write equations in an intuitive form and avoid the heavy notation of (4). The downside is that such inconsistent notation could be a source of confusion for the reader. We hope this note is enough to avoid that.

Variational integrators
A variational integrator is a numerical integrator for Lagrangian differential equations, obtained by discretizing the Lagrange function. The action integral L(q,q) dt is replaced by a sum j L disc (q j , q j+1 , h). The sequence (q j ) j∈Z is a critical point of the action sum if and only if it satisfies the discrete Euler-Lagrange equation where D 1 and D 2 denote partial derivatives with respect to the first and second variable. Assuming this difference equation can be solved for q j+1 , it provides a numerical approximation of the Euler-Lagrange equations of L. An excellent overview of the subject of variational integrators is given by Marsden and West [9]. For Lagrangians that are linear in the velocities, the continuous Euler-Lagrange equation (3) is of first order, but the discrete Euler-Lagrange equation (5) involves three points, i.e. it is of second order. This means that we are dealing with two-step methods.
We will discuss two examples of variational integrators in detail. Both are obtained by using a simple quadrature rule to approximate the exact discrete Lagrangian where q(jh) = q j , q((j + 1)h) = q j+1 , and q(t) solves the continuous Euler-Lagrange equation.

Midpoint rule
Using q j+1 −q j 2 to approximateq and the average q j +q j+1 2 to approximate q in the integrand, we find the discrete Lagrangian with discrete Euler-Lagrange equation In case α is linear, i.e. α(q) = Aq, this simplifies to where A skew is defined in Equation (2). In the case of a non-degenerate Lagrangian this discretization would lead to a variational integrator that is equivalent to the implicit midpoint rule applied to the corresponding symplectic system. Also in the present context we will refer to it as the midpoint rule.

Trapezoidal rule
To obtain the second discretization we use the trapezoidal quadrature rule to approximate the exact discrete Lagrangian: we take the average of the integrand evaluated with q = q j and with q = q j+1 , while still using q j+1 −q j 2 to approximate the derivativeq. We find the discrete Lagrangian with discrete Euler-Lagrange equation In case α is linear this simplifies to This discretization is sometimes called the explicit midpoint rule, but we will not use this name to avoid confusion with the previous method. Instead we call this method the trapezoidal rule. In the case of a non-degenerate Lagrangian the trapezoidal rule would lead to the Störmer-Verlet method.

Modified equations for multistep methods
The classical theory of modified equations does not capture parasitic solutions of multistep methods. An extension of this theory for linear multistep methods was developed by Hairer [6]. (See also [7,Chapter XV].) Here we mention some of the main results, restricted to the case of two-step methods. For a first order ODEq = f (q), consider the linear two-step method We call the method (8) symmetric if a 0 = −a 2 , a 1 = 0, and b 0 = b 2 . We say that it is stable if all roots of the polynomial ρ(ζ) = a 0 + a 1 ζ + a 2 ζ 2 satisfy |ζ| ≤ 1, and the roots with |ζ| = 1 are simple. A method is stable if and only if the numerical solution forq = 0 is bounded for any initial condition. Note that the trapezoidal rule is a stable symmetric linear two-step method, but that the midpoint rule is not of the form (8).
The theory of modified equations for one-step methods is easily extended to yield the following.
. Consider a consistent method of the form (8). Then there exist unique functions (f n (q)) n∈N such that for every truncation index k, every solution oḟ In general the right hand side of Equation (9) will not converge as k → ∞. Nevertheless, we will call the formal differential equatioṅ Proposition 3.2 (Special case of Theorem XV.3.5 from [7]). Assume that the method (8) is stable, consistent, and symmetric. Then there exist functions (f n (x, y)) n∈N and (g n (x, y)) n∈N such that for every truncation index k, for every solution oḟ with for every choice of t.
We will call the corresponding system of formal differential equationṡ the full system of modified equations. We call Equation (14) the parasitic modified equation.
If y = 0, then Equation (13) reduces to the principal modified equation (10) and Equation (14) readsẏ = 0. Hence to determine whether parasitic solutions become dominant over time we need to determine the stability of the invariant manifold {y = 0} of the system (13)- (14).
In general, even if the difference equation is not of the form (8), we have the following definition.
(a) Equation (10) is the principal modified equation for the difference equation if for every truncation index k, every solution of the truncated equation (9) satisfies at all times t.
(b) The system of equations (13)-(14) is the full system of modified equations for the Equation (15) if for every truncation index k, for every solution (x, y) of the truncated system (11)-(12), the discrete curve q j = x(t + jh) for all choices of t.

A Lagrangian for the principal modified equation
In [15] we constructed a modified Lagrangian for variational integrators in the case of non-degenerate Lagrangian systems. A straightforward adaptation of this construction will give us a Lagrangian for the principal modified equation. Here we present the construction and a rough sketch of the proof. The details of the proof are perfectly analogous to the non-degenerate case, so we refer to [15] for their discussion. We identify points q j of a numerical solution with step size h with evaluations q(jh) of an interpolating curve. Using a Taylor expansion we can write the discrete Lagrangian L disc (q j−1 , q j , h) as a function of the interpolating curve q and its derivatives, all evaluated at the point jh − h 2 , where the square brackets denote dependence on q and any number of its derivatives. We want to write the discrete action as an integral. This can be done using the Euler-Maclaurin formula. We obtain the meshed modified Lagrangian where B 2i are the Bernoulli numbers. The power series defining L mesh generally does not converge. Formally, it satisfies Note that L mesh depends on higher derivatives of q. Below we will construct a modified Lagrangian that only depends on q andq. The word meshed refers to the fact that the discrete system provides additional structure for the continuous variational problem. In the meshed variational problem, nondifferentiable curves are admissible as long as their singular points are consistent with the mesh, i.e. if they occur at times that are an integer multiple of h away from each other. This imposes additional conditions on critical curves, related to the natural boundary conditions and to the Weierstrass-Erdmann Corner conditions (see e.g. [4,Sec. 6 and 13] for these concepts). These conditions are We will call them the natural interior conditions. Because the action integral of L mesh equals the discrete action, variations supported on a single mesh interval (i.e. in between consecutive points of the discrete curve) do not change the action integral of L mesh . This implies that the natural interior conditions are automatically satisfied on solutions of the Euler-Lagrange equation (for the particular Lagrangian L mesh , but not in general).
Because the natural interior conditions (16) are automatically satisfied on critical curves, it is equivalent to This equation is of the form E 0 (q,q) + hE 1 (q,q,q) + h 2 E 2 q,q,q, q (3) + . . . = 0.
In the leading order we find a first order differential equation, which we can use to eliminate higher derivatives in the next order (assuming that the derivatives of q are bounded as h → 0). This can be applied recursively up to any order. Hence we can write the Euler-Lagrange equation formally as a first order differential equation, saẏ Then expressions for all higher derivatives follow by differentiation and substitution, 8 The assumption that the derivatives of q are bounded as h → 0 is not restrictive in practice. The same assumption is necessary to state many other results regarding modified equations rigorously. Families of curves (q h ) h∈(0,∞) that satisfy this condition are called admissible families in [15]. In particular there holds for admissible families that if the functions are small, q h = O(h k ), then so are their derivatives, q . We will use this implicitly later on.
Using (17) and (18) we can replace second and higher derivatives in the meshed Lagrangian to find a first order modified Lagrangian, Or, avoiding formal power series, a truncated modified Lagrangian where T k denotes truncation of the power series after order k. In general the replacements q (j) = F j (q, h) would change the Euler-Lagrange equations, but because of the natural interior conditions (16) this is not the case here. Indeed, one finds It follows that so up to a truncation error, both Lagrangians yield the same Euler-Lagrange equations. Note that the natural interior conditions do not imply that ∂L mesh /∂q = 0, so replacing first derivatives usingq = F (q, h) is not allowed! The details presented in [15] carry over to the degenerate case and yield the following result.
Theorem 4.1. Consider a discrete Lagrangian that is a consistent discretization of a Lagrangian of the form (1). Let L be either L mesh or L mod,k , derived from this discrete Lagrangian. Solve the equation forq, and truncate the resulting power series after order k. The result, is a truncation of the principal modified equation.
From the discrete Lagrangian (6) we find It follows that where the argument q of A skew , α, H, and their derivatives is omitted. From this expression we obtain L mod,3 by replacing all second derivatives of q using the derivative of the leading order equation, In case that α is linear we haveq and we find the following expression for the modified Lagrangian (truncated after h 3 ):

Trapezoidal rule
From the discrete Lagrangian (7) we find and Again we assume that α is linear. Using Equation (19) we find the modified Lagrangian

The full system of modified equations
For linear symmetric two-step methods, Proposition 3.2 describes the full system of modified equations. Here we will show that for variational integrators, without assuming linearity, the full system of modified equations is of the same form. In order to construct the system of modified equations, we split the variable q j of the discrete system into two parts, The motivation for this is that we want to use one variable, x j , to encode the principal behavior and the other, y j , for the parasitic behavior. This is inspired by the formula q j = x(t + jh) + (−1) j y(t + jh) from Proposition 3.2 and Definition 3.3.

The Lagrangian approach
A key property of the doubling of variables is that the extended system is still variational.
Proposition 5.1. The discrete curve (x j , y j ) j∈Z is critical for if and only if the discrete curves (q + j ) j∈Z and (q − j ) j∈Z , defined by q ± j = x j ± (−1) j y j , are critical for L(q j , q j+1 , h).
Proof. The discrete Euler-Lagrange equations for L(x j , y j , x j+1 , y j+1 , h) are Taking the sum resp. the difference of these equations we find Depending on the parity of j, either the first or the second of those equations is Hence (x j , y j ) j∈Z satisfies the Euler-Lagrange equations for L(x j , y j , x j+1 , y j+1 , h) if and only if (q + j ) j∈Z and (q − j ) j∈Z satisfy the Euler-Lagrange equation for L(q j , q j+1 , h).
be the k-th truncation of the principal modified equation for the difference equation described by the discrete Lagrangian L from Proposition 5.1. Then (20) is the k-th truncation of the full system of modified equations for the variational integrator described by L.
Proof. Let (x(t), y(t)) be a solution of the system (20). By definition of the principal modified equation, the discrete curve x(t + jh), y(t + jh) j∈Z satisfies the discrete Euler-Lagrange equations for L up to a truncation error for any choice of t. Hence, by Proposition 5.1, the discrete curve satisfies the discrete Euler-Lagrange equations for L up to a truncation error. This is the defining property of the system of modified equations, see Definition 3.3(b).

Corollary 5.3. Up to a truncation error of arbitrarily high order, the full system of modified equations (20) for a variational integrator is Lagrangian.
Let us illustrate this construction by applying it to our two methods.

Midpoint rule
We have Hence This is also the leading order term of the modified Lagrangian, L mod,0 (x, y,ẋ,ẏ, h). If α is linear, its Euler-Lagrange equations arė Since y is constant in leading order, we need to look at higher order terms to determine whether parasitic solutions occur. No higher order terms of the modified Lagrangian contain y itself, and those terms that contain derivatives of y are at least quadratic in the derivatives of y. From these observations one can deduce that the parasitic modified equation isẏ = 0 to any order of accuracy. It follows that the parasitic oscillations are of constant magnitude. Hence if the initialization of the discrete system is close to a solution of the principal modified equation, then the discrete solution will remain close to it.

Trapezoidal rule
We have Hence . Separating the alternating terms from the rest, we finḋ Unsurprisingly, we find the same system of modified equations as with the Lagrangian method.

Trapezoidal rule
Now we consider the difference equation and make the same identifications as before. We finḋ . If we assume that y = O(h), then the system of modified equations iṡ

Examples
To illustrate the theory above, we apply our two integrators to two examples. Since the calculations tend to be quite long in real-world problems, we start with a minimal toy problem. After that, we discuss the dynamics of point vortices in the plane.

Toy Problem
Consider the Lagrangian on T R 2 . Its Euler-Lagrange equations arė As a concrete example, the choice V (q) = − cos(q) and U (p) = 1 2 p 2 describes the pendulum.
We have L disc (p j , q j , p j+1 , q j+1 , h) = 1 2 This corresponds to the following system of difference equations: By Taylor expansion we obtain It follows that Its Euler-Lagrange equations are Solving for q andq we find the principal modified equationṡ Eliminating higher derivatives in L mesh we find As discussed in the previous section we do not expect parasitic solutions with this method (see Figure 1).

Trapezoidal rule
We have L disc (p j , q j , p j+1 , q j+1 , h) = 1 2 The corresponding discrete Euler-Lagrange equations are By Taylor expansion we obtain It follows that Its Euler-Lagrange equations are Solving for q andq we find the principal modified equationṡ Eliminating higher derivatives in L mesh we find For the pendulum, V (q) = − cos(q) and U (p) = 1 2 p 2 , we have hence the matrix in Equation (21) is This matrix has a pair of real eigenvalues if cos(q) < 0 and a pair of purely imaginary eigenvalues if cos(q) > 0. This suggests (but does not prove; q is not constant) that exponentially growing parasites occur in the regions where cos(q) < 0.
In the top right image of Figure 1 one clearly observes parasitic solutions for this method. Note the parasites only seem to grow where |q| > π 2 , i.e. where cos(q) < 0. In the region where |q| < π 2 there is no noticeable growth in the amplitude of the oscillations. Instead we observe a rotation in the direction of the oscillations, as expected when the eigenvalues are purely imaginary. This is visualized in Figure 1 by line segments connecting the points of the discrete solution with the corresponding points on the solution of the principal modified equation.
When the initial conditions are chosen such that q remains in the stable region |q| < π 2 no parasites are observed (bottom right image of Figure 1), even if the simulation is continued for many periods (not pictured).

Point vortices
Our second example involves vortices on a planar surface. If all vorticity is contained in a finite number of points, then the movement of those points is described by first order ODEs [10,12]. To be precise, the dynamics of N point vortices in the (complex) plane is described by the Lagrangian where z j and Γ j are the position and circulation of the j-th vortex, and the bar denotes the complex conjugate. The equations of motion arė It follows thatz

Midpoint rule
We have To obtain the modified Lagrangian we evaluate the second derivatives in L mesh using the leading order equation (22). We find Therefore, L mod (z, z,ż,ż, h) = L(z, z,ż,ż) + h 2 24

Trapezoidal rule
For the Trapezoidal rule, we find in the same way that

Legend:
Dashed curves exact solution.
Bullets discrete solution. Solid curves solution of the truncated principal modified equation.

Conclusion
We have described a Lagrangian algorithm to calculate the modified equation of a variational integrator applied to a degenerate continuous Lagrangian that is linear in the velocities. To obtain the principal modified equation this was a straightforward adaptation of the procedure developed for non-degenerate Lagrangians. To obtain the full system of modified equations we doubled the dimension of the discrete system in a suitable way. As a consequence, we proved that the system of modified equations is variational. We have illustrated the construction of modified Lagrangians and the possible issue of parasitic solutions with examples. Our construction is potentially useful to create more accurate variational integrators, for example in the spirit of [2,3].