On the Derivation of Multisymplectic Variational Integrators for Hyperbolic PDEs Using Exponential Functions

: We investigated the derivation of numerical methods for solving partial differential equations, focusing on those that preserve physical properties of Hamiltonian systems. The formulation of these properties via symplectic forms gives rise to multisymplectic variational schemes. By using analogy with the smooth case, we deﬁned a discrete Lagrangian density through the use of exponential functions, and derived its Hamiltonian by Legendre transform. This led to a discrete Hamiltonian system, the symplectic forms of which obey the conservation laws. The integration schemes derived in this work were tested on hyperbolic-type PDEs, such as the linear wave equations and the non-linear seismic wave equations, and were assessed for their accuracy and the effective-ness by comparing them with those of standard multisymplectic ones. Our error analysis and the convergence plots show signiﬁcant improvements over the standard schemes.


Introduction
The theory of variations has been used so far to derive integrators that preserve the physical properties of the system in question, e.g., a symplectic form [1][2][3][4]. When tested on several Hamiltonian ordinary differential equations (ODEs), symplectic integrators have shown robustness, efficiency and accuracy even for long integrations [5][6][7][8][9][10]. Hamiltonian systems that lead to partial differential equations (PDEs) allow a similar description, referred to as the "Veselov discretisation for PDEs in variational form" [3,4]. Due to the presence of more preserved quantities, and hence more symplectic forms, these integrators (named multisymplectic schemes) have been derived on the basis of the structure of Hamiltonian PDEs from Lagrangian or Hamiltonian perspectives [11][12][13][14].
In our recent works [15,16] we investigated the weaknesses (e.g., numerical instabilities) of standard multisymplectic integrators derived via simple quadratic expressions adopted for the spatial and time derivatives [11]. We also explored several schemes that use complex expressions based on the non-standard finite differences of [17]. The latter proved to be more accurate, and their performances are better than those the standard ones. They use trigonometric functions in the derivative expressions and preserve certain properties and structures dictated by the physical system [16,18,19].
In this work, we extend our previous multisymplectic schemes by exploiting the advantages of exponential functions to express the spatial and time derivatives. The focus is on the solution of the wave equation, a hyperbolic-type equation within the theory of linear PDEs. Unlike other PDEs where the solutions are smooth when the initial conditions are smooth, the hyperbolic equations are challenging due to their sharp solutions behaviour [20,21].
In addition, we consider a non-linear version of the standard wave equation, namely, the non-linear seismic wave equation [22]. To that end, we obtain its multisymplectic form directly from the variational principle, i.e., by using the Lagrangian approach. For 2 of 11 simplicity, in this article we work in local coordinates, although a more general coordinate system can be considered. Our attempt leads to multi-symplectic integrators with very good energy performance with respect to conservation of the nearby Hamiltonian, up to an exponentially small error [11,15,16,[23][24][25][26].
The paper is organised as follows. Firstly, we introduce some notions from the smooth multisymplectic Hamiltonian systems in Section 2. These are used in Section 3 to formulate a discrete version of the Lagrangian density relying on exponential functions and triangular space-time discretisation (a square space-time discretisation is additionally presented in the Appendix A). The resulting discrete Lagrangian is then applied in the discrete version of the Euler-Lagrange field equations in order to deduce the corresponding discrete equations. The proposed scheme is tested in the linear wave equation in Section 4 and in the non-linear seismic wave equation in Section 5. Subsequently, error analysis and convergence plots are presented to demonstrate the improved numerical properties and the behaviour of the proposed scheme as compared to the standard ones. Finally, the main conclusions extracted from our study are summarised in Section 6.

Multisymplectic Hamiltonian Systems
A large class of partial differential equations can be cast in the form where z ∈ R n (n ≥ 3), M and K are skew-symmetric operators in R n×n , and S : R n → R is a smooth function [12][13][14]. Equation (1) involves the spatial derivative with respect to only one spatial dimension, the x-coordinate, but it can be extended to arbitrary number of dimensions if required.
For this system one can define the pre-symplectic forms ω associated with time, and κ associated with space, as follows [12].
Through the pre-symplectic forms of the Hamiltonian system, one may write down a multisymplectic conservation law of the form Specifically, for the energy density E and the energy flux F, the corresponding energy conservation law should be written as For the momentum density I and the momentum flux G, the corresponding momentum conservation law takes the form In general, the aim of deriving multisymplectic integrators is to strictly enforce the conservation laws when solving the Hamiltonian system (1); see [12][13][14]. Such integrators can be deduced via symplectic (in both space and time) Runge-Kutta schemes [23], but in this work we use the Lagrangian density of the system, L(u, u x , u t ), referred from now on simply as the Lagrangian.
Subsequently, we employ a Legendre transform to define the Hamiltonian H as a function of the Lagrangian density as Appl. Sci. 2021, 11, 7837 3 of 11 Following these steps, the Hamiltonian system (1) can afterwards be written in terms of L as Furthermore, equivalent expressions for the energy conservation law (4) and the momentum conservation law (5) can also be obtained; see [18,19].

Discrete Lagrangian Density
By mimicking the smooth case and following [11,15,16], the discrete Lagrangian density can readily be figured out. Towards this purpose, we consider fields over a higherdimensional oriented manifold M, together with its tangent bundle TM and cotangent bundle T * M. Then, the Lagrangian density L, is a smooth bundle map over X; i.e., where J 1 (Y) is the first jet bundle over Y defined via the fibre bundle π XY : Y → X over the manifold X. In the latter equation, Λ n+1 (X) represents the bundle of (n + 1)-forms on X (for more details, see [11] and references therein). Next, we consider the space-time discretisation of X when dim(X) = 2, which generalises the Veselov discretisation [3,4] in the multisymplectic field theory. In such a case, the discrete representation of the manifold can be considered as X = Z × Z = (i, j). On the other hand, the fibre bundle Y will obey the standard definition of the smooth case, and hence it will come out of the X × F for some smooth manifold F [11,15,16].
In this work, we restrict the discretisation to a uniform quadrangular mesh in the base space, with mesh lengths ∆x and ∆t. The nodes in this mesh, which are objects in Z × Z, can be denoted by (i, j) and correspond to points (x i , t j ) := (i∆x, j∆t) ∈ R 2 . Using the symbol u for the field values, we denote by u j i its value at the node (i, j). Then, the triplet ((i, j), (i + 1, j), (i, j + 1)) forms a triangle at the (i, j) which in our notation is denoted by ij . Finally, the X will denote the set of all such triangles; see Figure 1.
Multisymplectic Geometry, Variational Integrators, and Nonlinear PDEs 377 Given φ ∈ C U and a vector field V , there is the 1-parameter family of sections for all vector fields V .
The discrete Euler-Lagrange equations. The variational principle gives certain field equations, the discrete Euler-Lagrange field equations (DELF equations), as follows. Focus upon some (i, j) ∈ int U , and abuse notation by writing φ(i, j) ≡ y ij . The action, written with its summands containing y ij explicitly, is (see Fig. 5.2) so by differentiating in y ij , the DELF equations are After the above explainations, the discrete jet bundle can be defined as [11,15,16] which is equivalent to X × R 3 . We remind the reader that an arbitrary value of u is assumed to be an averaging of the field over the triangle vertices (see Figure 1); i.e., and the derivatives of u are expressed by (∆t and ∆x are the time step and spatial step sizes, respectively) where φ(∆t) and ψ(∆x) involve the type of the numerical scheme. Specifically, the selection φ(∆t) = ∆t and ψ(∆x) = ∆x leads to the midpoint rule and the generated schemes can be found in [15,16]. Further, the selection of φ(∆t) = 2 sin ∆t 2 and ψ(∆x) = 2 sin ∆x 2 leads to the nonstandard finite differences, and the schemes obtained can be found in [17].
In this work, in order to truncate the high oscillations, we introduce new expressions for φ(∆t) and ψ(∆x) based on exponential functions as The latter expressions are, furthermore, exploited to obtain the discrete Lagrangian at any triangle, which depends on the chosen edges, as and also to write the discrete Euler-Lagrange field equations as (see Figure 1 right) It is worth noting that the assumed domain can also be discretised by using squares (instead of triangles), and then the appropriate discrete Lagrangian can be defined for each one of the four neighbouring squares around the node (i, j). Obviously, the resulting Euler-Lagrange equations are more complicated this way (see Appendix A).
In the following two Sections, in order to assess the integrators proposed in this section, we test their performances on two concrete examples: (i) the linear wave equation, and (ii) the non-linear seismic wave equation. Both of them are hyperbolic-type PDEs.

The Linear Wave Equation
It is well known that the linear wave equation contains second order partial derivatives of the unknown wave function u(x, t) with respect to time and space, respectively (see, e.g., [20,21] where u t and u x denote the derivatives ∂u/∂t and ∂u/∂x, respectively. This equation describes one-dimensional (1D) waves-e.g., the longitudinal string oscillations propagating along the x-axis. The corresponding Lagrangian is [15,16] L(u, u t , By adopting the triangle discretisation described in Section 3, one gets the discrete Lagrangian where ∆t and ∆x can be thought of as the mesh lengths for time and space, respectively, as stated in Section 3. The insertion of the discrete Lagrangian L d into the discrete Euler-Lagrange field Equation (14) yields the discrete version of (15); i.e., The present example was tested on the space domain bounded as −10 < x < 10, with initial conditions and boundary conditions The chosen grid discretisation is ∆t = 0.1 and ∆x = 0.2. The resulting solution u(x, t) is shown in Figure 2 by a 3D diagram. It can be seen that the time evolution of the solution u(x = const., t) is a continuous function in which the periodicity is preserved. As a second test and for the sake of comparison with the standard multisyplectic methods, we examine the energy behaviour in both x-directions. Firstly, we considered the temporal energy error (deviation from the initial value) at a given spatial point x via Secondly, we calculated the spatial energy error at a given time instant t via The obtained results are plotted in Figure 3 and are compared with those of the standard multisymplectic scheme of [11]. The graphs of this figure represent relative errors, i.e., differences from the initial values. While during the numerical simulation, no energy loss or energy blow-up appeared to occur in the system in either case, the proposed scheme improved both energy errors dramatically-the spatial and the temporal ones.  [11] (blue) versus the proposed one (red). Top: temporal evolution of the discrete energy error (20). Bottom: spatial evolution of the discrete energy error (21).

The Non-Linear Seismic Wave Equation
In the second concrete application with our method, we seek solutions of the nonlinear seismic wave equation. Its general form, which describes the propagation of waves in elastic media, is [22] − 1 [c(x, y, z)] 2 u tt + u xx + u yy + u zz = 0, where the function c(x, y, z) is the wave velocity reflecting spatially variable elastic properties of the medium. It should be noted that the non-linearity of the seismic wave Equation (22) enters through the square of the spatially dependent velocity c(x, y, z). Equation (22) has two types of solutions that correspond to compressional and shear waves, which share many common characteristics with the solutions of the 1D wave Equation (15). Moreover, an expression equivalent to (22) can be obtained by starting from the Hamiltonian system (1) where [14] M = for z = [u v w] T and S(z) = (v 2 − w 2 )/2. The pre-symplectic forms ω and κ given by (2) can also be defined so that the Hamiltonian system satisfies the multisymplectic conservation law (3). In addition, Equation (22) can also be derived from the Lagrangian inserted into the Euler-Lagrange field Equations (7). For simplicity, we consider a plane wave propagating in x-direction; i.e., we seek solutions of (22) for compressional waves travelling in x-direction. The resulting equation is written as By discretising the corresponding Lagrangian of (24) and using a simple quadrature rule for the discrete version of the wave velocity, we obtain the integrator Furthermore, in order to compare the latter integrator with previously presented schemes, we follow the setup of [27]. The wave velocity is assumed to be given by The solution is sought in the domain −10 < x < 10, with the initial conditions [27] u(x, 0) = sechx, u t (x, 0) = 0, and with the periodic boundary condition u(−10, t) = u(10, t).
(28) Figure 4 shows the wave-forms of the seismic wave obtained via the discrete Equation (26) with ∆x = 0.1 and ∆t = 0.01. In addition, Figure 5 illustrates the waveforms at four distinct time instances: t = 0, t = 4, t = 10 and t = 15. These have been selected following [27] to visualise the split of energy of the initial wave to two smaller but equal waves travelling in different (opposite) directions, which, due to periodic boundary conditions, reform the initial wave. In a similar way to that followed in the previous Section, we calculated the evolution of the relative energy error in the time interval t ∈ [0, 400], as the difference between the current and the initial energy. The results are plotted in Figure 6, where the blue dashed line corresponds to the error obtained by the standard multisymplectic scheme of [11], which is identical to the one presented in [27], and the red line represents the error evolution in the new scheme. Although both results show the same behaviour qualitatively-i.e., the resulting peaks appear at the same time frames-the ones obtained by using the exponential functions (12) are much smaller. Finally, Figure 7 demonstrates the comparison of the convergence plots of the numerical solution of the non-linear seismic wave equation from the proposed method and the standard multisymplectic scheme of [11]. For the different choices of ∆x employed, both schemes show the same order of convergence, which is indicated by the slopes of the curves. However, the proposed scheme shows better behaviour, with the error being an order of magnitude smaller compared to that obtained in [27]. We mention that about the same result was found in the application of the square space-time discretisation presented in the Appendix A. Before closing, it is worth noting that this work can be extended in several different directions, and here we briefly summarise some specific cases. The first application may involve non-linear phenomena described through the Schrödinger and Klein-Gordon equation [28]. Furthermore, the magneto-hydro-dynamical description of the astrophysical plasma, for example, inside the astrophysical outflows (jets), may be treated in a similar way to that followed for the non-linear seismic wave. In [29,30], Runge-Kutta numerical integrators were included that could be replaced with an extended and improved version of the present work. The latter two applications are of great interest.

Conclusions
The interest in deriving numerical schemes that preserve the physical properties of a system is growing as computational power increases. When the principle was first introduced, the geometric properties were well proven, but the computational costs were significant. The present computers, however, allow us to take full advantage of the multisymplectic integrators, which have been proven to have very good energy behaviour in terms of conservation of the nearby Hamiltonian and the symplectic forms, especially for long term integration processes.
In this work, we have proposed an advanced integration scheme that uses exponential functions at the level of the discrete Lagrangian density. The great advantages of this scheme have been demonstrated through accurate solutions of the linear wave equations and the non-linear seismic wave equations, i.e., of hyperbolic-type equations within the context of the linear PDE theory. In general, the solutions of such equations, when the initial conditions are smooth, are not always smooth, and thus the presence of the exponential functions in the derived integrators significantly improves the performance of multisymplectic-type schemes.
Future work may involve extension and refinement of the present method with applications in various non-linear phenomena and descriptions of the known shock waves inside astrophysical plasmas.

Acknowledgments:
The authors acknowledge gratefully the support of the Engineering and Physical Sciences Research Council (EPSRC), UK, via grand EP/N026136/1.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A. Square Discretization
For square discretisation we denote a square at (i, j) with ordered quaternion ((i, j), (i + 1, j), (i + 1, j + 1), (i, j + 1)) by i j , and the set of all such squares by X ; see Figure A1. The discrete jet bundle is defined as which is equal to X × R 4 . For more details, see [11,16] and references therein.
(i, j) (i, j + 1) (i + 1, j + 1) (i + 1, j) 8 y i j+1 = k 2 (y i+1 j − 2 y ij + y i−1 j ) + 2 y i which is stable whenever the Courant stability condition is Extensions: Jets from rectangles and other polygons. Our is obviously not restricted to triangles, and can be extended polygons (left of Fig. 5.4). A rectangle is a quadruple of t = (i, j), (i, j + 1), (i + 1, j + 1), (i a point is an interior point of a subset U of rectangles if U touching that point, the discrete Lagrangian depends on v DELF equations become ∂L ∂y 1 (y ij , y i j+1 , y i+1 j+1 , y i+1 j ) + ∂L ∂y 2 (y i j−1 , y ij , y i+1 j + ∂L ∂y 3 (y i−1 j−1 , y i−1 j , y ij , y i j−1 ) + ∂L ∂y 4 (y i−1 j The extension to polygons with even higher numbers of example is illustrated on the right of Fig. 5   By averaging the fields over all vertices of the square, the field u can be obtained as (see Figure A1) As above, the expressions for the derivatives can be taken from [16] for the discrete Lagrangian, which now depends on the edges of the square; i.e., see Figure A1 (right).