Variational Integrators for Nonvariational Partial Differential Equations
VVariational Integrators for NonvariationalPartial Differential Equations
Michael Kraus ([email protected])
Max-Planck-Institut für PlasmaphysikBoltzmannstraße 2, 85748 Garching, DeutschlandTechnische Universität München, Zentrum MathematikBoltzmannstraße 3, 85748 Garching, Deutschland
Omar Maj ([email protected])
Max-Planck-Institut für PlasmaphysikBoltzmannstraße 2, 85748 Garching, Deutschland
November 23, 2015
Abstract
Variational integrators for Lagrangian dynamical systems provide a systematicway to derive geometric numerical methods. These methods preserve a discrete mul-tisymplectic form as well as momenta associated to symmetries of the Lagrangianvia Noether’s theorem. An inevitable prerequisite for the derivation of variationalintegrators is the existence of a variational formulation for the considered problem.Even though for a large class of systems this requirement is fulfilled, there are manyinteresting examples which do not belong to this class, e.g., equations of advection-diffusion type frequently encountered in fluid dynamics or plasma physics.On the other hand, it is always possible to embed an arbitrary dynamical systeminto a larger Lagrangian system using the method of formal (or adjoint) Lagrangians.We investigate the application of the variational integrator method to formal La-grangians, and thereby extend the application domain of variational integrators toinclude potentially all dynamical systems.The theory is supported by physically relevant examples, such as the advec-tion equation and the vorticity equation, and numerically verified. Remarkably,the integrator for the vorticity equation combines Arakawa’s discretisation of thePoisson brackets with a symplectic time stepping scheme in a fully covariant waysuch that the discrete energy is exactly preserved. In the presentation of the results,we try to make the geometric framework of variational integrators accessible to nonspecialists.
Keywords:
Conservation Laws, Geometric Discretization, Lagrangian Field Theory, Linear andNonlinear PDEs, Noether Theorem, Variational Integrators, Variational Methods for PDEs a r X i v : . [ m a t h . NA ] N ov ontents In recent years, the field of structure-preserving or geometric discretisation [17, 22, 14] hasbecome a flourishing discipline of numerical analysis and scientific computing. One particularfamily of geometric discretisation methods is that of variational integrators [64, 46, 48, 34,47, 42, 40], which are based on the discretisation of Hamilton’s principle of stationary action[31, 6, 23, 45, 3]. Variational integrators preserve a discrete multisymplectic form and havegood longtime energy behaviour. As we will see, they can be designed to preserve energy evenexactly, which in practice means up to machine precision. Furthermore, they preserve momentaassociated to symmetries of the discrete equations of motion via a discrete version of Noether’stheorem [53, 33]. hile in most standard discretisation techniques for dynamical systems the equations ofmotion are directly discretised, the basic idea of variational integrators is to construct a discretecounterpart to the considered system. This means that the fundamental building blocks ofclassical mechanics and field theory, namely the action functional, the Lagrangian, the variationalprinciple, and the Noether theorem, all have discrete equivalents. The application of the discretevariational principle to the discrete action then leads to discrete Euler-Lagrange equations. Theevolution map that corresponds to the discrete Euler-Lagrange equations is what is called avariational integrator. The discrete Noether theorem can be used to relate symmetries of thediscretised system to discrete momenta that are exactly preserved by this integrator. Whereasmost standard techniques put emphasis on the minimisation of local errors, for variationalintegrators the focus is rather on the preservation of global or geometric properties of the system.An obvious limitation of the variational integrator method is its applicability to Lagrangiansystems only. This excludes a large class of interesting systems, for example the problems ofadvection-diffusion type often found in fluid dynamics and plasma physics. We propose herethat the method of formal (or adjoint) Lagrangians [8] can be used as an expedient to avoidthis limitation. More specifically, formal Lagrangians allow us to embed any given system into alarger system which, in turn, admits a Lagrangian formulation. To obtain a formal Lagrangian L , the equation at hand, say F [ u ] = 0, is multiplied by an adjoint variable v , giving L = v · F [ u ].Variation of the resulting action functional, A = R L d n +1 x, with respect to v gives the originalequation F [ u ] = 0. Variation of the action functional with respect to the physical variable u gives an additional equation that determines the evolution of the adjoint variable v .At first sight one might be tempted to regard the formal Lagrangian formalism as merelya method for obtaining a weak formulation of the problem at hand. Then, if our goal is toobtain an integrator, the details of the dynamics of the adjoint variable v would seem irrelevant.However, it turns out that the dynamics of v play a key role in relating symmetries of the formalLagrangian to conservation laws satisfied by u . Ibragimov [25, 26, 27] developed a theory forthe analysis of conservation laws of arbitrary differential equations by extending the Noethertheorem to formal Lagrangians. This leads to conservation laws for the extended system ( u, v ),which can be restricted to the original system provided that it is possible to express the solutionof the adjoint variable v in terms of u .In this work, we propose the combination of the discrete variational principle with Ibragi-mov’s theory in order to derive variational integrators for systems without a natural Lagrangianformulation and to determine the associated discrete conservation laws. Thereby we extendsignificantly the applicability of the variational integrator method. The goal of this approachis to design numerical schemes which respect certain conservation laws of a given system in arather systematic way.We proceed as follows. In section 2, we present the theory of variational integrators in simpleterminology. To set the stage and fix notation we review the continuous action principle for fieldtheories and the corresponding Noether theorem before passing over to the discrete theory,which is extended to account for discrete divergence symmetries. The style of presentation ischosen to make the theory accessible to a wide audience without extensive background in moderndifferential geometry. This implies some loss of generality, but hopefully not too much of thegeometric beauty of the original work is lost. In section 3, we recall the inverse problem ofthe calculus of variations, review the theory of formal Lagrangians and explain the derivation ofconserved quantities in this setting. We also provide a geometric formulation of the theory, whichto our knowledge has not been presented, yet. Finally, in section 4, we apply the method to someprototypical examples, including the advection equation and the vorticity equation, and verifythe theoretical properties in numerical experiments. More elaborate numerical examples for theVlasov-Poisson system as well as ideal and reduced magnetohydrodynamics will be presentedelsewhere [37, 36, 38]. The examples we provide here can be seen as building blocks for thesemore complicated systems. Variational Integrators
In this work, we are concerned with the discretisation of partial differential equations (PDEs) ofevolution type. A field is a map u : X → F from a bounded domain X ⊂ R n +1 taking values inan open set F ⊆ R m . Most often, X corresponds to some region of spacetime with coordinatesx = (x µ ) = (x , x i ) = ( t, x, y, z ) with 0 ≤ µ ≤ n, ≤ i ≤ n, and n = dim X − X being a subset of an Euclidean space.Thus X can be replaced by a differentiable manifold, in which case the results are valid locallyin a coordinate chart. Points on F are denoted by y = (y a ) with 1 ≤ a ≤ m and m = dim F being the number of field components. We make use of the Einstein summation convention onrepeated indices, both for coordinates and fields.The configuration of a field u is geometrically represented by its graphgraph( u ) = (cid:8) (x , y) (cid:12)(cid:12) y = u (x) (cid:9) , which is a subset of the Cartesian product Y = X × F = (cid:8) (x , y) (cid:12)(cid:12) x ∈ X, y ∈ F (cid:9) . With the projection onto the first factor, π : Y → X, (x , y) x , we can define a geometrical structure ( Y, X, π ) called (trivial) fibre bundle. Here, X is calledthe base space, Y the configuration space, and F the fibre. Local coordinates on Y are given by(x µ , y a ) with 0 ≤ µ ≤ n, ≤ a ≤ m. A field u can be identified with a section of the bundle, i.e., a map ϕ : X → Y , satisfying thecondition π ◦ ϕ = id X , where id X is the identity map on X . This condition ensures that the image of a section ϕ corresponds to the graph of a field u , ϕ ( X ) = (cid:8) (x , y) = ϕ (x) (cid:12)(cid:12) x ∈ X (cid:9) = (cid:8) (x , u (x)) (cid:12)(cid:12) x ∈ X (cid:9) . In local coordinates, a section ϕ : X → Y can be written as ϕ (x) = (x µ , ϕ a (x)), that is ϕ (x) isgiven by functions y a = ϕ a (x) = u a (x).Partial derivatives of a field component u a with respect to x µ or x i are denoted u aµ or u ai ,respectively. The collection of all partial derivatives of a given order k is denoted u ( k ) . Theappropriate geometric setting for partial derivatives of a field u is the theory of jet bundles[58, 32, 55, 29, 15]. Jets combine into a single object the values of a field at a point and thevalues of its derivatives. More specifically, the jet prolongation of a section ϕ is a mapj ϕ : x (cid:0) x µ , ϕ a (x) , ϕ aµ (x) (cid:1) , aking values in Y × R m ( n +1) . This space is identified with the first jet bundle of Y , denotedJ Y , so that the first jet prolongation of a field is a section of J Y . Local coordinates on J Y are given by (x µ , y a , z aµ ) with 0 ≤ µ ≤ n, ≤ a ≤ m, where z aµ represents the possible values of partial derivatives. The Lagrangian of a first orderfield theory for example will be a function defined on the first jet bundle J Y . The jet bundleJ Y has a natural projection on both Y and X , denoted by π Y and π X , respectively. The latterdefines a fibre bundle (J Y, X, π X ). It is worth noting that not every section ψ of the bundleJ Y → X is the jet prolongation of a section ϕ : X → Y . When that happens we say that ψ isholonomic.Higher order jet bundles J k Y of Y can be defined analogously by considering the k th jetprolongation of a section ϕ , j k ϕ : x (cid:0) x µ , ϕ a (x) , ϕ aµ (x) , ϕ aµν (x) , ... (cid:1) , which includes derivatives of ϕ up to k th order. Higher-order jet bundles J k Y can also bedefined as submanifolds of iterated first-order jet bundles, e.g., J Y is a submanifold of J (J Y ).If coordinates on J (J Y ) are denoted by(x µ , y a , z aµ , ¯z aµ , w aµν ) with 0 ≤ µ, ν ≤ n, ≤ a ≤ m, then J Y is obtained by requiring that z aµ = ¯z aµ for all a and all µ , and accordingly for larger k (for more details see e.g. [32, Chapter 32] or [15, Section 4.1.4]).Jet bundles have become the standard framework for PDE analysis and variational calculus.They provide a natural setting for the formulation of field theories and the analysis of conser-vation laws which are central to this work. This framework generalises in a standard way tonon-trivial cases in which X is a manifold and Y cannot be written globally as a Cartesianproduct. For a thorough introduction into the geometric framework we refer to Gotay et al.[21]. For definiteness, let us consider only first order Lagrangian field theories, i.e., the Lagrangian L shall depend only on the coordinates, the fields and their first derivatives. Therefore, theLagrangian is a function defined on the first jet bundle, L : J Y → R . (1)The corresponding action functional is A [ ϕ ] = Z X L (cid:0) j ϕ (cid:1) d n +1 x . (2)Hamilton’s principle of stationary action [31, 6, 3] states that among all possible field configu-rations ϕ : X → Y , the one chosen by nature makes the action functional (2) stationary. Asusual in geometric mechanics [23, 45, 21, 48], stationary points of the action are meant in theformal sense. The variations of a field ϕ are defined in terms of a geometric transformation ofthe underlying bundle Y , that is a sufficiently regular one-parameter group of transformations σ (cid:15) (x , y) of Y into itself, defined for (cid:15) in a neighbourhood of zero. In order to have near identitytransformations, it is required that σ (cid:15) reduces to the identity map at (cid:15) = 0, i.e., σ (cid:15) | (cid:15) =0 = id. urthermore, variations should vanish at boundary points, so that σ (cid:15) = id for x ∈ ∂X . Themap σ (cid:15) (x , y) is the flow of the vector field V over Y given by V (x , y) = dd(cid:15) σ (cid:15) (x , y) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 , (3)and written in components as V (x , y) = η a (x , y) ∂∂ y a , (4)where we have η a (x , y) = 0 for x ∈ ∂X . Here, vector fields are identified with first-orderdifferential operators. The variation of a field ϕ in the direction of V is then defined by ϕ (cid:15) = ( σ (cid:15) ◦ ϕ )(x) = σ (cid:15) ( ϕ (x)) , (5)and we have ϕ (cid:15) | (cid:15) =0 = ϕ and ϕ (cid:15) = ϕ on boundary points x ∈ ∂X . Loosely speaking, ϕ (cid:15) isobtained by dragging ϕ along the flow of the vector field V . For sake of simplicity, we considertransformations in σ (cid:15) that leave the point x unchanged, thus, preserving the fibres of the bundle Y . Those are referred to as vertical transformations . We say that ϕ is a stationary point of A [ ϕ ] if for every σ (cid:15) we find that (cid:15) = 0 is a stationary point of A [ ϕ (cid:15) ] viewed as a function of (cid:15) inan open interval about (cid:15) = 0. Then Hamilton’s principle amounts to dd(cid:15) A [ σ (cid:15) ◦ ϕ ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = 0 for every σ (cid:15) . (6)Explicitly, dd(cid:15) A [ σ (cid:15) ◦ ϕ ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = Z X dd(cid:15) h L (cid:0) j ( σ (cid:15) ◦ ϕ ) (cid:1)i(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 d n +1 x (7)= Z X (cid:20) ∂L∂ y a (cid:0) j ϕ (cid:1) · d ( σ (cid:15) ) a d(cid:15) + ∂L∂ z aµ (cid:0) j ϕ (cid:1) · dd(cid:15) ∂ ( σ (cid:15) ) a ∂ x µ (cid:21)(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 d n +1 x . (8)The second term in (8) can be integrated by parts, dd(cid:15) A [ σ (cid:15) ◦ ϕ ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = Z X (cid:20) ∂L∂ y a (cid:0) j ϕ (cid:1) − ∂∂ x µ ∂L∂ z aµ (cid:0) j ϕ (cid:1)! (cid:21) · η a ( ϕ ) d n +1 x , (9)where we used (3) and (4). As the variation of the action has to vanish for arbitrary transfor-mations σ (cid:15) and therefore for arbitrary vector fields V , the term in square brackets has to vanishidentically, that is (cid:0) j ϕ (cid:1) ∗ (cid:0) D EL ( L ) (cid:1) a = ∂L∂ y a (cid:0) j ϕ (cid:1) − ∂∂ x µ ∂L∂ z aµ (cid:0) j ϕ (cid:1)! = 0 . (10)These are the Euler-Lagrange field equations, i.e., the equations of motion for a first orderLagrangian field theory, which for a non-degenerate Lagrangian L live in the second-order jetbundle J Y . By D EL we denote the Euler-Lagrange operator acting on the Lagrangian L .The pull-back notation (j ϕ ) ∗ states that the result is evaluated on the second jet prolongationj ϕ . The theory is fully covariant, that is for time-dependent problems, time is regarded as acomponent of x. As noted in [48, 49], one really should consider general transformations that might affect the pointsof the base space as well as the fields. Otherwise, the Cartan form and thus the multisymplectic form willnot be correctly obtained from the variational principle. However, in this work we are only concernedwith the Euler-Lagrange equations and the Noether theorem. Both will be obtained correctly even ifonly transformations of the fields (vertical transformations in the language of jet bundles) are consid-ered. Transformations of the base space (horizontal transformations) would substantially complicate thederivations and are neglected for sake of simplicity. .3 Continuous Noether Theorem The Noether theorem [53, 33] states that each Lie point symmetry of a Lagrangian correspondsto a conservation law of the associated Euler-Lagrange equations.We restrict our attention to conservation laws that are generated by vertical transformationsof the configuration bundle Y , i.e., transformations which leave the base space X invariant.In the framework of formal Lagrangians addressed below, this is often sufficient to uncoverinteresting conservation laws (including conservation of momentum and energy). Our derivationof the Noether theorem is essentially based on reference [48].As in section 2.2, a transformation is generated by a map σ (cid:15) (x , y), that is ϕ (cid:15) = σ (cid:15) ◦ ϕ with σ (cid:15) | (cid:15) =0 = id and V = dσ (cid:15) d(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 , (11)but it is not required that σ (cid:15) reduces to the identity at boundary points. In our analysis we willusually just prescribe the generating vector field V instead of the actual transformation σ (cid:15) . Incomponents, V can be written as V (x , y) = η a (x , y) ∂∂ y a with η a (x , y) = dd(cid:15) ( σ (cid:15) ) a (x , y) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 , (12)where ( σ (cid:15) ) a is the a -th component of the transformation map. Since Lagrangians are functionsdefined on J Y , we need to compute the first prolongation of the generating vector field. Theprolongation of V is defined via the jet prolongation of the transformation map σ (cid:15) ,j V = dd(cid:15) h j σ (cid:15) (x , y) i(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 , (13)and is given by j V = η a ∂∂ y a + (cid:18) ∂η a ∂ x µ + z bµ ∂η a ∂ y b (cid:19) ∂∂ z aµ = η a ∂∂ y a + η aµ ∂∂ z aµ . (14)A vertical transformation σ (cid:15) is a symmetry transformation for the Lagrangian (1) if the invari-ance condition, L (cid:0) j ( σ (cid:15) ◦ ϕ ) (cid:1) = L (cid:0) j ϕ (cid:1) , (15)is satisfied . Taking the (cid:15) derivative of (15), we obtain an infinitesimal invariance condition, dd(cid:15) L (cid:0) j ( σ (cid:15) ◦ ϕ ) (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = 0 , (16)which is equivalent to (15). Explicitly computing the (cid:15) derivative, we obtainj V (cid:0) L (cid:1)(cid:0) j ϕ (cid:1) = ∂L∂ y a (cid:0) j ϕ (cid:1) · η a ( ϕ ) + ∂L∂ z aµ (cid:0) j ϕ (cid:1) · (cid:20) ∂η a ∂ x µ + ∂η a ∂ y b ∂ϕ b ∂ x µ (cid:21) = 0 . (17) Here, invariance of the Lagrangian suffices, but for general transformations, i.e., transformationsthat transform the base space in addition to the configuration space, invariance of the Lagrangian doesnot guarantee that ϕ (cid:15) is a solution whenever ϕ is a solution. Invariance of solutions requires invarianceof the action, and invariance of the action requires equivariance of the Lagrangian, which also takes thedeformation of the integration domain X into account. For a detailed discussion see, e.g., Lew et al. [41],section 6.5. f ϕ solves the Euler-Lagrange field equations (10), we can replace the first term on the right-hand side of (17) to obtain (cid:20) ∂∂ x µ ∂L∂ z aµ (cid:0) j ϕ (cid:1)(cid:21) · η a ( ϕ ) + ∂L∂ z aµ (cid:0) j ϕ (cid:1) · (cid:20) ∂η a ∂ x µ + ∂η a ∂ y b ∂ϕ b ∂ x µ (cid:21) = 0 . (18)This, at last, amounts to a total divergence,div (cid:2) J (j ϕ ) (cid:3) = 0 , (19)with the Noether current J given by J µ (cid:0) j ϕ (cid:1) = ∂L∂ z aµ (cid:0) j ϕ (cid:1) · η a ( ϕ ) . (20)The fact that the Noether current is divergence-free, expresses the conservation law satisfied bysolutions ϕ of the Euler-Lagrange field equations (10). The flux of J through the boundary ofany domain Ω ⊆ X is zero. Global Form of Conservation Laws
Let us consider an alternative point of view that is better suited for actual computations on thediscrete level. If the Lagrangian is invariant under the vertical transformation (11), the actionis also invariant under this transformation, dd(cid:15) A [ σ (cid:15) ◦ ϕ ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = Z X j V (cid:0) L (cid:1)(cid:0) j ϕ (cid:1) d n +1 x = 0 . (21)Repeating the steps that lead from (17) to (19) this becomes dd(cid:15) A [ σ (cid:15) ◦ ϕ ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = Z X ∂∂ x µ (cid:20) ∂L∂ z aµ (cid:0) j ϕ (cid:1) · η a ( ϕ ) (cid:21) d n +1 x = 0 . (22)With appropriate boundary conditions, one has Z ∂∂ x i (cid:20) ∂L∂ z ai (cid:0) j ϕ (cid:1) · η a ( ϕ ) (cid:21) d n x = 0 , (23)and thus, integrating (22) for x = t ∈ [ t , t ], dd(cid:15) A [ σ (cid:15) ◦ ϕ ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = (cid:20) Z ∂L∂ z at (cid:0) j ϕ (cid:1) · η a ( ϕ ) d n x (cid:21) t t = 0 . (24)This is equivalent to integrating the divergence of the Noether current (19) over the spatialdimensions and using (23). Since t and t are arbitrary, this implies the conservation of J = Z ∂L∂ z at (cid:0) j ϕ (cid:1) · η a ( ϕ ) d n x = Z J t (j ϕ ) d n x , (25)which is called Noether charge. ivergence Symmetries Sometimes it is necessary to consider a slightly more general version of Noether’s theorem[54, 26]. The invariance condition (17) can be weakened toj V ( L )(j ϕ ) = div B (j ϕ ) , (26)with B being a vector field over Y of the form B = B µ (x , y) ∂/∂ x µ , and div B is a scalar fielddefined over J Y given by div B (x , y , z) = ∂B µ (x , y) ∂ x µ + ∂B µ (x , y) ∂ y a z aµ . (27)This form of the invariance condition is tightly related to the gauge freedom of the Lagrangian.We can add any divergence to the Lagrangian without changing the equations of motion as A [ ϕ ] = Z X L (j ϕ ) d n +1 x = Z X (cid:0) L + div H (cid:1) (j ϕ ) d n +1 x = Z X L (j ϕ ) d n +1 x , (28)assuming appropriate boundary conditions. Here, H is a vector field over Y of the form H = H µ (x , y) ∂/∂ x µ , div H is defined as in (27), and L (x , y , z) = (cid:0) L + div H (cid:1) (x , y , z) . (29)One can check that the Euler-Lagrange equations for L and L are the same. Requiring theinvariance condition (17) to be satisfied for L ,j V ( L )(j ϕ ) = j V ( L )(j ϕ ) + div ˜ H (j ϕ ) = 0 , with ˜ H = V ( H µ ) ∂∂x µ , (30)we find the following relation between the vector fields B and H , B = − ˜ H. (31)In summary, if V is a divergence symmetry for L in the sense of (26) and there exists a vectorfield H satisfying (31), then V is a Lie-point symmetry (17) of the equivalent Lagrangian L defined in (29). We can therefore apply the Noether theorem to L , obtaining the Noethercurrent J (j ϕ ) = J (j ϕ ) − B ( ϕ ) , (32)with J the conserved Noether current associated to the divergence symmetry (26) and J asdefined in (20), which is not conserved in this case. In order to derive a discrete version of Hamilton’s action principle, we proceed in three steps.In this section, we define discrete analogues of the base space X , the configuration space Y , andthe jet bundles J Y and J Y . In the next section, we discretise the Lagrangian and the actionfunctional, and finally, we work out the discrete action principle.The continuous base space X is replaced by its discrete analogue X d , a bounded subset X d ⊂ Z n +1 which corresponds to a grid of points in X . In what follows we assume X tobe two-dimensional and choose an equidistant rectangular discretisation like it is depicted inFigure 1(a). The theory is easily applicable to other types of grids as well, e.g., triangular as inFigure 1(b), hexagonal as in Figure 1(c), or even staggered or irregular ones. a) Rectangular Grid (b) Triangular Grid (c) Hexagonal Grid Figure 1: Different regular grids.
The coordinates on the discrete base space X d ⊂ Z × Z are denoted ( i, j ) and assumed totake values i ∈ { , ..., N } and j ∈ { , ..., N } , respectively, where N and N are the numberof points for each dimension. This corresponds to a grid of points x i,j = ( ih, jh ) in X with, forsimplicity, the same step size h in both directions. Then X d ∼ = (cid:8) ( ih, jh ) ∈ X (cid:12)(cid:12) i = 1 , ..., N , j = 1 , ..., N (cid:9) . (33)The discrete configuration space is defined as the cartesian product Y d = X d × F, (34)where F is the same as in the continuous case and we have an analogous projection π d : Y d → X d , (y i,j ) ( i, j ) . (35)Coordinates of Y d are denoted y ai,j with 1 ≤ a ≤ dim F . The coordinates of the base point ( i, j )are already implied and therefore not specified separately. While coordinates on X d and Y d aredefined point-wise, coordinates on J Y d will be defined grid-cell-wise. We therefore introducethe following abstract but convenient notation. A square (cid:3) on X d is an ordered quadruplet (cid:3) = (cid:0) ( i, j ) , ( i, j + 1) , ( i + 1 , j + 1) , ( i + 1 , j ) (cid:1) , (36)defining a primal grid cell. The set of such cells on X d is denoted X (cid:3) . Vertices (cid:3) l of a square (cid:3) with 1 ≤ l ≤ (cid:3) = ( i, j ) , (cid:3) = ( i, j + 1) , (cid:3) = ( i + 1 , j + 1) , (cid:3) = ( i + 1 , j ) . (37)On a quadrilateral grid, the first jet bundle of Y d becomes [48]J Y d = X (cid:3) × F , (38)with coordinates given by (cid:8) ( (cid:3) , y a (cid:3) l ) (cid:12)(cid:12) (cid:3) ∈ X (cid:3) , y a (cid:3) l ∈ R , ≤ l ≤ , ≤ a ≤ dim F (cid:9) . (39)This implies that if (cid:3) = ( i, j ), then we havey (cid:3) = y i,j , y (cid:3) = y i,j +1 , y (cid:3) = y i +1 ,j +1 , y (cid:3) = y i +1 ,j . (40)A section of Y d , representing a discrete field, is a map ϕ d : X d → Y d such that π ◦ ϕ d = id X d .The components of a discrete field ϕ d at the vertices (37) of a specific square are ϕ (cid:3) = ϕ d ( (cid:3) ) , ϕ (cid:3) = ϕ d ( (cid:3) ) , ϕ (cid:3) = ϕ d ( (cid:3) ) , ϕ (cid:3) = ϕ d ( (cid:3) ) , (41) i, j ) ( i, j + 1)( i + 1 , j + 1)( i + 1 , j ) (a) Grid Points y (cid:3) y (cid:3) y (cid:3) y (cid:3) (b) Field Components Figure 2: Vertices of a primal grid cell in a two-dimensional rectangular grid and fieldcomponents at those vertices. as depicted in Figure 2(b). The discrete first jet prolongation of a discrete field ϕ d is a mapj ϕ d : X (cid:3) → J Y d , defined as j ϕ d ( (cid:3) ) = (cid:0) (cid:3) , ϕ (cid:3) , ϕ (cid:3) , ϕ (cid:3) , ϕ (cid:3) (cid:1) . (42)This contains all the information necessary to define discrete first-order derivatives.We will see that the discrete Euler-Lagrange equations are living on J Y d , i.e., the second jetbundle of Y d , similar to the continuous case (see [34] for another discussion of discrete secondorder jet bundles). Coordinates on J Y d will be defined on a quadruple of cells (cid:1) , namely thosecells (cid:3) in X (cid:3) which share a common vertex. The cells (cid:3) m of (cid:1) with 1 ≤ m ≤ (cid:1) = (cid:8) (cid:3) , (cid:3) , (cid:3) , (cid:3) ∈ X (cid:3) (cid:12)(cid:12) (cid:3) = (cid:3) = (cid:3) = (cid:3) (cid:9) . (43)The set of all such quadruples on X d is denoted by X (cid:1) . Vertices (cid:1) l of a quadruple (cid:1) with1 ≤ l ≤ (cid:1) = ( i − , j − , (cid:1) = ( i − , j ) , (cid:1) = ( i − , j + 1) , (cid:1) = ( i, j − , (cid:1) = ( i, j ) , (cid:1) = ( i, j + 1) , (cid:1) = ( i + 1 , j − , (cid:1) = ( i + 1 , j ) , (cid:1) = ( i + 1 , j + 1) . On a quadrilateral grid, the second jet bundle of Y d can be identified withJ Y d = X (cid:1) × F , (44)with coordinates given by (cid:8) ( (cid:1) , y a (cid:1) l ) (cid:12)(cid:12) (cid:1) ∈ X (cid:1) , y a (cid:1) l ∈ R , ≤ l ≤ , ≤ a ≤ dim F (cid:9) , (45)in analogy with (39). The second jet prolongation of a discrete field ϕ d is a map j ϕ d : X (cid:1) → J Y d , defined as j ϕ d ( (cid:1) ) = (cid:0) (cid:1) , ϕ (cid:1) , ϕ (cid:1) , ϕ (cid:1) , ϕ (cid:1) , ϕ (cid:1) , ϕ (cid:1) , ϕ (cid:1) , ϕ (cid:1) , ϕ (cid:1) (cid:1) . (46)Similar to the continuous case, J Y d can also be defined via iteration as the first jet bundle ofJ Y d . If we write coordinates on J (J Y d ) as (cid:8) ( (cid:1) , (cid:3) m , y a (cid:3) lm ) (cid:12)(cid:12) (cid:1) ∈ X (cid:1) , (cid:3) m ∈ (cid:1) , y a (cid:3) lm ∈ R , ≤ l, m ≤ , ≤ a ≤ dim F (cid:9) , (47)then J Y d consists of all elements of J (J Y d ), which satisfyy a (cid:3) = y a (cid:3) = y a (cid:3) = y a (cid:3) , y a (cid:3) = y a (cid:3) , y a (cid:3) = y a (cid:3) , y a (cid:3) = y a (cid:3) , y a (cid:3) = y a (cid:3) , (48)for all a , which in our construction is always guaranteed due to (40) and (43). Therefore, J Y d can be identified with J (J Y d ), unlike the continuous case, where J Y is strictly embedded intoJ (J Y ). .5 Discrete Action Principle The discretisation of the Lagrangian is based on the observation that the continuous actionfunctional can be written as A [ ϕ ] = X (cid:3) ∈ X (cid:3) Z Vol( (cid:3) ) L (cid:0) j ϕ (cid:1) d n +1 x , (49)where Vol( (cid:3) ) ⊂ X is the physical domain enclosed by (cid:3) . The integral in (49) is approximatedby a function of values of ϕ d in four different points in the spacetime grid, which corresponds tothe discrete jet prolongation defined in (42). This is the discrete Lagrangian L d , i.e., Z Vol( (cid:3) ) L (cid:0) j ϕ (cid:1) d n +1 x ≈ L d (cid:0) j ϕ d ( (cid:3) ) (cid:1) = L d (cid:0) (cid:3) , ϕ (cid:3) , ϕ (cid:3) , ϕ (cid:3) , ϕ (cid:3) (cid:1) . (50)The action functional (49) is then approximated by A d [ ϕ d ] = X (cid:3) ∈ X (cid:3) L d (cid:0) j ϕ d ( (cid:3) ) (cid:1) , (51)which can also be written explicitly as A d [ ϕ d ] = N − X i =1 N − X j =1 L d (cid:0) ϕ i,j , ϕ i,j +1 , ϕ i +1 ,j +1 , ϕ i +1 ,j (cid:1) , (52)but for most of our derivations the abstract notation is more practical.Specifically, the discrete Lagrangian is obtained by introducing a quadrature rule (e.g., trape-zoidal, midpoint, Simpson) to approximate the integral in (50) as well as approximations of thefields and their derivatives. In the spirit of the Veselov discretisation [62, 63, 51, 48], we will usethe midpoint rule and first-order finite differences. This entails that the continuous fields in theLagrangian are replaced with averages of field values at the four vertices of the grid cell (cid:3) , i.e.,y a → (cid:16) y a (cid:3) + y a (cid:3) + y a (cid:3) + y a (cid:3) (cid:17) ≡ y a ( (cid:3) ) . (53)There are two possibilities for the definition of each of the derivatives. With reference to Fig-ure 2(b), let the vertical dimension correspond to x (time) and the horizontal dimension to x (space). Then ∂ can be defined along the left as well as along the right edge of the grid cell.Similarly, ∂ can be defined along the upper as well as along the lower edge. Best results areusually obtained for the most symmetric discretisation of the Lagrangian. Therefore we use theaverage of the respective options and replace the derivatives in the Lagrangian according toz a → (cid:18) y a (cid:3) − y a (cid:3) h + y a (cid:3) − y a (cid:3) h (cid:19) ≡ z a ( (cid:3) ) , (54)z a → (cid:18) y a (cid:3) − y a (cid:3) h + y a (cid:3) − y a (cid:3) h (cid:19) ≡ z a ( (cid:3) ) . (55)At last, in order to obtain the discrete equations of motion, we only need to apply a discreteversion of Hamilton’s principle of stationary action. The vertical transformation σ (cid:15) is discretisedin the same way as the fields ϕ , by considering its values over points ( i, j ) of the discrete basespace X d , i.e., we replace σ (cid:15) (x , y) with σ (cid:15)d = (cid:0) σ (cid:15)i,j (y i,j ) (cid:1) . The generating vector field is computedas in (3), but on grid points only, η ai,j = d ( σ (cid:15)i,j ) a d(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 . (56) i − ,j − y i,j − y i +1 ,j − y i − ,j y i,j y i +1 ,j y i − ,j +1 y i,j +1 y i +1 ,j +1 Figure 3: Contributions to the discrete Euler-Lagrange field equations. The derivativesof the discrete Lagrangian at the highlighted grid cells with respect to y i,j correspond tothe terms in equation (59).
The discrete version of the action principle (6) then reads dd(cid:15) A d [ σ (cid:15)d ◦ ϕ d ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = 0 for every σ (cid:15)d . (57)The explicit computation of (57) leads to dd(cid:15) A d [ σ (cid:15)d ◦ ϕ d ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = X (cid:3) ∈ X (cid:3) dd(cid:15) L d (cid:0) j ϕ (cid:15)d ( (cid:3) ) (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = X (cid:3) ∈ X (cid:3) X l =1 ∂L d ∂ y a (cid:3) l (cid:0) j ϕ d ( (cid:3) ) (cid:1) · η a (cid:3) l ( ϕ (cid:3) l ) . (58)As the variation of the action has to vanish for each η ai,j on the spacetime grid independently,it is sufficient to consider only those contributions that are multiplied with the vector field at afixed grid point ( i, j ). In total there are four such contributions (see Figure 3), dd(cid:15) A d [ σ (cid:15)d ◦ ϕ d ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = ... + ∂L d ∂ y a (cid:3) (cid:16) ϕ i,j , ϕ i,j +1 , ϕ i +1 ,j +1 , ϕ i +1 ,j (cid:17) · η ai,j ( ϕ i,j ) + ...... + ∂L d ∂ y a (cid:3) (cid:16) ϕ i,j − , ϕ i,j , ϕ i +1 ,j , ϕ i +1 ,j − (cid:17) · η ai,j ( ϕ i,j ) + ...... + ∂L d ∂ y a (cid:3) (cid:16) ϕ i − ,j − , ϕ i − ,j , ϕ i,j , ϕ i,j − (cid:17) · η ai,j ( ϕ i,j ) + ...... + ∂L d ∂ y a (cid:3) (cid:16) ϕ i − ,j , ϕ i − ,j +1 , ϕ i,j +1 , ϕ i,j (cid:17) · η ai,j ( ϕ i,j ) + ... = 0 . The variation of the discrete action vanishes, if the sum of all expressions multiplying η ai,j vanishesidentically for all ( i, j ) (since the ( σ (cid:15)i,j ) a and therefore the η ai,j are arbitrary). For each a , this equirement yields the discrete Euler-Lagrange field equations at ( i, j ),0 = ∂L d ∂ y a (cid:3) (cid:16) ϕ i,j , ϕ i,j +1 , ϕ i +1 ,j +1 , ϕ i +1 ,j (cid:17) + ∂L d ∂ y a (cid:3) (cid:16) ϕ i,j − , ϕ i,j , ϕ i +1 ,j , ϕ i +1 ,j − (cid:17) + ∂L d ∂ y a (cid:3) (cid:16) ϕ i − ,j − , ϕ i − ,j , ϕ i,j , ϕ i,j − (cid:17) + ∂L d ∂ y a (cid:3) (cid:16) ϕ i − ,j , ϕ i − ,j +1 , ϕ i,j +1 , ϕ i,j (cid:17) , (59)which can be compactly written as an equation on the discrete second jet bundle J Y d , namely (cid:0) j ϕ d ( (cid:1) ) (cid:1) ∗ (cid:0) D EL ( L d ) (cid:1) a = X l =1 (cid:3) l = (cid:1) ∂L d ∂ y a (cid:3) l (cid:0) j ϕ d ( (cid:3) ) (cid:1) = 0 for all a and all (cid:1) . (60)Here, D EL ( L d ) is the discrete Euler-Lagrange operator acting on the discrete Lagrangian L d ,which is evaluated on the prolongation j ϕ d ( (cid:1) ) as indicated by the pull-back notation. Theserelations define the variational integrator for a first-order Lagrangian field theory according tothe Veselov discretisation of the Lagrangian as it was described above. Following the derivation of the continuous theory from the previous section, we consider verticaltransformations σ (cid:15)i,j (y i,j ) and define ϕ (cid:15)i,j = σ (cid:15)i,j ◦ ϕ i,j with σ (cid:15)i,j | (cid:15) =0 = id and η ai,j = d ( σ (cid:15)i,j ) a d(cid:15) (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 . (61)The transformation σ (cid:15)i,j is a symmetry for the discrete Lagrangian (50) if L d (cid:0) j ϕ (cid:15)d ( (cid:3) ) (cid:1) = L d (cid:0) j ϕ d ( (cid:3) ) (cid:1) , (62)which is equivalent to the infinitesimal symmetry condition dd(cid:15) L d (cid:0) j ϕ (cid:15)d ( (cid:3) ) (cid:1)(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = X l =1 ∂L d ∂ y a (cid:3) l (cid:0) j ϕ d ( (cid:3) ) (cid:1) · η a (cid:3) l (cid:0) ϕ (cid:3) l (cid:1) = 0 . (63)On each grid cell (cid:3) we can define four discrete momentum maps in analogy to (20), (cid:0) j ϕ d ( (cid:3) ) (cid:1) ∗ J (cid:3) l = ∂L d ∂ y a (cid:3) l (cid:0) j ϕ d ( (cid:3) ) (cid:1) · η a (cid:3) l (cid:0) ϕ (cid:3) l (cid:1) , l ∈ { , , , } . (64)Instead of two components of J , corresponding to the coordinates of the base space, we nowhave four contributions, corresponding to the vertices of a grid cell. With that, we can writethe discrete symmetry condition (63) as J (cid:3) + J (cid:3) + J (cid:3) + J (cid:3) = 0 , for every (cid:3) ∈ X (cid:3) . (65)In addition to (64), we define some explicit shorthand notation that will become useful in thefollowing derivation, J i,j = ∂L d ∂ y a (cid:3) (cid:16) ϕ i,j , ϕ i,j +1 , ϕ i +1 ,j +1 , ϕ i +1 ,j (cid:17) · η ai,j ( ϕ i,j ) , (66a) J i,j = ∂L d ∂ y a (cid:3) (cid:16) ϕ i,j , ϕ i,j +1 , ϕ i +1 ,j +1 , ϕ i +1 ,j (cid:17) · η ai,j +1 ( ϕ i,j +1 ) , (66b) J i,j = ∂L d ∂ y a (cid:3) (cid:16) ϕ i,j , ϕ i,j +1 , ϕ i +1 ,j +1 , ϕ i +1 ,j (cid:17) · η ai +1 ,j +1 ( ϕ i +1 ,j +1 ) , (66c) = 1 i = N j − j j + 1 j + 2 N − Figure 4: Contributions to the discrete Noether theorem. J i,j = ∂L d ∂ y a (cid:3) (cid:16) ϕ i,j , ϕ i,j +1 , ϕ i +1 ,j +1 , ϕ i +1 ,j (cid:17) · η ai +1 ,j ( ϕ i +1 ,j ) . (66d)In order to derive the discrete conservation law, we follow the argument from section 2.3.If the Lagrangian is invariant under a given vertical transformation, then the action is alsoinvariant under this transformation, namely, dd(cid:15) A d (cid:2) ϕ (cid:15)d (cid:3)(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = 0 , (67)which becomes dd(cid:15) A d (cid:2) ϕ (cid:15)d (cid:3)(cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 = X (cid:3) ∈ X (cid:3) X l =1 ∂L d ∂ y a (cid:3) l (cid:0) j ϕ d ( (cid:3) ) (cid:1) · η a (cid:3) l (cid:0) ϕ (cid:3) l (cid:1) = X (cid:3) ∈ X (cid:3) X l =1 (cid:0) j ϕ d ( (cid:3) ) (cid:1) ∗ J (cid:3) l = 0 . (68)Those contributions to the sum that originate from interior points vanish in virtue of the discreteEuler-Lagrange field equations (59). Specifically, J i,j + J i,j − + J i − ,j − + J i − ,j = (cid:0) j ϕ d ( (cid:1) ) (cid:1) ∗ (cid:0) D EL ( L d ) (cid:1) ai,j · η ai,j ( ϕ i,j ) = 0 , (69)so that only boundary cells contribute to (68), c.f. Figure (4). For simplicity, we will assumeperiodic boundary conditions in the space-like variable, then the only contribution comes fromthe boundary of the time-like variable,0 = N − X j =1 h J ,j + J ,j − + J N − ,j − + J N − ,j i . (70)Using the discrete symmetry condition (65) to replace the first two terms as well as the periodicityof the spatial domain, we obtain the discrete counterpart of (24), that is the discrete conservationlaw N − X j =1 h J ,j + J ,j i = N − X j =1 h J N − ,j + J N − ,j i , (71) ee Figure 4. Since the number of timesteps N is arbitrary, this is equivalent to J d = N − X j =1 (cid:3) =( i,j ) h J (cid:3) (cid:0) j ϕ d ( (cid:3) ) (cid:1) + J (cid:3) (cid:0) j ϕ d ( (cid:3) ) (cid:1)i = const. for all i, (72)in analogy to (25). This proves that variational integrators preserve discrete invariants J d tomachine accuracy . Discrete Divergence Symmetries
We can generalise the discrete Noether theorem in a similar way as we did with the contin-uous Noether theorem for divergence symmetries. In analogy to (26), we rewrite the discretesymmetry condition (65) as J (cid:3) + J (cid:3) + J (cid:3) + J (cid:3) = B (cid:3) + B (cid:3) + B (cid:3) + B (cid:3) . (73)Making no assumption on the B (cid:3) l , equation (70) takes the form N − X j =1 h J ,j + J ,j − − B ,j − B ,j − + J N − ,j − + J N − ,j − B N − ,j − − B N − ,j i == N − X i =2 N − X j =1 h B i,j + B i,j − + B i − ,j − + B i − ,j i . (74)We conclude, that we have the discrete equivalent of a divergence symmetry, if the sum of all B li,j in the interior of X d vanishes, N − X i =2 N − X j =1 h B i,j + B i,j − + B i − ,j − + B i − ,j i = 0 , (75)such that only the fluxes at the boundaries of X d contribute to (74). Then, the generalisedNoether charge (72) becomes J d = N − X j =1 (cid:3) =( i,j ) h J (cid:3) (cid:0) j ϕ d ( (cid:3) ) (cid:1) + J (cid:3) (cid:0) j ϕ d ( (cid:3) ) (cid:1) − B (cid:3) (cid:0) j ϕ d ( (cid:3) ) (cid:1) − B (cid:3) (cid:0) j ϕ d ( (cid:3) ) (cid:1)i = const. for all i, (76)in correspondence with (32). In this section we review the idea of formal (or adjoint) Lagrangians [8, 25] and the correspondingNoether theorem as introduced by Ibragimov [26]. We provide a geometric view of this theoryin line with the formalism adopted in the rest of the paper, so that the application of the resultsof section 2 to formal Lagrangians becomes natural. At last, we propose a discrete version ofIbragimov’s extension of Noether’s theorem, which allows us to prove discrete conservation lawsfor variational integrators obtained from formal Lagrangians. In practice, many variational integrators require the solution of nonlinear algebraic equations. Inthat case, the accuracy to which the invariants are preserved will generally depend on the accuracy ofthe nonlinear solver. .1 Inverse Problem of the Calculus of Variations Consider a generic system of nonlinear partial differential equations, F (x , u, u (1) , ..., u ( k ) ) = 0 , (77)for a field u : X → F , where u ( k ) denotes all derivatives of order k . Here, F is a (sufficientlyregular) function of the point x ∈ X , the field u (x), and its derivatives up to order k , takingvalues in R m , where m = dim F , so that we have as many equations as variables. The componentsof F are denoted F a with 1 ≤ a ≤ m .Let B be the space of functions u : X → F where the solution of (77) is sought. Weassume that B is a Banach space and denote by B ∗ its topological dual, i.e., the space of linearcontinuous functionals from B to R . The specific choice of the space B depends on the problemat hand. For definiteness, we require B , → (cid:2) L ( X ) (cid:3) m , i.e., B is continuously embedded in (cid:2) L ( X ) (cid:3) m . If the function x
7→ F (x , u, u (1) , ..., u ( k ) ) defined in (77) belongs to (cid:2) L ( X ) (cid:3) m , thenit can be identified with the mapping F : B → B ∗ defined by hF [ u ] , v i = Z X v (x) · F (x , u, u (1) , ..., u ( k ) ) d n +1 x , (78)where h· , ·i : B ∗ × B → R is the duality pairing between elements of B ∗ and B . With some abuseof notation we denote by F both the function in (77), which is defined on a finite-dimensionalspace, and the mapping F : B → B ∗ , which is defined on a Banach space.The nonlinear operator F admits a natural variational formulation if there exists a functional A : B → R of class C ( B ), such that F = D A , where D denotes the Fréchet derivative operator.We recall that the Fréchet derivative of A : B → R is a map D A : B → B ∗ , such that for every u, v ∈ B , h D A [ u ] , v i coincides with the Gateaux derivative of A along v [54, 50], that is h D A [ u ] , v i = lim (cid:15) → A [ u + (cid:15)v ] − A [ u ] (cid:15) = dd(cid:15) A [ u + (cid:15)v ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 . (79)In this framework, the functional derivative δ A [ u ] /δu , if it exists, is defined as the L -represen-tative of the Fréchet derivative. Specifically, δ A [ u ] /δu is the unique element of (cid:2) L ( X ) (cid:3) m , if itexists, such that h D A [ u ] , v i = Z X δ A [ u ] δu · v d x , (80)for all v ∈ B [4]. Here, we are implicitly assuming that boundary conditions on the field u areencoded in the definition of the linear space B ; this is possible in particular for both homogeneousDirichlet and periodic boundary conditions.Analogously to (79), the Fréchet derivative at u ∈ B of a functional F : B → B ∗ of class C is the linear continuous operator D F [ u ] : B → B ∗ , such that for every v ∈ B , D F [ u ] v coincideswith the Gateaux derivative of F along v , i.e., D F [ u ] v = dd(cid:15) F [ u + (cid:15)v ] (cid:12)(cid:12)(cid:12)(cid:12) (cid:15) =0 . (81)Therefore, D F : B → L ( B, B ∗ ) ∼ = B ∗ ⊗ B ∗ , where L ( B, B ∗ ) is the space of linear continuousoperators from B to B ∗ . If F admits a variational formulation, then the second-order Fréchetderivative D F = D A : B → B ∗ ⊗ B ∗ should be symmetric at each u , hence a necessarycondition for the existence of the action A is h D F [ u ] v, w i = h D F [ u ] w, v i ∀ u, v, w ∈ B. (82) t turns out that this condition is also sufficient [54, 11, 60].Unfortunately, for many interesting systems this condition is not fulfilled. Some of thesesystems admit a variational principle after a variable transformation. This is the case for manyequations from fluid dynamics [59]. But such transformations can be inconvenient since the newvariables might suffer from problems with respect to non-uniqueness, boundary conditions orregularity. In addition, the resulting variational principle might be subject to constraints on thevariations which are not easily dealt with at the discrete level. Instead, they call for extensionsof the theory that complicate its application [56, 20]. For a still larger class of systems, avariational formulation is not known, even after a change of coordinates. Nevertheless, we canderive variational integrators for these systems by considering the following construction. A generic system of differential equations (77) can be treated as part of a Lagrangian system,constructed by doubling the number of variables. The action of such a Lagrangian system is A [ u, v ] = hF [ u ] , v i , (83)with u, v ∈ B and the pairing as defined in (78). Then, the Lagrangian is L (x , u, u (1) , ..., u ( k ) , v, v (1) , ..., v ( k ) ) = v (x) · F (x , u, u (1) , ..., u ( k ) ) . (84)This is referred to as formal (or adjoint) Lagrangian [8, 25]. The variational principle (6) appliedto (84) gives the original equation, δ A [ u, v ] δv = F [ u ] = 0 , (85)as well as the adjoint equation , which is defined by δ A [ u, v ] δu = F ∗ [ u, v ] = 0 , (86)assuming that the functional derivative exists.Summarising, if we consider a dynamical system described by the variables ( u, v ), whosedynamics are governed by the equations F [ u ] = 0 and F ∗ [ u, v ] = 0, then this system has avariational formulation and contains the original system (77) as a subsystem. In the following,we refer to the system (85-86) as extended system describing the dynamics of the variables ( u, v ). Geometry of Formal Lagrangians
As before, we identify a field u with a section ϕ : X → Y of the configuration bundle Y . Withj k ϕ the k th jet prolongation of ϕ , we can write the system of equations (77) as a function onthe jet bundle J k Y , F (j k ϕ ) = 0 . (87)The configuration bundle ˜ Y of the extended system is given by the Cartesian product˜ Y = X × F × F, (88)together with the projection ˜ π : ˜ Y → X, (x , y , ˜y) x . (89) ocal coordinates on ˜ Y are denoted by (x µ , y a , ˜y a ) with 1 ≤ a ≤ m = dim F . Sections ˜ ϕ : X → ˜ Y are at the same time associated to the graph of u and v ,˜ ϕ : x (cid:0) x µ , ϕ a (x) , ˜ ϕ a (x) (cid:1) = (cid:0) x µ , u a (x) , v a (x) (cid:1) , (90)where the ϕ a are associated with the physical variables u a and the ˜ ϕ a are associated with theadjoint variables v a . With j k ˜ ϕ the k th jet prolongation of ˜ ϕ , we can write the formal Lagrangian(84) as a function on the jet bundle J k ˜ Y , L (cid:0) j k ˜ ϕ (cid:1) = v · F (x , u, u (1) , ..., u ( k ) ) , (91)and the formal action functional (83) as A [ ˜ ϕ ] = Z X L (cid:0) j k ˜ ϕ (cid:1) d n +1 x . (92)In the following we shall consider first-order Lagrangians only, i.e., k = 1. Ibragimov [26] has shown that the adjoint equations (86) inherit all symmetries of the originalequations (85). Therefore it is possible to determine conservation laws of any system of differen-tial equations, even without a natural Lagrangian, by application of the Noether theorem fromsection 2.3 to the corresponding formal Lagrangian (91). Again, we restrict our attention to thecase of first order field theories.More specifically, if the the original system (77) admits a symmetry related to the infinites-imal vector field V = η a (x , y) ∂∂ y a , (93)then the extended system (85-86) admits a symmetry related to the vector field˜ V = η a (x , y) ∂∂ y a + ˜ η a (x , y , ˜y) ∂∂ ˜y a , (94)with appropriately chosen coefficients ˜ η a [26]. By definition, the vector field (93) describes a Liepoint symmetry of (77), if there exists a matrix-valued function λ = ( λ ba (x , y , z)), such thatj V ( F a ) = λ ba F b , (95)where F (x , y , z) is the function defining the system of first order equations F ( x, u, u (1) ). Weapply the jet prolongation of the extended vector field (94) to the Lagrangian (91) of the extendedsystem (85-86), treating ˜ η a as unknown coefficients,j ˜ V ( L ) = ˜ η a F a + ˜y a j V ( F a ) = ˜ η a F a + ˜y a λ ba F b . (96)We want to choose the ˜ η a so that the symmetry condition (17) is satisfied. We see that it issufficient to set ˜ η a = − λ ab ˜y b (97)for (96) to vanish. The extended vector field (94) therefore becomes˜ V = η a ∂∂ y a − λ ab ˜y b ∂∂ ˜y a , (98) ith λ given in the symmetry condition (95). In the line of the standard Noether theorem, c.f.section 2.3, if (94) is a symmetry of the formal Lagrangian (91), we obtain the local conservationlaw div ˜ J (j ˜ ϕ ) = 0 , (99)for all solutions ˜ ϕ of the Euler-Lagrange equations on ˜ Y , with the Noether current ˜ J given by˜ J µ (cid:0) j ˜ ϕ (cid:1) = ∂L∂ z aµ (cid:0) j ˜ ϕ (cid:1) · η a ( ˜ ϕ ) + ∂L∂ ˜z aµ (cid:0) j ˜ ϕ (cid:1) · ˜ η a ( ˜ ϕ ) . (100)For formal Lagrangians of the form (91), the second terms on the right-hand side are zero.Sometimes, however, we might want to symmetrise the Lagrangian (see e.g. section 4.3), inwhich case these terms do contribute to the Noether current. From equation (100), we obtainthe (global) Noether charge (25),˜ J = Z (cid:20) ∂L∂ z at (cid:0) j ˜ ϕ (cid:1) · η a ( ˜ ϕ ) + ∂L∂ ˜z at (cid:0) j ˜ ϕ (cid:1) · ˜ η a ( ˜ ϕ ) (cid:21) d n x , (101)which is constant for all times t . Usually, formal Lagrangians of the form (91) do not feature atime derivative of the adjoint variables, so that the Noether charge simplifies accordingly. Restriction of Conservation Laws
The conservation law thus obtained is of course a conservation law of the extended system (85-86), and therefore generally depends on both u and v . In order to restrict conservation lawsof the extended systems to conservation laws of the physical system, we have to find a suitableway of restricting solutions of the extended system to solutions of the physical system. For thatpurpose, the following two definitions will become useful.If the adjoint equation (86), restricted to v = u , becomes equivalent to the original equation(77), i.e., F ∗ [ u, u ] = λ F [ u ] , (102)for some matrix λ , possibly depending on the fields and their derivatives, the system (77) iscalled self-adjoint in the sense of Ibragimov [25]. If F [ u ] is a linear operator and λ the identitymatrix m × m , the above definition coincides with the standard definition of formally self-adjointoperators. In general, however, the systems we are considering are not self-adjoint. Ibragimovrelaxed the requirement of self-adjointness by introducing the concept of quasi-self-adjointness[27, 28]. This is a generalisation of self-adjointness where v = u is generalised to v = φ ( u ) forsome function φ : F → F . The advantage of self-adjointness or quasi-self-adjointness in thesense of Ibragimov is the possibility to build a solution of the full extended system (85-86) givena solution of the original problem (85).Given a diffeomorphism φ : F → F of F into itself, we can build the embeddingΦ : Y , → ˜ Y , (cid:0) x , y (cid:1) (cid:0) x , y , φ (y) (cid:1) . (103)In order to use this embedding to restrict the Noether current to the physical system, we needto lift Φ to a map from J Y to J ˜ Y . With this aim, let us consider a generic section ϕ of Y forwhich the composed map ˜ ϕ = Φ ◦ ϕ : X → ˜ Y (104) mounts to ˜ ϕ : x (cid:0) x , u (x) , φ ( u (x)) (cid:1) = (cid:0) x , u (x) , v (x) (cid:1) , (105)where u (x) is the field component corresponding to ϕ , and we have defined the second field v = φ ◦ u : X → F . It follows that the condition ˜ ϕ ◦ ˜ π = id X is satisfied and ˜ ϕ = Φ ◦ ϕ is indeeda section of ˜ Y , i.e., the composition with Φ maps sections of Y into sections of ˜ Y . By the chainrule, we compute the first jet of ˜ ϕ , that isj ˜ ϕ (x) = (cid:0) x , u (x) , v (x) , Du (x) , Dv (x) (cid:1) = (cid:0) x , u (x) , φ ( u (x)) , Du (x) , Dφ ( u (x)) · Du (x) (cid:1) , (106)where Df denotes the Jacobian matrix of a function f . We now can define the liftj Φ : J Y → J ˜ Y (107)by j Φ (x , y , z) = (cid:0) x , y , φ (y) , z , Dφ (y) · z (cid:1) , (108)which by construction satisfies the identityj (cid:0) Φ ◦ ϕ (cid:1) = j Φ ◦ j ϕ, (109)as desired. With an abuse of notation, the symbol j Φ is not used here in the usual sense of jetprolongation, instead j Φ is a lift of Φ up to the first jet bundle. We can now use (108) to pullback the Noether current in (99) if we assume that ˜ ϕ can be realised in the form ˜ ϕ = Φ ◦ ϕ ,˜ J (cid:0) j ˜ ϕ (cid:1) = ˜ J (cid:0) j (Φ ◦ ϕ ) (cid:1) = ˜ J (cid:0) j Φ ◦ j ϕ (cid:1) = (cid:0) ˜ J ◦ j Φ (cid:1)(cid:0) j ϕ (cid:1) . (110)If ˜ ϕ = Φ ◦ ϕ solves the equations of the extended system (85-86), then div ˜ J (cid:0) j ˜ ϕ (cid:1) = 0. Upondefining the restricted Noether current by J = ˜ J ◦ j Φ , (111)the Noether theorem (99) takes the formdiv J (j ϕ ) = 0 , (112)which expresses the local conservation law for the physical system.This result is crucial for the application of variational integrators to formal Lagrangians.Without the construction of a solution ( u, φ ( u )) of the extended system from a solution u ofthe physical system, it is in general not possible to determine the discrete momenta that areconserved by the variational integrator due to symmetries of the physical system. Generalisations
To simplify the derivations and the analysis of the conservation laws, it is sometimes usefulto add a term G (j k ϕ ) to the Lagrangian, that is a function of the coordinates, the physicalvariables and their derivatives, but not of the adjoint variables and their derivatives. Themodified Lagrangian, L (cid:0) j k ˜ ϕ (cid:1) = L (cid:0) j k ˜ ϕ (cid:1) + G (cid:0) j k ˜ ϕ (cid:1) with ∂G∂ ˜y a (cid:0) j k ˜ ϕ (cid:1) = 0 and ∂G∂ ˜z aµ (cid:0) j k ˜ ϕ (cid:1) = 0 for all µ, a, (113) ields the same physical equations of motion as L , but in general will lead to different adjointequations. This freedom can be used to simplify the search for the embedding Φ used to restrictthe conservation laws of the extended system to the physical system.For example, for equations of advection type, u t + A ( u, u x ) = 0 , (114)where A is some possibly nonlinear function of the field u and its spatial derivatives, it is possibleto construct an extended system which is guaranteed to be always self-adjoint. It suffices tochoose G = − u · A , so that the Lagrangian becomes L (cid:0) j k ˜ ϕ (cid:1) = v · u t + v · A ( u, u x ) − u · A ( u, u x ) . (115)While this simplifies the construction of the embedding Φ, it complicates the extension of sym-metries. In general it will not be possible to determine the components of the generating vectorfield in an algorithmic way as in (97). Discrete Embedding
At this point, the discretisation of formal Lagrangians can be carried out straightforwardly asan application of the theory reviewed in section 2. The only missing element is the restrictionof discrete conservation laws of the extended system to the physical system.Since the fibres of Y d are copies of F , we can consider a diffeomorphism φ : F → F , as inthe continuous case and define the discrete embeddingΦ d : Y d , → ˜ Y d , (cid:0) y i,j (cid:1) (cid:0) y i,j , φ (y i,j ) (cid:1) , (116)and its lift to the discrete first jet bundle,j Φ d : J Y d → J ˜ Y d , (117)which is given byj Φ d ( (cid:3) ) = (cid:0) (cid:3) , y (cid:3) , y (cid:3) , y (cid:3) , y (cid:3) , φ (y (cid:3) ) , φ (y (cid:3) ) , φ (y (cid:3) ) , φ (y (cid:3) ) (cid:1) . (118)With this we can pullback discrete conservation laws (72) of the extended system,˜ J d = N − X j =1 (cid:3) =( i,j ) h ˜ J (cid:3) (cid:0) j ˜ ϕ d ( (cid:3) ) (cid:1) + ˜ J (cid:3) (cid:0) j ˜ ϕ d ( (cid:3) ) (cid:1)i = const. for all i, (119)to the physical system by defining the restricted discrete momentum map J (cid:3) l = ˜ J (cid:3) l ◦ j Φ d ( (cid:3) ) , (120)such that the restricted Noether charge becomes J d = N − X j =1 (cid:3) =( i,j ) h J (cid:3) (cid:0) j ϕ d ( (cid:3) ) (cid:1) + J (cid:3) (cid:0) j ϕ d ( (cid:3) ) (cid:1)i = const. for all i. (121)It is worth noticing that we do not pull back the Noether charge ˜ J but the momentum maps˜ J (cid:3) l , just as in the continuous case (111), where we pull back the Noether current ˜ J . Applications
In this section we present two applications, the linear advection equation and the vorticityequation. The linear advection equation in one spatial dimension is a prototypical examplethat shares many characteristics with more complicated systems from plasma physics and fluiddynamics. The vorticity equation in two spatial dimensions is an important equation of fluiddynamics which is widely used for example in atmospheric sciences but also in plasma physics.It consists of an advection equation which is coupled to a Poisson equation.On the continuous level, we construct the formal Lagrangian and compute variations of thecorresponding action functional. We check the resulting equations of motion for self-adjointnessor quasi-self-adjointness in order to define the restriction map φ . Then we check the equationsfor symmetries and compute the extended generating vector field. This vector field is applied tothe Lagrangian in order to compute the corresponding conservation laws. After discretising theformal Lagrangian, we derive the actual variational integrator and check the discrete equationsof motion for discrete self-adjointness or quasi-self-adjointness. Then we check the discreteLagrangian for symmetries and compute the discrete conservation laws. We will discuss differentdiscretisations of the Lagrangian, especially with respect to different quadrature rules, and theirconsequences for properties of the resulting variational integrator (e.g., explicit vs. implicit) andthe corresponding conservation laws. Discretisation of Formal Lagrangians
Before we proceed, we want to give a short comment on the discretisation of formal Lagrangians.A straight forward application of the Veselov discretisation, i.e., using the midpoint quadraturerule for both time and space integration, as it was presented in section 2.5 and as it is used inmany works on variational integrators, does not appear to be suitable for formal Lagrangiansdue to their particular structure. As we will see in section 4.2 on the linear advection equation,it results in numerical schemes, which are prone to unphysical oscillations in space. Instead,when discretising the derivatives with first-order finite differences, the spatial integral can beapproximated by the trapezoidal or Simpson rule. Another option to avoid this issue is to moveto a slightly more abstract Galerkin framework and use e.g. Lagrange polynomials [16] or splines[35] as basis functions.
So far, we have only considered the midpoint rule (Veselov discretisation) for approximating theLagrangian L ( u, u t , u x ), namely, L mp d (cid:0) j ϕ d ( (cid:3) ) (cid:1) = h t h x L (cid:18) (cid:16) ϕ (cid:3) + ϕ (cid:3) + ϕ (cid:3) + ϕ (cid:3) (cid:17) , (cid:16) ϕ (cid:3) − ϕ (cid:3) h t + ϕ (cid:3) − ϕ (cid:3) h t (cid:17) , (cid:16) ϕ (cid:3) − ϕ (cid:3) h x + ϕ (cid:3) − ϕ (cid:3) h x (cid:17)(cid:19) . (122)As will be discussed below, for formal Lagrangians the midpoint rule does not appear to lead tostable variational integrators. Therefore we will also consider some alternative discretisations of he Lagrangian. First, the trapezoidal rule, L tr d (cid:0) j ϕ d ( (cid:3) ) (cid:1) = h t h x (cid:18) L (cid:16) ϕ (cid:3) , ϕ (cid:3) − ϕ (cid:3) h t , ϕ (cid:3) − ϕ (cid:3) h x (cid:17) + L (cid:16) ϕ (cid:3) , ϕ (cid:3) − ϕ (cid:3) h t , ϕ (cid:3) − ϕ (cid:3) h x (cid:17) + L (cid:16) ϕ (cid:3) , ϕ (cid:3) − ϕ (cid:3) h t , ϕ (cid:3) − ϕ (cid:3) h x (cid:17) + L (cid:16) ϕ (cid:3) , ϕ (cid:3) − ϕ (cid:3) h t , ϕ (cid:3) − ϕ (cid:3) h x (cid:17)(cid:19) , (123)which leads to the (explicit) Leapfrog scheme. Second, a combination of both, the midpoint rulefor time integration and the trapezoidal rule for spatial integration, L mt d (cid:0) j ϕ d ( (cid:3) ) (cid:1) = h t h x L (cid:18) ϕ (cid:3) + ϕ (cid:3) , ϕ (cid:3) − ϕ (cid:3) h t , (cid:16) ϕ (cid:3) − ϕ (cid:3) h x + ϕ (cid:3) − ϕ (cid:3) h x (cid:17)(cid:19) + h t h x L (cid:18) ϕ (cid:3) + ϕ (cid:3) , ϕ (cid:3) − ϕ (cid:3) h t , (cid:16) ϕ (cid:3) − ϕ (cid:3) h x + ϕ (cid:3) − ϕ (cid:3) h x (cid:17)(cid:19) . (124)This leads to a simplified implicit scheme with the same conservation properties as the midpointscheme derived from (122) but without the aforementioned instability.We want to define some shorthand notation similar to (53-55), u → (cid:16) u (cid:3) + u (cid:3) + u (cid:3) + u (cid:3) (cid:17) ≡ u ( (cid:3) ) (125) ∂u∂t → (cid:18) u (cid:3) − u (cid:3) h t + u (cid:3) − u (cid:3) h t (cid:19) ≡ u t ( (cid:3) ) , (126) ∂u∂x → (cid:18) u (cid:3) − u (cid:3) h x + u (cid:3) − u (cid:3) h x (cid:19) ≡ u x ( (cid:3) ) . (127)For the product of two fields as well as the product of a field and a derivative, according to themixed use of trapezoidal and midpoint rule, we obtain vu → (cid:20) v (cid:3) + v (cid:3) u (cid:3) + u (cid:3) v (cid:3) + v (cid:3) u (cid:3) + u (cid:3) (cid:21) ≡ vu ( (cid:3) ) , (128) v ∂u∂t → (cid:20) v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h t + v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h t (cid:21) ≡ vu t ( (cid:3) ) , (129) v ∂u∂x → (cid:20) v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h x + v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h x (cid:21) ≡ vu x ( (cid:3) ) . (130)This will prove handy in writing the discrete Lagrangians. The initial value problem for the linear advection equation is an instructive example for thederivation of numerical schemes. Its analytic solution is known and its structure is similar tothe systems we will consider in subsequent publications [37, 36, 38]. It is given by u t + cu x = 0 , (131)with initial condition u (0 , x ) = u ( x ) , (132) nd describes the advection of a field u ( t, x ) with a constant velocity c that is given as aparameter. The analytic solution is u ( t, x ) = u ( x − ct ) . (133)Here, x denotes the spatial variable only, while the base space X comprises both space andtime, i.e., ( t, x ) ∈ X . The field u can be interpreted as a density, such that R cu dx is the totalmomentum and R c u dx is the total energy of the system. Both are proportional to the “mass” R u dx which is a conserved quantity. In this case, we have F [ u ] = u t + cu x , (134)and the Fréchet derivative (81), D F [ u ] v = v t + cv x , (135)is not self-adjoint, as under the assumption of appropriate boundary conditions (e.g., periodicor homogeneous Dirichlet), integration by parts gives h D F [ u ] v, w i = − h D F [ u ] w, v i . (136)It follows that there is no functional A [ u ] such that D A = F , i.e., F cannot be written as avariational problem. The formal Lagrangian for the advection equation is obtained by multiplying (131) with theauxiliary variable v ( t, x ). The solution vector of the extended system is denoted ( u, v ) with thecorresponding section of the configuration bundle ˜ Y denoted by ˜ ϕ . In coordinates, ˜ ϕ is givenexplicitly as ˜ ϕ : ( t, x ) ( t, x, u, v ) , (137)and its jet prolongation as j ˜ ϕ : ( t, x ) ( t, x, u, v, u t , u x , v t , v x ) . (138)With that, the Lagrangian can be written as L = ˜y (cid:0) z t + c z x (cid:1) , (139)with action functional A [ ˜ ϕ ] = Z L (j ˜ ϕ ) dt dx = Z v (cid:0) u t + cu x (cid:1) dt dx. (140)The variation of the action with respect to the adjoint variable v gives the advection equation, δ A [ ˜ ϕ ] δv = u t + cu x = 0 , (141)while the variation with respect to the original variable u yields the adjoint equation, δ A [ ˜ ϕ ] δu = − v t − cv x = 0 , (142)always assuming vanishing variations at the boundaries. It is immediately observed that theadjoint equation has the same solution as the original equation, such that if u is a solution of theadvection equation, then ( u, u ) solves the Euler-Lagrange equations of the extended Lagrangian(139). This means, that the advection equation is actually self-adjoint in the sense of Ibragimov,cf. equation (102), with λ = − .2.2 Continuous Conservation Laws We will consider the conservation of mass and the L norm of u . The general form of the Noethercharge (101) is˜ J = Z (cid:20) ∂L∂ z at (cid:0) j ˜ ϕ (cid:1) · η a ( ˜ ϕ ) + ∂L∂ ˜z at (cid:0) j ˜ ϕ (cid:1) · ˜ η a ( ˜ ϕ ) (cid:21) d n x = Z ˜ J t (j ˜ ϕ ) d n x = Z vη d n x , (143)with ˜ J = const. for all t . Since the advection equation is self-adjoint, we can use the identityas embedding map, i.e., φ ( u ) = u , such that v = u , and the restriction of conservation laws forthe extended system to the original system is straight forward.The infinitesimal generator of the transformation that leads to mass conservation and itsfirst jet prolongation are V = ∂∂ y , j V = ∂∂ y . (144)Applying j V to (134), we obtain j V ( F ) = 0 , (145)such that by (95) and (97) we find ˜ η = 0 and the extended vector field (94) is˜ V = η ∂∂ y + ˜ η ∂∂ ˜y = ∂∂ y , (146)and its prolongation is j ˜ V = ∂∂ y . (147)The invariance of the Lagrangian (139) is easily confirmed, asj ˜ V ( L ) = ∂L∂ y = 0 . (148)The resulting conservation law (143) reads˜ J = Z ˜ J t (j ˜ ϕ ) dx = Z v dx = const. for all t. (149)Upon restricting the Noether current ˜ J with v = φ ( u ) = u , this becomes conservation of thetotal mass in the system, J = Z (cid:0) ˜ J t ◦ j Φ (cid:1) (j ϕ ) dx = Z u dx = const. for all t. (150)The conservation of momentum and energy follows exactly in the same way for η = c and η = c , respectively.The infinitesimal generator corresponding to conservation of the L norm is given by V = y ∂∂ y . (151)Applying its prolongation, j V = y ∂∂ y + z t ∂∂ z t + z x ∂∂ z x , (152) o (131), we obtain j V ( F ) = z t + c z x = λ F , (153)with λ = 1, such that by (95) and (97) the extended vector field (94) becomes˜ V = y ∂∂ y − ˜y ∂∂ ˜y , (154)and its prolongation isj ˜ V = y ∂∂ y + z t ∂∂ y t + z x ∂∂ y x − ˜y ∂∂ ˜y − ˜z t ∂∂ ˜y t − ˜z x ∂∂ ˜y x . (155)The invariance condition of the Lagrangian (139),j ˜ V ( L )(j ˜ ϕ ) = − v (cid:0) u t + cu x (cid:1) + v (cid:0) u t + cu x (cid:1) = 0 , (156)is again fulfilled. The corresponding Noether charge (143) is˜ J = Z ˜ J t (j ˜ ϕ ) dx = Z vu dx. (157)Upon restricting the Noether current ˜ J with v = φ ( u ) = u , this becomes conservation of the L norm, J = Z (cid:0) ˜ J t ◦ j Φ (cid:1) (j ϕ ) dx = Z u dx = const. for all t. (158)In the next sections we will derive the same conservation laws on the discrete level. We discretise the Lagrangian (139) as described in section 2.5 according to the midpoint rule(122), obtaining L mp d (cid:0) j ˜ ϕ ( (cid:3) ) (cid:1) = h t h x (cid:16) v (cid:3) + v (cid:3) + v (cid:3) + v (cid:3) (cid:17) ×× (cid:20) (cid:18) u (cid:3) − u (cid:3) h t + u (cid:3) − u (cid:3) h t (cid:19) + c (cid:18) u (cid:3) − u (cid:3) h x + u (cid:3) − u (cid:3) h x (cid:19)(cid:21) . (159)In the shorthand notation defined in (125-127), this becomes L d (cid:0) j ˜ ϕ ( (cid:3) ) (cid:1) = h t h x v ( (cid:3) ) (cid:2) u t ( (cid:3) ) + c u x ( (cid:3) ) (cid:3) , (160)which resembles the continuous Lagrangian. Variational Integrator
The discrete Euler-Lagrange field equations (59), corresponding to the variation of v d , are com-puted as 0 = 14 (cid:20) u i +1 ,j − − u i − ,j − h t + 2 u i +1 ,j − u i − ,j h t + u i +1 ,j +1 − u i − ,j +1 h t (cid:21) + c (cid:20) u i − ,j +1 − u i − ,j − h x + 2 u i,j +1 − u i,j − h x + u i +1 ,j +1 − u i +1 ,j − h x (cid:21) . (161) s in the continuous case, the discrete adjoint equation has the exact same form as the discreteadvection equation. We see that the derivatives are approximated by second-order centred finitedifferences. This means that we need initial data at two consecutive points in time, even thoughthe advection equation is first order in time, so that initial data at one point in time shouldsuffice. This problem is typical for applying the Veselov discretisation to formal Lagrangians.A simple solution will be provided in section 4.3.4.We also see that the time derivative is averaged in space and the spatial derivative is averagedin time. Under certain conditions, the spatial average of the time derivative can lead to grid-scale oscillations. This can be seen as follows. The time derivative in (161) features a spatialaverage of the form [1 2 1]. So if on top of the actual solution there is some oscillation, e.g.,of the form [ − − h x [ − x . This is very similar to the phenomenon of checker-boarding often observedin simulations of the incompressible Euler equations.In the next section we will prove some important discrete conservation laws of the resultingscheme. This exemplifies that even though a variational integrator has favourable conservationproperties it might be unstable, and that not every discretisation of the Lagrangian leads to agood scheme. Discrete Conservation Laws
For a formal Lagrangian of a single equation, the discrete symmetry condition (65) reads0 = ∂L d ∂ y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · η i,j + ∂L d ∂ ˜y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · ˜ η i,j + ∂L d ∂ y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · η i,j +1 + ∂L d ∂ ˜y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · ˜ η i,j +1 + ∂L d ∂ y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · η i +1 ,j +1 + ∂L d ∂ ˜y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · ˜ η i +1 ,j +1 + ∂L d ∂ y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · η i +1 ,j + ∂L d ∂ ˜y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · ˜ η i +1 ,j . (162)Here and in the following, evaluation of η i,j and ˜ η i,j at the fields ϕ i,j and ˜ ϕ i,j is implied. Forthe Lagrangian (159), this becomes0 = h t h x (cid:20) v i,j + v i,j +1 + v i +1 ,j +1 + v i +1 ,j (cid:21)(cid:20) η i +1 ,j − η i,j h t + η i +1 ,j +1 − η i,j +1 h t (cid:21) + h t h x c (cid:20) v i,j + v i,j +1 + v i +1 ,j +1 + v i +1 ,j (cid:21)(cid:20) η i,j +1 − η i,j h x + η i +1 ,j +1 − η i +1 ,j h x (cid:21) + h t h x (cid:20) ˜ η i,j + ˜ η i,j +1 + ˜ η i +1 ,j +1 + ˜ η i +1 ,j (cid:21)(cid:20) u i +1 ,j − u i,j h t + u i +1 ,j +1 − u i,j +1 h t (cid:21) + h t h x c (cid:20) ˜ η i,j + ˜ η i,j +1 + ˜ η i +1 ,j +1 + ˜ η i +1 ,j (cid:21)(cid:20) u i,j +1 − u i,j h x + u i +1 ,j +1 − u i +1 ,j h x (cid:21) . (163)The corresponding discrete conservation law (72) is˜ J d = N − X j =1 (cid:20) ∂L d ∂ y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · η i +1 ,j +1 + ∂L d ∂ y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · η i +1 ,j ∂L d ∂ ˜y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · ˜ η i +1 ,j +1 + ∂L d ∂ ˜y (cid:3) (cid:16) ˜ ϕ i,j , ˜ ϕ i,j +1 , ˜ ϕ i +1 ,j +1 , ˜ ϕ i +1 ,j (cid:17) · ˜ η i +1 ,j (cid:21) = const. for all i. (164)In contrast to the continuous case (143), here we also find contributions due to derivatives ofthe discrete Lagrangian L d with respect to ˜y. For the Lagrangian (159), the Noether chargebecomes˜ J d = h t h x N − X j =1 (cid:20) v i,j + v i,j +1 + v i +1 ,j +1 + v i +1 ,j (cid:21)(cid:20) η i +1 ,j +1 + η i +1 ,j h t + c η i +1 ,j +1 − η i +1 ,j h x (cid:21) + h t h x N − X j =1 (cid:20) ˜ η i +1 ,j +1 + ˜ η i +1 ,j (cid:21)(cid:20) (cid:18) u i +1 ,j − u i,j h t + u i +1 ,j +1 − u i,j +1 h t (cid:19) + c (cid:18) u i,j +1 − u i,j h x + u i +1 ,j +1 − u i +1 ,j h x (cid:19)(cid:21) . (165)We begin with mass conservation. The discrete generator (c.f. (146)) is η i,j = 1 , ˜ η i,j = 0 . (166)The discrete symmetry condition (163) is fulfilled for (166). The second summation in (165) iszero, hence ˜ J d = h x N − X j =1 h v i,j + v i,j +1 + v i +1 ,j +1 + v i +1 ,j i = const. for all i. (167)Identifying v i,j = u i,j , we obtain the discrete conservation law for the total mass in the system J d = h x N − X j =1 h u i,j + u i,j +1 + u i +1 ,j +1 + u i +1 ,j i = const. for all i. (168)This expression supports grid oscillations. For the L norm, the discrete generator (c.f. (154))is η i,j = u i,j , ˜ η i,j = − v i,j . (169)The discrete symmetry condition (163) is again easily checked to be fulfilled. The correspondingconservation law (72) is ˜ J d = const. for all i with˜ J d = h t h x N − X j =1 (cid:20) v i,j + v i,j +1 (cid:21)(cid:20) u i +1 ,j +1 + u i +1 ,j h t + c u i +1 ,j +1 − u i +1 ,j h x (cid:21) + h t h x N − X j =1 (cid:20) v i +1 ,j +1 + v i +1 ,j (cid:21)(cid:20) u i,j + u i,j +1 h t − c u i,j +1 − u i,j h x (cid:21) . (170)Upon identifying v i,j = u i,j and assuming periodic boundary conditions, the discrete Noethercharge becomes J d = h x N − X j =1 u i,j (cid:20) u i +1 ,j − + 2 u i +1 ,j + u i +1 ,j +1 ch t u i +1 ,j +1 − u i +1 ,j − h x (cid:21) . (171)Just as the discrete mass (168), this expression also supports grid oscillations. It is somewhatunexpected that the discrete conservation law contains terms corresponding to the spatial com-ponent J x of the Noether current. We therefore found an example where the discrete conservedquantity differs from the continuous conserved quantity. Nevertheless, J d is consistent with the L norm of u , as in the limit of h t →
0, holding c fixed, the additional term vanishes. .2.4 Trapezoidal Discretisation Writing out the discrete Lagrangian (123) for the linear advection equation, we obtain L tr d (cid:0) j ˜ ϕ ( (cid:3) ) (cid:1) = h t h x (cid:20) (cid:18) v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h t + v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h t (cid:19) + c (cid:18) v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h x + v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h x (cid:19)(cid:21) , or in the compact notation of (129-130), L tr d (cid:0) j ˜ ϕ ( (cid:3) ) (cid:1) = h t h x (cid:2) vu t ( (cid:3) ) + c vu x ( (cid:3) ) (cid:3) . (172) Variational Integrator
The discrete Euler-Lagrange equation (59) for the physical variable u d reads0 = u i +1 ,j − u i − ,j h t + c u i,j +1 − u i,j − h x , (173)which is the leapfrog scheme. The discrete Euler-Lagrange equation for the adjoint variable v d takes exactly the same form. Discrete Conservation Laws
The discrete symmetry condition (162) applied to (172) becomes0 = h t h x (cid:20) (cid:18) v i,j + v i +1 ,j η i +1 ,j − η i,j h t + v i,j +1 + v i +1 ,j +1 η i +1 ,j +1 − η i,j +1 h t (cid:19) + c (cid:18) v i,j + v i,j +1 η i,j +1 − η i,j h x + v i +1 ,j +1 + v i +1 ,j η i +1 ,j +1 − η i +1 ,j h x (cid:19) + 12 (cid:18) ˜ η i,j + ˜ η i +1 ,j u i +1 ,j − u i,j h t + ˜ η i,j +1 + ˜ η i +1 ,j +1 u i +1 ,j +1 − u i,j +1 h t (cid:19) + c (cid:18) ˜ η i,j + ˜ η i,j +1 u i,j +1 − u i,j h x + ˜ η i +1 ,j +1 + ˜ η i +1 ,j u i +1 ,j +1 − u i +1 ,j h x (cid:19)(cid:21) . (174)The corresponding discrete Noether charge (164) is˜ J d = h t h x N − X j =1 (cid:20)(cid:18) v i,j +1 + v i +1 ,j +1 h t + c v i +1 ,j + v i +1 ,j +1 h x (cid:19) · η i +1 ,j +1 + (cid:18) v i,j + v i +1 ,j h t − c v i +1 ,j + v i +1 ,j +1 h x (cid:19) · η i +1 ,j + (cid:18) u i +1 ,j +1 − u i,j +1 h t + c u i +1 ,j +1 − u i +1 ,j h x (cid:19) · ˜ η i +1 ,j +1 + (cid:18) u i +1 ,j − u i,j h t + c u i +1 ,j +1 − u i +1 ,j h x (cid:19) · ˜ η i +1 ,j (cid:21) . (175)For the discrete generator of mass conservation (166), the symmetry condition (174) is fulfilledand the discrete Noether charge (175) amounts to (168). Similarly, the symmetry condition(174) is fulfilled for the discrete generator (169). The restriction of the corresponding discreteNoether charge (175) becomes J d = h x N − X j =1 (cid:20) u i,j +1 u i +1 ,j +1 + u i,j u i +1 ,j (cid:21) = const. for all i, (176)which is indeed different from (171) and an immediate discretisation of the L norm. .2.5 Simplified Implicit Scheme The discrete Lagrangian (124) reads L mt d (cid:0) j ˜ ϕ ( (cid:3) ) (cid:1) = h t h x (cid:20) (cid:18) v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h t + v (cid:3) + v (cid:3) u (cid:3) − u (cid:3) h t (cid:19) + c (cid:18) v (cid:3) + v (cid:3) + v (cid:3) + v (cid:3) (cid:19)(cid:18) u (cid:3) − u (cid:3) h x + u (cid:3) − u (cid:3) h x (cid:19)(cid:21) , or equivalently L mt d (cid:0) j ˜ ϕ ( (cid:3) ) (cid:1) = h t h x (cid:2) vu t ( (cid:3) ) + c v ( (cid:3) ) u x ( (cid:3) ) (cid:3) . (177) Variational Integrator
The discrete Euler-Lagrange equations (59) for the physical variable are obtained in the form0 = u i +1 ,j − u i − ,j h t + c (cid:20) u i − ,j +1 − u i − ,j − h x + 2 u i,j +1 − u i,j − h x + u i +1 ,j +1 − u i +1 ,j − h x (cid:21) . (178)As becomes already apparent by comparing the discrete Lagrangian (177) with (172) and (159),the time derivative is discretised in the same way as in the trapezoidal scheme (173) and thespatial derivative is discretised in the same way as in the midpoint scheme (161). Again, thediscrete Euler-Lagrange equation for the adjoint variable takes the same form so that the systemis self-adjoint. Discrete Conservation Laws
In this case, the discrete symmetry condition (162) can be written as0 = h t h x (cid:20) (cid:18) v i,j + v i +1 ,j η i +1 ,j − η i,j h t + v i,j +1 + v i +1 ,j +1 η i +1 ,j +1 − η i,j +1 h t (cid:19) + c (cid:18) v i,j + v i,j +1 + v i +1 ,j +1 + v i +1 ,j (cid:19)(cid:18) η i,j +1 − η i,j h x + η i +1 ,j +1 − η i +1 ,j h x (cid:19) + 12 (cid:18) ˜ η i,j + ˜ η i +1 ,j u i +1 ,j − u i,j h t + ˜ η i,j +1 + ˜ η i +1 ,j +1 u i +1 ,j +1 − u i,j +1 h t (cid:19) + c (cid:18) ˜ η i,j + ˜ η i,j +1 + ˜ η i +1 ,j +1 + ˜ η i +1 ,j (cid:19)(cid:18) u i,j +1 − u i,j h x + u i +1 ,j +1 − u i +1 ,j h x (cid:19)(cid:21) , (179)which yields the discrete Noether charge, (c.f., (164)),˜ J d = h t h x N − X j =1 (cid:20)(cid:18) v i,j +1 + v i +1 ,j +1 h t + c v i,j + v i,j +1 + v i +1 ,j + v i +1 ,j +1 h x (cid:19) · η i +1 ,j +1 + (cid:18) v i,j + v i +1 ,j h t − c v i,j + v i,j +1 + v i +1 ,j + v i +1 ,j +1 h x (cid:19) · η i +1 ,j + (cid:18) u i +1 ,j +1 − u i,j +1 h t + c u i,j +1 − u i,j h x + c u i +1 ,j +1 − u i +1 ,j h x (cid:19) · ˜ η i +1 ,j +1 + (cid:18) u i +1 ,j − u i,j h t + c u i,j +1 − u i,j h x + c u i +1 ,j +1 − u i +1 ,j h x (cid:19) · ˜ η i +1 ,j (cid:21) . (180)The symmetry condition (179) is satisfied for both, mass conservation and the L norm. Thecorresponding discrete conservation laws reduce to (168) and (176), respectively. c/c grid = 0 .
75 (left) and c/c grid = 1 . In this section we want to carry out a simple von Neumann analysis of some of the variationalintegrators we have derived in order to analyse their dissipation and dispersion behaviour.The oscillatory function u ( t, x ) = e − ıωt + ıkx is a solution of the advection equation (131) ifand only if ( ω, k ) satisfies the exact dispersion relation ω = ck . When collocated at grid points t i = ih t , x j = jh x , this amounts to u ij = e − ı ( τi − ξj ) , (181)with τ = ωh t and ξ = kh x . In analogy with the continuous case, the collocated plane wave u ij isa solution of the numerical scheme if and only if ( τ, ξ ) satisfies the numerical dispersion relation.For the case of the Veselov scheme (161), the dispersion relation can be written assin( τ ) (cid:2) ξ ) (cid:3) − cc grid (cid:2) τ ) (cid:3) sin( ξ ) = 0 , (182)where c grid = h x /h t . This is the same dispersion relation as the one obtained by Ascher andMcLachlan [7] for the box scheme. For comparison, we consider the simplified implicit scheme(178), for which the dispersion relation issin( τ ) − cc grid (cid:2) τ ) (cid:3) sin( ξ ) = 0 , (183)as well as the leap frog scheme (173), for which the dispersion relation issin( τ ) − cc grid sin( ξ ) = 0 . (184)One can observe that in all cases the left-hand side is real valued, i.e., numerical dissipationis exactly zero, but there is numerical dispersion. Figure 5 shows the solution of the numericaldispersion relation in comparison with the exact linear dispersion of the advection equation n the ( ξ, τ )-plane for different values of the ratio c/c grid . As usual, for ( ξ, τ ) near the origin,the numerical dispersion is a good approximation of the exact linear relation, and thus thecorresponding harmonics have a phase velocity close to the exact value c . On the other hand,when the wave number is larger than half of the grid wave number 1 /h x , the frequency of theharmonic deviates from the exact dispersion relation. The corresponding harmonics have aphase velocity different from the exact value c , leading to numerical dispersion.For the simplified implicit scheme, one can also observe local maxima of the numericaldispersion, which correspond to a zero numerical group velocity, while the exact value for theadvection equation is equal to the phase velocity c .The left-hand-side panel of Figure 5 shows the dispersion relation for the case c/c grid < ξ ∈ [ − π, + π ] thereare two values of τ . The high-frequency branches are called parasitic modes. The Veselov andthe simplified implicit scheme also feature such modes at the τ = − π and τ = + π lines. As wewill see in the numerical experiments, these branches will be populated only if the spectrum ofthe solution is sufficiently broad to cover the nonlinear region of the dispersion relation, that isfor waves which are not well represented by the chosen grid parameters.For didactical purposes, we have also considered the case c/c grid > ξ = ± sin − (0 . c/c grid > We consider two test cases, a sum of cosines and a Gaussian. The first is purely meant toverify the dispersion relation, whereas the second is used to judge the quality of the solution. Inboth cases, the domain is Ω = [ − . , +0 .
5) with periodic boundary conditions, the timestep is h t = 2 . × − , the phase velocity is set to c = 1, and the number of points in space and timeis n x = 255 and n t = 4000, respectively. Hence, the time interval is [0 , Experimental Dispersion Relation
In this example, we initialise u with a sum of cosines, u C, ( x ) = n x / X i =1 cos( i πx ) , (185)which excites all the modes supported by the grid in order to experimentally verify the dispersionrelations (182-184) from the previous section. Figure 6 shows the theoretical and experimentaldispersion relations for the Veselov scheme, the simplified implicit scheme and the leapfrogscheme. The theoretical dispersion relations are well matched by the simulation, including someof the parasitic modes. Gaussian Wave
In this section, we consider a Gaussian, u G, ( x ) = 1 σ √ π exp (cid:18) − (cid:20) xσ (cid:21) (cid:19) , (186) a) Veselov Scheme (b) Simplified Implicit Scheme (c) Leapfrog Scheme Figure 6: Theoretical and experimental dispersion relation for the sum of cosines. Theexperimental values lie right on top of the theoretical curve. Red dashed line: Veselovscheme. Blue dash-dotted line: Simplified implicit scheme. Green dotted line: Leapfrogscheme. Black solid line: Analytical dispersion relation. with σ = 0 .
1, so that the spectrum is well within the linear region of the dispersion relation.Figure 7 shows that after 10 passes the form of the Gaussian is close to the initial one for both,the simplified implicit variational integrator (178) and the leapfrog scheme (173), and seems toexactly match the initial conditions for the Veselov scheme (161). The leapfrog scheme showsslightly less distortion than the simplified implicit scheme, which is explained by their dispersionrelations, which is closer to the analytic dispersion relation for the leapfrog scheme and evenmore so for the Veselov scheme. For all three integrators, conservation of energy and the L norm is also shown in Figure 7. Next we consider an example with more than one equation and more than one spatial dimension,the vorticity equation, ω t + { ψ, ω } = 0 , ω = ∆ ψ. (187)Here, ψ is the streaming function and ω is the vorticity. The vorticity equation arises bycomputing the curl of the incompressible Euler equation in two dimensions, u t + ( u · ∇ ) u + ∇ p = 0 , (188)where u is the fluid velocity and p the pressure. Specifically, ψ and ω are related to the velocity u by u = ∇ ⊥ ψ and ω = ∇ ⊥ · u with ∇ ⊥ = ( − ∂ y , ∂ x ) T . The Poisson bracket {· , ·} is defined by { ψ, ω } = ψ x ω y − ψ y ω x , (189)which coincides with the determinant of the Jacobian of the transformation ( x, y ) ( ψ, ω ). Forthis reason, the 2-dimensional Poisson bracket is also referred to as Jacobian. The variationaldiscretisation of (187) leads to a very interesting result, namely Arakawa’s famous discretisationof the Poisson bracket [5]. The formal Lagrangian of the vorticity equation (187) is L (j ˜ ϕ ) = ζ ( ω t + { ψ, ω } ) + χ ( ω − ∆ ψ ) , (190) .5 0.0 0.5 x u ( x ) Distribution Function u ( x,t =10 . t ( E ( t ) − E ) / E
1e 14
Error in the Total Energy ∆ E t ( L ( t ) − L ) / L
1e 14
Error in the L Norm ∆ L (a) Veselov Scheme x u ( x ) Distribution Function u ( x,t =10 . t ( E ( t ) − E ) / E
1e 13
Error in the Total Energy ∆ E t ( L ( t ) − L ) / L
1e 13
Error in the L Norm ∆ L (b) Simplified Implicit Scheme x u ( x ) Distribution Function u ( x,t =10 . t ( E ( t ) − E ) / E
1e 15
Error in the Total Energy ∆ E t ( L ( t ) − L ) / L
1e 15
Error in the L Norm ∆ L (c) Leapfrog Scheme Figure 7: Gaussian after ten passes with the Veselov scheme, the simplified implicitscheme and the leapfrog scheme. Comparing the solution (green dashed curve) to theinitial condition (blue curve) the effect of dispersion becomes visible. Energy and the L norm are well preserved. For the Veselov scheme as well as the simplified implicit schemethere is some drift in the errors, but this is extremely small.35 ith solution vector ( ω, ψ, ζ, χ ). Here, the formal Lagrangian L is defined on the second jet-bundle J ˜ Y due to the Laplace operator. In order to avoid second order derivatives, we apply asymmetrisation (corresponding to an integration by parts in the action functional) in the termthat produces the Laplacian. This yields an equivalent Lagrangian, which however is defined onthe first jet-bundle J ˜ Y .Some care has to be taken when discretising the Poisson bracket (see Salmon and Talley[57]). In order to retain the antisymmetry property of the continuous bracket at the discretelevel, a symmetrisation has to be introduced in the Lagrangian. Using integration by partswhile assuming appropriate boundary conditions, it is seen that the cyclic permutations of thefunctions in the integrand are all identical , Z ζ { ψ, ω } dx dy = Z ψ { ω, ζ } dx dy = Z ω { ζ, ψ } dx dy. (191)Instead of selecting one of those equivalent forms, a convex combination can be considered,namely, Z ζ { ψ, ω } dx dy = Z h α ζ { ψ, ω } + β ψ { ω, ζ } + γ ω { ζ, ψ } i dx dy, (192)with α + β + γ = 1. One finds that the symmetric case, α = β = γ = 1 /
3, is the one that retainsthe properties of the bracket at the discrete level. We therefore use the modified Lagrangian L (j ˜ ϕ ) = ζω t + 13 h ζ { ψ, ω } + ψ { ω, ζ } + ω { ζ, ψ } i + χω + ∇ χ · ∇ ψ. (193)Computing the variational derivatives of the action with respect to the adjoint variables, weobtain (187). The variational derivatives with respect to the physical variables give ζ t + { ψ, ζ } = χ, ∆ χ = { ω, ζ } . (194)This extended system of equations is obviously not self-adjoint, but we can still find a simplecompatible solution of the adjoint equations, namely ( ζ, χ ) = φ ( ω, ψ ) = ( ω, ζ = ω , then the Poisson bracket in the adjoint Poisson equation vanishes, which implies that χ is a harmonic function. In particular, we can choose χ = 0, so that the right-hand side of theequation for ζ becomes zero and ζ fulfils the same equation as ω . This justifies the choice of φ . In order to determine the conservation laws for the vorticity equation, it is not enough to justconsider vertical transformations as in (2.3). We have to consider also horizontal transforma-tions, as the interesting symmetries of the vorticity equation are not generated by purely verticalvector fields. To do so, the discrete Noether theorem as it is presented in section 2.6 is not suffi-ciently general. Unfortunately, it is not straight forward to consider symmetry generators withhorizontal components in the framework of finite difference discretisations of the Lagrangian aswe use it here. For that reason, we follow a simplified approach to determine the continuousconservation laws for the vorticity equation which is amenable to discretisation. A rigorousanalysis of the continuous and discrete conservation laws for the vorticity equation is postponedto a subsequent exposition.In the following, we will consider the vorticity equation for ω alone, F [ ω ] = ω t + { ψ, ω } = 0 , (195) This holds assuming boundary conditions such that the surface terms of the integration vanish, e.g.,periodic or homogeneous Dirichlet. ithout the Poisson equation, assuming a constant-in-time streaming function ψ . The Lagran-gian for this system reads ¯ L = ˜y (z t + ψ x z y − ψ y z x ) . (196)It is easily seen that Euler-Lagrange equations (10) for (196) amount to (195) and ζ t + { ψ, ζ } = 0 , (197)such that the extended system is self-adjoint and we can construct a solution for ζ simply by φ ( ω ) = ω . The generating vector field for vertical transformations of solutions of the vorticityequation (195) and its first jet prolongation are written as V = η ∂∂ y and j V = η ∂∂ y + η µ ∂∂ z µ . (198)The action of a general vector field V on F is given by j V ( F ) = η t + ψ x η y − ψ y η x . (199)This is used to determine the component ˜ η of the corresponding extended vector field ˜ V and itsfirst jet prolongation j ˜ V ,˜ V = η ∂∂ y + ˜ η ∂∂ ˜y and j ˜ V = η ∂∂ y + η µ ∂∂ z µ + ˜ η ∂∂ ˜y + ˜ η µ ∂∂ ˜z µ . (200)Its action on the Lagrangian (196) isj ˜ V ( ¯ L ) = ˜y ( η t + ψ x η y − ψ y η x ) + ˜ η (z t + ψ x z y − ψ y z x ) . (201)For conservation of circulation (“mass”) we have η = 1, and therefore j V = ∂∂ y . (202)The action of the prolonged vector field on (195) becomes j V ( F ) = 0 , (203)so that ˜ η = 0 and j ˜ V = ∂∂ y . (204)The action of the prolonged vector field on the Lagrangian (196) isj ˜ V ( ¯ L ) = ∂ ¯ L∂ y = 0 , (205)that is we have a symmetry, and the corresponding conservation law (Noether charge) reads˜ J = Z ˜ J t dx = Z ∂L∂ z t (cid:0) j ˜ ϕ (cid:1) · η ( ˜ ϕ ) dx = Z ζ dx = const. for all t. (206)Upon restricting the Noether current ˜ J with ζ = φ ( ω ) = ω , this becomes conservation of thetotal circulation in the system, J = Z (cid:0) ˜ J t ◦ j Φ (cid:1) (j ϕ ) dx = Z ω dx = const. for all t. (207) s here, ψ is treated as constant in time, energy conservation is obtained in the same way for η = ψ .Next we consider enstrophy, for which η = y, so that the first jet prolongation of thegenerating vector field, j V = y ∂∂ y + z t ∂∂ z t + z x ∂∂ z x + z y ∂∂ z y , (208)acts on the vorticity equation (196) as j V ( F ) = z t + ψ x z y − ψ y z x = F , (209)hence ˜ η = − ˜y and j ˜ V = y ∂∂ y + z t ∂∂ z t + z x ∂∂ z x + z y ∂∂ z y − ˜y ∂∂ ˜y − ˜z t ∂∂ ˜z t − ˜z x ∂∂ ˜z x − ˜z y ∂∂ ˜z y . (210)The prolongation of the extended vector field acts on the Lagrangian (196) asj ˜ V ( ¯ L ) = ˜y (z t + ψ x z y − ψ y z x ) − ˜y (z t + ψ x z y − ψ y z x ) = 0 , (211)so that we have a symmetry whose Noether charge is computed as˜ J = Z ˜ J t dx = Z ∂L∂ z t (cid:0) j ˜ ϕ (cid:1) · η ( ˜ ϕ ) dx = Z ζω dx = const. for all t. (212)Restricting the Noether current ˜ J with ζ = φ ( ω ) = ω , this becomes conservation of enstrophy, J = Z (cid:0) ˜ J t ◦ j Φ (cid:1) (j ϕ ) dx = Z ω dx = const. for all t. (213)In summary, we obtain conservation of(a) circulation (“mass”) Z ω dx dy, (214a)(b) enstrophy ( L norm) Z ω dx dy, (214b)(c) and kinetic energy 12 Z ψω dx dy. (214c)For the symmetrised Lagrangian, c.f. equation (193),¯ L (j ˜ ϕ ) = ζω t + 13 h ζ { ψ, ω } + ψ { ω, ζ } + ω { ζ, ψ } i , (215)we obtain the same conservation laws, although in that case they correspond to divergencesymmetries. t : y x : y y : xyt Figure 8: For a given spacetime grid cell, there are four possible ways of defining deriva-tives in the different coordinate directions ( t, x, y ), highlighted along the black lines.
To discretise the Lagrangian L ( ˜ ϕ, ˜ ϕ t , ˜ ϕ x , ˜ ϕ y ) from (193) we adopt the same strategy as in section4.2.5. We use simple finite differences to approximate the derivatives (see Figure 8), the midpointrule for the time integral and the trapezoidal rule for the two spatial integrals, L d (cid:0) j ˜ ϕ ( (cid:3) ) (cid:1) = h t h x h y L (cid:16) ˜ ϕ (cid:3) + ˜ ϕ (cid:3) , ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h t , (cid:16) ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h x + ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h x (cid:17) , (cid:16) ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h y + ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h y (cid:17)(cid:17) + h t h x h y L (cid:16) ˜ ϕ (cid:3) + ˜ ϕ (cid:3) , ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h t , (cid:16) ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h x + ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h x (cid:17) , (cid:16) ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h y + ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h y (cid:17)(cid:17) + h t h x h y L (cid:16) ˜ ϕ (cid:3) + ˜ ϕ (cid:3) , ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h t , (cid:16) ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h x + ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h x (cid:17) , (cid:16) ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h y + ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h y (cid:17)(cid:17) + h t h x h y L (cid:16) ˜ ϕ (cid:3) + ˜ ϕ (cid:3) , ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h t , (cid:16) ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h x + ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h x (cid:17) , (cid:16) ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h y + ˜ ϕ (cid:3) − ˜ ϕ (cid:3) h y (cid:17)(cid:17) . (216)In 2 + 1 dimensions, the discrete Euler-Lagrange field equations have eight contributions insteadof four,0 = ∂L d ∂ ˜ ϕ a (cid:3) (cid:16) ˜ ϕ i,j,k , ˜ ϕ i,j +1 ,k , ˜ ϕ i,j +1 ,k +1 , ˜ ϕ i,j,k +1 , ˜ ϕ i +1 ,j,k , ˜ ϕ i +1 ,j +1 ,k , ˜ ϕ i +1 ,j +1 ,k +1 , ˜ ϕ i +1 ,j,k +1 (cid:17) + ∂L d ∂ ˜ ϕ a (cid:3) (cid:16) ˜ ϕ i,j − ,k , ˜ ϕ i,j,k , ˜ ϕ i,j,k +1 , ˜ ϕ i,j − ,k +1 , ˜ ϕ i +1 ,j − ,k , ˜ ϕ i +1 ,j,k , ˜ ϕ i +1 ,j,k +1 , ˜ ϕ i +1 ,j − ,k +1 (cid:17) + ∂L d ∂ ˜ ϕ a (cid:3) (cid:16) ˜ ϕ i,j − ,k − , ˜ ϕ i,j,k − , ˜ ϕ i,j,k , ˜ ϕ i,j − ,k , ˜ ϕ i +1 ,j − ,k − , ˜ ϕ i +1 ,j,k − , ˜ ϕ i +1 ,j,k , ˜ ϕ i +1 ,j − ,k (cid:17) + ∂L d ∂ ˜ ϕ a (cid:3) (cid:16) ˜ ϕ i,j,k − , ˜ ϕ i,j +1 ,k − , ˜ ϕ i,j +1 ,k , ˜ ϕ i,j,k , ˜ ϕ i +1 ,j,k − , ˜ ϕ i +1 ,j +1 ,k − , ˜ ϕ i +1 ,j +1 ,k , ˜ ϕ i +1 ,j,k (cid:17) + ∂L d ∂ ˜ ϕ a (cid:3) (cid:16) ˜ ϕ i − ,j,k , ˜ ϕ i − ,j +1 ,k , ˜ ϕ i − ,j +1 ,k +1 , ˜ ϕ i − ,j,k +1 , ˜ ϕ i,j,k , ˜ ϕ i,j +1 ,k , ˜ ϕ i,j +1 ,k +1 , ˜ ϕ i,j,k +1 (cid:17) + ∂L d ∂ ˜ ϕ a (cid:3) (cid:16) ˜ ϕ i − ,j − ,k , ˜ ϕ i − ,j,k , ˜ ϕ i − ,j,k +1 , ˜ ϕ i − ,j − ,k +1 , ˜ ϕ i,j − ,k , ˜ ϕ i,j,k , ˜ ϕ i,j,k +1 , ˜ ϕ i,j − ,k +1 (cid:17) + ∂L d ∂ ˜ ϕ a (cid:3) (cid:16) ˜ ϕ i − ,j − ,k − , ˜ ϕ i − ,j,k − , ˜ ϕ i − ,j,k , ˜ ϕ i − ,j − ,k , ˜ ϕ i,j − ,k − , ˜ ϕ i,j,k − , ˜ ϕ i,j,k , ˜ ϕ i,j − ,k (cid:17) + ∂L d ∂ ˜ ϕ a (cid:3) (cid:16) ˜ ϕ i − ,j,k − , ˜ ϕ i − ,j +1 ,k − , ˜ ϕ i − ,j +1 ,k , ˜ ϕ i − ,j,k , ˜ ϕ i,j,k − , ˜ ϕ i,j +1 ,k − , ˜ ϕ i,j +1 ,k , ˜ ϕ i,j,k (cid:17) . (217)The resulting integrator for the vorticity equation takes the form ω i +1 − ω i − h t + 18 h A ( ψ i +1 , ω i +1 ) + A ( ψ i , ω i +1 ) + A ( ψ i +1 , ω i )+ 2 A ( ψ i , ω i ) + A ( ψ i − , ω i ) + A ( ψ i , ω i − ) + A ( ψ i − , ω i − ) i , (218) here ω i = ( ω ij ) j and ψ i = ( ψ ij ) j are the rows of the solution matrix corresponding to constant-time slices, whereas A denotes Arakawa’s discretisation of the Poisson bracket [5], given by A ( ψ, ω ) = 13 (cid:16) A ++ ( ψ, ω ) + A + × ( ψ, ω ) + A × + ( ψ, ω ) (cid:17) , (219)with ω and ψ being row vectors and A ++ ( ψ, ω ) = 14 h x h y h(cid:0) ψ +0 − ψ − (cid:1)(cid:0) ω − ω − (cid:1) − (cid:0) ψ − ψ − (cid:1)(cid:0) ω +0 − ω − (cid:1)i , (220a) A + × ( ψ, ω ) = 14 h x h y h ψ +0 (cid:0) ω + − − ω ++ (cid:1) − ψ − (cid:0) ω −− − ω − + (cid:1) − ψ (cid:0) ω − + − ω ++ (cid:1) + ψ − (cid:0) ω −− − ω + − (cid:1)i , (220b) A × + ( ψ, ω ) = 14 h x h y h ψ ++ (cid:0) ω +0 − ω (cid:1) − ψ −− (cid:0) ω − − ω − (cid:1) − ψ − + (cid:0) ω − − ω (cid:1) + ψ + − (cid:0) ω − − ω +0 (cid:1)i . (220c)The indices ( − , , +) indicate the increment of the corresponding index on the grid relative tothe point where the bracket is computed. The integrator for the Poisson equation is h ∆ x ψ i,j,k i t + h ∆ y ψ i,j,k i t = h ω i,j,k i t . (221)Here, ∆ x and ∆ y denote the standard centred finite difference second-order derivative withrespect to x and y , i.e., ∆ x ψ i,j,k = ψ i,j − ,k − ψ i,j,k + ψ i,j +1 ,k h x , ∆ y ψ i,j,k = ψ i,j,k − − ψ i,j,k + ψ i,j,k +1 h y , and the angle brackets h·i t denote an average in time of the form (cid:2) (cid:3) , namely, h ω i,j,k i t = ω i − ,j,k + 2 ω i,j,k + ω i +1 ,j,k . It is remarkable that Arakawa’s discretisation of the Poisson bracket is recovered. Indeed,a similar derivation was proposed by Salmon and Talley [57]. Our approach, however, is fullycovariant, leading to a complete spacetime discretisation, that is a combination of Arakawa’sscheme in space with a symplectic integrator in time. Whereas the Arakawa bracket aloneguarantees energy conservation only for the spatial discretisation, i.e., up to errors due to thediscretisation of the time derivative, the variational integrator is exactly energy preserving.
The integrator (218) for the vorticity equation has one drawback. Even though we are consid-ering a partial differential equation that is first order in time, we need to initialise the fieldsat two consecutive points in the discrete time domain, that is (218) is a multistep integrator.Nevertheless, we can reduce it to a single-step integrator. Observe that we can rewrite (218) inthe following way,0 = ψ i +1 − ψ i h t + 18 h A ( ψ i +1 , ω i +1 ) + A ( ψ i +1 , ω i ) + A ( ψ i , ω i +1 ) + A ( ψ i , ω i ) i + ψ i − ψ i − h t + 18 h A ( ψ i , ω i ) + A ( ψ i , ω i − ) + A ( ψ i − , ω i ) + A ( ψ i − , ω i − ) i , (222) hich is the average over two points in time of the single-step integrator0 = ψ i − ψ i − h t + 14 h A ( ψ i , ω i ) + A ( ψ i , ω i − ) + A ( ψ i − , ω i ) + A ( ψ i − , ω i − ) i . (223)Analogously, the Poisson equation (221) is the average over three points in time of∆ x ψ i,j,k + ∆ y ψ i,j,k = ω i,j,k . (224)The solution of the single-step integrator (223-224) is fully determined by specifying the vorticity ω = ω at t = t on the spatial grid. Then ψ can be obtained by (224) and the vorticity canbe advected by solving (223). We observe that if the sequence ( ω i , ψ i ) solves (223-224) withinitial conditions ω , then it is also a solution of (218-221) with initialisation ( ω , ψ , ω , ψ ).Vice versa, if the sequence ( ω i , ψ i ) is a solution of (223-224) initialised by using (223-224) withthe initial data ω , then it is also a solution of (223-224). Equations (223-224) are called theunderlying one-step method of (218-221).The simplified equations (223) and (224) can also be obtained directly from a discrete La-grangian as is shown in [37]. This is of course preferable as it allows us to check for symmetriesand to compute the discrete conserved quantities more easily. As solutions of the simplifiedintegrator (223-224) will also be solutions of (218-221), they will satisfy the same conservationlaws. Nevertheless, it is possible that (223-224) will have additional symmetries and thereforeadditional invariants. The discrete conservation laws of the variational integrator (218) for the vorticity equation arecomputed in the same way as for the linear advection equation (section 4.2) except that we haveto account for discrete divergence symmetries as introduced in section 2.6 and that we are usingthe generators discussed in section 4.3.2. As in the continuous case, we restrict our analysis tothe linear case, that is the Lagrangian (215) discretised according to (216).Because the vorticity equation is defined on three-dimensional spacetime, the actual compu-tations are tediously long and hardly possible to carry out by hand. We therefore provide onlythe results, obtained via computer aided calculations, namely(a) circulation (“mass”) h x h y X j,k ω j,k = const ., (225a)(b) enstrophy ( L norm) h x h y X j,k ω j,k = const ., (225b)(c) and kinetic energy h x h y X j,k ω j,k φ j,k = const .. (225c)It is well known that these quantities are preserved by Arakawa’s discrete bracket [5]. In nu-merical simulations we indeed find that our integrator conserves them to machine accuracy orat least the tolerance of the nonlinear solver, not only in the linear case (when ψ is constant intime), but even in the nonlinear case (with self-consistent streaming function ψ ). .3.6 Numerical Examples We implemented the simplified variational integrator (223-224) using Python [2, 39], Cython[13, 12], PETSc [9, 10] and petsc4py [18]. Visualisation was done using NumPy [61], SciPy[30] and matplotlib [1, 24]. The nonlinear system is solved via Picard iteration. Within eachnonlinear step the two linear systems corresponding to the vorticity equation and the Poissonequation are solved separately. The vorticity equation is solved with GMRES and the Poissonequation via LU decomposition with SuperLU [44, 43]. The tolerance of the nonlinear solver isset to 10 − or smaller, which is reached after 5 −
10 iterations. Usually, the GMRES solverneeds between 5 and 25 iterations to converge with a relative tolerance of 10 − or an absolutetolerance of 10 − . The Linear Case
At first we consider the linear case, where we prescribe the streaming function ψ and only solvethe vorticity equation. The streaming potential is set to ψ ( x, y ) = y + 1 − cos( x ) , (226)while the vorticity is initialised with a localised, symmetric Gaussian, ω ( x, y ) = f ( x, x , σ x ) f ( y, y , σ y ) , (227)with f ( z, z , σ ) = 1 σ √ π exp (cid:18) − (cid:20) z − z σ (cid:21) (cid:19) , z ∈ R . (228)The domain is Ω = [ − π, +2 π ) × [ − π, +2 π ) with periodic boundaries, the spatial resolution is1024 × − .The parameters are set to x = 0, y = 2 and σ x = σ y = 0 .
2, so that the Gaussian is placedon the separatrix of the streaming function (c.f., Figure 9). As the Gaussian is moving alongthe separatrix, it is stretched while the contours of the vorticity are preserved. In fact, if a field ω is advected by a smooth flow, the topology of the contours of ω should be preserved. Thisbehaviour appears to be respected by the integrator. In Figure 10a, the contours for differentvalues of the vorticity are shown. Until t = 3 all three contours stay intact. At about t = 6the stretching of the contours is so strong that resolution becomes insufficient and only theoutermost contour is still preserved. Towards the centre and the boundaries in x the effect ofdispersion becomes visible. Figure 10a shows the evolution of the errors of circulation, enstrophyand energy. Lamb Dipole
The lamb dipole [52] is a stationary solution of the vorticity equation which is at rest in itsframe of reference. The vorticity is initialised as ω ( x, y ) = λU cos θ J ( λr ) J ( λR ) r ≤ R, r > R, (229)which leads to a dipole of radius R moving in the y -direction with velocity U . Here, r and θ denote polar coordinates, i.e., r = q x + y and θ = arctan( y/x ) , (230) i is the i -th Bessel function of first kind and λR is the first root of J . In the example, weuse R = 0 . U = 1. The domain is Ω = [ − , +1) × [ − , +1) with periodic boundaries, thespatial resolution is 1024 × − . The tolerance of the nonlinear solveris set to 10 − .After 10 ,
000 timesteps (at t = 10) no distortion of the lamb dipole is visible (comparefigures 11(a) and 11(b)). The dipole keeps its shape just as it is supposed to. It does not end upexactly at the initial position due to the finite size of the box which reduces the velocity of thereference frame [52, section IV.B]. Figure 11(c) shows that circulation, enstrophy and energyare preserved to about machine accuracy. There is a small drift in the errors of enstrophy andenergy, but this is of O (10 − ) and therefore comparable to the tolerance of the nonlinear solver. Gaussian Vortex
Vorticity is initialised by a nonsymmetric Gaussian, ω ( x, y ) = f ( x, x , σ x ) f ( y, y , σ y ) , (231)with f ( x, x , σ ) = 1 σ √ π exp (cid:18) − (cid:20) x − x σ (cid:21) (cid:19) , (232)which develops into a spiral structure. The domain is Ω = [ − , +1) × [ − , +1) with periodicboundaries, the spatial resolution is 1024 × − . The tolerance of thenonlinear solver is 10 − . We set σ x = 0 . σ y = 0 . x = y = 0. The solution followsthe expected dynamics (c.f., Figure 12) and conservation of circulation, enstrophy and energyis respected to machine accuracy (c.f., Figure 14(a)). Vortex Sheet Rollup
The vortex sheet is a popular model in fluid dynamics used to approximate shear layers. Fol-lowing [19], we use initial conditions u = (cid:18) tanh( ρ ( y − . y ≤ . , tanh( ρ (0 . − y )) for y > . , , .
05 sin(2 πx ) (cid:19) T , (233)which corresponds to the initial vorticity ω ( x, y ) = ∇ ⊥ · u = 0 . π cos(2 πx ) + − ρ cosh ( ρ ( y − . y ≤ r, + ρ cosh ( ρ (0 . − y )) y > r. (234)The domain is Ω = [0 , × [0 ,
1) with periodic boundaries, discretised by a grid of 1024 × h t = 10 − . The tolerance of the nonlinear solver is 10 − and theparameter ρ = 30. We observe good qualitative agreement of the solution in Figure 13 with [19,p. 328] together with excellent conservation properties (c.f., Figure 14(b)). a) t = 0 . t = 0 . t = 1 . t = 2 . t = 3 . t = 6 . Figure 9: Linear case with a separatrix. Vorticity ω at t = 0 .
0, 0 .
5, 1 .
0, 2 .
0, 3 . . ψ are also shown.44 a) ( C ( t ) − C ) / C
1e 14
Total Circulation Error ∆ C ( t ) ( W ( t ) − W ) / W
1e 13
Total Enstrophy Error ∆ W ( t ) t ( E ( t ) − E ) / E
1e 14
Total Energy Error ∆ E ( t ) (b) Figure 10: Linear case with a separatrix. (a) Contours of the streaming function ψ (black)and of the vorticity ω (at 0 . ω , 0 . ω and 0 . ω ) at times t = 0 .
0, 0 . .
0, 2 .
0, 3 . .
0. (b) Conservation laws.45 a) t = 0 . t = 10 . C ( t )
1e 13
Total Circulation C ( t ) ( W ( t ) − W ) / W
1e 13
Total Enstrophy Error ∆ W ( t ) t ( E ( t ) − E ) / E
1e 13
Total Energy Error ∆ E ( t ) (c) Figure 11: Lamb Dipole. Top: Vorticity ω at t = 0 . .
0. Bottom: ConservationLaws. 46 a) t = 0 (b) t = 2(c) t = 4 (d) t = 6(e) t = 8 (f) t = 10 Figure 12: Gaussian Vortex. Vorticity ω at t = 0, 2, 4, 6, 8, 10.47 a) t = 0 .
00 (b) t = 1 . t = 1 .
25 (d) t = 1 . t = 1 .
75 (f) t = 2 . Figure 13: Vortex sheet rollup. Vorticity ω at t = 0 .
00, 1 .
00, 1 .
25, 1 .
50, 1 .
75, 2 . .80.60.40.20.00.20.40.60.81.0 ( C ( t ) − C ) / C
1e 13
Total Circulation Error ∆ C ( t ) ( W ( t ) − W ) / W
1e 14
Total Enstrophy Error ∆ W ( t ) t ( E ( t ) − E ) / E
1e 13
Total Energy Error ∆ E ( t ) (a) Gaussian Vortex C ( t )
1e 13 7.16957e 8
Total Circulation C ( t ) ( W ( t ) − W ) / W
1e 13
Total Enstrophy Error ∆ W ( t ) t ( E ( t ) − E ) / E
1e 13
Total Energy Error ∆ E ( t ) (b) Vortex Sheet Rollup Figure 14: Conservation Laws. The errors in the circulation, enstrophy and energy areall of the order of the machine accuracy. 49
Summary and Outlook
We created a link between variational integrators, formal Lagrangians and Ibragimov’s theory ofconservation laws for arbitrary differential equations. The proposed method allows us to derivevariational discretisation schemes for potentially arbitrary systems of differential equations, eventhe ones that do not possess a classical variational formulation. Thereby we were able to extendthe applicability of the variational integrator method to a class of systems much larger thanoriginally envisaged. The main strength of our method is that it allows for the straight forwarddesign of numerical schemes that respect certain conservation laws of the system at hand.We presented an introduction to the theory of variational integrators that tries to make thegeometric framework of the variational integrators available also to non specialists. Thereby wehope to make this class of methods accessible to a wider audience. We extended the discreteNoether theorem to include discrete divergence symmetries.The power of the method was demonstrated by several examples that are prototypical forproblems arising in fluid dynamics and plasma physics. We emphasised the analysis of discreteconservation laws, which is seldom found in the variational integrator literature, and verifiedthese theoretical properties in numerical experiments. In follow-up papers we will present nu-merical results for the Vlasov-Poisson system [37] as well as ideal and reduced magnetohydro-dynamics [36, 38]. There, explicit computations will demonstrate the favourable properties ofthe variational integrators for more elaborate and challenging applications.Remarkably, we recovered Arakawa’s discretisation of the Poisson bracket combined with asymplectic integrator in time. That is, we constructed a spacetime generalisation of Arakawa’sdiscretisation. With our method, it is also straight-forward to generalise Arakawa’s method tohigher spatial dimensions and to higher order schemes.An open question within our framework is the meaning of the discrete multisymplectic formarising from the boundary terms of the action sum (see [48] for further details) and its restrictionto the original system. We observed that the multisymplectic form vanishes identically if theextended system is self-adjoint in the sense of Ibragimov. This is a topic we have not includedin our discussion but one that certainly deserves attention.The limits of the finite-difference approach to the discretisation of the Lagrangian becameobvious in several places in this work. The treatment of discrete divergence symmetries is notstraight-forward so that we are limited to global conservation laws. We cannot treat horizontaltransformations in the symmetry generator but only vertical transformations. This is a problemfor more complicated systems like the vorticity equations or the nonlinear advection equation(inviscid Burgers equation). For nonlinear systems there is the additional complication thateven if the system is self-adjoint on the continuous level, it will in general not be self-adjointon the discrete level due to the absence of a discrete Leibniz rule. Moreover, we are limited tolow-order schemes. While it is in principle possible to design higher-order methods, this quicklybecomes very cumbersome and confusing, especially in the analysis of the discrete conservationlaws. All those issues appear easier to deal with in a Galerkin framework with finite elements orsplines as basis functions. The former has already been considered to a certain extend [16], thatis for the discretisation of the Lagrangian and the approximation of the action integral, but notfor the analysis of conservation laws. The latter is a topic under active development [35].So far, we only considered examples of advection-diffusion type. A detailed analysis of theapplicability of the proposed method to other types of equations would be most interesting.
Acknowledgements
The first author would like to thank Bruce D. Scott for his guidance and the freedom to pursuethe path that seemed most interesting, as well as Eric Sonnendrücker for supporting this project.We are thankful to Emanuele Tassi and Jonathan Squire for raising and discussing important uestions that added a lot to the understanding of the matters presented as well as readingvarious versions of the manuscript. Furthermore, we wish to thank Claudio Dappiaggi for thediscussions on the inverse problem of calculus of variations and for suggesting the paper byBampi and Morro [11]. References [1] The matplotlib Plotting Library. .[2] The Python Programming Language. .[3] Ralph Abraham and Jerrold E. Marsden.
Foundations of Mechanics . American Mathemat-ical Society, 1978. URL http://authors.library.caltech.edu/25029/ .[4] Ralph Abraham, Jerrold E. Marsden, and Tudor S. Ratiu.
Manifolds, Tensor Analysis andApplication . Springer, 1988.[5] Akio Arakawa. Computational design for long-term numerical integration of the equationsof fluid motion: Two-dimensional incompressible flow. part i.
Journal of ComputationalPhysics , 1:119–143, 1966. doi: 10.1016/0021-9991(66)90015-5.[6] Vladimir I. Arnold.
Mathematical Methods of Classical Mechanics . Springer, 1989.[7] Uri M. Ascher and Robert I. McLachlan. Multisymplectic box schemes and the korteweg–devries equation.
Applied Numerical Mathematics , 48:255–269, 2004. doi: http://dx.doi.org/10.1016/j.apnum.2003.09.002.[8] Robert W. Atherton and George M. Homsy. On the existence and formulation of variationalprinciples for nonlinear differential equations.
Studies in Applied Mathematics , 54:31–60,1975.[9] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschel-man, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, LoisCurfman McInnes, Karl Rupp, Barry F. Smith, and Hong Zhang. PETSc Web page. .[10] Satish Balay, Shrirang Abhyankar, Mark F. Adams, Jed Brown, Peter Brune, Kris Buschel-man, Victor Eijkhout, William D. Gropp, Dinesh Kaushik, Matthew G. Knepley, LoisCurfman McInnes, Karl Rupp, Barry F. Smith, and Hong Zhang. PETSc users manual.Technical Report ANL-95/11 - Revision 3.5, Argonne National Laboratory, 2014. URL .[11] Franco Bampi and Angelo Morro. The inverse problem of the calculus of variations appliedto continuum physics.
Journal of Mathematical Physics , 23:2312–2321, 1982. doi: 10.1063/1.525322.[12] Stefan Behnel, Robert Bradshaw, Craig Citro, Lisandro Dalcin, Dag Sverre Seljebotn, andKurt W. Smith. Cython: The best of both worlds.
Computing in Science Engineering , 13:31–39, 2011. doi: 10.1109/MCSE.2010.118.[13] Robert Bradshaw, Stefan Behnel, Dag Sverre Seljebotn, Greg Ewing, and et al. The CythonCompiler. .[14] Chris J. Budd and Matthew D. Piggott. Geometric integration and its applications. In inHandbook of numerical analysis , pages 35–139. North-Holland, 2000.
15] Cédric M. Campos.
Geometric Methods in Classical Field Theory and Continuous Media .PhD thesis, Universidad Autónoma de Madrid, 2010.[16] Jing-Bo Chen. Variational integrators and the finite element method.
Applied Mathematicsand Computation , 196:941–958, 2008.[17] Snorre H. Christiansen, Hans Z. Munthe-Kaas, and Brynjulf Owren. Topics instructure-preserving discretization.
Acta Numerica , 20:1–119, 2011. doi: 10.1017/S096249291100002X.[18] Lisandro Dalcin, Rodrigo R. Paz, Pablo A. Kler, and Alejandro Cosimo. Parallel distributedcomputing using python.
Advances in Water Resources , 34:1124–1139, 2011. doi: 10.1016/j.advwatres.2011.04.013.[19] Michel O. Deville, Paul F. Fischer, and Ernest H. Mund.
High-Order Methods for Incom-pressible Fluid Flow . Cambridge University Press, 2002.[20] Evan S. Gawlik, Patrick Mullen, Dmitry Pavlov, Jerrold E. Marsden, and Mathieu Des-brun. Geometric, variational discretization of continuum theories.
Physica D: NonlinearPhenomena , 240:1724–1760, 2011. doi: 10.1016/j.physd.2011.07.011. arXiv:1010.4851.[21] Mark J. Gotay, James Isenberg, and Jerrold E. Marsden. Momentum maps and classicalfields, 1998.[22] Ernst Hairer, Christian Lubich, and Gerhard Wanner.
Geometric Numerical Integration .Springer, 2006.[23] Darryl D. Holm, Tanya Schmah, and Cristina Stoica.
Geometric Mechanics and Symmetry .Oxford University Press, 2009.[24] John D. Hunter. Matplotlib: A 2d graphics environment.
Computing In Science & Engi-neering , 9:90–95, 2007.[25] Nail H. Ibragimov. Integrating factors, adjoint equations and lagrangians.
Journal ofMathematical Analysis and Applications , 318:742–757, 2006. doi: 10.1016/j.jmaa.2005.11.012.[26] Nail H. Ibragimov. A new conservation theorem.
Journal of Mathematical Analysis andApplications , 333:311–328, 2007. doi: 10.1016/j.jmaa.2006.10.078.[27] Nail H. Ibragimov. Quasi-self-adjoint differential equations.
Archives of ALGA , 4:55–60,2007.[28] Nail H. Ibragimov, Mariano Torrisi, and Rita Tracinà. Quasi self-adjoint nonlinear waveequations.
Journal of Physics A: Mathematical and Theoretical , 43:442001, 2010. doi:10.1088/1751-8113/43/44/442001.[29] Vladimir G. Ivancevic and Tijana T. Ivancevic.
Applied Differential Geometry . WorldScientific, 2007.[30] Eric Jones, Travis Oliphant, Pearu Peterson, and et al. SciPy: Open source scientific toolsfor Python. .[31] Jorge V. Jose and Eugene J. Saletan.
Classical Dynamics: A Contemporary Approach .Cambridge University Press, 1998.
32] Ivan Kolar, Jan Slovak, and Peter W. Michor.
Natural operations in differential geometry .Springer, 1993. URL .[33] Yvette Kosmann-Schwarzbach.
The Noether Theorems: Invariance and Conservation Lawsin the Twentieth Century . Springer, 2010.[34] Shinar O. Kouranbaeva and Steve Shkoller. A Variational Approach to Second-OrderMultisymplectic Field Theory.
Journal of Geometry and Physics , 25:333–366, 2000. doi:10.1016/S0393-0440(00)00012-7. arXiv:math/9909100.[35] Michael Kraus. Spline Variational Integrators. In preparation, 2015.[36] Michael Kraus and Omar Maj. Variational Integrators for Ideal Magnetohydrodynamics.In preparation, 2015.[37] Michael Kraus, Omar Maj, and Eric Sonnendrücker. Variational Integrators for the Vlasov-Poisson System. In preparation, 2015.[38] Michael Kraus, Emanuele Tassi, and Daniela Grasso. Variational Integrators for ReducedMagnetohydrodynamics. In preparation, 2015.[39] Hans Petter Langtangen.
A Primer on Scientific Programming with Python . Springer, 4thedition, 2014.[40] A. Lew, J. E. Marsden, M. Ortiz, and M. West. Variational time integrators.
InternationalJournal for Numerical Methods in Engineering , 60(1):153–212, 2004. doi: 10.1002/nme.958.[41] Adrian Lew, Jerrold E. Marsden, Michael Ortiz, and Matthew West. Asynchronous Vari-ational Integrators.
Archive for Rational Mechanics and Analysis , 167:85–146, 2003. doi:10.1007/s00205-002-0212-y.[42] Adrian Lew, Jerrold E. Marsden, Michael Ortiz, and Matthew West.
An Overview ofVariational Integrators , volume Finite Element Methods: 1970’s and Beyond. Theory andengineering applications of computational methods, pages 98–115. International Center forNumerical Methods in Engineering (CIMNE), Barcelona, 2004.[43] Xiaoye S. Li. An overview of SuperLU: Algorithms, implementation, and user interface.
ACM Transactions on Mathematical Software , 31:302–325, 2005.[44] Xiaoye S. Li, James W. Demmel, John R. Gilbert, Laura Grigori, Piyush Sao, Meiyue Shao,and Ichitaro Yamazaki. SuperLU Users’ Guide. Technical Report LBNL-44289, LawrenceBerkeley National Laboratory, 1999.[45] Jerrold E. Marsden and Tudor S. Ratiu.
Introduction to Mechanics and Symmetry . Springer,2002.[46] Jerrold E. Marsden and Jeffrey M. Wendlandt. Mechanical systems with symmetry, vari-ational principles, and integration algorithms. In Mark Alber, Bei Hu, and JoachimRosenthal, editors,
Current and Future Directions in Applied Mathematics , pages 219–261.Birkhäuser Boston, 1997. ISBN 978-1-4612-7380-6. doi: 10.1007/978-1-4612-2012-1.[47] Jerrold E. Marsden and Matthew West. Discrete mechanics and variational integrators.
Acta Numerica , 10:357–514, 2001. doi: 10.1017/S096249290100006X.[48] Jerrold E. Marsden, George W. Patrick, and Steve Shkoller. Multisymplectic Geometry,Variational Integrators, and Nonlinear PDEs.
Communications in Mathematical Physics ,199:351 – 395, 1998. doi: 10.1007/s002200050505. arXiv:math/9807080.
49] Jerrold E. Marsden, Sergey Pekarsky, Steve Shkoller, and Matthew West. VariationalMethods, Multisymplectic Geometry and Continuum Mechanics.
Journal of Geometry andPhysics , 38:253–284, 2001. doi: 10.1016/S0393-0440(00)00066-8. arXiv:math/0005034.[50] Philip J. Morrison. Hamiltonian Description of the Ideal Fluid.
Reviews of Modern Physics ,70:467–521, 1998. doi: 10.1103/RevModPhys.70.467.[51] Jürgen Moser and Alexander P Veselov. Discrete versions of some classical integrablesystems and factorization of matrix polynomials.
Communications in Mathematical Physics ,139:217–243, 1991.[52] Anders Henry Nielsen and Jens Juul Rasmussen. Formation and temporal evolution of thelamb-dipole.
Physics of Fluids , 9:982–991, 1997. doi: 10.1063/1.869193.[53] Emmy Noether. Invariante variationsprobleme.
Nachrichten der Königlichen Gesellschaftder Wissenschaften zu Göttingen , pages 235–257, 1918. arXiv:physics/0503066.[54] Peter J. Olver.
Applications of Lie Groups to Differential Equations . Springer, 1993.[55] Peter J. Olver.
Equivalence, Invariants and Symmetry . Cambridge University Press, 1995.[56] Dmitry Pavlov, Patrick Mullen, Yiying Tong, Eva Kanso, Jerrold E. Marsden, and MathieuDesbrun. Structure-preserving discretization of incompressible fluids.
Physica D: NonlinearPhenomena , 240:443–458, 2011. doi: 10.1016/j.physd.2010.10.012. arXiv:0912.3989.[57] Rick Salmon and Lynne D. Talley. Generalizations of arakawa’s jacobian.
Journal ofComputational Physics , 83:247–259, 1989. doi: 10.1016/0021-9991(89)90118-6.[58] David J. Saunders.
The Geometry of Jet Bundles . Cambridge University Press, 1989.[59] Robert L. Seliger and Gerald B. Whitham. Variational Principles in Continuum Mechanics.
Proceedings of the Royal Society A , 305:1–25, 1968. doi: 10.1098/rspa.1968.0103.[60] Enzo Tonti. Variational formulations of nonlinear differential equations.
Bulletins del’Academie Royale des Sciences, des Lettres et des Beaux Arts de Belgique , LV:137–165,262–278, 1969.[61] Stefan van der Walt, S. Chris Colbert, and Gaël Varoquaux. The numpy array: A structurefor efficient numerical computation.
Computing in Science Engineering , 13:22–30, 2011.ISSN 1521-9615. doi: 10.1109/MCSE.2011.37. arXiv:1102.1523.[62] Alexander P. Veselov. Integrable discrete-time systems and difference operators.
FunctionalAnalysis and Its Applications , 22:83–93, 1988. doi: 10.1007/BF01077598.[63] Alexander P. Veselov. Integrable lagrangian correspondences and the factorization of matrixpolynomials.
Functional Analysis and Its Applications , 25:112–122, 1991. doi: 10.1007/BF01079590.[64] Jeffrey M. Wendlandt and Jerrold E. Marsden. Mechanical integrators derived from adiscrete variational principle.
Physica D Nonlinear Phenomena , 106:223–246, 1997. doi:10.1016/S0167-2789(97)00051-1., 106:223–246, 1997. doi:10.1016/S0167-2789(97)00051-1.