Mimetic Finite Difference methods for Hamiltonian wave equations in 2D
Abstract
In this paper we consider the numerical solution of the Hamiltonian wave equation in two spatial dimensions. We construct a two step procedure in which we first discretize the space by the Mimetic Finite Difference (MFD) method and then we employ a standard symplectic scheme to integrate the semidiscrete Hamiltonian system derived. The main characteristic of the MFD methods, when applied to stationary problems, is to mimic important properties of the continuous system. This approach yields a full numerical procedure suitable to integrate Hamiltonian problems. A complete theoretical analysis of the method and some numerical simulations are developed in the paper.
1 Introduction
Because of the symplectic structures, Hamiltonian partial differential equations (PDEs) are used to give a mathematical representation of many physical systems and are of interest to various applicative fields, see for instance quantum field theory, meteorology, nonlinear optics, weather forecast.
An important requirement that any numerical method for Hamiltonian PDEs has to satisfy is the preservation of the intrinsic geometric properties of the original continuous problem. In particular, the numerical procedure should preserve the symplectic structure of the Hamiltonian system during numerical simulations. A standard procedure to derive a suitable method for an infinitedimensional Hamiltonian PDE consists into two steps: in the first one the system is discretized in space in order to obtain a finitedimensional Hamiltonian system, and then the semidiscretized system is solved in time by a symplectic integrator [MR1904823, MR2657217, MR2132573, gawlik2011, gawlik2014, demoures2015]. There exists also a recent approach in which the space and time are considered on equal footing, this approach requires a multisymplectic formulation of the system and leads to the multisymplectic numerical schemes for the numerical solution of the PDEs (see [MR1854689, MR2220764, MR2222808, marsden1998]).
The effectiveness of this approach is ensured by the property that the derived semidiscrete system is a finitedimensional Hamiltonian system of ordinary differential equations (ODEs). The space discretization of a Hamiltonian system is usually performed by one of the following techniques: finite difference methods, finite element methods, spectral methods, pseudospectral methods, Fourier expansion, wavelet based methods (see for instance [MR1472194, MR1997061, MR2211034, MR2425152, MR2586202]). However, these semidiscretization approaches could become very expensive or could not be applicable when the space dimension is greater than .
Instead, in this paper we consider the Mimetic Finite Difference (MFD) method to approximate the continuous problem combined with a standard symplectic integration in time to integrate the derived semidiscrete Hamiltonian system.
The main results about MFD methods, for stationary problems, can be found in the recent book [BLM11book] and papers [MR3133437, MR2192322] where, in particular, the theoretical framework of the mimetic spaces and the discretization of the operators are introduced. Significative applications of MFD methods may be found for instance in [MR2339994, MR2512497, MR2534128, MR2568590, MR2655949, MR2737630]. Among the first publication in this field it is worth mentioning [MR1383592, MR1426277] where a first approach to mimetic discretization of the continuous operators can be found and the fundamental papers [MR2168945, MR3133438] where the modern approach to MFD was introduced. A generalization of the MFD methods has been recently proposed, the virtual element methods (VEMs); we cite [MR3073346, MR2997471, VEMhelmholtz, vaccabis, vaccahyper, vaccadivfree] as a very short representative list.
Recently in [vacca], MDF methods has been applied to the space discretization of PDEs of parabolic type in two dimension, showing how this technique preserves invariants of the solution better than classical space discretizations such as finite difference methods.
The main characteristic of the MFD methods is to mimic important properties of the continuous system, e.g., conservation laws, symmetry and positivity of the solutions, and the most important properties of the continuous differential operators, including duality and selfadjointness relations. Furthermore MFD methods can be applied for general polygonal and polyhedral meshes of the space domain instead of more standard triangular/quadrilateral grids.
The main novelty of this paper is the use of MFD methods for the space discretization of the nonlinear wave equation in 2D coupled with a standard symplectic method (the implicit midpoint scheme) for the time integration. We derive a full numerical discretization procedure which will exploit the conservative properties of the MFD approach associated to the symplectic features of the time integrator. We show that the mimetic semidiscrete Hamiltonian is preserved in time and we derive the conservation law for the mimetic semidiscrete energy. Furthermore we give a bound for the conservation of the full discretized Hamiltonian and for the conservation of the full discretized energy. We also prove the convergence of the semidiscrete and fully discrete solutions to the solution of the original problem
The paper is organized in the following way. In Section 2 we recall the basic elements of the MFD approach. In Section 3 we recall the mathematical form of the Hamiltonian PDE we wish to study. In Section 4 we apply the MFD method to the continuous problem and we give a result of the convergence of the semidiscrete solution to the continuous solution of the original problem; we define the semidiscrete Hamiltonian and energy density, show their conservation laws. In Section 5 we discretize the semidiscrete system by using a symplectic time integrator, the implicit midpoint rule, of the second order in time. We will prove the convergence of the full discrete numerical solution by providing an error estimate of the second order in space and time. Hence we give a result about the conservation of the discrete Hamiltonian and of the discrete energy of the system. Section 6 is devoted to show some numerical results.
2 Background on Mimetic Finite Differences Methods
In this section, for ease of reading, we recall the basic concepts and notations on MFD methods which will be used to discretize PDEs in the spatial domain where we assume bounded polygon. For more details on this subject we refer the interested reader to the recent book [BLM11book] or to the papers [MR2192322, MR3133437, vacca]. Let a measurable subset of the domain and let a full symmetric positive definite tensor. By making use of standard notation, we consider the following scalar products:
(1)  
(2) 
It is clear that, in the sense of distribution,
thus we get the duality relation with respect to the scalar product (1) and (2)
(3) 
Let be an unstructured mesh of into nonoverlapping simplyconnected polygons with flat faces, where
Let be the set of edges of the polygons in . We use the following notations for the mesh objects: denotes a general cell in the mesh with measure and centroid ; denotes a general edge of the cell with measure and centroid ; indicates the unit normal vector to the edge with preassigned direction; represents the mutual orientation of the vector and the outward normal vector to with respect to the cell .
Moreover, let , and let , then we denote with the subset of of all the elements that are related with , and we indicate with the cardinality of this set. For example denotes all cells sharing face and denotes all faces forming the boundary of cell .
In the following we take on the element the shape regularity assumptions listed, for instance, in [BLM11book, MR2192322]. A possibility is to assume that for all , each element in satisfies:

is starshaped with respect to a ball of radius greater then ,

any two vertexes in are at least apart,
where is the diameter of . The constants and are positive and uniform with respect to the mesh family.
The mesh objects will define the degrees of freedom of the discrete system, that is these will define the space of the discrete pressures and discrete fluxes.
Let
Let be the set of the pressures that are piecewise constant on , i.e.
Given a pressure , we define the interpolant discrete pressure with
The space of the discrete velocities is defined as follows. For all edge we associate a real number and we denote with the vector with components given by the collection of all the . The symbol will represent the vector space of all . Let a vector function, and let us assume that all faceintegrals
exist. Then the interpolant discrete flux of in the space is defined by with
Remark 2.1.
The discrete spaces , and the interpolation operators are defined starting from the degrees of freedom:

, for all , and ,

, for all , and .
Remark 2.2.
There are obvious correspondences:
With a slight abuse of notation we can refer to a function in the discrete functional spaces as a vector and vice versa.
The definition of the mimetic scheme carries on with the discretisation of the differential operators. Let with , then the Divergence Theorem states that
where is the unit outward normal to . Therefore, the continuous operator admits the immediate discretisation with
The operator is called discrete primary operator.
The next step in the construction of the MFD method is the definition of suitable inner products on the discrete functional spaces and that allow to construct the derived operators imposing the duality relations for the discrete operators.
We assume, for the moment, the following scalar products on the vector spaces and :
(4)  
(5) 
where , are suitable symmetric positive definite matrices. These matrices are locally constructed in such a way, on each cell, the corresponding local discrete inner products have to “mimic” the scalar products defined in (1)) and (2). Therefore we would like that
where, in general, with the notation we denote the vector with the degrees of freedom of the function relative to the cell .
As regards the first local inner products, we observe that the vector has a single component, representing the (constant) value of in the cell . Then the only possible quadrature formula is
therefore and . It is clear that the discrete inner products gives the exact value of the continuous one whenever .
The definition of the local scalar product for the fluxes requires a different approach. The key idea is to define suitable consistency and stability constraints in order to introduce algebraic conditions on the elements of the matrix . Without spelling things out, we requires that the following properties are satisfied

consistency: let two vector fields and let their interpolant functions. If is constant in and for each edge in , is constant, then

stability: there exist two positive independent constants and such that
The last preliminary step in the construction of the MFD method is the definition of the derived discrete operators, which are obtained through a duality relation from the primary operators. Let us consider the spaces , equipped respectively with the scalar products (4), (5). From continuous duality relations (3), we can introduce the discrete operator
and impose the duality relation:
for all , from which it follows that
Finally we can introduce the discrete counterpart of the continuous operator , by defining the operator
given by
(6) 
3 The continuous problem
Let be a bounded polygon and let us consider the nonlinear wave equation with homogeneous boundary value problem
(7) 
where is a full symmetric positive definite tensor, and the source term is the derivative of a smooth function . We would observe that no particularly restrictive assumptions on are required, for instance in the sineGordon equation or the ones of polynomial type with respect to may be considered. For seek of simplicity we consider in the proof global Lipschitz, however the convergence results are still valid for local Lipschitz (see Remark 4.2)
(7) admits the equivalent formulation
(8) 
where the initial and boundary conditions are given by
(8) is said Hamiltonian formulation of (7) for which the Hamiltonian
(9) 
is invariant with respect to time along the solution, that is
(10) 
The energy density of the system is defined by
(11) 
The total derivative of with respect to , along the solution of (8), is given by
Let the energy flux, then we have the energy conservation law
(12) 
which is more general than the global conservation of the Hamiltonian. Indeed if the energy conservation law holds, then it is easy to prove that .
4 The semidiscrete problem
By using the MFD approach we can approximate the continuous operators by discrete ones, in order to derive the semidiscrete problem for the wave (7). Then the resulting semidiscrete wave equation reads:
(13) 
where and are the interpolant functions in of the initial data. In the same way, (8) can be discretized in the following form
(14) 
We observe that the semidiscrete (14) preserves the Hamiltonian structure of (8). In light of the definition in Section 2, the Hamiltonian functional in (9) admits the natural mimetic semidiscretization:
(15) 
that will be called mimetic semidiscrete Hamiltonian functional.
We can observe now that, if we denote with the gradient with respect to the variable and with the gradient with respect to the variable , then
and
Hence (13) may be written as a Hamiltonian system of ordinary differential equations (ODEs), that is as:
where is the canonical symplectic matrix while denotes the gradient with respect the variables and .
We can conclude that the MFD approach gives a finitedimensional system of ODEs that retains the Hamiltonian character of the given PDE. Therefore MFD methods can be considered powerful scheme for the spatial discretization of Hamiltonian PDEs.
4.1 Convergence for the semidiscrete problem
Now we will investigate the convergence of the solution of the semidiscrete wave equation (13) to the solution of (8) in norm. Before analysing the error between the solution, we have to show some preliminary technical results.
Let us introduce the energy projection , with defined as the solution of the diffusion problem
(16) 
In particular, for the duality relation between the operators, the projection satisfies
(17) 
In the following we use to denote the norm induced by the scalar product (that is equivalent to norm on ). We will denote with a generic constant, possibly different at each occurrence, independent from the mesh size and the time step size . In order to prove the convergence results, we need the following Lemma (see [MR2192322] for the proof).
Lemma 4.1.
Let us assume the convexity of the domain . Let and let the energy projection of . Then the following estimate holds:
While the next Lemma shows the spectral properties of the operator (see [vacca] for more details).
Lemma 4.2.
The spectrum of satisfies
(18) 
where and are positive and independent constants.
For the treatment of the nonlinear term we have the following lemma.
Lemma 4.3.
Let be a smooth function and let . Then
Proof.
Let for every . Then the nonlinear term may be treated in the following way. Using the Taylor expansion, since
(19) 
for suitable and for every . Then, setting and using (19), we have:
Now, since is constant and, by definition, , we obtain
with
(20) 
Now we observe that for classic Sobolev embedding theory, implies that , then, since is bounded by and the constant value we obtain that for all element . Therefore being a smooth function, . Now, using the Hölder Theorem
(21) 
Now, using standard polynomial approximation results [scott], we have
(22) 
∎
Now we have the instruments for proving the following convergence theorem.
Theorem 4.1.
Proof.
The proof follows the guidelines of Theorem 1 in [nonlinear] for given for the finite element approximation. Let us set
(23) 
We study separately the two terms. The second term represents the error generated by the energy projection; using Lemma 4.1, we obtain
(24) 
For the first term, from (13) and (17), we get
and, since is the solution of (7), we obtain
and in particular
(25) 
for all . For let use define
(26) 
and let in (25) be a function of . Then it is straightforward to see that
(27) 
Let us fix and we set in (27)
in particular we can observe that . Now the duality relation among discrete operators and simple computations yield
thus
(28) 
Integrating (28) with respect to from 0 to , observing that , and by definition , we get
and then
(29) 
Now by definition (26), from Lipschitz assumption on the load , and Lemma 4.3 we have
Therefore, from CauchySwartz inequality and (24)
and always from CauchySwartz and (24)
By collecting the previous estimates in (29), from (24) we get
(30) 
It is straightforward to check that
then by Gronwall inequality it holds that
from which follows the thesis.
∎
Remark 4.1.
The use of the projection in the proof of the theorem seems to be necessary. Indeed if we compute directly as done for example in [MR2586202], we obtain a term of the form
and does not converge to zero. For instance in Figure 1 we plot the asymptotic behaviour of as a function of for , tensor and domain