R-adaptive multisymplectic and variational integrators

Moving mesh methods (also called r-adaptive methods) are space-adaptive strategies used for the numerical simulation of time-dependent partial differential equations. These methods keep the total number of mesh points fixed during the simulation, but redistribute them over time to follow the areas where a higher mesh point density is required. There are a very limited number of moving mesh methods designed for solving field-theoretic partial differential equations, and the numerical analysis of the resulting schemes is challenging. In this paper we present two ways to construct r-adaptive variational and multisymplectic integrators for (1+1)-dimensional Lagrangian field theories. The first method uses a variational discretization of the physical equations and the mesh equations are then coupled in a way typical of the existing r-adaptive schemes. The second method treats the mesh points as pseudo-particles and incorporates their dynamics directly into the variational principle. A user-specified adaptation strategy is then enforced through Lagrange multipliers as a constraint on the dynamics of both the physical field and the mesh points. We discuss the advantages and limitations of our methods. Numerical results for the Sine-Gordon equation are also presented.


Introduction
The purpose of this work is to design, analyze, and implement variational and multisymplectic integrators for Lagrangian partial differential equations with space-adaptive meshes. In this paper, we combine geometric numerical integration and r-adaptive methods for the numerical solution of Partial Differential Equations (PDEs). We show that these two fields are compatible, mostly due to the fact that, in r-adaptation, the number of mesh points remains constant and we can treat them as additional pseudo-particles of which the dynamics are coupled to the dynamics of the physical field of interest.
Geometric (or structure-preserving) integrators are numerical methods that preserve geometric properties of the flow of a differential equation (see Reference [1]). This encompasses symplectic integrators for Hamiltonian systems, variational integrators for Lagrangian systems, and numerical methods on manifolds, including Lie group methods and integrators for constrained mechanical systems. Geometric integrators proved to be extremely useful for numerical computations in astronomy, molecular dynamics, mechanics, and theoretical physics. The main motivation for developing structure-preserving algorithms lies in the fact that they show excellent numerical behavior, especially for long-time integration of equations possessing geometric properties.
An important class of structure-preserving integrators are variational integrators for Lagrangian systems [1,2]. This type of integrator is based on discrete variational principles. The variational approach provides a unified framework for the analysis of many symplectic algorithms and is characterized by a natural treatment of the discrete Noether theorem, as well as forced, dissipative, and constrained systems. Variational integrators were first introduced in the context of finite-dimensional mechanical systems, but later, Marsden, Patrick, and Shkoller [3] generalized this idea to field theories. Variational integrators have since been successfully applied in many computations, for example in elasticity [4], electrodynamics [5], or fluid dynamics [6]. Existing variational integrators so far have been developed on static, mostly uniform spatial meshes. The main goal of this paper is to design and analyze variational integrators that allow for the use of space-adaptive meshes.
Adaptive meshes used for the numerical solution of partial differential equations fall into three main categories: h-adaptive, p-adaptive, and r-adaptive. R-adaptive methods, which are also known as moving mesh methods [7,8], keep the total number of mesh points fixed during the simulation but relocate them over time. These methods are designed to minimize the error of the computations by optimally distributing the mesh points, contrasting with h-adaptive methods for which the accuracy of the computations is obtained via insertion and deletion of mesh points. Moving mesh methods are a large and interesting research field of applied mathematics, and their role in modern computational modeling is growing. Despite the increasing interest in these methods in recent years, they are still in a relatively early stage of their development compared to the more matured h-adaptive methods.

Overview
There are three logical steps to r-adaptation: • Discretization of the physical PDE The key ideas of this paper regard the first and the last step. Following the general spirit of variational integrators, we discretize the underlying action functional rather than the PDE itself and then derive the discrete equations of motion. We base our adaptation strategies on the equidistribution principle and the resulting moving mesh partial differential equations (MMPDEs). We interpret MMPDEs as constraints, which allows us to consider novel ways of coupling them to the physical equations. Note that we will restrict our explanations to one time and one space dimension for the sake of simplicity.
Let us consider a (1+1)-dimensional scalar field theory with the action functional where φ : [0, X max ] × [0, T max ] −→ R is the field and L : R × R × R −→ R its Lagrangian density. For simplicity, we assume the following fixed boundary conditions In order to further consider moving meshes let us perform a change of variables X = X(x, t) such that for all t the map X(., t) : [0, X max ] −→ [0, X max ] is a "diffeomorphism"-more precisely, we only require that X(., t) is a homeomorphism such that both X(., t) and X(., t) −1 are piecewise C 1 . In the context of mesh adaptation, the map X(x, t) represents the spatial position at time t of the mesh point labeled by x. Define ϕ(x, t) = φ(X(x, t), t). Then, the partial derivatives of φ are φ X (X(x, t), t) = ϕ x /X x and φ t (X(x, t), t) = ϕ t − ϕ x X t /X x . Plugging these equations in Equation (1), we get where the last equality defines two modified, or "reparametrized", action functionals. For the first one, S is considered as a functional of ϕ only, whereas in the second one we also treat it as a functional of X. This leads to two different approaches to mesh adaptation, which we dub the control-theoretic strategy and the Lagrange multiplier strategy, respectively. The "reparametrized" field theories defined byS [ϕ] andS[ϕ, X] are both intrinsically covariant; however, it is convenient for computational purposes to work with a space-time split and to formulate the field dynamics as an initial value problem.

Outline
This paper is organized as follows. In Sections 2 and 3 we take the view of infinite dimensional manifolds of fields as configuration spaces and develop the control-theoretic and Lagrange multiplier strategies in that setting. It allows us to discretize our system in space first and consider time discretization later on. It is clear from our exposition that the resulting integrators are variational. In Section 4, we show how similar integrators can be constructed using the covariant formalism of multisymplectic field theory. We also show how the integrators from the previous sections can be interpreted as multisymplectic. In Section 5, we apply our integrators to the Sine-Gordon equation and we present our numerical results. We summarize our work in Section 6 and discuss several directions in which it can be extended.

Control-Theoretic Approach to r-Adaptation
At first glance, it appears that the simplest and most straightforward way to construct an r-adaptive variational integrator would be to discretize the physical system in a similar manner to the general approach to variational integration, i.e., discretize the underlying variational principle, then derive the mesh equations, and couple them to the physical equations in a way typical of the existing r-adaptive algorithms. We explore this idea in this section and show that it indeed leads to space adaptive integrators that are variational in nature. However, we also show that those integrators do not exhibit the behavior expected of geometric integrators, such as good energy conservation. We will refer to this strategy as control-theoretic, since, in this description, the field ϕ represents the physical state of the system, while X can be interpreted as a control variable and the mesh equations can be interpreted as feedback (see, e.g., Reference [9]).

Reparametrized Lagrangian
For the moment, let us assume that X(x, t) is a known function. We denote by ξ(X, t) the function such that ξ(., t) = X(., t) −1 , that is ξ(X(x, t), t) = x (We allow a little abuse of notation here: X denotes both the argument of ξ and the change of variables X(x, t). If we wanted to be more precise, we would write X = h(x, t).). We thus haveS[ϕ] = S[ϕ(ξ(X, t), t)].
The corresponding instantaneous LagrangianL : (ϕ, ϕ x , ϕ t , t) dx (5) with the Lagrangian densityL The function spaces Q and W must be chosen appropriately for the problem at hand so that the Lagrangian in Equation (5) makes sense. For instance, for a free field, we will have Q = H 1 ([0, X max ]) and W = L 2 ([0, X max ]). Since X(x, t) is a function of t, we are looking at a time-dependent system. Even though the energy associated with the Lagrangian in Equation (5) is not conserved, the energy of the original theory associated with the action functional in Equation (1) is conserved. To see this, note that if φ(X, t) extremizes S[φ], then dE/dt = 0 (computed from Equation (7)). Trivially, this means that dE/dt = 0 when Equation (8) is invoked as well. Moreover, as we have noted earlier, φ(X, t) extremizes S[φ] iff ϕ(x, t) extremizesS [ϕ]. This means that the energy (Equation (8)) is constant on solutions of the reparametrized theory.

Spatial Finite Element Discretization
We begin with a discretization of the spatial dimension only, thus turning the original infinite-dimensional problem into a time-continuous finite-dimensional Lagrangian system. Let ∆x = X max /(N + 1), and define the reference uniform mesh x i = i · ∆x for i = 0, 1, ..., N + 1 and the corresponding piecewise linear finite elements We now restrict X(x, t) to be of the form with X 0 (t) = 0, X N+1 (t) = X max , and arbitrary X i (t), i = 1, 2, ..., N as long as X(., t) is a homeomorphism for all t. In our context of numerical computations, the functions X i (t) represent the current position of the ith mesh point. Define the finite element spaces Q N = W N = span(η 0 , ..., η N+1 ) (11) and assume that Q N ⊂ Q, W N ⊂ W. Let us denote a generic element of Q N by ϕ and a generic element of W N byφ. We have the decompositions The numbers (y i ,ẏ i ) thus form natural (global) coordinates on Q N × W N . We can now approximate the dynamics of the system described by Equation (5) in the finite-dimensional space Q N × W N . Let us consider the restrictionL N =L| Q N ×W N ×R of the Lagrangian of Equation (5) to Q N × W N × R. In the chosen coordinates, we havẽ Note that, given the boundary conditions (Equation (2)), y 0 , y N+1 ,ẏ 0 , andẏ N+1 are fixed. We will thus no longer write them as arguments ofL N .
The advantage of using a finite element discretization lies in the fact that the symplectic structure induced on Q N × W N byL N is strictly a restriction (i.e., a pull-back) of the (pre-)symplectic structure (In most cases, the symplectic structure of (Q × W,L) is only weakly nondegenerate; see Reference [10]) on Q × W. This establishes a direct link between symplectic integration of the finite-dimensional mechanical system (Q N × W N ,L N ) and the infinite-dimensional field theory (Q × W,L).

Differential-Algebraic Formulation and Time Integration
We now consider time integration of the Lagrangian system (Q N × W N ,L N ). If the functions X i (t) are known, then one can perform variational integration in the standard way, that is, define the discrete LagrangianL d : R × Q N × R × Q N → R and solve the corresponding discrete Euler-Lagrange equations (see References [1,2]). Let t n = n · ∆t for n = 0, 1, 2, . . . be an increasing sequence of times and {y 0 , y 1 , . . .} be the corresponding discrete path of the system in Q N . The discrete Lagrangian L d is an approximation to the exact discrete Lagrangian L E d , such that where y n = (y n 1 , ..., y n N ), y n+1 = (y n+1 1 , ..., y n+1 N ), and y(t) are the solutions of the Euler-Lagrange equations corresponding toL N with boundary values y(t n ) = y n and y(t n+1 ) = y n+1 . Depending on the quadrature we use to approximate the integral in Equation (14), we obtain different types of variational integrators. As will be discussed below, in r-adaptation, one has to deal with stiff differential equations or differential-algebraic equations; therefore, higher order implicit integration in time is advisable (see References [11,12]). We will employ variational partitioned Runge-Kutta methods. An s-stage Runge-Kutta method is constructed by choosing where t i = t n + c i (t n+1 − t n ); the right-hand side is extremized under the constraint It can be shown that the variational integrator with the discrete Lagrangian defined by Equation (15) is equivalent to an appropriately chosen symplectic partitioned Runge-Kutta method applied to the Hamiltonian system corresponding toL N (see Reference [1,2]). With this in mind, we turn our semi-discrete Lagrangian system (Q N × W N ,L N ) into the Hamiltonian system (Q N × W * N ,H N ) via the standard Legendre transform H N (y 1 , ..., y N , p 1 , ..., p N ; X 1 , ..., X N ,Ẋ 1 , ..., where p i = ∂L N /∂ẏ i and we explicitly stated the dependence on the positions X i and velocitiesẊ i of the mesh points. The Hamiltonian equations take the form (It is computationally more convenient to directly integrate the implicit Hamiltonian system p i = ∂L N /∂ẏ i ,ṗ i = ∂L N /∂y i , but as long as the system described by Equation (1) is at least weakly nondegenerate, there is no theoretical issue with passing to the Hamiltonian formulation, which we do for the clarity of our exposition.) Suppose that the function X i (t) is C 1 and that H N is smooth as a function of the y i 's, p i 's, X i 's, andẊ i 's (note that these assumptions are used for simplicity and can be easily relaxed, if necessary, depending on the regularity of the considered Lagrangian system). Then, the assumptions of Picard's theorem are satisfied and there exists a unique (17). This flow is symplectic.
However, in practice, we do not know the X i 's and we in fact would like to be able to adjust them "on the fly", based on the current behavior of the system. We will do that by introducing additional constraint functions g i (y 1 , ..., y N , X 1 , ..., X N ) and demanding that the conditions g i = 0 be satisfied at all times (In the context of Control Theory, the constraints g i = 0 are called strict static state feedback. See Reference [9].). The choice of these functions will be discussed in Section 2.4. This leads to the following system of differential-algebraic equations (DAEs) of index 1 (see References [11][12][13]) for i = 1, ..., N. Note that an initial condition for X is fixed by the constraints. This system is of index 1 because one has to differentiate the algebraic equations with respect to time once in order to reduce it to an implicit system of ordinary differential equations (ODEs). In fact, the implicit system will take the formẏ = ∂H N ∂p y, p; X,Ẋ , where X (0) is a vector of arbitrary initial condition for the X i 's. Suppose again that H N is a smooth function of y, p, X, andẊ. Futhermore, suppose that g is a C 1 function of y and X and that ∂g ∂X − ∂g ∂y ∂ 2 H N ∂Ẋ∂p is invertible with its inverse bounded in a neighborhood of the exact solution. (Again, these assumptions can be relaxed if necessary.) Then, by the Implicit Function Theorem, Equation (19) can be solved explicitly forẏ,ṗ, andẊ and the resulting explicit ODE system will satisfy the assumptions of Picard's theorem. Let (y(t), p(t), X(t)) be the unique C 1 solution to this ODE system (and hence to Equation (19)). We have the trivial result Proposition 2. If g(y (0) , X (0) ) = 0, then (y(t), p(t), X(t)) is a solution to Equation (18). (Note that there might be other solutions, as for any given y (0) , there might be more than one X (0) that solves the constraint equations.) In practice, we would like to integrate Equation (18). A question arises in what sense is the system defined by this equation symplectic and in what sense a numerical integration scheme for this system can be regarded as variational. Let us address these issues. Proposition 3. Let (y(t), p(t), X(t)) be a solution to Equation (18) and use this X(t) to form the Hamiltonian system (17). Then, we have where F t 0 ,t (ŷ,p) is the symplectic flow for Equation (17).
Suppose we now would like to find a numerical approximation of the solution to Equation (17) using an s-stage partitioned Runge-Kutta method with coefficients a ij , b i ,ā ij ,b i , and c i [1,14]. The numerical scheme will take the forṁ where Y i ,Ẏ i , P i , andṖ i are the internal stages and ∆t is the integration timestep. Let us apply the same partitioned Runge-Kutta method to Equation (18). In order to compute the internal stages Q i ,Q i of the X variable, we use the state-space form approach, that is, we demand that the constraints and their time derivatives be satisfied (see Reference [12]). The new step value X n+1 is computed by solving the constraints as well. The resulting numerical scheme is thuṡ We have the following trivial observation.

Proposition 4.
If X(t) is defined to be a C 1 interpolation of the internal stages Q i ,Q i at times t n + c i ∆t (that is, if the values X(t n + c i ∆t) andẊ(t n + c i ∆t) coincide with Q i andQ i ), then the schemes defined by Equations (20) and (21) give the same numerical approximations y n and p n to the exact solution y(t) and p(t).
Intuitively, Proposition 4 states that we can apply a symplectic partitioned Runge-Kutta method to the DAE system in Equation (18), which solves both for X(t) and (y(t), p(t)), and the result will be the same as if we performed a symplectic integration of the Hamiltonian system in Equation (17) for (y(t), p(t)) with a known X(t).

Moving Mesh Partial Differential Equations
The concept of equidistribution is the most popular paradigm of r-adaptation (see References [7,8]). Given a continuous mesh density function ρ(X), the equidistribution principle seeks to find a mesh 0 = X 0 < X 1 < ... < X N+1 = X max such that the following holds that is, the quantity represented by the density function is equidistributed among all cells. In the continuous setting, we will say that the reparametrization where σ = X max 0 ρ(X) dX is the total amount of the equidistributed quantity. Differentiate this equation with respect to x to obtain It is still a global condition in the sense that σ has to be known. For computational purposes, it is convenient to differentiate this relation again and to consider the following partial differential equation with the boundary conditions X(0) = 0, X(X max ) = X max . The choice of the mesh density function ρ(X) is typically problem-dependent and the subject of much research. A popular example is the generalized solution arclength given by It is often used to construct meshes that can follow moving fronts with locally high gradients [7,8]. With this choice, Equation (25) is equivalent to assuming X x > 0, which we demand anyway. A finite difference discretization on the mesh x i = i · ∆x gives us the set of constraints with the previously defined y i 's and X i 's. This set of constraints can be used in Equation (18).

Example
To illustrate these ideas, let us consider the Lagrangian density The reparametrized Lagrangian in Equation (5) takes the form Let N = 1 and φ L = φ R = 0. Then The semi-discrete Lagrangian is The Legendre transform gives p 1 = ∂L N /∂ẏ 1 = X maxẏ1 /3; hence, the semi-discrete Hamiltonian isH The corresponding DAE system iṡ This system is to be solved for the unknown functions y 1 (t), p 1 (t), and X 1 (t). It is of index 1 because we have three unknown functions and only two differential equations-the algebraic equation has to be differentiated once in order to obtain a missing ODE.

Backward Error Analysis
The true power of symplectic integration of Hamiltonian equations is revealed through backward error analysis: It can be shown that a symplectic integrator for a Hamiltonian system with the Hamiltonian H(q, p) defines the exact flow for a nearby Hamiltonian system, of which the Hamiltonian can be expressed as the asymptotic series H (q, p) = H(q, p) + ∆tH 2 (q, p) + ∆t 2 H 3 (q, p) + . . . (35) Owing to this fact, under some additional assumptions, symplectic numerical schemes nearly conserve the original Hamiltonian H(q, p) over exponentially long time intervals. See Reference [1] for details.
Let us briefly review the results of backward error analysis for the integrator defined by Equation (21). Suppose g(y, X) satisfies the assumptions of the Implicit Function Theorem. Then, at least locally, we can solve the constraint X = h(y). The Hamiltonian DAE system in Equation (18) can be then written as the following (implicit) ODE system for y and ṗ Since we used the state-space formulation, the numerical scheme of Equation (21) is equivalent to applying the same partitioned Runge-Kutta method to Equations (36), that is, we have We computed the corresponding modified equation for several symplectic methods, namely Gauss and Lobatto IIIA-IIIB quadratures. Unfortunately, none of the quadratures resulted in a form akin to Equation (36) for some modified Hamiltonian functionH N related toH N by a series similar to Equation (35). This hints at the fact that we should not expect this integrator to show excellent energy conservation over long integration times. One could also consider the implicit ODE system of Equation (19), which has an obvious triple partitioned structure and could apply a different Runge-Kutta method to each variable y, p, and X. Although we did not pursue this idea further, it seems unlikely it would bring a desirable result.
We therefore conclude that the control-theoretic strategy, while yielding a perfectly legitimate numerical method, does not take the full advantage of the underlying geometric structures. Let us point out that, while we used a variational discretization of the governing physical PDE, the mesh equations were coupled in a manner that is typical of the existing r-adaptive methods (see References [7,8]). We now turn our attention to a second approach, which offers a novel way of coupling the mesh equations to the physical equations.

Lagrange Multiplier Approach to r-Adaptation
As we saw in Section 2, discretization of the variational principle alone is not sufficient if we would like to accurately capture the geometric properties of the physical system described by the action functional in Equation (1). In this section, we propose a new technique of coupling the mesh equations to the physical equations. Our idea is based on the observation that, in r-adaptation, the number of mesh points is constant; therefore, we can treat them as pseudo-particles and we can incorporate their dynamics into the variational principle. We show that this strategy results in integrators that much better preserve the energy of the considered system.

Reparametrized Lagrangian
In this approach, we treat X(x, t) as an independent field, that is, another degree of freedom, and we will treat the "modified" action defined by Equation (3) as a functional of both ϕ and X: S =S[ϕ, X]. For the purpose of the derivations below, we assume that ϕ(., t) and X(., t) are continuous and piecewise C 1 . One could consider the closure of this space in the topology of either the Hilbert or Banach spaces of sufficiently integrable functions and interpret differentiation in a sufficiently weak sense, but this functional-analytic aspect is of little importance for the developments in this section. We refer the interested reader to References [15,16]. As in Section 2.1, let ξ(X, t) be the function such that ξ(., t) = X(., t) −1 , that is ξ(X(x, t), t) = x. Then,S[ϕ, X] = S[ϕ(ξ(X, t), t)]. We begin with two propositions and one corollary which will be important for the rest of our exposition.

Proposition 5.
Extremizing S[φ] with respect to φ is equivalent to extremizingS[ϕ, X] with respect to both ϕ and X.

Spatial Finite Element Discretization
The Lagrangian of the "reparametrized" theoryL : has the same form as in Equation (5) (we only treat it as a functional of X and X t as well), where Q, G, W, and Z are spaces of continuous and piecewise C 1 functions, as mentioned before. We again let ∆x = X max /(N + 1) and define the uniform mesh x i = i · ∆x for i = 0, 1, ..., N + 1. Define the finite element spaces where we used the finite elements introduced in Equation (9). We have Q N ⊂ Q, G N ⊂ G, W N ⊂ W, and Z N ⊂ Z. In addition to Equation (12), we also consider The numbers (y i , We again consider the restricted LagrangianL N =L| Q N ×G N ×W N ×Z N . In the chosen coordinates where ϕ(x), X(x),φ(x),Ẋ(x) are defined by Equations (12) and (40). Once again, we refrain from writing y 0 , y N+1 ,ẏ 0 ,ẏ N+1 , X 0 , X N+1 ,Ẋ 0 , andẊ N+1 as arguments ofL N in the remainder of this section, as those are not actual degrees of freedom.

Invertibility of the Legendre Transform
For simplicity, let us restrict our considerations to Lagrangian densities of the form We chose a kinetic term that is most common in applications. The corresponding "reparametrized" Lagrangian isL where we kept only the terms that involve the velocities ϕ t and X t . The semi-discrete Lagrangian becomes Let us define the conjugate momenta via the Legendre Transform This can be written as  where the 2N × 2N mass matrixM N (y, X) has the following block tridiagonal structurẽ with the 2 × 2 blocks where From now on, we will always assume δ i > 0, as we demand that Proposition 7. The mass matrixM N (y, X) is non-singular almost everywhere (as a function of the y i 's and X i 's) and singular iff γ i−1 = γ i for some i.
Proof. We will compute the determinant ofM N (y, X) by transforming it into a block upper triangular form by zeroing the blocks B i below the diagonal. Let us start with the block B 1 . We use linear combinations of the first two rows of the mass matrix to zero the elements of the block B 1 below the diagonal. Suppose γ 0 = γ 1 . Then, it is easy to see that the first two rows of the mass matrix are not linearly independent, so the determinant of the mass matrix is zero. Assume γ 0 = γ 1 . Then by Equation (50), the block A 1 is invertible. We multiply the first two rows of the mass matrix by B 1 A −1 1 and subtract the result from the third and fourth rows. This zeroes the block B 1 below the diagonal and replaces the block A 2 by We now zero the block B 2 below the diagonal in a similar fashion. After n − 1 steps of this procedure, the mass matrix is transformed into In a moment we will see that C n is singular iff γ n−1 = γ n , and in that case, the two rows of the matrix above that contain C n and B n are linearly dependent, thus making the mass matrix singular. Suppose γ n−1 = γ n , so that C n is invertible. In the next step of our procedure the block A n+1 is replaced by Together with the condition C 1 = A 1 , this gives us a recurrence. By induction on n, we find that and which justifies our assumptions on the invertibility of the blocks C i . We can now express the determinant of the mass matrix as det We see that the mass matrix becomes singular iff γ i−1 = γ i for some i, and this condition defines a measure zero subset of R 2N .
Remark 1. This result shows that the finite-dimensional system described by the semi-discrete Lagrangian of Equation (44) is non-degenerate almost everywhere. This means that, unlike in the continuous case, the Euler-Lagrange equations corresponding to the variations of the y i 's and X i 's are independent of each other (almost everywhere) and the equations corresponding to the X i 's are in fact necessary for the correct description of the dynamics. This can also be seen in a more general way. Owing to the fact we are considering a finite element approximation, the semi-discrete action functionalS N is simply a restriction ofS, and therefore, Equation (37) still holds. The corresponding Euler-Lagrange equations take the form which must hold for all variations δϕ( . Since we are working in a finite dimensional subspace, the second equation now does not follow from the first equation. To see this, consider a particular variation δX( which is discontinuous at x = x k and cannot be expressed as ∑ N i=1 δy i (t)η i (x) for any δy i (t) unless γ k−1 = γ k . Therefore, we cannot invoke the first equation to show that δ 2S [ϕ, X] · δX(x, t) = 0. The second equation becomes independent.

Remark 2.
It is also instructive to realize what exactly happens when γ k−1 = γ k . This means that, locally in the interval [X k−1 , X k+1 ], the field φ(X, t) is a straight line with the slope γ k . It also means that there are infinitely many values (X k , y k ) that reproduce the same local shape of φ(X, t). This reflects the arbitrariness of X(x, t) in the infinite-dimensional setting. In the finite element setting, however, this holds only when the points (X k−1 , y k−1 ), (X k , y k ), and (X k+1 , y k+1 ) line up. Otherwise any change to the middle point changes the shape of φ(X, t). See Figure 1.
then any change to the middle point changes the local shape of φ(X, t).
(Right) If γ k−1 = γ k , then there are infinitely many possible positions for (X k , y k ) that reproduce the local linear shape of φ(X, t).

Existence and Uniqueness of Solutions
Since the Legendre Transform in Equation (46) becomes singular at some points, this raises a question about the existence and uniqueness of the solutions to the Euler-Lagrange Equation (57). In this section, we provide a partial answer to this problem. We will begin by computing the Lagrangian symplectic formΩ where p i and S i are given by Equation (45). For notational convenience, we will collectively denote q = (y 1 , , the symplectic form can be represented by the matrix where the 2N × 2N block∆ N (q,q) has the further block tridiagonal structure In this form, it is easy to see that so the symplectic form is singular whenever the mass matrix is.
The energy corresponding to the Lagrangian in Equation (44) can be written as In the chosen coordinates, dẼ N can be represented by the row vector dẼ N = (∂Ẽ N /∂q 1 , ..., ∂Ẽ N /∂q 2N ). It turns out that where the vector ξ has the following block structure Each of these blocks has the form ξ k = (ξ k,1 , ξ k,2 ) T . Through basic algebraic manipulations and integration by parts, one finds that and We are now ready to consider the generalized Hamiltonian equation which we solve for the vector field Equations of this form are called (quasilinear) implicit ODEs (see References [17,18]). If the symplectic form is non-singular in a neighborhood of (q (0) ,q (0) ), then the equation can be solved directly via Z = [Ω T N (q,q)] −1 dẼ T N (q,q) to obtain the standard explicit ODE form, and standard existence/uniqueness theorems (Picard's, Peano's, etc.) of ODE theory can be invoked to show local existence and uniqueness of the flow of Z in a neighborhood of (q (0) ,q (0) ). If, however, the symplectic form is singular at (q (0) ,q (0) ), then there are two possibilities. The first case is and it means there is no solution for Z at (q (0) ,q (0) ). This type of singularity is called an algebraic one, and it leads to so-called impasse points (see References [17][18][19]). The other case is and it means that there exists a nonunique solution Z at (q (0) ,q (0) ). This type of singularity is called a geometric one. If (q (0) ,q (0) ) is a limit of regular points of Equation (70) (i.e., points where the symplectic form is nonsingular), then there might exist an integral curve of Z passing through (q (0) ,q (0) ). See References [17][18][19][20][21][22][23] for more details.

Proposition 8.
The singularities of the symplectic formΩ N (q,q) are geometric.
Proof. Suppose that the mass matrix (and thus the symplectic form) is singular at (q (0) ,q (0) ). Using the block structures described by Equations (60) and (65), we can write Equation (70) as the system The second equation implies that there exists a solution α =q (0) . In fact, this is the only solution we are interested in, since it satisfies the second order condition: The Euler-Lagrange equations underlying the variationl principle are second order, so we are only interested in solutions of the form The first equation can be rewritten as Since the mass matrix is singular, we must have γ k−1 = γ k for some k. As we saw in Section 3.3, this means that the two rows of the kth "block row" of the mass matrix (i.e., the rows containing the blocks B k−1 , A k , and B k ) are not linearly independent. In fact, we have where a m * denotes the mth row of the matrix a. Equation (74) will have a solution for β iff the right-hand side satisfies a similar scaling condition in the the kth "block element". Using Equations (62), (67) and (68), we show that −ξ −∆ Nq (0) indeed has this property. Hence, 0) ) is a limit of regular points. (70)) has to deal with the singularity points of the symplectic form. While there are some numerical algorithms allowing one to get past singular hypersurfaces (see Reference [17]), it might not be very practical from the application point of view. Note that, unlike in the continuous case, the time evolution of the meshpoints X i 's is governed by the equations of motion, so the user does not have any influence on how the mesh is adapted. More importantly, there is no built-in mechanism that would prevent mesh tangling. Some preliminary numerical experiments show that the mesh points eventually collapse when started with nonzero initial velocities. (47) bear some similarities to the singularities of the mass matrices encountered in the Moving Finite Element method. In References [24,25], the authors proposed introducing a small "internodal" viscosity which penalizes the method for relative motion between the nodes and thus regularizes the mass matrix. A similar idea could be applied in our case: One could add some small ε kinetic terms to the Lagrangian in Equation (44) in order to regularize the Legendre Transform. In light of the remark made above, we did not follow this idea further and decided to take a different route instead, as described in the following sections. However, investigating further similarities between our variational approach and the Moving Finite Element method might be worthwhile. There also might be some connection to the r-adaptive method presented in Reference [26]: The evolution of the mesh in that method is also set by the equations of motion, although the authors considered a different variational principle and different theoretical reasoning to justify the validity of their approach.

Constraints and Adaptation Strategy
As we saw in Section 3.4, upon discretization, we lose the arbitrariness of X(x, t) and the evolution of X i (t) is governed by the equations of motion, while we still want to be able to select a desired mesh adaptation strategy, like Equation (28). This could be done by augmenting the Lagrangian in Equation (44) with Lagrange multipliers corresponding to each constraint g i . However, it is not obvious that the dynamics of the constrained system as defined would reflect in any way the behavior of the approximated system with the Lagrangian density given in Equation (42). We will show that the constraints can be added via Lagrange multipliers already at the continuous level (Equation (42)) and the continuous system as defined can be then discretized to arrive at the Lagrangian of Equation (44) with the desired adaptation constraints.

Global Constraint
As mentioned before, eventually we would like to impose the constraints on the semi-discrete system defined by the Lagrangian in Equation (44). Let us assume that g : .., g N ) T is C 1 and 0 is a regular value of g, so that Equation (76) defines a submanifold. To see how these constraints can be introduced at the continuous level, let us select uniformly distributed points x i = i · ∆x, i = 0, ..., N + 1, and ∆x = X max /(N + 1) and demand that the constraints be satisfied by ϕ(x, t) and X(x, t). One way of imposing these constraints is solving the system This system consists of one Euler-Lagrange equation that corresponds to extremizingS with respect to ϕ (we saw in Section 3.1 that the other Euler-Lagrange equation is not independent) and a set of constraints enforced at some preselected point x i . Note, that upon finite element discretization on a mesh coinciding with the preselected points, this system reduces to the approach presented in Section 2: We minimize the discrete action with respect to the y i 's only and supplement the resulting equations with the constraints of Equation (76).
Another way that we want to explore consists in using Lagrange multipliers. Define the auxiliary action functional as We will assume that the Lagrange multipliers λ i (t) are at least continuous in time. According to the method of Lagrange multipliers, we seek the stationary points ofS C . This leads to the following system of equations: where, for clarity, we suppressed writing the arguments of Equation (78) is more intuitive because we directly use the arbitrariness of X(x, t) and simply restrict it further by imposing constraints. It is not immediately obvious how the solutions of Equations (78) and (80) relate to each other. We would like both systems to be "equivalent" in some sense or at least their solution sets to overlap. Let us investigate this issue in more detail.
Suppose (ϕ, X) satisfies Equation (78). Then, it is quite trivial to see that (ϕ, X, λ 1 , ..., λ N ) such that λ k ≡ 0 satisfies Equation (80): The second equation is implied by the first one and the other equations coincide with those of Equation (78). At this point, it should be obvious that Equation (80) may have more solutions for ϕ and X than Equation (78).
Proof. Suppose (ϕ, X, λ 1 , ..., λ N ) satisfy both Equations (78) and (80). Equation (78) implies that δ 1S · δϕ = 0 and δ 2S · δX = 0. Using this in Equation (80) gives In particular, this has to hold for variations δϕ and δX such that are continuous and we get Since we assumed that 0 is a regular value of g and the constraint g = 0 is satisfied by ϕ and X, we have that for all t the matrix Dg has full rank-that is, there exists a non-singular We see that considering Lagrange multipliers in Equation (79) makes sense at the continuous level. We can now perform a finite element discretization. The auxiliary LagrangianL C : Q × G × W × Z × R N −→ R corresponding to Equation (79) can be written as whereL is the Lagrangian of the unconstrained theory and has been defined by Equation (38). Let us choose a uniform mesh coinciding with the preselected points x i . As in Section 3.2, we consider the restrictionL CN =L C | Q N ×G N ×W N ×Z N ×R N and we get We see that the semi-discrete LagrangianL CN is obtained from the semi-discrete LagrangianL N by adding the constraints g i directly at the semi-discrete level, which is exactly what we set out to do at the beginning of this section. However, in the semi-discrete setting, we cannot expect the Lagrange multipliers to vanish for solutions of interest. This is because there is no semi-discrete counterpart of Proposition 9. On one hand, the semi-discrete version of Equation (78) (that is, the approach presented in Section 2) does not imply that δ 2S · δX = 0, so the above proof will not work. On the other hand, if we supplement Equation (78) with the equation corresponding to variations of X, then the finite element discretization will not have solutions, unless the constraint functions are integrals of motion of the system described byL N (y i , X j ,ẏ k ,Ẋ l ), which generally is not the case. Nonetheless, it is reasonable to expect that if the continuous system given by Equation (78) has a solution, then the Lagrange multipliers of the semi-discrete system defined by the Lagrangian in Equation (84) should remain small.
Defining constraints by Equation (77) allowed us to use the same finite element discretization for bothL and the constraints and to prove some correspondence between the solutions of Equations (78) and (80). However, the constraints of Equation (77) are global in the sense that they depend on the values of the fields ϕ and X at different points in space. Moreover, these constraints do not determine unique solutions to Equations (78) and (80), which is a little cumbersome when discussing multisymplecticity (see Section 4).

Local Constraint
In Section 2.4, we discussed how some adaptation constraints of interest can be derived from certain partial differential equations based on the equidistribution principle, for instance Equation (27). We can view these PDEs as local constraints that only depend on pointwise values of the fields ϕ, X and their spatial derivatives. Let G = G(ϕ, X, ϕ x , X x , ϕ xx , X xx , ...) represent such a local constraint. Then, similarly to Equation (78), we can write our control-theoretic strategy from Section 2 as Note that higher-order derivatives of the fields may require the use of higher degree basis functions than the ones in Equation (9) or of finite differences instead.
The Lagrange multiplier approach consists in defining the auxiliary Lagrangian: Suppose that the pair (ϕ, X) satisfies Equation (85). Then, much like in Section 3.5.1, one can easily check that the triple (ϕ, X, λ ≡ 0) satisfies the Euler-Lagrange equations associated with Equation (86). However, an analog of Proposition 9 does not seem to be very interesting in this case; therefore, we are not proving it here.
Introducing the constraints this way is convenient because the Lagrangian given in Equation (86) then represents a constrained multisymplectic field theory with a local constraint, which makes the analysis of multisymplecticity easier (see Section 4). The disadvantage is that discretization of the Lagrangian in Equation (86) requires mixed methods. We will use the linear finite elements of Equation (9) to discretizeL[ϕ, X, ϕ t , X t ], but the constraint term will be approximated via finite differences. This way, we again obtain the semi-discrete Lagrangian of Equation (84), where g i represents the discretization of G at the point x = x i .
In summary, the methods presented in Sections 3.5.1 and 3.5.2 both lead to the same semi-discrete Lagrangian but have different theoretical advantages.

DAE Formulation of the Equations of Motion
The Lagrangian (84) can be written as where The Euler-Lagrange equations thus take the forṁ where Equation (89) is to be solved for the unknown functions q(t), u(t) and λ(t). This is a DAE system of index 3, since we are lacking a differential equation for λ(t) and the constraint equation has to be differentiated three times in order to expressλ as a function of q, u, and λ, provided that certain regularity conditions are satisfied. Let us determine these conditions. Differentiate the constraint equation with respect to time twice to obtain the acceleration level constraint. where We can then write Equation (91) and the second equation of Equation (89) together as If we could solve this equation foru and λ in terms of q and u, then we could simply differentiate the expression for λ one more time to obtain the missing differential equation, thus showing Equation (89) is of index 3. Equation (93) is solvable if its matrix is invertible. Hence, for Equation (89) to be of index 3, the following condition has to be satisfied for all q or at least in a neighborhood of the points satisfying g(q) = 0. Note that, with suitably chosen constraints, this condition allows the mass matrix to be singular. We would like to perform time integration of this mechanical system using the symplectic (variational) Lobatto IIIA-IIIB quadratures for constrained systems (see References [1,2,12,[27][28][29][30][31]). However, due to the singularity of the Runge-Kutta coefficient matrices (a ij ) and (ā ij ) for the Lobatto IIIA and IIIB schemes, the assumption stated in Equation (94) does not guarantee that these quadratures define a unique numerical solution: The mass matrix would need to be invertible. To circumvent this numerical obstacle, we resort to a trick described in Reference [28]. We embed our mechanical system in a higher dimensional configuration space by adding slack degrees of freedom r andṙ and form the augmented LagrangianL A N by modifying the kinetic term ofL N to read Assuming Equation (94) holds, the augmented system has a non-singular mass matrix. If we multiply out the terms we obtain simplỹ This formula in fact holds for general Lagrangians, not only for Equation (44). In addition to g(q) = 0, we further impose the constraint r = 0. Then, the augmented constrained Lagrangian takes the formL The corresponding Euler-Lagrange equations arė It is straightforward to verify that r(t) = 0, w(t) = 0, and µ(t) = 0 are the exact solution and that the remaining equations reduce to Equation (89), that is, the evolution of the augmented system coincides with the evolution of the original system by construction. The advantage is that the augmented system is now regular and we can readily apply the Lobatto IIIA-IIIB method for constrained systems to compute a numerical solution. It should be intuitively clear that this numerical solution will approximate the solution of Equation (89) as well. What is not immediately obvious is whether a variational integrator based on the Lagrangian in Equation (96) can be interpreted as a variational integrator based onL N . This can be elegantly justified with the help of exact constrained discrete Lagrangians. Let N ⊂ Q N × G N be the constraint submanifold defined by g(q) = 0. The exact constrained discrete LagrangianL C,E N : N × N −→ R is defined bỹ where q(t) is the solution to the constrained Euler-Lagrange Equation (89) such that it satisfies the boundary conditions q(0) = q (1) and q(∆t) = q (2) . Note that N × {0} ⊂ (Q N × G N ) × R N is the constraint submanifold defined by g(q) = 0 and r = 0. Since necessarily r (1) = r (2) = 0, we can define the exact augmented constrained discrete LagrangianL A,C,E where q(t) and r(t) are the solutions to the augmented constrained Euler-Lagrange Equation (98) such that the boundary conditions q(0) = q (1) , q(∆t) = q (2) , and r(0) = r(∆t) = 0 are satisfied. Proof. Let q(t) and r(t) be the solutions to Equation (98) such that the boundary conditions q(0) = q (1) , q(∆t) = q (2) , and r(0) = r(∆t) = 0 are satisfied. As argued before, we in fact have r(t) = 0 and q(t) satisfies Equation (89) as well. By Equation (96), we havẽ This means that any discrete LagrangianL d : to order s also approximatesL C,E N to the same order, that is, a variational integrator for Equation (98), in particular our Lobatto IIIA-IIIB scheme, is also a variational integrator for Equation (89).

Backward error analysis
The advantage of the Lagrange multiplier approach is the fact that upon spatial discretization we deal with a constrained mechanical system. Backward error analysis of symplectic/variational numerical schemes for such systems shows that the modified equations also describe a constrained mechanical system for a nearby Hamiltonian (see Theorem 5.6 in Section IX.5.2 of Reference [1]). Therefore, we expect the Lagrange multiplier strategy to demonstrate better performance in terms of energy conservation than the control-theoretic strategy. The Lagrange multiplier approach makes better use of the geometry underlying the field theory we consider, the key idea being to treat the reparametrization field X(x, t) as an additional dynamical degree of freedom on equal footing with ϕ(x, t).

Multisymplectic Field Theory Formalism
In Sections 2 and 3, we took the view of infinite dimensional manifolds of fields as configuration spaces and presented a way to construct space-adaptive variational integrators in that formalism. We essentially applied symplectic integrators to semi-discretized Lagrangian field theories. In this section, we show how r-adaptive integrators can be described in the more general framework of multisymplectic geometry. In particular, we show that some of the integrators obtained in the previous sections can be interpreted as multisymplectic variational integrators. Multisymplectic geometry provides a covariant formalism for the study of field theories in which time and space are treated on equal footing, as a conseqence of which multisymplectic variational integrators allow for more general discretizations of spacetime, such that, for instance, each element of space may be integrated with a different timestep (see Reference [4]). For the convenience of the reader, below, we briefly review some background material and provide relevant references for further details. We then proceed to reformulate our adaptation strategies in the language of multisymplectic field theory.

Lagrangian Mechanics and Veselov-Type Discretizations
Let Q be the configuration manifold of a certain mechanical system and TQ be its tangent bundle. Denote the coordinates on Q by q i and on TQ by (q i ,q i ), where i = 1, 2, ..., n. The system is described by defining the Lagrangian L : TQ −→ R and the corresponding action functional S[q(t)] = b a L q i (t),q i (t) dt. The dynamics are obtained through Hamilton's principle, which seeks the curves q(t) for which the functional S[q(t)] is stationary under variations of q(t) with fixed endpoints, i.e., we seek q(t) such that for all δq(t) with δq(a) = δq(b) = 0, where q (t) is a smooth family of curves satisfying q 0 = q and d d =0 q = δq. By using integration by parts, the Euler-Lagrange equations follow as The canonical symplectic form Ω on T * Q, the 2n-dimensional cotangent bundle of Q, is given by Ω = dq i ∧ dp i , where summation over i is implied and (q i , p i ) is the canonical coordinates on T * Q. The Lagrangian defines the Legendre transformation FL : TQ −→ T * Q, which in coordinates is given by (q i , p i ) = (q i , ∂L ∂q i ). We then define the Lagrange 2-form on TQ by pulling back the canonical symplectic form, i.e., Ω L = FL * Ω. If the Legendre transformation is a local diffeomorphism, then Ω L is a symplectic form. The Lagrange vector field is a vector field X E on TQ that satisfies X E Ω L = dE, where the energy E is defined by E(v q ) = FL(v q ) · v q − L(v q ) and denotes the interior product, i.e., the contraction of a differential form with a vector field. It can be shown that the flow F t of this vector field preserves the symplectic form, that is, F * t Ω L = Ω L . The flow F t is obtained by solving the Euler-Lagrange Equation (102).
For a Veselov-type discretization, we essentially replace TQ with Q × Q, which serves as a discrete approximation of the tangent bundle. We define a discrete Lagrangian L d as a smooth map L d : Q × Q −→ R and the corresponding discrete action S = ∑ N−1 k=0 L d (q k , q k+1 ). The variational principle now seeks a sequence q 0 , q 1 , ..., q N that extremizes S for variations holding the endpoints q 0 and q N fixed. The discrete Euler-Lagrange equations follow This implicitly defines a discrete flow F : ) denotes the coordinates on Q × Q. It then follows that the discrete flow F is symplectic, i.e., F * ω L = ω L .
Given a continuous Lagrangian system with L : TQ −→ R. one chooses a corresponding discrete Lagrangian as an approximation L d (q k , q k+1 ) ≈ t k+1 t k L q(t),q(t) dt, where q(t) is the solution of the Euler-Lagrange equations corresponding to L with the boundary values q(t k ) = q k and q(t k+1 ) = q k+1 .
For more details regarding Lagrangian mechanics, variational principles, and symplectic geometry, see Reference [32]. Discrete Mechanics and variational integrators are discussed in Reference [2].

Multisymplectic Geometry and Lagrangian Field Theory
Let X be an oriented manifold representing the (n + 1)-dimensional spacetime with local coordinates (x 0 , x 1 , . . . , x n ) ≡ (t, x), where x 0 ≡ t is time and (x 1 , . . . , x n ) ≡ x are space coordinates. Physical fields are sections of a configuration fiber bundle π X Y : Y −→ X , that is, continuous maps φ : X −→ Y such that π X Y • φ = id X . This means that for every (t, x) ∈ X , φ(t, x) is in the fiber over (t, x), which is Y (t,x) = π −1 X Y ((t, x)). The evolution of the field takes place on the first jet bundle J 1 Y, which is the analog of TQ for mechanical systems. J 1 Y is defined as the affine bundle over Y such that, for y ∈ Y (t,x) , the fiber J 1 y Y consists of linear maps ϑ : Intuitively, the first jet bundle consists of the configuration bundle Y and of the first partial derivatives of the field variables with respect to the independent variables. Let φ(x 0 , . . . , x n ) = (x 0 , . . . , x n , y 1 , . . . , y m ) in coordinates and let v a µ = y a ,µ = ∂y a /∂x µ denote the partial derivatives. We can think of J 1 Y as a fiber bundle over X . Given a section φ : X −→ Y, we can define its first jet prolongation j 1 φ : X −→ J 1 Y, in coordinates given by j 1 φ(x 0 , x 1 , . . . , x n ) = (x 0 , x 1 , . . . , x n , y 1 , . . . , y m , y 1 ,0 , . . . , y m ,n ), which is a section of the fiber bundle J 1 Y over X . For higher-order field theories, we consider higher order jet bundles, defined iteratively by J 2 Y = J 1 (J 1 Y) and so on. The local coordinates on J 2 Y are denoted (x µ , y a , v a µ , w a µ , κ a µν ). The second jet prolongation j 2 φ : X −→ J 2 Y is given in coordinates by j 2 φ(x µ ) = (x µ , y a , y a ,µ , y a ,µ , y a ,µ,ν ).
Lagrangian density for first-order field theories is defined as a map L : J 1 Y −→ R. The corresponding action functional is S[φ] = U L(j 1 φ) d n+1 x, where U ⊂ X . Hamilton's principle seeks fields φ(t, x) that extremize S, that is for all η λ Y that keep the boundary conditions on ∂U fixed, where η λ Y : Y −→ Y is the flow of a vertical vector field V on Y. This leads to the Euler-Lagrange equations Given the Lagrangian density L, one can define the Cartan (n + 1)-form Θ L on J 1 Y in local coordinates given by x. The multisymplectic (n + 2)-form is then defined by Ω L = −dΘ L . Let P be the set of solutions of the Euler-Lagrange equations, that is, the set of sections φ satisfying Equation (104) or Equation (105). For a given φ ∈ P, let F be the set of first variations, that is, the set of vector fields is also a solution, where η Y is the flow of V. The multisymplectic form formula states that if φ ∈ P then for all V and W in F , where j 1 V is the jet prolongation of V, that is, the vector field on J 1 Y in local coordinates given by The multisymplectic form formula is the multisymplectic counterpart of the fact that, in finite-dimensional mechanics, the flow of a mechanical system consists of symplectic maps.
For a kth-order Lagrangian field theory with the Lagrangian density L : J k Y −→ R, analogous geometric structures are defined on J 2k−1 Y. In particular, for a second-order field theory the multisymplectic (n + 2)-form Ω L is defined on J 3 Y and a similar multisymplectic form formula can be proven. If the Lagrangian density does not depend on the second-order time derivatives of the field, it is convenient to define the subbundle J 2 0 Y ⊂ J 2 Y such that J 2 0 Y = {ϑ ∈ J 2 Y | κ a 00 = 0}. For more information about the geometry of jet bundles, see Reference [33]. The multisymplectic formalism in field theory is discussed in Reference [34]. The multisymplectic form formula for first-order field theories is derived in Reference [3] and generalized for second-order field theories in Reference [35]. Higher-order field theory is considered in Reference [36].

Multisymplectic Variational Integrators
Veselov-type discretization can be generalized to multisymplectic field theory. We take where for simplicity we consider dim X = 2, i.e., n = 1. The configuration fiber bundle is Y = X × F for some smooth manifold F . The fiber over (j, i) ∈ X is denoted Y ji and its elements are y ji . A rectangle of X is an ordered 4-tuple of the form = (j, i), (j, i + 1), (j + 1, i + 1), (j + 1, i) = ( 1 , 2 , 3 , 4 ). The set of all rectangles in X is denoted X . A point (j, i) is touched by a rectangle if it is a vertex of that rectangle. Let U ⊂ X . Then (j, i) ∈ U is an interior point of U if U contains all four rectangles that touch (j, i). The interior int U is the set of all interior points of U . The closure cl U is the union of all rectangles touching interior points of U . The boundary of U is defined by ∂U = (U ∩ cl U )\int U . A section of Y is a map φ : U ⊂ X → Y such that φ(j, i) ∈ Y ji . We can now define the discrete first jet bundle of Y as J 1 Y = (y ji , y j i+1 , y j+1 i+1 , y j+1 i ) (j, i) ∈ X , y ji , y j i+1 , y j+1 i+1 , y j+1 i ∈ F = X × F 4 . (107) Intuitively, the discrete first jet bundle is the set of all rectangles together with four values assigned to their vertices. Those four values are enough to approximate the first derivatives of a smooth section with respect to time and space using, for instance, finite differences. The first jet prolongation of a section φ of Y is the map j 1 φ : For a vector field V on Y, let V ji be its restriction to Y ji . Define a discrete Lagrangian L : J 1 Y → R, L = L(y 1 , y 2 , y 3 , y 4 ), where for convenience we omit writing the base rectangle. The associated discrete action is given by The discrete variational principle seeks sections that extremize the discrete action, that is, for all vector fields V on Y that keep the boundary conditions on ∂U fixed, where φ λ (j, i) = F V ji λ (φ(j, i)) and F V ji λ is the flow of V ji on F . This is equivalent to the discrete Euler-Lagrange equations ∂L ∂y 1 (y ji , y j i+1 , y j+1 i+1 , y j+1 i ) for all (j, i) ∈ int U , where we adopt the convention φ(j, i) = y ji . In analogy to the Veselov discretization of mechanics, we can define four 2-forms Ω l L on J 1 Y, where l = 1, 2, 3, 4 and Ω 1 L + Ω 2 L + Ω 3 L + Ω 4 L = 0, that is, only three 2-forms of these forms are independent. The 4-tuple (Ω 1 L , Ω 2 L , Ω 3 L , Ω 4 L ) is the discrete analog of the multisymplectic form Ω L . We refer the reader to the literature for details, e.g., Reference [3]. By analogy to the continuous case, let P be the set of solutions of the discrete Euler-Lagrange Equation (109). For a given φ ∈ P, let F be the set of first variations, that is, the set of vector fields V on J 1 Y defined similarly as in the continuous case. The discrete multisymplectic form formula then states that if φ ∈ P then for all V and W in F , where the jet prolongations are defined to be The discrete multisymplectic form formula given in Equation (110) is in direct analogy to the multisymplectic form formula that holds in the continuous case (Equation (106)).
Given a continuous Lagrangian density L, one chooses a corresponding discrete Lagrangian as an approximation L(y 1 , y 2 , y 3 , y 4 ) ≈ L • j 1φ dx dt, where is the rectangular region of the continuous spacetime that contains andφ(t, x) is the solution of the Euler-Lagrange equations corresponding to L with the boundary values at the vertices of corresponding to y 1 , y 2 , y 3 , and y 4 .
Similar constructions then follow, and a similar discrete multisymplectic form formula can be derived for a second order-field theory.
Multisymplectic variational integrators for first-order field theories are introduced in Reference [3] and generalized for second-order field theories in Reference [35].

Continuous Setting
We now discuss a multisymplectic setting for the approach presented in Section 2. Let the computational spacetime be X = R × R with coordinates (t, x) and consider the trivial configuration bundle Y = X × R with coordinates (t, x, y). Let U = [0, T max ] × [0, X max ] and let our scalar field be represented by a sectionφ : U −→ Y with the coordinate representationφ(t, x) = (t, x, ϕ(t, x)). Let (t, x, y, v t , v x ) denote local coordinates on J 1 Y. In these coordinates, the first jet prolongation ofφ is represented by j 1φ (t, x) = (t, x, ϕ(t, x), ϕ t (t, x), ϕ x (t, x)). Then, the Lagrangian density of Equation (6) can be viewed as a mappingL : J 1 Y −→ R. The corresponding action of Equation (3) can now be expressed asS Just like in Section 2, let us for the moment assume that the function X : U −→ [0, X max ] is known, so that we can viewL as being time-and space-dependent. The dynamics is obtained by extremizing S with respect toφ, that is, by solving forφ such that for all η λ Y that keep the boundary conditions on ∂U fixed, where η λ Y : Y −→ Y is the flow of a vertical vector field V on Y. Therefore, for an a priori known X(t, x), the multisymplectic form formula (Equation (106)) is satisfied for solutions of Equation (115).
Consider the additional bundle π X B : B = X × [0, X max ] −→ X of which the sectionsX : U −→ B represent our diffeomorphisms. LetX(t, x) = (t, x, X(t, x)) denote a local coordinate representation and assume X(t, .) is a diffeomorphism. Then, defineỸ = Y ⊕ B. We have J kỸ ∼ = J k Y ⊕ J k B. In Section 3.5.2, we argued that Equation (25) can be interpreted as a local constraint on the fieldsφ,X and their spatial derivatives. This constraint can be represented by a function G : J kỸ −→ R. Sections ϕ andX satisfy the constraint if G(j kφ , j kX ) = 0. Therefore our control-theoretic strategy expressed in Equation (85) can be rewritten as for all η λ Y , similarly as above. Let us argue how to interpret the notion of multisymplecticity for this problem. Intuitively, multisymplecticity should be understood in a sense similar to Proposition 3. We first solve Equation (116) forφ andX, given some initial and boundary conditions.
Then, we substitute thisX into Equation (115). Let P be the set of solutions to this problem. Naturally, ϕ ∈ P. The multisymplectic form formula (Equation (106)) will be satisfied for all fields in P, but the constraint G = 0 will be satisfied only forφ.

Discretization
Discretize the computational spacetime R × R by picking the discrete set of points t j = j · ∆t, x i = i · ∆x, and define X = {(j, i) | j, i ∈ Z}. Let X and X be the set of rectangles and 6-tuples in X , respectively. The discrete configuration bundle is Y = X × R, and for convenience of notation, let the elements of the fiber Y ji be denoted by y j i . Let U = {(j, i) | j = 0, 1, . . . , M + 1, i = 0, 1, . . . , N + 1}, where ∆x = X max /(N + 1) and ∆t = T max /(M + 1). Suppose we have a discrete LagrangianL : J 1 Y −→ R and the corresponding discrete actionS that approximates the action in Equation (114), where we assume that X(t, x) is known and of the form (10). A variational integrator is obtained by solving for a discrete sectionφ : U −→ Y, as described in Section 4.1. This integrator is multisymplectic, i.e., the discrete multisymplectic form formula (Equation (110)) is satisfied.
In Equation (20), consider the 1-stage symplectic partitioned Runge-Kutta method with the coefficients a 11 =ā 11 = c 1 = 1/2 and b 1 =b 1 = 1. This method is often called the midpoint rule and is a 2nd order member of the Gauss family of quadratures. It can be easily shown (see References [1,2]) that the discrete Lagrangian of Equation (15) for this method is given bỹ where ∆t = t j+1 − t j and y j = (y j 1 , . . . , y j N ). Using Equations (5) and (13), we can writẽ where we defined the discrete LagrangianL : J 1 Y −→ R by the formulã Given the Lagrangian densityL as in Equation (6) and assuming X(t, x) is known, one can evaluate the integral in Equation (120) explicitly. It is now a straightforward calculation to show that the discrete variational principle of Equation (117) for the discrete LagrangianL as defined is equivalent to the discrete Euler-Lagrange Equation (103) forL d and consequently to Equation (20).
This shows that the 2nd order Gauss method applied to Equation (20) (15) as a sum similar to Equation (119). The resulting integrators are still variational, since they are derived by applying the discrete variational principle of Equation (117) to some discrete actionS, but this action cannot be expressed as the sum ofL over all rectangles. Therefore, these integrators are not multisymplectic, at least not in the sense of Equation (110).

Constraints.
Let the additional bundle be B = X × [0, X max ], and denote by X n j the elements of the fiber B ji . DefineỸ = Y ⊕ B. We have J kỸ ∼ = J k Y ⊕ J k B. Suppose G : J kỸ −→ R represents a discretization of the continuous constraint. For instance, one can enforce a uniform mesh by defining G : J 1Ỹ → R, G(j 1φ , j 1X ) = X x − 1 at the continuous level. The discrete counterpart will be defined on the discrete jet bundle J 1Ỹ by the following formula: Arclength equidistribution can be realized by enforcing Equation (27), that is, G : The discrete counterpart will be defined on the discrete sub-bundle J 2 0Ỹ by the following formula: where, for convenience, we used the notation introduced in Equation (113) and l, r = 1, . . . , 6. Note that Equation (123) coincides with Equation (28). In fact, g i in Equation (28) is nothing else but G computed on an element of J 2 0Ỹ over the base 6-tuple such that 2 = (j, i). The only difference is that, in Equation (28), we assumed g i might depend on all the field values at a given time step, while G only takes arguments locally, i.e., it depends on at most 6 field values on a given 6-tuple.
A numerical scheme is now obtained by simultaneously solving the discrete Euler-Lagrange Equation (109) resulting from Equation (117)

Continuous Setting
We now turn to describing the Lagrange multiplier approach in a multisymplectic setting. Similarly as in Section 4.2, let the computational spacetime be X = R × [0, X max ] with coordinates (t, x) and consider the trivial configuration bundles π X Y : Y = X × R −→ X and π X B : B = X × [0, X max ] −→ X . Let our scalar field be represented by a sectionφ : X −→ Y with the coordinate representationφ(t, x) = (t, x, ϕ(t, x)) and our diffeomorphism be represented by a sectionX : X −→ B with the local representationX(t, x) = (t, x, X(t, x)). Let the total configuration bundle beỸ = Y ⊕ B. Then, the Lagrangian density of Equation (6) can be viewed as a mapping

corresponding action of Equation (3) can now be expressed as
where U = [0, T max ] × [0, X max ]. As before, the MMPDE constraint can be represented by a function G : J kỸ −→ R. Two sectionsφ andX satisfy the constraint if Vakonomic formulation.
We now face the problem of finding the right equations of motion. We want to extremize the action functional of Equation (124) in some sense, subject to the constraint in Equation (125). Note that the constraint is essentially nonholonomic, as it depends on the derivatives of the fields. Assuming G is a submersion, G = 0 defines a submanifold of J kỸ , but this submanifold will not in general be the kth jet of any subbundle ofỸ. Two distinct approaches are possible here. One could follow the Lagrange-d'Alembert principle and take variations ofS first but choose variations V (vertical vector fields onỸ) such that the jet prolongations j k V are tangent to the submanifold G = 0 and then enforce the constraint G = 0. On the other hand, one could consider the variational nonholonomic problem (also called vakonomic), and minimizeS over the set of all sections (φ,X) that satisfy the constraint G = 0, that is, enforce the constraint before taking the variations. If the constraint is holonomic, both approaches yield the same equations of motion. However, if the constraint is nonholonomic, the resulting equations are in general different. Which equations are correct is really a matter of experimental verification. It has been established that the Lagrange-d'Alembert principle gives the right equations of motion for nonholonomic mechanical systems, whereas the vakonomic setting is appropriate for optimal control problems (see References [37][38][39][40]).
We will argue that the vakonomic approach is the right one in our case. In Proposition 5, we showed that in the unconstrained case extremizing S[φ] with respect to φ was equivalent to extremizingS[φ,X] with respect toφ, and in Proposition 6, we showed that extremizing with respect toX did not yield new information. This is because there was no restriction on the fieldsφ and X and, for any givenX, there was a one-to-one correspondence between φ andφ given by the formula ϕ(t, x) = φ(t, X(t, x)), so extremizing over all possibleφ was equivalent to extremizing over all possible φ. Now, let N be the set of all smooth sections (φ,X) that satisfy Equation (125) such that X(t, .) is a diffeomorphism for all t. It should be intuitively clear that, under appropriate assumptions on the mesh density function ρ, for any given smooth function φ(t, X), Equation (25) together with ϕ(t, x) = φ(t, X(t, x)) define a unique pair (φ,X) ∈ N (since our main purpose here is to only justify the application of the vakonomic approach, we do not attempt to specify those analytic assumptions precisely). Conversely, any given pair (φ,X) ∈ N defines a unique function φ through the formula φ(t, X) = ϕ(t, ξ(t, X)), where ξ(t, .) = X(t, .) −1 , as in Section 3.1. Given this one-to-one correspondence and the fact that S[φ] =S[φ,X] by definition, we see that extremizing S with respect to all smooth φ is equivalent to extremizingS over all smooth sections (φ,X) ∈ N . We conclude that the vakonomic approach is appropriate in our case, since it follows from Hamilton's principle for the original, physically meaningful, action functional S.
Let us also note that our constraint depends on spatial derivatives only. Therefore, in the setting presented in Sections 2 and 3 it can be considered holonomic, as it restricts the infinite-dimensional configuration manifold of fields that we used as our configuration space. In that case, it is valid to use Hamilton's principle and to minimize the action functional over the set of all allowable fields, i.e., those that satisfy the constraint G = 0. We did that by considering the augmented instantaneous Lagrangian (86).
In order to minimizeS over the set of sections satisfying Equation (125), we will use the bundle-theoretic version of the Lagrange multiplier theorem, which we cite below after Reference [41].
Theorem 1 (Lagrange multiplier theorem). Let π M,E : E −→ M be an inner product bundle over a smooth manifold M, Ψ be a smooth section of π M,E , and h : M −→ R be a smooth function. Setting N = Ψ −1 (0), the following are equivalent: 1.
σ ∈ N is an extremum of h| N , 2.
There exists an extremumσ ∈ E ofh : E −→ R such that π M,E (σ) = σ, Let us briefly review the ideas presented in Reference [41], adjusting the notation to our problem and generalizing when necessary. Let be the set of smooth sections of π XỸ on U . Then,S : and let C ∞ U (V ) be the set of smooth sectionsλ : U −→ V, which represent our Lagrange multipliers and in local coordinates have the representationλ(t, x) = (t, x, λ(t, x)). The set C ∞ U (V ) is an inner product space with λ 1 ,λ 2 = U λ 1 λ 2 dt ∧ dx. Take This is an inner product bundle over C ∞ U (Ỹ) with the inner product defined by (σ,λ 1 ), (σ,λ 2 ) We now have to construct a smooth section Ψ : C ∞ U (Ỹ) −→ E that will realize our constraint of Equation (125). Define the fiber-preserving mappingG : J kỸ −→ V such that for ϑ ∈ J kỸ G(ϑ) = π X ,J kỸ (ϑ), G(ϑ) .
The augmented action functionalS C : E −→ R is now given bỹ or denotingσ = (φ,X,λ) Theorem 1 states that, if (φ,X,λ) is an extremum ofS C , then (φ,X) extremizesS over the set N of sections satisfying the constraint G = 0. Note that using the multisymplectic formalism we obtained the same result as Equation (86) in the instantaneous formulation, where we could treat G as a holonomic constraint. The dynamics is obtained by solving for a triple (φ,X,λ) such that for all η Y , η B , η V that keep the boundary conditions on ∂U fixed, where η denotes the flow of vertical vector fields on respective bundles. Note that we can defineỸ C = Y ⊕ B ⊕ V andL C : J kỸ C −→ R by settingL C =L − λ · G, i.e., we can consider a kth-order field theory. If k = 1, 2 then an appropriate multisymplectic form formula in terms of the fieldsφ,X, andλ will hold. Presumably, this can be generalized for k > 2 using the techniques put forth in Reference [35]. However, it is an interesting question whether there exists any multisymplectic form formula defined in terms ofφ,X, and objects on J kỸ only. It appears to be an open problem. This would be the multisymplectic analog of the fact that the flow of a constrained mechanical system is symplectic on the constraint submanifold of the configuration space.

Discretization
Let us use the same discretization as discussed in Section 4.2. Assume we have a discrete LagrangianL : J 1Ỹ −→ R, the corresponding discrete actionS[φ,X], and a discrete constraint G : J 1Ỹ −→ R or G : J 2 0Ỹ −→ R. Note thatS is essentially a function of 2MN variables and that we want to extremize it subject to the set of algebraic constraints G = 0. The standard Lagrange multiplier theorem proved in basic calculus textbooks applies here. However, let us work out a discrete counterpart of the formalism introduced at the continuous level. This will facilitate the discussion of the discrete notion of multisymplecticity. Let be the set of discrete sections of π XỸ :Ỹ −→ X . Similarly, define the discrete bundle V = X × R, and let C U 0 (V ) be the set of discrete sectionsλ : U 0 −→ V representing the Lagrange multipliers, where U 0 ⊂ U is defined below. Letλ(j, i) = (j, i, λ(j, i)) with λ j i ≡ λ(j, i) be the local representation. The set C U 0 (V ) is an inner product space with λ ,μ = ∑ (j,i)∈U 0 λ j i µ j i . Take E = C U (Ỹ) × C U 0 (V ). Just like at the continuous level, E is an inner product bundle. However, at the discrete level, it is more convenient to define the inner product on E in a slightly modified way. Since there are some nuances in the notation, let us consider the cases k = 1 and k = 2 separately.
Define the trivial bundleV = X × R and let C U (V ) be the set of all sections ofV defined on U . For a given sectionλ ∈ C U 0 (V ), we define its extension that is,λ assigns to the square the value thatλ takes on the first vertex of that square. Note that this operation is invertible: Given a section of C U (V ), we can uniquely determine a section of C U 0 (V ). We can define the inner product One can easily see that we have λ ,μ = λ ,μ , so by a slight abuse of notation, we can use the same symbol ., . for both inner products. It will be clear from the context which definition should be invoked. We can now define an inner product on the fibers of E as (σ,λ), (σ,μ) Let us now construct a section Ψ : C U (Ỹ) −→ E that will realize our discrete constraint G. First, in analogy to Equation (130), define the fiber-preserving mappingG : J 1Ỹ −→V such that where l, r = 1, 2, 3, 4. We now define Ψ by requiring that for σ ∈ C U (Ỹ) the extension of Ψ(σ), defined similar to Equation (136), is given byΨ The set of allowable sections N ⊂ C U (Ỹ) is now defined by N = Ψ −1 (0)-that is, (φ,X) ∈ N , provided that G(j 1φ , j 1X ) = 0 for all ∈ U . The augmented discrete actionS C : E −→ R is thereforẽ By the standard Lagrange multiplier theorem, if (φ,X,λ) is an extremum ofS C , then (φ,X) is an extremum ofS over the set N of sections satisfying the constraint G = 0. The discrete Hamilton principle can be expressed as for all vector fields V on Y, W on B, and Z on V that keep the boundary conditions on ∂U fixed, The discrete Euler-Lagrange equations are obtained by differentiating with respect to y j i , X we can consider an unconstrained field theory in terms of the fieldsφ,X, andλ. Then, the solutions of Equation (144) satisfy the multisymplectic form formula (Equation (110)) in terms of objects defined on J 1Ỹ C .
Define the trivial bundleV = X × R, and let C U (V ) be the set of all sections ofV defined on U . For a given sectionλ ∈ C U 0 (V ), we define its extension that is,λ assigns to the 6-tuple the value thatλ takes on the second vertex of that 6-tuple. Like before, this operation is invertible. We can define the inner product and the inner product on E as in Equation (138). Define the fiber-preserving mappingG : where l, r = 1, . . . , 6. We now define Ψ by requiring that, for σ ∈ C U (Ỹ), the extension of Ψ(σ), defined similar to Equation (146), is given byΨ Again, the set of allowable sections is N = Ψ −1 (0). That is, (φ,X) ∈ N , provided that G(j 2 0φ , j 2 0X ) = 0 for all ∈ U . The augmented discrete actionS C : E −→ R is therefore Writing out the terms involving y j i , X j i , and λ j i explicitly, as in Equation (143), and invoking the discrete Hamilton principle (142), one obtains the discrete Euler-Lagrange equations, which can be compactly expressed as Let us also set G(y 1 , . . . , X 4 ) = 0 if 2 = (j, 0), (j, N + 1). Define A = { | 2 , 5 ∈ U }. Then, Equation (150) can be written as where the last equality defines the augmented LagrangianL C : J 2 0Ỹ C −→ R forỸ C = Y ⊕ B ⊕ V. Therefore, we can consider an unconstrained second-order field theory in terms of the fieldsφ,X, andλ, and the solutions of Equation (151) will satisfy a discrete multisymplectic form formula very similar to the one proved in Reference [35]. The only difference is the fact that the authors analyzed a discretization of the Camassa-Holm equation and were able to consider an even smaller sub-bundle of the second jet of the configuration bundle. As a result, it was sufficient for them to consider a discretization based on squares rather than 6-tuples . In our case, there will be six discrete 2-forms Ω lL C for l = 1, . . . , 6 instead of just four.

Remark 5.
In both cases, we showed that our discretization leads to integrators that are multisymplectic on the augmented jets J kỸ C . However, just like in the continuous setting, it is an interesting problem whether there exists a discrete multisymplectic form formula in terms of objects defined on J kỸ only.
Consider the semi-discrete Lagrangian in Equation (44). We can use the trapezoidal rule to define the discrete Lagrangian in Equation (14) as where y j = (y where for brevity q j = (y . . , λ j N ) and g is an adaptation constraint, for instance Equation (28). If q j−1 , q j are known, then Equation (155) can be used to compute q j+1 and λ j . It is easy to verify that the condition of Equation (94) is enough to ensure solvability of Equation (155), assuming the time step ∆t is sufficiently small, so there is no need to introduce slack degrees of freedom as in Equation (95). If the mass matrix of Equation (47) was constant and non-singular, then Equation (155) would result in the SHAKE algorithm, or in the RATTLE algorithm if one passes to the position-momentum formulation (see References [1,2]).
Using Equations (38) and (41), we can writẽ where we defined the discrete LagrangianL : J 1Ỹ −→ R by the formulã and similarly forX(x). Given the Lagrangian densityL as in Equation (43), one can compute the integrals in Equation (157) explicitly. Suppose that the adaptation constraint g has a "local" structure, for instance as in Equation (122) or as in Equation (123). It is straightforward to show that Equation (144) or Equation (151) are equivalent to Equation (155), that is, the variational integrator defined by Equation (155) is also multisymplectic.
For reasons similar to the ones pointed out in Section 4.2, the 2nd and 4th order Lobatto IIIA-IIIB methods that we used for our numerical computations are not multisymplectic.

The Sine-Gordon Equation
We applied the methods discussed in the previous sections to the Sine-Gordon equation This equation results from the (1+1)-dimensional scalar field theory with the Lagrangian density The Sine-Gordon equation arises in many physical applications. For instance, it governs the propagation of dislocations in crystals, the evolution of magnetic flux in a long Josephson-junction transmission line, or the modulation of a weakly unstable baroclinic wave packet in a two-layer fluid. It also has applications in the description of one-dimensional organic conductors, one-dimensional ferromagnets, and liquid crystals or in particle physics as a model for baryons (see References [42,43]).
The Sine-Gordon equation has interesting soliton solutions. A single soliton traveling at the speed v is given by It is depicted in Figure 2. The backscattering of two solitons, each traveling with the velocity v, is described by the formula It is depicted in Figure 3. Note that if we restrict X ≥ 0, then this formula also gives a single soliton solution satisfying the boundary condition φ(0, t) = 0, that is, a soliton bouncing from a rigid wall.

Generating Consistent Initial Conditions
Suppose we specify the following initial conditions and assume they are consistent with the boundary conditions stated in Equation (2). In order to determine appropriate consistent initial conditions for Equations (18) and (98), we need to solve several equations. First, we solve for the y i 's and X i 's. We have y 0 = φ L , y N+1 = φ R , X 0 = 0, and X N+1 = X max . The rest are determined by solving the system y i = a(X i ), 0 = g i (y 1 , . . . , y N , X 1 , . . . , X N ), for i = 1, . . . , N. This is a system of 2N nonlinear equations for 2N unknowns. We solve it using Newton's method. Note, however, that we do not a priori know good starting points for Newton's iterations. If our initial guesses are not close enough to the desired solution, the iterations may converge to the wrong solution or may not converge at all. In our computations, we used the constraints given in Equation (28). We found that a very simple variant of a homotopy continuation method worked very well in our case. Note that for α = 0 Equation (28) generates a uniform mesh. In order to solve Equation (166) for some α > 0, we split [0, α] into d subintervals by picking α k = (k/d) · α for k = 1, . . . , d. We then solved Equation (166) with α 1 using the uniformly spaced mesh points X (0) i = (i/(N + 1)) · X max as our initial guess, resulting in X (1) i and y (1) i . Then, we solved Equation (166) with α 2 using X (1) i and y (1) i as the initial guesses, resulting in X (2) i and y (2) i . Continuing in this fashion, we got X i as the numerical solution to Equation (166) for the original value of α. Note that for more complicated initial conditions and constraint functions, predictor-corrector methods should be used-see Reference [44] for more information. Another approach to solving Equation (166) could be based on relaxation methods (see References [7,8]).
Next, we solve for the initial values of the velocitiesẏ i andẊ i . Since ϕ( We also require that the velocities be consistent with the constraints. Hence, the linear system iṡ This is a system of 2N linear equations for the 2N unknownsẏ i andẊ i , where y = (y 1 , . . . , y N ) and X = (X 1 , . . . , X N ). We can use those velocities to compute the initial values of the conjugate momenta. For the control-theoretic approach, we use p i = ∂L N /∂ẏ i , as in Section 2.3, and for the Lagrange multiplier approach, we use Equation (46). In addition, for the Lagrange multiplier approach, we also have the initial values for the slack variables r i = 0 and their conjugate momenta B i = ∂L A N /∂ṙ i = 0. It is also useful to use Equation (93) to compute the initial values of the Lagrange multipliers λ i that can be used as initial guesses in the first iteration of the Lobatto IIIA-IIIB algorithm. The initial guesses for the slack Lagrange multipliers are trivially µ i = 0. Note that both λ and µ are algebraic variables, so their values at each time step are completely determined by the Lobatto IIIA-IIIB algorithm (see References [1,27,28] for details), and therefore no further initial or boundary conditions are necessary.

Convergence
In order to test the convergence of our methods as the number of mesh points N is increased, we considered a single soliton bouncing from two rigid walls at X = 0 and X = X max = 25. We imposed the boundary conditions φ L = 0 and φ R = 2π, and as initial conditions we used Equation (163) with X 0 = 12.5 and v = 0.9. It is possible to obtain the exact solution to this problem by considering a multi-soliton solution to Equation (161) on the whole real line. Such a solution can be obtained using a Bäcklund transformation (see References [42,43]). However, the formulas quickly become complicated and, technically, one would have to consider an infinite number of solitons. Instead, we constructed a nearly exact solution by approximating the boundary interactions with (164): where n is an integer number and T satisfies φ SS (X max /2, T) = π (we numerically found T ≈ 13.84). Given how fast the functions given in Equations (163) and (164) approach their asymptotic values, one may check that Equation (168) can be considered exact to machine precision. We performed numerical integration with the constant time step ∆t = 0.01 up to the time T max = 50. For the control-theoretic strategy, we used the 1-stage and 2-stage Gauss method (2nd and 4th order, respectively), and the 2-stage and 3-stage Lobatto IIIA-IIIB method (also 2nd/4th order).
For the Lagrange multiplier strategy, we used the 2-stage and 3-stage Lobatto IIIA-IIIB method for constrained mechanical systems (2nd/4th order). See References [1,12,14] for more information about the mentioned symplectic Runge-Kutta methods. We used the constraints of Equation (28) based on the generalized arclength density of Equation (26). We chose the scaling parameter to be α = 2.5, so that approximately half of the available mesh points were concentrated in the area of high gradient. A few example solutions are presented in Figures 4-7. Note that the Lagrange multiplier strategy was able to accurately capture the motion of the soliton with merely 17 mesh points (that is, N = 15). The trajectories of the mesh points for several simulations are depicted in Figures 8 and 9. An example solution computed on a uniform mesh is depicted in Figure 10.  For the convergence test, we performed simulations for several N in the range 15-127. For comparison, we also computed solutions on a uniform mesh for N in the range 15-361. The numerical solutions were compared against the solution in Equation (168). The L ∞ errors are depicted in Figure 11. The L ∞ norms were evaluated over all nodes and over all time steps. Note that, in case of a uniform mesh, the spacing between the nodes is ∆x = X max /(N + 1); therefore, the errors are plotted versus (N + 1). The Lagrange multiplier strategy proved to be more accurate than the control-theoretic strategy. As the number of mesh points is increased, the uniform mesh solution becomes quadratically convergent, as expected, since we used linear finite elements for spatial discretization. The control-theoretic strategy also shows near quadratic convergence, whereas the Lagrange multiplier method seems to converge slightly slower. While there are very few analytical results regarding the convergence of r-adaptive methods, it has been observed that the rate of convergence depends on several factors, including the chosen mesh density function. Our results are consistent with the convergence rates reported in References [45,46]. Both papers deal with the viscous Burgers' equation but consider different initial conditions. Computations with the arc-length density function converged only linearly in Reference [45] but quadratically in Reference [46].

Energy Conservation
As we pointed out in Section 2.6, the true power of variational and symplectic integrators for mechanical systems lies in their excellent conservation of energy and other integrals of motion, even when a big time step is used. In order to test the energy behavior of our methods, we performed simulations of the Sine-Gordon equation over longer time intervals. We considered two solitons bouncing from each other and from two rigid walls at X = 0 and X max = 25. We imposed the boundary conditions φ L = −2π and φ R = 2π, and as initial conditions we used φ(X, 0) = φ SS (X − 12.5, −5) with v = 0.9. We ran our computations on a mesh consisting of 27 nodes (N=25). Integration was performed with the time step ∆t = 0.05, which is rather large for this type of simulations. The scaling parameter in Equation (28) was set to α = 1.5, so that approximately half of the available mesh points were concentrated in the areas of high gradient. An example solution is presented in Figure 12.   Figure 11. Comparison of the convergence rates of the discussed methods: Integration in time was performed using the 4th-order Lobatto IIIA-IIIB method for constrained systems in case of the Lagrange multiplier strategy and the 4th-order Gauss scheme in case of both the control-theoretic strategy and the uniform mesh simulation. The 4th-order Lobatto IIIA-IIIB scheme for the control-theoretic strategy and the uniform mesh simulation yield a very similar level of accuracy. Also, using 2nd-order integrators gives very similar error plots. The exact energy of the two-soliton solution can be computed using Equation (7). It is possible to compute that integral explicitly to obtain E = 16/ √ 1 − v 2 ≈ 36.71. The energy associated with the semi-discrete Lagrangian in Equation (44) can be expressed by the formula where R N was defined in Equation (88) and for our Sine-Gordon system is given by and M N is the mass matrix in Equation (47). The energy E N is an approximation to Equation (7) if the integrand is sampled at the nodes X 0 ,. . .,X N+1 and then piecewise linearly approximated. Therefore, we used E N to compute the energy of our numerical solutions.
The energy plots for the Lagrange multiplier strategy are depicted in Figure 13. We can see that the energy stays nearly constant in the presented time interval, showing only mild oscillations, which are reduced as a higher order of integration in time is used. The energy plots for the control-theoretic strategy are depicted in Figure 14. In this case, the discrete energy is more erratic and not as nearly preserved. Moreover, the symplectic Gauss and Lobatto methods show virtually the same energy behavior as the non-symplectic Radau IIA method, which is known for its excellent stability properties when applied to stiff differential equations (see Reference [12]). It seems that we do not gain much by performing symplectic integration in this case. It is consistent with our observations in Section 2.6 and shows that the control-theoretic strategy does not take the full advantage of the underlying geometry. t Energy Figure 13. The discrete energy E N for the Lagrange multiplier strategy: Integration in time was performed with the 2nd (top) and 4th (bottom) order Lobatto IIIA-IIIB method for constrained mechanical systems. The spikes correspond to the times when the solitons bounce off of each other or of the walls.
As we did not use adaptive time-stepping and did not implement any mesh smoothing techniques, the quality of the mesh deteriorated with time in all the simulations, eventually leading to mesh crossing, i.e., two mesh points collapsing or crossing each other. The control-theoretic strategy, even though less accurate, retained good mesh quality longer, with the breakdown time T break > 1000, as opposed to T break ∼ 600 in case of the Lagrange multiplier approach (both using a rather large constant time step). We discuss extensions to our approach for increased robustness in Section 6.

Summary and Future Work
We have proposed two general ideas on how r-adaptive meshes can be applied in geometric numerical integration of Lagrangian partial differential equations. We have constructed several variational and multisymplectic integrators and discussed their properties. We have used the Sine-Gordon model and its solitonic solutions to test our integrators numerically.
Our work can be extended in many directions. Interestingly, it also opens many questions in geometric mechanics and multisymplectic field theory. Addressing those questions may have a broad impact on the field of geometric numerical integration.

Non-Hyperbolic Equations
The special form of the Lagrangian density considered in Equation (42) leads to a hyperbolic PDE, which poses a challenge to r-adaptive methods, as at each time step the mesh is adapted globally in response to local changes in the solution. Causality and the structure of the characteristic lines of hyperbolic systems make r-adaptation prone to instabilities, and integration in time has to be performed carefully. The literature on r-adaptation almost entirely focuses on parabolic problems (see References [7,8] and references therein). Therefore, it would be interesting to apply our methods to PDEs that are first-order in time, for instance the Korteweg-de Vries, Nonlinear Schrödinger, or Camassa-Holm equations. All three equations are first-order in time and are not hyperbolic in nature. Moreover, all can be derived as Lagrangian field theories (see References [35,42,[47][48][49][50][51]). The Nonlinear Schrödinger equation has applications to optics and water waves, whereas the Korteweg-de Vries and Camassa-Holm equations were introduced as models for waves in shallow water. All equations possess interesting solitonic solutions. The purpose of r-adaptation would be to improve resolution, for instance, to track the motion of solitons by placing more mesh points near their centers and making the mesh less dense in the asymptotically flat areas.

Hamiltonian Field Theories
Variational multisymplectic integrators for field theories have been developed in the Lagrangian setting [3,35]. However, many interesting field theories are formulated in the Hamiltonian setting. They may not even possess a Lagrangian formulation. It would be interesting to construct Hamiltonian variational integrators for multisymplectic PDEs by generalizing the variational characterization of discrete Hamiltonian mechanics. This would allow to handle Hamiltonian PDEs without the need for converting them to the Lagrangian framework. Recently Leok and Zhang [52] and Vankerschaver, Ciao, and Leok [53] have laid foundations for such integrators. It would also be interesting to see if the techniques we used in our work could be applied in order to construct r-adaptive Hamiltonian integrators.

Time Adaptation Based on Local Error Estimates
One of the challenges of r-adaptation is that it requires solving differential-algebraic or stiff ordinary differential equations. This is because there are two different time scales present: one defined by the physics of the problem and one following from the strategy we use to adapt the mesh. Stiff ODEs and DAEs are known to require time integration with an adaptive step size control based on local error estimates (see References [11,12]). In our work, we used constant time-stepping, as adaptive step size control is difficult to combine with geometric numerical integration. Classical step size control is based on past information only; time symmetry is destroyed and with it the qualitative properties of the method. Hairer and Söderlind [54] developed explicit, reversible, symmetry-preserving, adaptive step size selection algorithms for geometric integrators, but their method is not based on local error estimation, thus it is not useful for r-adaptation. Symmetric error estimators are considered in Reference [28] and some promising results are discussed. Hopefully, the ideas presented in those papers could be combined and generalized. The idea of Asynchronous Variational Integrators (see Reference [4]) could also be useful here, as this would allow the use of a different time step for each cell of the mesh.

Constrained Multisymplectic Field Theories
The multisymplectic form formula stated in Equation (106) was first introduced in Reference [3]. The authors, however, consider only unconstrained field theories. In our work, we start with the unconstrained field theory of Equation (1), but upon choosing an adaptation strategy represented by the constraint G = 0, we obtain a constrained theory, as described in Sections 3 and 4.3. Moreover, this constraint is essentially non-holonomic, as it contains derivatives of the fields and the equations of motion are obtained using the vakonomic approach (also called variational nonholonomic) rather than the Lagrange-d'Alembert principle. All that gives rise to many very interesting and general questions. Is there a multisymplectic form formula for such theories? Is it derived in a similar fashion? Do variational integrators obtained this way satisfy some discrete multisymplectic form formula? These issues have been touched upon in Reference [41] but are by no means resolved.

Mesh Smoothing and Variational Nonholonomic Integrators
The major challenge of r-adaptive methods is mesh crossing, which occurs when two mesh points collapse or cross each other. In order to avoid mesh crossing and to retain good mesh quality, mesh smoothing techniques were developed [7,8]. They essentially attempt to regularize the exact equidistribution constraint G = 0 by replacing it with the condition ∂X/∂t = G, where is a small parameter. This can be interpreted as adding some attraction and repulsion pseudoforces between mesh points. If one applies the Lagrange multiplier approach to r-adaptation as described in Section 3, then upon finite element discretization, one obtains a finite dimensional Lagrangian system with a non-holonomic constraint. This constraint is enforced using the vakonomic (non-holonomic variational) formulation. Variational integrators for systems with non-holonomic constraints have been developed mostly in the Lagrange-d'Alembert setting, but there have also been some results regarding discrete vakonomic mechanics. The ideas presented in References [55][56][57] may be used to design structure-preserving mesh smoothing techniques.