Mimetic Finite Difference methods for Hamiltonian wave equations in 2D

# Mimetic Finite Difference methods for Hamiltonian wave equations in 2D

L. Beirão da Veiga lourenco.beirao@unimib.it Dipartimento di Matematica e Applicazioni, Università degli Studi di Milan-Bicocca, Via Roberto Cozzi, 55 - 20125 L. Lopez luciano.lopez@uniba.it Dipartimento di Matematica, Università degli Studi di Bari, Via Edoardo Orabona, 4 - 70125 Bari G. Vacca giuseppe.vacca@unimib.it Dipartimento di Matematica e Applicazioni, Università degli Studi di Milan-Bicocca, Via Roberto Cozzi, 55 - 20125
July 11, 2019
###### 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 semi-discrete 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 infinite-dimensional Hamiltonian PDE consists into two steps: in the first one the system is discretized in space in order to obtain a finite-dimensional Hamiltonian system, and then the semi-discretized 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 multi-symplectic formulation of the system and leads to the multi-symplectic 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 semi-discrete system is a finite-dimensional 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 semi-discretization 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 semi-discrete 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, VEM-helmholtz, 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 self-adjointness 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 semi-discrete Hamiltonian is preserved in time and we derive the conservation law for the mimetic semi-discrete 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 semi-discrete 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 semi-discrete solution to the continuous solution of the original problem; we define the semi-discrete Hamiltonian and energy density, show their conservation laws. In Section 5 we discretize the semi-discrete 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:

 (u,v)L2(ω) :=∫ωuvdxfor all u,v∈L2(ω), (1) (ω,σ)K,ω :=∫ωK−1ω⋅σdxfor all ω,σ∈(L2(ω))2. (2)

It is clear that, in the sense of distribution,

 (K∇u,σ)K,Ω=−(u,divσ)L2(Ω)for all u∈L2(Ω), % σ∈H(div,Ω)

thus we get the duality relation with respect to the scalar product (1) and (2)

 K∇=−(div)∗. (3)

Let be an unstructured mesh of into nonoverlapping simply-connected polygons with flat faces, where

 h:=supc∈Thdiameter(c).

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 star-shaped 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

 Nc:=|Th|,Nf:=|Eh|,N∗:=maxc|Eh(c)|.

Let be the set of the pressures that are piecewise constant on , i.e.

 Ch:=\setu∈L2(Ω)u|c=const,∀c∈Th.

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 face-integrals

 ∫fω⋅nfdS,for all f∈Eh

exist. Then the interpolant discrete flux of in the space is defined by with

 ωf=1|f|∫fω⋅nfdS,for all f∈Eh.
###### 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:

 Ch≅RNcu↦(uc)c∈Th,andFh≅RNfω↦(ωf)f∈Eh.

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

 ∫cdivωdx=∫∂cω⋅ndS ,

where is the unit outward normal to . Therefore, the continuous operator admits the immediate discretisation with

 (DIVωh)c=1|c|∑f∈Eh(c)αc,f|f|ωf∀ωh∈Fh.

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 :

 [uh,vh]Ch:=uThMChvhfor all uh,vh∈Ch, (4) [ωh,σh]Fh:=ωThMFhσhfor all ωh,σh∈Fh, (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

 [uh,c,vh,c]Ch,c=:(uh,c)TMCh,cvh,c≈(uh,vh)L2(c),for % all uh,vh∈Ch, Fh,c=:(ωh,c)TMFh,cσh,c≈(ωh,σh)K,c,for all ωh,σh∈Fh,

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

 [uh,c,vh,c]Ch,c=(uh,c)TMCh,cvh,c=|c|ucvc,

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

 [ωh,c,σh,c]Fh,c=∫cK−1ω⋅σdc;
• stability: there exist two positive -independent constants and such that

 C∗|c|(ωh,c)Tωh,c≤(ωTh,c)MFh,cωh,c≤C∗|c|(ωTh,c)ωh,c∀ωh,c∈Fh.

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

 Δh:Ch→Ch

given by

## 3 The continuous problem

Let be a bounded polygon and let us consider the nonlinear wave equation with homogeneous boundary value problem

 ⎧⎪⎨⎪⎩utt(x,t)=divK∇u(x,t)−f′(u(x,t))in Ω×(0,T)u(x,0)=u0(x) ,ut(x,0)=v0(x)in Ωu(x,t)=0on ∂Ω×(0,T)% (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 sine-Gordon 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)

 {ut(x,t)=v(x,t)in % Ω×(0,T)vt(x,t)=divK∇u(x,t)−f′(u(x,t))in Ω×(0,T) (8)

where the initial and boundary conditions are given by

 u(x,0)=u0(x),v(x,0)=v0(x)in Ωu(x,t)=0,v(x,t)=0 ,on ∂Ω×(0,T) .

(8) is said Hamiltonian formulation of (7) for which the Hamiltonian

 H[u,v]:=∫Ω(12v2+12∇u⋅K∇u+f(u))dx (9)

is invariant with respect to time along the solution, that is

 ddtH[u,v]=0 . (10)

The energy density of the system is defined by

 E(u,v):=12v2+12∇u⋅K∇u+f(u). (11)

The total derivative of with respect to , along the solution of (8), is given by

 Et=(divK∇u)v+∇u⋅K∇v=div(vK∇u).

Let the energy flux, then we have the energy conservation law

 Et(u,v)+divω(u,v)=0 , (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 semi-discrete problem

By using the MFD approach we can approximate the continuous operators by discrete ones, in order to derive the semi-discrete problem for the wave (7). Then the resulting semi-discrete wave equation reads:

 {uh,tt(t)=Δhuh(t)−f′(uh(t))for t∈(0,T),uh(0)=uh,0,uh,t(0)=vh,0 , (13)

where and are the interpolant functions in of the initial data. In the same way, (8) can be discretized in the following form

 ⎧⎪⎨⎪⎩uh,t(t)=vh(t)% for t∈(0,T),vh,t(t)=Δhuh(t)−f′(uh(t))for t∈(0,T),uh(0)=uh,0,vh(0)=vh,0 . (14)

We observe that the semi-discrete (14) preserves the Hamiltonian structure of (8). In light of the definition in Section 2, the Hamiltonian functional in (9) admits the natural mimetic semi-discretization:

that will be called mimetic semi-discrete 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

 M−1Ch∇vhHh[uh,vh]=M−1ChMChvh=vh ,

and

Hence (13) may be written as a Hamiltonian system of ordinary differential equations (ODEs), that is as:

 (uh,tvh,t)=JNcM−1Ch∇H[uh,vh] ,

where is the canonical symplectic matrix while denotes the gradient with respect the variables and .

We can conclude that the MFD approach gives a finite-dimensional 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 semi-discrete problem

Now we will investigate the convergence of the solution of the semi-discrete 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

 Δh(Phu)=(divK∇u)I. (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:

 ∥uI−Phu∥Ch≤Ch2|u|H2(Ω).

While the next Lemma shows the spectral properties of the operator (see [vacca] for more details).

###### Lemma 4.2.

The spectrum of satisfies

 σ(−Δh)⊆[s∗,s∗h−2], (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

 ∥f′(u)I−f′(uI)∥Ch≤C(u)h2.
###### Proof.

Let for every . Then the nonlinear term may be treated in the following way. Using the Taylor expansion, since

 f′(u(x))=f′(uc)+f′′(uc)(u(x)−uc)+12f′′′(ˆuc(x))(u(x)−uc)2for % a.e. x∈c (19)

for suitable and for every . Then, setting and using (19), we have:

 f′(u)c−f′(uc)=1|c|∫c(f′(u(x))−f′(uc))dx=1|c|∫c(f′(uc)+f′′(uc)(u(x)−uc)+12f′′′(ˆuc(x))(u(x)−uc)2−f′(uc))dx.

Now, since is constant and, by definition, , we obtain

 f′(u)c−f′(uc)=σ(u)c

with

 σ(u)c:=12|c|∫cf′′′(ˆuc(x))(u(x)−uc)2dx. (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

 ∥σ(u,t)∥Ch=12∑c∈Th∣∣∣∫cf′′′(ˆuc(x))(u(x)−uc)2dx∣∣∣≤12∑c∈Th∫c|f′′′(ˆuc(x))|(u(x)−uc)2dx≤12∑c∈Th∥f′′′(ˆuc)∥L∞(c)∥(u−uc)2∥L1(c)≤C∑c∈Th∥(u−uc)2∥L1(c). (21)

Now, using standard polynomial approximation results [scott], we have

 ∥σ(u,t)∥Ch≤C∑c∈Th∥(u−uc)2∥L1(c)=C∑c∈Th∥(u−uc)∥2L2(c)≤C∑c∈Thh2|u|2H1(c)=Ch2|u|2H1(Ω). (22)

Now we have the instruments for proving the following convergence theorem.

###### Theorem 4.1.

Under the assumptions of Lemma 4.1 and Lemma 4.2, let be the solution of (7) and be the solution of (13). Let us assume that for all . Moreover let us assume that is globally Lipschitz. Then, for all , it follows that:

 ∥u(t)I−uh(t)∥Ch≤Cψ(T)h2(|u0|H2(Ω)+|v0|H2(Ω)+|ut|L1(0,t,H2(Ω))++|utt(t)|L2(0,t,H2(Ω))+|u(t)|L2(0,t,H2(Ω))+|u(t)|2L2(0,t,H1(Ω))),

where denotes the interpolant of in and the scalar function is bounded for all .

###### Proof.

The proof follows the guidelines of Theorem 1 in [nonlinear] for given for the finite element approximation. Let us set

 uh(t)−u(t)I=(uh(t)−Phu(t))+(Phu(t)−u(t)I)=:ϑ(t)+ϱ(t). (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

 ϑtt(t)−Δhϑ(t)=−f′(uh(t))−(Phu(t))tt+(divK∇u(t))I

and, since is the solution of (7), we obtain

 ϑtt(t)−Δhϑ(t)=−f′(uh(t))−Phutt(t)+utt(t)I+(f′(u(t)))I=−ϱtt(t)−(f′(uh(t))−f′(u(t))I)

and in particular

 [ϑtt(t),χ]Ch−[Δhϑ(t),χ]Ch=−[ϱtt(t),χ]Ch−[f′(uh(t))−f′(u(t))I,χ]Ch (25)

for all . For let use define

 G(t):=∫t0−(f′(uh(s))−f′(u(s))I)ds, (26)

and let in (25) be a function of . Then it is straightforward to see that

 −[ϑt(t),χt(t)]Ch−[Δhϑ(t),χ(t)]Ch=ddt[(uIt−uh,t)(t)+G(t),χ(t)]Ch++[ϱt(t),χt(t)]Ch−[G(t),χt(t)]Ch. (27)

Let us fix and we set in (27)

 χ(t):=∫τtϑ(s)ds,for t∈[0,T],

in particular we can observe that . Now the duality relation among discrete operators and simple computations yield

thus

Integrating (28) with respect to from 0 to , observing that , and by definition , we get

and then

 ∥ϑ(τ)∥2Ch≤∥ϑ(0)∥2Ch+2∫τ0∥G(t)∥Ch∥ϑ(t)∥Chdt+2∫τ0∥ϱ(t)∥Ch∥ϑ(t)∥Chdt. (29)

Now by definition (26), from Lipschitz assumption on the load , and Lemma 4.3 we have

 ∥G(t)∥Ch≤∫t0∥f′(uh(s))−f′(u(s))I∥Chds≤∫t0∥f′(uh(s))−f′(u(s)I)∥Chds+∫t0∥f′(u(s))I−f′(u(s)I)∥Chds≤C∫t0∥ϑ(s)∥Chds+C∫t0∥ϱ(s)∥Chds+Ch2∫t0|u(s)|2H1(Ω)ds.

Therefore, from Cauchy-Swartz inequality and (24)

 2∫τ0∥G(t)∥Ch∥ϑ(t)∥Chdt≤2C∫τ0(∫t0∥ϑ(s)∥Chds+∫t0∥ϱ(s)∥Chds+h2∫t0|u(s)|2H1(Ω)ds)∥ϑ(t)∥Chdt≤C∫τ0(∫t0∥ϱ(s)∥Chds+h2∫t0|u(s)|2H1(Ω)ds)2dt+C∫τ0(1+τ)∥ϑ(t)∥2Chdt≤C(u)Th4+C(1+T)∫τ0∥ϑ(t)∥2Chdt

and always from Cauchy-Swartz and (24)

 2∫τ0∥ϱ(t)∥Ch∥ϑ(t)∥Chdt≤∫τ0∥ϱ(t)∥2Chdt+∫τ0∥ϑ(t)∥2Chdt≤C(u)Th4+∫τ0∥ϑ(t)∥2Chdt.

By collecting the previous estimates in (29), from (24) we get

 ∥ϑ(τ)∥2Ch≤∥ϑ(0)∥2Ch+C(u)Th4+(1+C+CT)∫τ0∥ϑ(t)∥2Chdt. (30)

It is straightforward to check that

 ∥ϑ(0)∥Ch=∥ϱ(0)∥Ch≤C(u0)h2,

then by Gronwall inequality it holds that

 ∥ϑ(τ)∥2Ch≤C(u,u0)h4TeLT.

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

 L(u):=∥∥(divK∇u)I−ΔhuI∥∥Ch ,

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