Spectral properties and conservation laws in Mimetic Finite Difference methods for PDEs

Accepted Manuscript Spectral properties and conservation laws in Mimetic Finite Difference methods for PDEs Luciano Lopez, Giuseppe Vacca PII: DOI: Reference:

S0377-0427(15)00036-9 http://dx.doi.org/10.1016/j.cam.2015.01.024 CAM 9978

To appear in:

Journal of Computational and Applied Mathematics

Received date: 17 September 2014 Revised date: 9 December 2014 Please cite this article as: L. Lopez, G. Vacca, Spectral properties and conservation laws in Mimetic Finite Difference methods for PDEs, Journal of Computational and Applied Mathematics (2015), http://dx.doi.org/10.1016/j.cam.2015.01.024 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Click here to view linked References

Spectral Properties and Conservation Laws in Mimetic Finite Difference Methods for PDEs Luciano Lopeza , Giuseppe Vaccaa,∗ a Dipartimento

di Matematica, Università degli Studi di Bari, Via Orabona, 4 - 70125 Bari

Abstract The Mimetic Finite Difference (MFD) methods for PDEs mimic crucial properties of mathematical systems: duality and self-adjointness of differential operators, conservation laws and properties of the solution on general polytopal meshes. In this article the structure and the spectral properties of the linear systems derived by the spatial discretization of diffusion problem are analysed. In addition, the numerical approximation of parabolic equations is discussed where the MFD approach is used in the space discretization while implicit ϑ-method and explicit Runge Kutta Chebyshev schemes are used in time discretization. Moreover, we will show how the numerical solution preserves certain conservation laws of the theoretical solution. Keywords: Mimetic Finite Difference, Runge-Kutta Chebyshev schemes, Conservation Laws.

1. Introduction The Mimetic Finite Difference (MFD) methods are numerical techniques to approximate the solutions of partial differential equations (PDEs). The main idea of MFD methods is to mimic important properties of the continuous models, 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 the MFD methods can be applied for general polygonal and polyhedral meshes instead of more standard triangular/quadrilateral (tetrahedral/hexahedral) grids. This feature has motived a recent growth of interest in the engineering and numerical analysis literature. The MFD methods (see for instance the recent book [7]) can be constructed in the following way. The first step is to define a mimetic discretization of scalar, vector and tensor fields through the choice of suitable degrees of freedom that are strictly connected with the mesh objects. We can associate the degrees of freedom with vertices, edges, faces and cells of the mesh and then define the nodal functions, edges functions, faces functions and cell functions. These discrete spaces, with the straightforward definition of addition and multiplication by a real number, set up a finite-dimensional linear spaces. The second step is to construct the mimetic discrete operators of the first-order operators that is gradient, divergence and curl. The idea is to define the discrete operators using the integral identity given by various formulation of the Stokes formula. The discrete operators came out of the Stokes formula are called primary operators. The third step is to define, on the spaces of the grid functions, suitable discrete scalar products. The definition of these scalar products has fundamental role in the construction of the mimetic scheme. The discrete scalar products have to satisfy, locally, suitable consistency and stability conditions. Finally we construct the other discrete operators, the derived operators, from the primary operators by imposing the discrete duality relations in the discrete spaces of the grid functions. In [19], the main results about MFD methods are shown, in particular the theoretical framework of the mimetic functions and the discretization of the operators are introduced and the most important applications of the MFD methods, e.g., the diffusion equation, the Maxwell’s ∗ Corresponding

author Email addresses: [email protected] (Luciano Lopez), [email protected] (Giuseppe Vacca)

Preprint submitted to Elsevier

December 9, 2014

equations, the magnetostatic problem, are considered. In [11] is proved a fundamental result about the convergence of MFD methods for diffusion problems on general polyhedral meshes. The MFD method has been successfully employed for solving magnetostatic field problem [18], electromagnetic problem [9], Stokes problem [6], linear elasticity problem [3] and diffusion problem [10, 13]. Among the first publication in this field it is worth mentioning [22] where can be found a first approach to mimetic discretization of the continuous operators. Recently a new approach that generalizes MFD methods has been proposed: the virtual element methods (VEMs) [1, 4, 5]. Instead, in [14, 2, 15] the authors analyse the structures of the matrices involved in MFD methods and suggest numerical methods for solving the resulting linear systems. In this paper we develop the analysis of the matrices involved in the MFD methods for linear diffusion problem with Dirichlet homogeneous boundary conditions. We study the sparsity and the spectral properties of the matrices and we derive the maximum and minimum eigenvalues and therefore the norm and the condition number of the matrices involved. Furthermore we apply the MFD methods to parabolic problems. We consider as a model problem the classical time-dependent diffusion equation. We discuss the stiffness property of the semi-discrete system. We consider two time integrator schemes for solving the semi-discrete system: the ϑ-method, that, when ϑ > 0, is an implicit method suitable for solving stiff ordinary differential equations (ODEs) and the Runge-Kutta Chebyshev (RKC) methods that are explicit methods but with large stability regions (see for instance [24, 26, 16, 27]). Finally we treat the parabolic problem coupled with Neumann boundary conditions, in particular we derive conservative properties for the semi-discrete and discrete problem. The paper is organized as follows. In Section 2 we consider the continuous linear diffusion problem and we recall the most important steps of the MFD method. In Section 3 we analyse the spectral properties of the matrices involved in the MFD methods. In Section 4 we treat the parabolic problem and we derive the discrete conservation laws for the system. Finally, in Section 5 we report some numerical tests. 2. Background Let Ω ⊆ Rd be a bounded Lipschitz polyhedron, when d = 3, or polygon, when d = 2, and let us consider the Darcy problem, e.g., the linear diffusion problem in mixed form: ( u = −K∇p in Ω constitutive equation, (1) ∇·u=f in Ω conservative equation, where K ∈ (W 1,∞ )d×d is a full symmetric positive definite tensor, f ∈ L2 (Ω) is a source term, while the boundary conditions, B(p, u) on ∂Ω, will be one of the following forms: p=ψ

on ∂Ω

Dirichlet boundary conditions,

(2a)

BN (u) :−u · n = ψ

on ∂Ω

Neumann boundary conditions,

(2b)

BD (p) :

with n unit outward normal vector to ∂Ω and ψ ∈ L2 (∂Ω). Let us consider the functional spaces H01 (Ω), H(Ω) := H 1 (Ω) ⊕ L2 (∂Ω), H(Ω) := H(div; Ω), L2 (∂Ω) that we endow with the following scalar products Z (p, q)H01 (Ω) := p q dV ∀p, q ∈ H01 (Ω), (3) Ω Z Z ((p, pb ), (q, qb ))H(Ω) := p q dV + pb qb dS ∀p, q ∈ H 1 (Ω), ∀pb , qb ∈ L2 (∂Ω), (4) ∂Ω ZΩ −1 (u, w)H(Ω) := K u · w dV ∀u, w ∈ H(Ω), (5) ZΩ (ψ, ϕ)L2 (∂Ω) := ψ ϕ dS ∀ψ, ϕ ∈ L2 (∂Ω). (6) ∂Ω

2

Let us define the operators DIV : H(Ω) → L2 (Ω) and GRAD : H01 (Ω) → L2 (Ω; Rd ) with DIVu := ∇ · u

in Ω

∀u ∈ H(Ω);

∀p ∈ H01 (Ω)

GRADp := K∇p

(7)

thus, using Green’s formula, we easily get the duality relation: GRAD = −(DIV)∗

(8)

where (DIV)∗ is the adjoint of DIV with respect to scalar products (3) and (5). Similarly, by introducing the operators DIV : H(Ω) → L2 (Ω) ⊕ L2 (∂Ω),

and GRAD : H 1 (Ω) → L2 (Ω; Rd )

with DIVu :=

(

∇·u −u · n

in Ω on ∂Ω

∀u ∈ H(Ω);

∀p ∈ H 1 (Ω)

GRADp := K∇p

(9)

and using Green’s Formula, we obtain the duality relation: GRAD = −(DIV)∗ . where DIV


∗

(10)

is the adjoint of DIV with respect to scalar products (4) and (5).

2.1. Meshes, grid functions and primary operators As mentioned above, we suppose that system (1) is posed in a bounded polyhedral(polygonal) domain Ω ⊆ Rd with a Lipschitz continuous boundary. Let Th be an unstructured mesh of Ω into nonoverlapping simply-connected polyhedra (polygons) with flat faces, where h := sup diameter(c). c∈Th

In the following we take on the element c ∈ Th the regularity assumptions listed, for instance, in [11]. Let Eh be the set of faces of the polyhedra (edges of the polygons) in Th and let Qh ⊆ Eh be the set of faces on the border of the domain Ω. We use the following notations for the mesh objects: c ∈ Th denotes a general cell in the mesh with measure |c| and centroid xc ; f ∈ Eh denotes a general face of the cell c with measure |f | and centroid xf ; nf indicates the unit normal vector to the face f with preassigned direction; αc,f = ±1 represents the mutual orientation of the vector nf and the outward normal vector to f respect to the cell c. Moreover, let Yh = Th , Eh , Qh , and let σ = c, f , then we denote with Yh (σ) the subset of Yh which elements are related with σ, and we indicate with |Yh (σ)| the cardinality of this set. The mesh objects will define the degrees of freedom of the discrete system, that is these will define the space of the discrete pressures, discrete fluxes and discrete boundary functions. Let Nc := |Th |, Nf := |Eh |, Nb := |Qh |, N ∗ := max |Eh (c)|. c

Let Ch be the set of the discrete pressures that are constant on each cell c ∈ Th : Ch := { p ∈ L2 (Ω) | p|c = const,

∀c ∈ Th } .

Given a pressure p ∈ L2 (Ω), we define the interpolant discrete pressure pI ∈ Ch with Z 1 pI|c = p dV, ∀c ∈ Th . (11) |c| c The space Fh of the discrete velocities is defined as follow. For all faces f ∈ Eh we associate a number uf and we denote with uh the vector of all uf and with Fh the vector space of all uh . 3

Let u ∈ H(Ω) a vector function and let us suppose that all face-integrals Z u · nf dS, ∀f ∈ Eh f

exist. Than the interpolant discrete flux of u in the space Fh is defined by uI := (uf )f ∈Eh with Z 1 u · nf dS, ∀f ∈ Eh . (12) uf = |f | f We denote with C h the set of the discrete pressures (with boundary values). In this case we associate a degree of freedom for every cell c ∈ Th and for every boundary face f ∈ Qh . Given a pressure p ∈ L2 (Ω) ⊕ L2 (∂Ω), we define the interpolant discrete pressure of p in C h with pI := (pc ; pf )c∈Th ,f ∈Qh where Z Z 1 1 pc = p dV, ∀c ∈ Th , and pf = p dS, ∀f ∈ Qh . (13) |c| c |f | f Finally, we define Bh the set of the discrete boundary functions that are constant on each face f ∈ Qh , i.e. Bh := { ψ ∈ L2 (∂Ω) | ψ|f = const, ∀f ∈ Qh } Given ψ ∈ L2 (∂Ω), the interpolant discrete boundary function ψ I ∈ Bh of ψ, is defined by Z 1 I ψ dS ∀f ∈ Qh . (14) ψ|f = |f | f Remark 2.1. The discrete spaces Ch , Fh , C h , Bh and the interpolation operators are defined starting from the degrees of freedom: •

1 |c|

R

•

1 |f |

R

•

1 |f |

R

c

p dV,

∀c ∈ Th ,

f

u · nf dS,

f

ψ dV,

∀p ∈ H 1 (Ω), ∀u ∈ H(Ω),

∀f ∈ Eh ,

∀ψ ∈ L2 (∂Ω).

∀f ∈ Qh ,

Remark 2.2. There are obvious correspondences: p 7→ (pc )c∈Th , with p|c = pc , Nf ∼ Fh = R u 7→ (uf )f ∈Eh ,

Ch ∼ = RN c

Ch ∼ p 7→ (pc ; pf )c∈Th ,f ∈Qh , = RNc ⊕ RNb = RNc +Nb Nb ∼ Bh = R ψ→ 7 (ψf )f ∈Qh , with ψ|f = ψf . With a slight abuse of notation we can refer to a function in the discrete functional spaces as a vector and viceversa. The definition of the mimetic scheme for the linear diffusion problem carries on with the discretization of the differential operators. Let u ∈ H(c) with c ∈ Th , then the Divergence Theorem states that Z Z ∇ · udV = u · ndS , c

∂c

where n is the unit outward normal to ∂c. Therefore, the continuous operator DIV admits the immediate discretization DIV : Fh → Ch , with (DIVuh )c =

1 |c|

X

αc,f |f |uf

f ∈Eh (c)

4

∀uh ∈ Fh .

(15)

Using similar arguments, the operator DIV has the natural discretization DIV : Fh → C h , with 1 X  αc,f |f |uf , if σ = c ∈ Th ,    |c| f ∈Eh (c) (DIVuh )σ := ∀uh ∈ Fh . (16)     −αc,f uf if σ = f ∈ Qh , and { c } = Th (f ), The operators DIV and DIV are called discrete primary operators. 2.2. Mimetic inner products and derived operators A crucial step in the construction of the MFD method is the definition of suitable inner products on the discrete functional spaces Ch , C h , Fh and Bh 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 Ch , C h , Fh and Bh : [ph , qh ]Ch := pTh MCh qh

for all ph , qh ∈ Ch ,

(17)

[ph , qh ]C h :=

pTh MC h qh

for all ph , qh ∈ C h ,

(18)

[uh , vh ]Fh :=

uhT MFh vh ψhT MBh ϕh

[ψh , ϕh ]Bh :=

for all uh , vh ∈ Fh ,

(19)

for all ψh , ϕh ∈ Bh ,

(20)

where MCh ∈ RNc ×Nc , MC h ∈ R(Nc +Nb )×(Nc +Nb ) , MFh ∈ RNf ×Nf , MBh ∈ RNb ×Nb are suitable symmetric positive definite matrices. The construction of these matrices is done locally, on each cell. The corresponding local inner products have to "mimic" the scalar products defined in (3), (4), (5). We would like that Z [ph,c , qh,c ]Ch,c =: (ph,c )T MCh,c qh,c ≈ p q dV, for all p, q ∈ H01 (c), (21) c Z Z p q dV, for all p, q ∈ H 1 (c), (22) [ph,c , qh,c ]C h,c =: (ph,c )T MC h,c qh,c ≈ pqdS + c ∂c Z [uh,c , vh,c ]Fh,c =: (uh,c )T MFh,c vh,c ≈ K−1 u · v dV, for all u, v ∈ H(c), (23) c

where, in general, by the notation rh,c we denote the vector with the degrees of freedom of r relative to the cell c. Let us start with (21). In this case the vector ph,c have a single component, i.e. ph,c = pc . Then the only possible quadrature formula for the integration in (21) is

therefore MCh,c = |c| and

[ph,c , qh,c ]Ch,c = (ph,c )T MCh,c qh,c = |c|pc qc ,

(24)

MCh := diag(|c1 |, . . . , |cNc |).

(25)

The scalar product (22) is slightly different from the previous one because of the integration along the border. Let c ∈ Th and { f1 , . . . , fk } = Qh (c), then, given ph ∈ C h , we have ph,c = pc , pf1 , . . . , pfk


T

and we set [ph,c , qh,c ]C h,c := (ph,c )T MC h,c qh,c = |c|pc qc +

k X j=1

Then we obtain

MC h,c = diag (|c|, |f1 |, . . . , |fk |) 5

|fj |pfj qfj .

(26)

and compute the global inner product on C h as: X X X [ph , qh ]C h = [ph,c , qh,c ]C h,c = |c|pc qc + |fb |pfb qfb . c∈Th

Therefore, we get MC h := where

fb ∈Qh

c∈Th



MCh 0

0 MBh



(27)

MBh := diag(|f1 |, . . . , |fNb |).

(28)

Remark 2.3. It is clear that (17) gives the exact value of the continuous scalar product for all p, q ∈ Ch . Furthermore (18) gives the exact value of the continuous inner product for all p, q ∈ C h . The definition of scalar product (23) requires a different approch. The key idea is to define suitable consistency and stability constraints in order to introduce algebraic conditions on the elements of the matrix MFh,c . These constraints can be formulated in the algebraic form as: MFh,c Nh,c = Rh,c ,

(29)

where Nh,c , Rh,c ∈ Rm×v , with m = |Eh (c)| and v number of constraints (see [19, 12]). Proposition 2.1. Let 

Nh,c and



nxf1  nxf  2 :=  .  .. nxfm

αc,f1 |f1 |(xf1 − xc ) αc,f2 |f2 |(xf2 − xc ) .. .

nyf1 nyf2 .. . nyfm

 nzf1 nzf2   ..  Kh,c , .  nzfm

αc,f1 |f1 |(yf1 − yc ) αc,f2 |f2 |(yf2 − yc ) .. .

αc,f1 |f1 |(zf1 − zc ) αc,f2 |f2 |(zf2 − zc ) .. .



    Rh,c :=     αc,fm |fm |(xfm − xc ) αc,fm |fm |(yfm − yc ) αc,fm |fm |(zfm − zc )  T T where fi ∈ Eh (c), and we set nfi = nxfi , nyfi , nzfi , xfi = (xfi , yfi , zfi ) for i = 1, . . . , m, and xc = (xc , yc , zc ) . Then the matrix T

T MFh,c := Rh,c Rh,c Nh,c

with


−1

T T T Rh,c + λh,c Im − Nh,c (Nh,c Nh,c )−1 Nh,c





−1 T  T λh,c := trace Rh,c Rh,c Nh,c Rh,c ,

(30) (31)

is symmetric and positive definite and verifies (29). Moreover there exist two positive hindependent constants C∗ and C ∗ such that T T C∗ |c|(vh,c )T vh,c ≤ (vh,c )MFh,c vh,c ≤ C ∗ |c|(vh,c )vh,c

∀vh,c ∈ Fh,c .

(32)

 −1 T T Remark 2.4. The term Rh,c Rh,c Nh,c Rh,c guarantees the consistency properties, while the   T T term λh,c Im − Nh,c (Nh,c Nh,c )−1 Nh,c gives the stability properties of the matrix MFh,c . Finally on the space of the discrete boundary functions Bh we define the global inner product: [ψh , ϕh ]Bh := ψhT MBh ϕh

for all ψh , ϕh ∈ Bh .

(33)

Remark 2.5. We can easily observe that (33) gives the exact value of correspondent continuous scalar product (6) for all ψ, ϕ ∈ Bh . 6

The last 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 Ch , C h , Fh equipped respectively with the scalar products (17), (18), (19). From continuous duality relations (8), we can introduce the discrete operator GRAD : Ch → Fh and impose the duality relation: [vh , GRADph ]Fh = −[DIVvh , ph ]Ch ⇔ vhT MFh GRADph = −vhT DIV T MCh ph , for all vh ∈ Fh , ph ∈ Ch , from which it follows that T GRAD := −M−1 Fh DIV MCh .

(34)

Similarly, from continuous relation (10), we have the discrete operator GRAD : C h → Fh defined by T [vh , GRADph ]Fh = −[DIVvh , ph ]C h ⇔ vhT MFh GRADph = −vhT DIV MC h ph , for all vh ∈ Fh , ph ∈ C h . Therefore we get (35)

GRAD = −M−1 Fh DIV MC h . T

Replacing the differential operators with the continuous operators, we can formulate the discrete mimetic system for the diffusion equation. Let us consider system (1) with homogeneous Dirichlet boundary conditions (2a). Let f I be the interpolant function of the source term f in Ch . Then the resulting linear system reads:   T uh = −GRADph , uh = M−1 Fh DIV MCh ph , ⇔ (36)   DIVuh = f I , DIVuh = f I . Similarly we can derive the discrete mimetic system for system (1) coupled with Neumann boundary conditions (2b), replacing DIV with DIV and GRAD with GRAD.

3. Spectral Properties In this section we study the spectral properties of the matrices arising from the discretization of spatial domains by the MFD approach. We consider the linear diffusion problem (1) with Dirichlet boundary conditions (2a), in particular, its discretization given by linear system (36). The discrete divergence operator DIV in (15) is clearly a linear function from Fh to Ch . In the light of Remark 2.2, we can represent DIV as a matrix in RNc ×Nf of the form DIV = M−1 Ch AF,

(37)

where F = diag(|f1 |, . . . , |fNf |) ∈ RNf ×Nf and A ∈ RNc ×Nf is defined by ( αci ,fj if fj ∈ Eh (ci ), for i = 1, . . . , Nc , and j = 1, . . . , Nf , ai,j = 0 otherwise. Thus, from (34) and (37), the discrete operator GRAD can be seen as a matrix in RNf ×Nc , that is T −1 T (38) GRAD = −M−1 Fh DIV MCh = −MFh F A and system (36) can be rewritten as  −1 uh − MFh F AT ph = 0 

M−1 Ch AF uh

=f

I



MFh

⇔ −AF 7

−F AT



0



uh ph





=

0 −MCh f

 I

 ,

(39)

the coefficient matrix of which is  MFh Rh := −AF

−F AT 0



∈ R(Nf +Nc )×(Nf +Nc ) .

(40)

The matrix Rh is symmetric and therefore it admits real eigenvalues. Moreover it is immediate e F , i.e. to observe that Rh is sparse. Indeed, MFh is the collection of local matrices M h,c MFh =

X

(41)

eF , M h,c

c∈Th

where F m e i,jh,c

=

if fi , fj ∈ Eh (c),

( F mσch,c (i),σc (j) 0

otherwise,

and σc (f ) denotes the local c-index of the global degree of freedom f ∈ Eh . Thus in each row of MFh there are at most 2N ∗ − 1 nonzero entries. Regarding the block AF , by definition, in the i-th row of the matrix A there are exactly |Eh (ci )| nonzero entries. Since F is a diagonal matrix we can conclude that the number nnz(Rh ) of nonzero entries of the matrix Rh satisfies X X nnz(Rh ) ≤ nnz(MFh ) + 2 nnz(A) ≤ (2N ∗ − 1) + 2 |Eh (c)| f ∈Eh

c∈Th

≤ 2N ∗ (Nf + Nc ) = 2N ∗ size(Rh ).

Figure 1: Sparsity of the matrix Rh for two polygonal decomposition of the domain [0, 1] × [0, 1]: on the left the classical squares mesh decomposition, on the right Voronoi mesh decomposition.

In order to explore the spectral properties of system (39), we first notice that Rh admits the following triangular factorization:     I 0 I −M−1 F AT MFh 0 F h Rh = −AF M−1 I 0 S 0 I Fh T where S := −AF M−1 Fh F A is the Schur complement of MFh . It is easy to observe that AF has full rank, and, being MFh positive definite, we can conclude that S is not singular, thus Rh is not singular also. The second property we can derive from the triangular factorization, is that Rh has spectrum in the following form

λNc +Nf (Rh ) ≤ λNc +Nf −1 (Rh ) ≤ · · · ≤ λNf +1 (Rh ) < 0 < λNf (Rh ) ≤ · · · ≤ λ2 (Rh ) ≤ λ1 (Rh ) that is Rh has Nc negative and Nf positive eigenvalues, in particular Rh is non definite. The next step is to study the boundary eigenvalues λNc +Nf (Rh ), λNf +1 (Rh ), λNf (Rh ), λ1 (Rh ), in order to estimate the condition number and the norm of the matrix Rh . For this goal we 8

have to introduce a basic assumption on the structure of the decomposition Th . (M1) Let us construct the graph with set of nodes Th and set of edges Lh ⊆ Th × Th defined by (c1 , c2 ) ∈ Lh

if and only if there exists f ∈ Eh , such that c1 , c2 ∈ Th (f ).

The mesh decomposition Th satisfies the following assumption: there exists a family of elementary paths Γ = { γi }i=1,...,Nγ such that – for every cell c ∈ Th there exists one and only one path γi ∈ Γ with c ∈ γi , i.e. the family of paths Γ defines a partition of Th , – every path γi ∈ Γ has as a start in cs,i and end in ce,i that have a face in Qh , – there exists a constant L∗ h-independent such that for i = 1, . . . , Nγ

Li ≤ L∗ h−1 , where Li is the length of the path γi .

Figure 2: Example of family of paths satisfying assumption (M1) for a Voronoi mesh decomposition. We denote with different colors the cells in different paths.

Remark 3.1. The matrix A can be seen as the incidence matrix of the graph defined by assumption (M1). Before we prove the main result on the spectral composition, we recall this result (see [8]): Lemma 3.1. Let A=



A B

BT 0

 ,

where A ∈ Rn×n is symmetric positive definite and B ∈ Rm×n has full rank. Let λ1 (A) and λn (A) denote respectively the largest and smallest eigenvalues of A, and let µ1 (B) and µm (B) denote respectively the largest and smallest singular values of B. Then the spectrum σ(A) of A satisfies σ(A) ⊆ I − ∪ I + , where I

−

   1  p p 1 2 2 2 2 := λn (A) − (λn (A)) + 4(µ1 (B)) , λ1 (A) − (λ1 (A)) + 4(µm (B)) , 2 2 9

and

  p 1 I + := λn (A), λ1 (A) + (λ1 (A))2 + 4(µ1 (B))2 . 2

Now, we can show the following result: Theorem 3.1. Let Rh be the matrix in (40). Then, under assumption (M1), the spectrum σ(Rh ) of Rh satisfies: σ(Rh ) ⊆ [−k1 hd−1 , −k2 hd ] ∪ [k3 hd , k4 hd−1 ]

(42)

where ki are positive constants h-independent. Therefore it holds: −d kR−1 , h k2 ≤ k∗ h

kRh k2 ≤ k ∗ hd−1 ,

(43)

where k ∗ = max{k1 , k4 } and k∗ = max{1/k2 , 1/k3 }. Moreover, the condition number of Rh in 2-norm is: µ2 (Rh ) ≤ k ∗ k∗ h−1 . Proof. In light of Lemma 3.1, since the matrix MFh is symmetric positive definite and −AF has full rank, we can bound the spectrum σ(Rh ) by using the previous estimates for I − and I +. Let us start with the eigenvalues of the matrix MFh . For stability conditions in (32), we have that the spectrum σ(MFh,c ) of the local matrices MFh,c satisfies σ(MFh,c ) ⊆ [C∗ |c|, C ∗ |c|]. Now, since MFh =

P

c∈Th

e F , using the Courant-Fischer formula, we have M h,c

P  P e T uT c∈Th MFh,c u uT MFh u c∈Th uc MFh,c uc = max = max λ1 (MFh ) = max u6=0 u6=0 u6=0 uT u uT u uT u P P T T c∈Th uc MFh,c uc c∈Th uc uc ∗  P P = max P ≤ 2C max |c| max P 2 T c∈Th u6=0 1 u6=0 2 Tu + c∈Th uc uc + f ∈Qh uf u u c c c∈T f ∈Q f 2 h h ≤ 2C ∗ max |c| ≤ 2C ∗ c∗ hd . c∈Th

Using similar arguments we can compute: P  P T e F ,c u T u M T c∈T u MFh u h c∈Th uc MFh,c uc λNf (MFh ) = min = min = min u6=0 u6=0 u6=0 uT u uT u uT u P P T T c∈Th uc MFh,c uc c∈Th uc uc P P P = min P ≥ C min |c| min ∗ 2 2 T T c∈Th u6=0 u6=0 c∈Th uc uc − f ∈Eh \Qh uf c∈Th uc uc − f ∈Eh \Qh uf ≥ C∗ min |c| ≥ C∗ c∗ hd . c∈Th

The next step is to bound the singular values of the matrix −AF . It is well known that (µi (−AF ))2 = λi (AF 2 AT ) for i = 1, . . . , Nc , thus we need: (µ1 (−AF ))2 = λ1 (AF 2 AT ),

(µNc (−AF ))2 = λNc (AF 2 AT ) .

For this goal, let us analyze the structure of the matrix AF 2 AT . We notice that P 2  if i = j,  P f ∈Eh (ci ) |f | k 2 T (AF A )i,j = − l=1 |fl |2 if i = 6 j, and { f1 , . . . , fk } = Eh (ci ) ∩ Eh (cj ),   0 otherwise. Indeed, ignoring the rescaling due to the multiplication by the diagonal matrix F 2 , if we multiply 10

the i-th row by the j-th row of A, by definition, we obtain respectively: • |Eh (ci )| if i = j, • −|Eh (ci ) ∩ Eh (cj )|: if the cells ci and cj have common faces f . Indeed, for each common face f , we have to compute αci ,f · αcj ,f that is necessarily −1 because the face f must have opposite orientation on two adjacent cells, • 0 otherwise. Moreover, we can note that AF 2 AT is diagonal dominant and it is strictly diagonal dominant in the rows relative to the cells with faces in the set Qh . Using the Gershgorin Circle Theorem we easily get   X X λ1 (AF 2 AT ) ≤ max  |f |2 + |f |2  ≤ 2N ∗ max |f |2 ≤ 2N ∗ f ∗ h2(d−1) . c∈Th

f ∈Eh (c)

f ∈Eh

f ∈Eh (c)\Qh (c)

The computation of the smallest singular value requires the additional assumption (M1). The idea is to use the Courant-Fischer formula for evaluating λNc (AF 2 AT ). Let us compute the entries of the vector F AT p for p ∈ RNc : X (F AT p)f = |f |αc,f pc , for f ∈ Eh , c∈Th (f )

then we have two different cases ( |f |αc,f pc T (F A p)f := |f |αc1 ,f (pc1 − pc2 )

for f ∈ Qh and { c } = Th (f ), for f ∈ Eh \ Qh and { c1 , c2 } = Th (f ).

Using the previous computation, we have λNc (AF 2 AT ) = min p6=0

(F AT p)T (F AT p) = min p6=0 pT p

P

f ∈Qh

|f |2 p2c +

P

f ∈Eh \Qh 2 c∈Th pc

|f |2 (pc1 − pc1 )2

P

,

where c, c1 and c2 are defined as above. From assumption (M1), it follows that P P 2 2 f ∈Qh pc + f ∈E \Q (pc1 − pc1 ) 2 T 2 P h 2h λNc (AF A ) ≥ min |f | min f ∈Eh p6=0 c∈Th pc  PNγ  2 P pcs,i + link of γi (pc1 ,i − pc2 ,i )2 + p2ce,i i=1 P ≥ f∗ h2(d−1) min . 2 p6=0 c∈Th pc Without loss of generality, it is possible to suppose that the cells of Th are ordered path by path, starting from pcs,i and ending with pce,i . The crucial observation is that the last row in the previous inequality is equivalent to the Courant-Fischer Formula for calculating the smallest eigenvalue of the block diagonal matrix P := f∗ h2(d−1) diag(Γ1 , Γ2 , . . . , ΓNγ −1 , ΓNγ ) ∈ RNc ×Nc where Γi = tridiag(−1, 2, −1) ∈ RLi ×Li . Now it is well known that the eigenvalues of Γi are   (Li + 1 − j)π , j = 1, . . . , Li , λj (Γi ) = 2 1 − cos Li + 1   and in particular λLi (Γi ) = 2 1 − cos Li1+1 . Therefore, we have 2(d−1)

λNc (P) = f∗ h

min

i=1,...,Nγ

2(d−1)

λLi (Γi ) = 2f∗ h 11

1 min 1 − cos Li Li + 1 

 .

From assumption (M1) on the maximum length of the paths, it is possible to conclude that     1 1 λNc (AF 2 AT ) ≥ λNc (P) = 2f∗ h2(d−1) min 1 − cos ≥ 2f∗ h2(d−1) 1 − cos ∗ −1 . Li Li + 1 L h +1 Then, asymptotically, the following relation holds λNc (AF 2 AT ) ≈ 2f∗ h2(d−1)

h2 = f∗ h2d . 2

By collecting the previous estimates we have λ1 (MFh ) ≤ 2C ∗ c∗ hd , and

λNf (MFh ) ≥ C∗ c∗ hd

(µ1 (−AF ))2 ≤ 2N ∗ f ∗ h2(d−1) ,

(µNc (−AF ))2 ≥ f∗ h2d .

To conclude the proof, it is sufficient to note that the function  p 1 λ − λ2 + 4µ2 (λ, µ) 7→ 2 is monotonically increasing in the first argument and monotonically decreasing in the second argument. Therefore, we obtain σ(Rh ) ⊆ I − ∪ I + with I− ⊆

    q q 1 C∗ c∗ hd − (C∗ c∗ hd )2 + 8N ∗ f ∗ h2(d−1) , C ∗ c∗ hd − (C ∗ c∗ hd )2 + f∗ h2d 2

and I

+

  q d ∗ ∗ d ∗ ∗ d 2 ∗ ∗ 2(d−1) ⊆ C∗ c∗ h , C c h + (C c h ) + 2N f h ,

from which immediately follows (42). Relations (43) can be easily derived from (42).

Figure 3: Asymptotic behaviour of the boundary eigenvalues as a functions of h for the sequence of Voronoi meshes in Figure 7.

Remark 3.2. Assumption (M1) is not too restrictive. In many general cases we can easily find the family of paths that satisfies the requirements described in (M1). It is important to observe 12

that the estimate of the maximum length of the paths is obtained in the worst mesh case. Indeed, if h is the mesh-size of the decomposition Th , then Nc = O(h−d ). If the number of cells having a face on the boundary is O(h−d+1 ), then we could find a constant L∗ h-independent, and a family of paths that satisfies (M1). If the number of cells having a face on the boundary is O(h−d ), then there exists a constant L∗ h-independent, and a family of paths that verifies (M1) but having maximum length less or equal to L∗ (not L∗ h−1 as set in (M1)). In this case the singular value µNc (−AF ) has order h(d−1) . Bounds (42) changes slightly, indeed I − has the form I − ⊆ [−c1 hd−1 , −c2 hd−1 ], while bounds (43) is asymptotically still valid.

Figure 4: Asymptotic behaviour of kRh k2 , kR−1 k2 , µ(Rh )2 as a functions of h for the sequence of meshes in h Figure 7.

Remark 3.3. In (31) we introduce the scalar factor λh,c for the stabilizing term in the construction of the matrix MFh . The choice in (31) is the classical choice, but we can change the parameter λh,c and still preserve the stability properties stated in Remark 2.4, in particular we can set 
−1 T  T λh,c := σ trace Rh,c Rh,c Nh,c Rh,c where σ > 0 is a constant h-independent. The new choice of the parameter λh,c does not change the asymptotic behaviour of the eigenvalues studied in Theorem 3.1, but, experimentally we can observe that a small values of σ get worse the condition number of the matrix Rh .

Figure 5: Behaviour of µ(Rh )2 as a functions of σ for the Voronoi mesh with h = 1/16 in Figure 7.

3.1. Convergence of the MFD methods Thanks to the spectral analysis of the matrix Rh we are able to study the order of convergence of the method. Indeed, as the mesh-size h becomes smaller, the condition number of Rh increases as h−1 and thus the error generated by the numerical solution of system (39) could grow up. We now recall the main results on the convergence of MFD method (see [11]): 13

Theorem 3.2. Let (u, p) be the solution of continuous system (1), and let (uh , ph ) be the discrete solution given by linear system (39). Moreover, let uI be the interpolant function of u in the set of discrete velocities Fh . Then it holds kuI − uh kFh ≤ γ1 hkpkH 2 (Ω) ,

(44)

where k·kFh is the induced norm by the inner product [·, ·]Fh and γ1 is a constant h-independent.

Theorem 3.3. Let f ∈ H 1 (Ω) and assumptions of Theorem 3.2 hold. Moreover, let pI be the interpolant of p in the set of discrete pressures Ch . Then it follows that:
kpI − ph kCh ≤ γ2 h2 kpkH 2 (Ω) + kf kH 1 (Ω) , (45) where k·kCh is the induced norm by the inner product [·, ·]Ch and γ2 is a constant h-independent.

Let us note that in the right-hand side of (39) we have the interpolant f I of the source term f in the set of the discrete pressures Ch , i.e. the function defined by Z 1 f dV, for all c ∈ Th . f|cI = |c| c In the applications we could not be able to compute the exact interpolant of the load term f I ∈ Ch . Actually an approximation feI of f I with a suitable degree of confidence, is achievable, in particular for a suitable γ ≥ 2 we can have fe|cI = f I |c + hγ τc where τc is a h-independent constant. Therefore the linear system in (39) has actually the form:      0 eh u MFh −F AT = (46) peh AF 0 −MCh feI Let (uh , ph ) be the solution of the linear system in (39), and let (e uh , peh ) be the solution of the linear system in (46). Then, if we set e h − uh , wh := u

εh := peh − ph ,

τh := (τc )c∈Th ,

we obtain that the errors wh , εh satisfy the linear system     wh 0 γ Rh =h . εh −MCh τh The next step is to estimate the errors kwh kFh and kεh kCh . For this goal, let us define the space Xh := Fh ⊕ Ch equipped with the norm k(vh , qh )k2Xh := kvh k2Fh + kqh k2Ch ,

for all vh ∈ Fh and for all qh ∈ Ch .

It is easy to see that k(vh , qh )k2Xh = vh with MXh :=



qh

MFh 0 14


T

MXh

0 MCh

 .

  vh , qh

Now, let ω := (wh , εh )T , δ := (0, MCh τh )T , thus, using the estimates of Theorem 3.1, we have −1 2 T γ −1 2γ kωk2Xh = (hγ R−1 h δ) MXh (h Rh δ) ≤ h λ1 (MXh )kRh δk2

2 2γ 2 ≤ h2γ max { λ1 (MCh ), λ1 (MFh ) } kR−1 max |c| max { 2C ∗ , 1 } kR−1 h δk2 ≤ h h δk2 c∈Th

2 2 2γ+d 2 −2d ≤ µh2γ+d kR−1 k∗ h Nc kMCh τ k2∞ h k2 kδk2 ≤ µh

≤ µh2γ−d nc h−d max |c|2 kτ k2∞ ≤ µh2γ−2d+2d kτ k2∞ c∈Th

≤ µhγ kτ k2∞ , where µ denotes a positive constant independent of h that may change at each occurrence, and it is set Nc = nc h−d . Then, we have: k(wh , εh )kXh ≤ µhγ kτh k∞ . Thus we have proved that the matrix Rh associated to a MFD method does not amplify the error due to the computation of the load term. In particular, we can prove the following theorem. Theorem 3.4. Let (e uh , pe) be the solution of linear system (46). Under assumption (M1) and assumptions of Theorem 3.2, Theorem 3.3, it holds that

and

e h kFh ≤ γ1 hkpkH 2 (Ω) + µhγ kτh k∞ kuI − u

(47)


kpI − peh kCh ≤ h2 γ2 kpkH 2 (Ω) + kf kH 1 (Ω) + µhγ kτh k∞

(48)

where τ is the vector above defined and µ, γ1 , γ2 are constants h-independent. We conclude that, if we chose γ ≥ 2, the spectral properties of the matrix Rh are consistent with the theoretical properties of the method: the order of accuracy of the solution as stated in Theorem 3.2 and Theorem 3.3 is preserved also in the solution of the linear systems. 4. MFD methods and parabolic equations This section deals with PDEs of parabolic type where a MFD method is used to discretize the space domain Ω while the ϑ-method and Runge-Kutta Chebyshev schemes will be considered for the time discretiziation. Let ΩT := Ω×(0, T ) be the cylinder with base Ω and ΣT := ∂Ω×(0, T ). Let us consider the initial-boundary value parabolic problem  in ΩT ,   pt = ∇ · K∇p + f B(p) on ΣT , (49)   p(·, 0) = p0 in Ω, where K is a full symmetric positive definite tensor with K ∈ (W 1,∞ )d×d , f ∈ L2 (ΩT ) is a source term, p0 ∈ L2 (Ω) is the initial value. For the boundary conditions we can refer to (2), ∂p ). For sake of simplicity we consider the case of homogeneous Dirichlet with B(p, u) := B(p, − ∂n boundary conditions (2a). Using the operators defined in Section 2, the PDE in (49) reads as: ( pt = DIV GRADp + f in ΩT , (50) p(·, 0) = p0 in Ω. Many numerical schemes to approximate problems of this type consist of two main steps: in the first step, often called method of lines (MOL), the PDE is discretized in the spatial variables, typically by a finite difference or finite element techniques, in order to obtain a system of ODEs while in the second step a time discretization is performed. Here we consider the MFD method 15

to approximate continuous system (50) in order to obtain the semi-discrete system:  ph,t (t) = DIV GRADph (t) + f I (t) for a.e. t ∈ (0, T ), ph (0) = pI0 ,

(51)

where ph,t denotes the derivative of ph with respect to t and f I (t) is the interpolant function of f (·, t) in the set of the discrete pressures Ch for every t ∈ (0, T ), pI0 is the interpolant function in Ch of p0 , DIV and GRAD are the discrete operators defined in Section 2. Let us recall that, with the notations of Section 3, it holds −1 T DIV GRAD = −M−1 Ch AF MFh F A .

The goal now is to investigate the spectral properties of the matrix DIV GRAD in order to discuss the stiffness nature of ODE (51). Theorem 4.1. Let DIV and GRAD be the matrices defined respectively in (37) and (38), for a mesh decomposition Th satisfying assumption (M1). Then, the spectrum σ(DIV GRAD) of the matrix DIV GRAD is such that: σ(DIV GRAD) ⊆ [−s∗ h−2 , −s∗ ],

(52)

where s∗ and s∗ are positive constant h-independent. Proof. Let us notice that the matrix MCh is non singular and it is definite, moreover the matrix −AF MFh F AT is symmetric. Therefore, we can apply the Courant-Fischer generalized formula to compute the boundary eigenvalues λNc (DIV GRAD) and λ1 (DIV GRAD). Using the Courant-Fischer generalized formula and the calculations of Theorem 3.1 we have λNc (DIV GRAD) = min p6=0

≥−

T T pT AF M−1 −pT AF M−1 Fh F A p Fh F A p = − max p6=0 pT MCh p pT MCh p

1

minc∈Th |c|

max p6=0

T pT AF M−1 Fh F A p pT p

pT AF 2 AT p 1 ≥− λ1 (AF 2 AT ) pT p C∗ (c∗ )2 h2d 1 2N ∗ f ∗ −2 ∗ ∗ 2(d−1) ≥− 2N f h = − h . C∗ (c∗ )2 h2d C∗ (c∗ )2

≥

kM−1 Fh k2 − max c∗ hd p6=0

Similarly λ1 (DIV GRAD) = max p6=0

T T −pT AF M−1 pT AF M−1 Fh F A p Fh F A p = − min T T p6=0 p MCh p p MCh p

1

T pT AF M−1 Fh F A p maxc∈Th |c| p6=0 pT p T 1 1 1 p AF 2 AT p ≤− ∗ d min ≤ − ∗ ∗ 2 2d λNc (AF 2 AT ) c h λN1 (MF ) p6=0 pT p 2C (c ) h f∗ f∗ ≤ − ∗ ∗ 2 2d h2d = − ∗ ∗ 2 . 2C (c ) h 2C (c )

≤−

min

Theorem 4.1 confirms a result that is well known when we discretize the Laplacian operator by finite difference or finite element technique, that is the stiff nature of semi-discrete system (51). In the case of mimetic schemes the stiffness ratio satisfies λNc (DIV GRAD) N ∗f ∗C ∗ ≤4 λ1 (DIV GRAD) f∗ C ∗ 16



c∗ c∗

2

h−2

Figure 6: Asymptotic behaviour of λ1 (DIV GRAD), λNc (DIV GRAD), µ(DIV GRAD)2 as a functions of h for the sequence of meshes in Figure 7.

that increases quadratically when the mesh-size h is made smaller. For stiff problems several classes of numerical integrators have been proposed and among these the class of IMEX (implicit-explicit) Runge-Kutta scheme is surely effective (see for instance [16], [17]). IMEX methods employ a pair of Runge-Kutta schemes, the first one is an implicit formula and integrates the stiff terms, the second one is an explicit formula and discretizes the non-stiff terms of the system of ODEs. In this way the stiff term is treated by a method with good stability properties, but it is necessary to solve linear systems of algebraic equations. Since the linear systems arising from the spatial discretization are of large size and ill-conditioned, it is necessary to find a compromise between the stability requirements and their computational costs. Instead, in this section, time discretizations by explicit schemes, such as Runge-Kutta Chebyshev schemes (RKC) will be compared with implicit methods as the well known ϑ-method with 0 < ϑ ≤ 1 (see, for instance, [21, 20]). 4.1. Time discretization: the ϑ-method Let us consider a partition of the time interval [0, T ] by tn = τ n, with n = 0, 1, . . . , N , and time step τ := T /N . The well known ϑ-method (0 ≤ ϑ ≤ 1) applied to semi-discrete problem (51) provides the approximate sequence:  I pn+1 = pn + τ DIV GRAD(ϑpn+1 + (1 − ϑ)pn ) + τ (ϑfn+1 − (1 − ϑ)fnI ), n = 0, . . . , N − 1, I p0 = p0 (53) with fnI := f I (·, tn ), from which I − (1 − ϑ)fnI ) for n = 0, . . . , N − 1, Rϑ pn+1 = τ (1 − ϑ)DIV GRADpn + τ (ϑfn+1

(54)

where Rϑ := (INc − τ ϑDIV GRAD). We can easily see that Rϑ is symmetric, positive definite and has the form −1 T Rϑ = INc − τ ϑM−1 Ch AF MFh F A . Thus, in order to solve linear system (54), we need the explicit form of the inverse matrix M−1 Fh , losing the sparsity of MFh . Therefore, (53) requires the solution of N full and large linear systems. Regarding the condition number of the matrix Rϑ we have µϑ (Rϑ ) :=

λ1 (Rϑ ) 1 + τ ϑs∗ h−2 ≤ ≤ µ(1 + τ ϑs∗ h−2 ) λNc (Rϑ ) 1 + τ ϑs∗

which depends on the time step τ and on the mesh-size h. Thus if τ = O(hδ ), then µϑ (Rϑ ) = O(hδ−2 ). 17

The right hand side of (54) requires the evaluation of the interpolant loading term f I (·) at tn+1 and the computation of −1 T DIV GRADpn = M−1 Ch AF MFh F A pn .

To avoid the inverse M−1 Fh we can perform the following steps: • first we compute

Yn := F AT pn ,

the matrix F is diagonal and the matrix AT is sparse with at most two non-zero entries per row; • then, we solve the linear system

MFh Xn = Yn ,

where the coefficient matrix MFh is sparse, symmetric, positive definite and well conditioned; • finally, we compute

M−1 Ch AF Xn ,

where the matrix MCh is diagonal and the matrix A is sparse with at most N ∗ non-zero entries per row. 4.2. Time discretization: RKC methods Usually, for a fixed number of stages s, the main feature of a standard Runge-Kutta formula for solving ODEs is to achieve the highest possible order. While the main feature of stabilized Runge-Kutta is to construct a formula with a region of absolute stability that is as large as possible. Runge-Kutta Chebyshev methods are stabilized methods that use explicit RungeKutta formulas to solve efficiently stiff systems of ODEs (see [24], [26]). All these formulas have order of accuracy two and stability regions that include narrow strips around the negative real axis, strips of length approximated by 0.653s2 . In order to recall these methods, consider the model initial value problem:  0 Y (t) = F (t, Y (t)), for t > 0 (55) Y (0) = Y0 where F : [0, T ] × Rn → Rn , Yn is an approximation of Y (t) at t = tn and τ is the step size at the current step. The explicit s-stage RKC formula of second-order has the following form: K0 = Yn , K1 = K0 + κ1 τ F0 , Kj = (1 − µj − νj )K0 + µj Kj−1 + νj Kj−2 + κj τ Fj−1 − aj−1 κj τ F0 ,

j = 2, . . . , s,

(56)

Yn+1 = Ks , where Fj = F (tn + cj τ, Kj ), for j = 0, 1, . . . , s, and the coefficients κj , µj , νj , aj , bj and cj , for j = 0, 1, . . . , s, are related to Chebyshev polynomials of first kind (see [16]). Let Rs (z) be the stability function of the RKC method in (56), then the stability boundary β(s), i.e. the maximum number such that the stability region R(s) := { z ∈ C | |Rs (z)| ≤ 1 } includes the interval [−β(s), 0], is given by: β(s) ≈

2 2 2 (s − 1)(1 − ε), 3 15

(57)

2 where ε = 13 is the dumping parameter (see [26, 24]). In order to preserve the stability properties of RKC methods and reduce the evaluations of the non stiff-terms a subclass of RKC methods has been introduced, the partitioned RungeKutta Chebyshev (PRKC) methods (see [27] for more details). Now, let us suppose that the

18

right hand side of (55) is in the following form: Y 0 (t) = F (t, Y (t)) + G(t, Y (t)),

for t ∈ [0, T ] ,

(58)

where F : [0, T ] × Rn → Rn , is the diffusion term and G : [0, T ] × Rn → Rn corresponds to non-stiff reaction or source terms. An arbitrary s-stage PRKC method is given by the following formula K−1 =Yn , K0 =K−1 + α0 τ G−1 K1 =K0 + κ1 τ F0 , Kj =(1 − µj − νj )K0 + µj Kj−1 + νj Kj−2 + κj τ Fj−1 − aj−1 κj τ F0 , Ks =(1 − µs − νs )K0 + µs Ks−1 + νs Ks−2 + κs τ Fs−1

j = 2, . . . , s − 1, (59)

− as−1 κs τ F0 + α1 τ G−1 + α2 τ G0 + α3 τ Gs−1 , Ks+1 =(1 − µs − νs )K0 + µs Ks−1 + νs Ks−2 + κs τ Fs−1 − as−1 κs τ F0 + α4 τ G−1 + α5 τ G0 + α6 τ Gs−1 + α7 τ Gs , Yn+1 =Ks+1 , where αi for i = 1, . . . 7, are suitable real constants, G−1 = G(tn , K−1 ), G0 = G(tn + α0 τ, K0 ), Gs−1 = G(tn + α0 τ, Ks−1 ), Gs = G(tn + τ, Ks ), and the other coefficients are the same as the RKC method (56). In [27] a possible choice of the parameters αi , for i = 1, . . . , 7, is suggested in order to achieve the order 3 for method in (59). Moreover, with this choice, method (59) has stability region R(s) that contains the rectangle [−β(s), 0] × [−1.7273, 1.7273], where β(s) is defined in (57). An important aspect of stabilized methods for advection-diffusion problem is how to choose the number of stages s and the step size τ . Let %F be the spectral radius of the Jacobian matrix JF (t, Y ) and let %G be the spectral radius of the jacobian matrix JG (t, Y ). For any given step size τ , it possible to choose the minimal stage number s satisfying the stability condition τ %G ≤ 1.7273.

τ %F ≤ β(s),

(60)

If one of these two relations is not satisfied, τ has to be reduced conveniently, or chose the minimal s such that the (60) is satisfied. We now discuss the time integration of semi-discrete problem (51) by means of PRKC methods (59). In this case, with the notation used above, we have F (t, Y ) = DIV GRADY, therefore %F =

2N ∗ f ∗ −2 h , C∗ (c∗ )2

G(t, Y (t)) = f I , %G = 0.

By using (60), the selection of the largest step size τ or the optimal number of stages s is set by 2N ∗ f ∗ −2 min { s ∈ N | βR (s) ≥ τh } , C∗ (c∗ )2 then, these may be chosen as s≈

p γτ h−2 + 1,

where γ=

τ≈

s2 − 1 2 h , γ

3N ∗ f ∗ 1 2 , C∗ (c∗ )2 1 − 15 ε

is a constant independent by h and τ . 19

(61)

The basic step simply requires the evaluation of the interpolant function f I at t = tn . For the first step we have K1 = K0 + κ1 τ DIV GRADK0 , and for the evaluation of the internal stages of an arbitrary PRKC method in (59), we have to compute Kj = (1 − µj − νj )K0 + µj Kj−1 + νj Kj−2 + κj τ DIV GRADKj−1 − aj−1 κj τ DIV GRADK0 ,

j = 2, . . . , s − 1 .

For the two final stages, since f I does not depend on Kj , the computation of Ks is not necessary, therefore we have to compute only Ks+1 = (1 − µs − νs )K0 + µj Ks−1 + νj Ks−2 + κs τ DIV GRADKs−1 + I − as−1 κs τ DIV GRADK0 + τ (α4 fnI + (α5 + α6 )f I (tn + τ α0 ) + α7 fn+1 ).

We conclude that the computation of Kj , for j = 0, . . . , s − 1, requires the solution of −1 T Fj−1 = DIV GRADKj−1 = −M−1 Ch AF MFh F A Kj−1 ,

for j = 1, . . . , s,

(62)

and the evaluation of f I at tn , tn + τ α0 and tn+1 . In (62) we may avoid the computation of M−1 Fh . Indeed, for the arguments of Section 2, the matrix MFh is sparse, symmetric, positive definite and well conditioned, therefore systems of the form M−1 Fh x = y can be easily solved by iterative methods. Then, the Fj−1 ’s can be calculated by using the method introduced in the previous subsection for the computation of the right hand side term of the ϑ-method. Each step of the PRKC method requires the solution of s linear systems of the form (62) and three evaluations of the interpolant function. Remark 4.1. Comparing standard RKC method in (56) with PRKC method in (59) we observe that both methods require the solution of s linear systems in (62), but RKC method requires s evaluations of the interpolant function, while PRKC method requires only three evaluations. We have to notice that if we use an arbitrary IMEX method, the internal stages will be given by Kj = τ ajj Fj + Gj , for j = 1, . . . , s, where Gj ’s are suitable vectors. In this case, the resulting linear systems will have the following form −1 T (INc + τ ajj M−1 for j = 1, . . . , s. Ch AF MFh F A )Kj = Gj , Let us observe that in the previous equation in order to compute the matrix −1 T M−1 Ch AF MFh F A

we have to solve Nc systems of the form MFh x = y and the resulting matrix is not sparse. Thus the linear systems above lose the structure of the matrix MFh , therefore these turn out to be very expensive. Thus, PRKC methods seem to be appropriate for parabolic problems where the spatial discretization is performed by mimetic MFD methods. The presence of the stiff terms is well treated by the stability properties of the methods. Moreover, the explicit nature of the formulas limits the computational cost of the resulting linear systems. 4.3. Discrete Conservation laws The duality property and the adjointness relations of Section 2 allow to derive conservative properties of the continuous and approximate problem (see also [23]). Let us consider system (49) with Neumann boundary conditions (2b) B(p) := K∇p(·, t) · n = ψ(·, t) 20

for t ∈ (0, T ).

Using the Divergence Theorem, we may easily prove that the solution p(x, t) of (49) satisfies the conservation law: Z Z Z d p(x, t) dx = f (x, t) dx + ψ(s, t) ds for each t ∈ (0, T ). (63) dt Ω Ω ∂Ω Let us set U(t) :=

Z Ω

p(x, t) dx,

Z

H(t) :=

Ω

f (x, t) dx +

Z

ψ(s, t) ds,

∂Ω

hence, the conservation law becomes:  d   U(t) = H(t) dt   U(0) = U0 .

U0 :=

Z Ω

p0 (x) dx,

for t ∈ (0, T ),

(64)

The solution of ODE (64) is given by U(t) = U0 + Y(t), where Y(t) :=

Z 0

t

H(σ) dσ

for t ∈ (0, T ).

Therefore the function V(t) := U(t)−Y(t) is constant with respect to t, that is V(t) ≡ V(0) = U0 is a first integral of dynamical system (64). Similar conservation law can be carried out for semi-discrete system. We have already observed that the space of the discrete pressures C h can be seen as the linear vector space C h = RNc ⊕ RNb , and every vector p ∈ C h can be written as p := (pc , pb ), where pc ∈ RNc and pb ∈ RNb . Moreover, from (16), we can also split the functional DIV : Fh → C h in the following way   DIV DIV := B where B : Fh → Bh is given by (Buh )f := −αc,f

for f ∈ Qh and c = Eh (f ).

By considering the matrix form of the functional B, we can factorize B ∈ RNb ×Nf as B := −M−1 Bh AB F, where (AB )fb ,f

( αc,fb := 0

when f = fb , with fb ∈ Qh and { c } = Th (fb ), otherwise.

If we discretize the continuous problem in (49), in the space of the discrete pressures C h with the prescribed Neumann boundary conditions, we derive the semi-discrete ODE:  c ph,t = DIV GRADph + f I , (65) −B GRADph = ψ I ,  c ph (0) = pI0 , where pch,t denotes the derivative with respect to t of the vector function pch , and f I , pI0 are defined as above and ψ I is the interpolant function of ψ in the space of discrete boundary functions Bh . Theorem 4.2. Let ph be the solution of system (65). Then the following semidiscrete conservation law holds: d c [p (t), 1]Ch = [f I (t), 1]Ch + [ψ I (t), 1]Bh dt h 21

Proof. From (65) it easy to see that d c [p (t), 1]Ch = [pch,t (t), 1]Ch = [DIV GRADph (t), 1]Ch + [f I (t), 1]Ch . dt h

(66)

Using the notation introduced above, we have that  −1   MCh 0 A DIV = F, −AB 0 M−1 Bh ∗

thus GRAD = −DIV is GRAD = −MFh F AT


−ATB .

Let 1 ∈ C h the interpolant function of the constant function f (x) ≡ 1. By computing   T A 1T GRAD = −1T F M−T Fh , −AB we obtain

(P    A c∈Th (f ) αc,f 1T = −AB α c,f − αc,f f

if f ∈ Eh \ Qh , if f ∈ Qh and { c } = Th (f ),

hence it follows that GRAD1 = 0. Thus [DIV GRADq, 1]C h = −[GRADq , GRAD1]C h = 0

for all q ∈ C h .

Some simple algebra yields 0 = [DIV GRADph , 1]C h

   DIV = GRADph , 1 B C

h

= [DIV GRADph , 1]Ch + [B GRADph , 1]Bh = [DIV GRADph , 1]Ch − [ψ I , 1]Bh , therefore

[DIV GRADph , 1]Ch = [ψ I , 1]Bh .

Thus by using (66), we get d c [p (t), 1]C = [f I (t), 1]C + [ψ I (t), 1]B , dt h

for every t ∈ (0, T ).

Now, we set Uh (t) := [pch (t), 1]Ch

Hh (t) := [f I (t), 1]Ch + [ψ I (t), 1]Bh

from Theorem 4.2 we have  d   Uh (t) = Hh (t), dt   Uh (0) = Uh,0 ,

for t ∈ (0, T ),

(67)

with Uh,0 := [pI0 , 1]Ch . Moreover from the definition of interpolant functions in the spaces Ch and Bh , and from the consistency properties of the local inner products [·, ·]Ch and [·, ·]Bh (see the Remark 2.3 and 22

2.5), it follows that: Hh (t) = [f I (t), 1]Ch + [ψ I (t), 1]Bh =

X

|c|fcI (t) +

=

c∈Th

f (x, t) dx +

c

X Z f ∈Qh

|f |ψfI (t) =

f ∈Qh

c∈Th

XZ

X

ψ(s, t) ds =

f

Z Ω

f (x, t) dx +

Z

ψ(s, t) ds = H(t).

∂Ω

Similarly Uh,0 = U0 , then ODE in (67) is actually equivalent to (64), in particular we obtain that Z c [ph (t), 1]Ch = p(x, t) dx, for all t ∈ (0, T ), Ω

and that the function V(t) is still a first integral of semi-discrete dynamical system (67). Remark 4.2. The semi-discrete conservation law stated in Theorem 4.2 can be seen as a consequence of the more general local semi-discrete conservation law. Indeed, let c ∈ Ch and let us consider the functions 1c ∈ Ch that is one in c and zero otherwise. By testing (66) with 1c ∈ Ch instead of the constant function 1, we obtain d c d [p (t), 1]Ch,c = [pch (t), 1c ]Ch = [pch,t (t), 1]Ch = [DIV GRADph (t), 1c ]Ch + [f I (t), 1c ]Ch dt h dtX (68) = αc,f |f |(GRADph (t))f + [f I (t), 1]Ch,c . f ∈Eh (c)

It easy to see that the previous formula mimic the local conservation law of the continuous system Z Z Z d p(x, t) dx = GRADp(s, t) · n ds + f (x, t) dx dt c ∂c c Z X Z GRADp(s, t) · n ds + f (x, t) dx = f ∈Eh (c)

=

X

f

c

αc,f |f |(GRADp(t))If +

Z

f (x, t) dx

c

f ∈Eh (c)

Moreover, by adding equation (68) on the cell c ∈ Th , we obtain the semi-discrete conservation law. Now, let us consider a full discrete scheme for (49) with the prescribed boundary conditions. The first time integrator method we use is the classical ϑ-method. In this case we derive the sequence { ph,n } ⊆ C h given by  c I ph,n+1 = pch,n + τ DIV GRAD(ϑpch,n+1 + (1 − ϑ)pch,n ) + τ (ϑfn+1 + (1 − ϑ)fnI ),     for n = 0, . . . , N − 1    (69) −B GRADph,n = ψnI for n = 0, . . . , N ,        c ph,0 = pI0 . Also in this case we can define a conservation law and a first integral of the system. Theorem 4.3. Let { ph,n } the sequence given by (69). Then the following discrete conservation law holds [pch,n+1 , 1]Ch = [pch,n , 1]Ch + I I τ [ϑfn+1 + (1 − ϑ)fnI , 1]Ch + τ [ϑψn+1 + (1 − ϑ)ψnI , 1]Bh

23

for n = 0, . . . , N − 1. (70)

Proof. Let τ > 0, from (69) and Theorem 4.2, we have  c  [pch,n+1 , 1]Ch − [pch,n , 1]Ch ph,n+1 − pch,n = ,1 = τ τ Ch I + (1 − ϑ)fnI , 1]Ch = = [DIV GRAD(ϑpch,n+1 + (1 − ϑ)pch,n ), 1]Ch + [ϑfn+1 I = −[B GRAD(ϑph,n+1 + (1 − ϑ)ph,n ), 1]Bh + [ϑfn+1 + (1 − ϑ)fnI , 1]Ch = I I = [ϑfn+1 + (1 − ϑ)fnI , 1]Ch + [ϑψn+1 + (1 − ϑ)ψnI , 1]Bh .

Now, if we define 1. Un := [pch,n , 1]Ch , for n = 0, . . . , N , 2. Hn,ϑ := ϑH(tn ) + (1 − ϑ)H(tn−1 ), for n = 1, . . . , N , Pn 3. Yn,ϑ := i=1 Hi,ϑ , for n = 1, . . . , N , then Therefore, if we set

Un = U0 + τ Yn,ϑ . Vn := Un − τ Yn,ϑ ,

for n = 1, . . . , N ,

(71)

the sequence { Vn } is constant, in particular Vn = U0 for n = 1, . . . , N . Finally we consider an explicit s-stage Runge-Kutta scheme defined by the Butcher tableau: c

A b

where A = (ai,j ) ∈ Rs×s , c = (ci ) ∈ Rs , b = (bi ) ∈ Rs . Applying this RK scheme to (51) we have: Ps Ps  c I ph,n+1 = pch,n + τ DIV GRAD( i=1 bi Yi ) + τ ( i=1 bi fn,i ) for n = 0, . . . , N − 1,         −B GRADph,n = ψnI for n = 0, . . . , N ,    Pi−1 Ps I Yic = pch,n + τ DIV GRAD( j=1 aij Yj ) + τ ( j=1 aij fn,j ) for i = 1, . . . , s, n = 0, . . . , N − 1,         I  −B GRADYi = ψn,i for i = 1, . . . , s, n = 0, . . . , N − 1,   c I ph,0 = p0 , (72) I where fn,i = f I (tn + τ ci ). In this case the conservation law has the following form: Theorem 4.4. Let { ph,n } the sequence given by system (72). Then the following discrete conservation law holds " s # " s # X X c c I I [ph,n+1 , 1]Ch = [ph,n , 1]Ch +τ bi fn,i , 1 +τ bi ψn,i , 1 for n = 0, . . . , N − 1. (73) i=1

i=1

Ch

24

Bh

Proof. Let τ > 0, using the usual simple algebra, we have  c  [pch,n+1 , 1]Ch − [pch,n , 1]Ch ph,n+1 − pch,n = ,1 = τ τ Ch " ! # s X = DIV GRAD bi Yi , 1 i=1

"

s X

= − B GRAD

bi Yi

,1

=

s X

#

+

I bi fn,i ,1

i=1

#

+

Ch

" s X i=1

Bh

" s X

s X

I bi fn,i ,1

#

i=1

Ch

!

i=1

"

+

"

=

Ch

#

=

I bi fn,i ,1

Ch

#

I bi ψn,i ,1

i=1

. Bh

In this case, we set 1. Un := [pch,n , 1]Ch , for i = 0, . . . , N , Ps 2. Hn,b := i=1 bi H(tn−1 + τ ci ), for n = 1, . . . , N , Ps 3. Yn,b := i=1 Hi,b , for n = 1, . . . , N , then Un = U0 + τ Yn,b . Therefore setting Vn := Un − τ Yn,b

for n = 1, . . . , N ,

(74)

the sequence { Vn } is a constant sequence, in particular Vn = U0 , for n = 1, . . . , N . Remark 4.3. If the functions f I and ψ I do not depend on t, then Yn = Y(tn ) for n = 0, . . . , N , and Un = U(tn ) , for n = 0, . . . , N + 1. Otherwise the error generated by the method satisfies |U(tn ) − Un | = O(τ γ ), where γ is the order of the chosen time integrator method. Remark 4.4. RKC methods are a class of explicit Runge-Kutta methods, hence Theorem 4.4 shows that the discrete conservation law in (73) is still satisfied for such methods. 5. Numerical Tests We implemented the fully discrete system with the implicit ϑ-method (ϑ = 0.5) and the PRKC method coupled with the MFD spatial discretization. The convergence of MFD for the parabolic problem with homogeneous boundary conditions has been evaluated in the discrete relative L2 (Ω) norm of the difference between the interpolant pI ∈ Ch of the exact solution p and the numerical solution ph at the final time T , i.e. Eh,τ :=

kpI (·, T ) − ph,N kCh , kpI (·, T )kCh

while for the parabolic problem with Neumann boundary conditions we considered E h,τ :=

kpI (·, T ) − ph,N kC h kpI (·, T )kC h

,

where pI is the interpolant of the solution in the space C h . We have considered the spatial domain Ω = [0, 1] × [0, 1] ⊆ R2 , and a general sequence of Voronoi meshes with h = 1/2, 1/4, 1/8, 1/16, and the time interval [0, 1] with τ = 25

1/5, 1/10, 1/20, 1/40. For the generation of the Voronoi meshes we used the code Polymesher [25].

Figure 7: Sequence of Voronoi meshes with h = 1/2, 1/4, 1/8, 1/16.

Example 5.1. Let us consider the parabolic equation (50) where the load term f and the initial data p0 are chosen in accordance with the anisotropic tensor K = I2 and the exact solution p(x, y, t) = e20t sin(2πx) sin(2πy),

for (x, y) ∈ [0, 1]2 , and t ∈ [0, 1].

We have considered the fully discrete system with the ϑ-method (ϑ = 0.5) and the PRKC method as time integrator schemes, coupled with the MFD discretization for the sequences of polygonal meshes in Figure 7. In Tables 1 we show the relative error Eh,τ for the ϑ-method. In Tables 2 we display the relative error for the PRKC method where we use the optimal number of stages (in accordance with (61)) shown in Table 3. Finally in Table 4 is displayed the optimal τ with number of stages s = 5 as stated in (61). Table 1: Eh,τ for MFD coupled with ϑ-method (ϑ = 0.5).

h = 1/2 h = 1/4 h = 1/8 h = 1/16

τ = 1/5

τ = 1/10

τ = 1/20

τ = 1/40

2.858924e − 00 1.068373e − 00 4.280854e − 01 2.048457e − 01

2.168115e − 00 8.504572e − 01 3.220719e − 01 1.282729e − 01

1.955095e − 00 7.753271e − 01 2.838127e − 01 1.001973e − 01

1.869308e − 00 7.438731e − 01 2.675311e − 01 8.818675e − 02

Table 2: Eh,τ for MFD coupled with PRKC method.

h = 1/2 h = 1/4 h = 1/8 h = 1/16

τ = 1/5

τ = 1/10

τ = 1/20

τ = 1/40

4.031130e − 00 5.622589e − 00 4.562020e − 00 3.424300e − 00

1.071179e − 00 1.047091e − 00 1.557681e − 01 7.652583e − 01

7.356420e − 01 6.763476e − 01 1.936138e − 01 1.218420e − 01

9.482819e − 01 2.931891e − 01 6.238054e − 02 2.440203e − 02

In this test we notice that the errors generated by the two methods are similar. In particular, for both methods, the error due to the spatial discretization and the error generated by the time discretization have the expected order of convergence i.e. O(h2 ) in h and O(t2 ) in τ . 26

Table 3: Optimal number of stages s for PRKC method.

τ = 1/5

τ = 1/10

τ = 1/20

τ = 1/40

2 3 5 10

2 2 4 7

2 2 3 6

2 2 3 5

h = 1/2 h = 1/4 h = 1/8 h = 1/16

Table 4: Optimal τ for PRKC method with number of stages s = 5.

Eh,τ τ

h = 1/2

h = 1/4

h = 1/8

h = 1/16

4.816850e − 00 9.792812e − 01

7.654074e − 01 2.114460e − 01

1.057269e − 01 5.625078e − 02

8.248277e − 02 1.488710e − 02

Example 5.2. Let us consider the parabolic problem (49) with Neumann boundary conditions where the load term f , the Neumann boundary data ϕ and the initial data p0 are chosen in accordance with the anisotropic tensor K = I2 and the exact solution p(x, y, t) = et sin(πx) sin(πy),

for (x, y) ∈ [0, 1]2 , and t ∈ [0, 1].

In this test we have evaluated the error E h,τ and, in order to check the conservation law, the parameter δh,τ := |VN − U0 | with VN defined in (71) for the ϑ-method and in (74) for the PRKC method, and U0 defined in Section 4. In Tables 5 and 6 we show the errors E h,τ and the check parameter δh,τ for the usual sequences of polygonal meshes for MFD coupled respectively with the ϑ-method and PRKC method. Table 5: E h,τ and δh,τ for ϑ-method (ϑ = 0.5).

τ = 1/5

τ = 1/10

τ = 1/20

τ = 1/40

E h,τ

h = 1/2 h = 1/4 h = 1/8 h = 1/16

5.612115e − 00 1.254160e − 00 3.897966e − 01 9.197202e − 02

5.610517e − 00 1.195797e − 00 3.595624e − 01 9.064013e − 02

5.610309e − 00 1.073417e − 00 3.598382e − 01 9.082012e − 02

5.610243e − 00 1.078023e − 00 3.595761e − 01 9.050926e − 02

δh,τ

h = 1/2 h = 1/4 h = 1/8 h = 1/16

1.221245e − 15 1.332267e − 14 4.840572e − 14 1.388889e − 13

4.773959e − 15 6.217248e − 15 1.052491e − 13 3.378408e − 13

2.253752e − 14 1.332267e − 14 2.708944e − 14 3.236300e − 13

8.326672e − 15 6.128431e − 14 1.016964e − 13 1.744160e − 13

We can observe that both methods, in accordance with our theoretical results, satisfy the conservation law, indeed for the argument explained in Section 4, we expect VN = U0 , and then δh,τ = 0 as the numerical tests confirm. In order to underline the peculiarity of the MFD method for the attainment of the conservation law, we tested the parabolic problem with prescribed Neumann boundary conditions with classical finite difference method on standard mesh on the square domain [0, 1]2 . In this case we have approximated Z pi−1,j−1 + pi,j−1 + pi,j + pi−1,j , (75) p dV ≈ h2 4 Qi,j where Qi,j is the sub-square with vertex (i − 1, j − 1), (i, j − 1), (i, j), (i − 1, j) and pi,j is the approximate value of p(xi , yj ). 27

Table 6: E h,τ and δh,τ for PRKC method.

τ = 1/5

τ = 1/10

τ = 1/20

τ = 1/40

E h,τ

h = 1/2 h = 1/4 h = 1/8 h = 1/16

4.393782e − 00 3.458312e − 00 3.884202e − 00 5.620932e − 00

4.565921e − 00 2.732892e − 00 9.104456e − 01 9.831745e − 01

3.129834e − 00 2.358713e − 00 5.229021e − 01 2.961378e − 01

3.729754e − 00 2.313776e − 00 4.873710e − 01 1.433209e − 01

δh,τ

h = 1/2 h = 1/4 h = 1/8 h = 1/16

9.436895e − 15 9.769962e − 15 5.551115e − 14 1.388889e − 13

1.221245e − 15 9.769962e − 15 5.906386e − 14 2.383648e − 13

1.543210e − 14 2.042810e − 14 4.884981e − 15 2.277067e − 13

1.221245e − 15 4.352074e − 14 7.682743e − 14 1.246780e − 13

In Table 7 we show the value of the check parameter δh,τ for the finite difference method with approximation (75) coupled with ϑ-method. Table 7: δh,τ for finite difference method coupled with ϑ-method (ϑ = 0.5).

h = 1/2 h = 1/4 h = 1/8 h = 1/16

τ = 1/5

τ = 1/10

τ = 1/20

τ = 1/40

1.598666e − 00 3.703703e − 01 1.355357e − 02 8.701903e − 03

1.598753e − 00 3.703706e − 01 1.302854e − 02 8.703313e − 03

1.598144e − 00 3.703703e − 01 1.341502e − 02 8.703744e − 03

1.598524e − 00 3.703703e − 01 1.335741e − 02 8.704179e − 03

We notice that the finite difference method generates values of δh,τ considerably greater than MFD method. Indeed the conservation law is strictly connected with the adjointness and duality properties of the discrete operators DIV and GRAD. 6. Conclusions In this paper we have analysed the structure and spectral properties of the matrices arising from the mimetic discretization of the linear diffusion problem. Under simple assumptions on the structure of the mesh we have derived the norm and the condition number of these matrices. Moreover we have discussed the parabolic problem and we have solved the semi-discrete problem using two time integrator schemes: the implicit ϑ-method, that is suitable for stiff problems but requires the solution of unstructured systems, and the explicit PRKC method, that has large stability regions and requires the solution of sparse and well conditioned systems. Finally using the duality properties of discrete mimetic operators we derived conservation laws for the numerical solution of the parabolic problem. In light of these results the authors think that the spatial discretization by MFD techniques may be useful in the context of different PDEs where the conservation of nonlinear laws is required. 7. Acknowledgements The authors wish to thank Lourenço Beirão da Veiga for several interesting discussions and suggestions which helped us to improve the quality of the paper. This work has been supported in part by the GNCS of INDAM. References [1] B. Ahmad, A. Alsaedi, F. Brezzi, L. D. Marini, and A. Russo. Equivalent projectors for virtual element methods. Comput. Math. Appl., 66(3):376–391, 2013. [2] E. D. Batista and J. E. Castillo. Mimetic schemes on non-uniform structured meshes. Electron. Trans. Numer. Anal., 34:152–162, 2009. 28

[3] L. Beirão da Veiga. A mimetic discretization method for linear elasticity. M2AN Math. Model. Numer. Anal., 44(2):231–250, 2010. [4] L. Beirão da Veiga, F. Brezzi, A. Cangiani, G. Manzini, L. D. Marini, and A. Russo. Basic principles of virtual element methods. Math. Models Methods Appl. Sci., 23(1):199–214, 2013. [5] L. Beirão da Veiga, F. Brezzi, L. D. Marini, and A. Russo. The Hitchhiker’s Guide to the Virtual Element Method. Math. Models Methods Appl. Sci., 24(8):1541–1573, 2014. [6] L. Beirão da Veiga and K. Lipnikov. A mimetic discretization of the Stokes problem with selected edge bubbles. SIAM J. Sci. Comput., 32(2):875–893, 2010. [7] L. Beirão da Veiga, K. Lipnikov, and G. Manzini. The mimetic finite difference method for elliptic problems, volume 11 of MS&A. Modeling, Simulation and Applications. Springer, 2014. [8] M. Benzi, G. H. Golub, and J. Liesen. Numerical solution of saddle point problems. Acta Numer., 14:1–137, 2005. [9] F. Brezzi and A. Buffa. Innovative mimetic discretizations for electromagnetic problems. J. Comput. Appl. Math., 234(6):1980–1987, 2010. [10] F. Brezzi, A. Buffa, and K. Lipnikov. Mimetic finite differences for elliptic problems. M2AN Math. Model. Numer. Anal., 43(2):277–295, 2009. [11] F. Brezzi, K. Lipnikov, and M. Shashkov. Convergence of the mimetic finite difference method for diffusion problems on polyhedral meshes. SIAM J. Numer. Anal., 43(5):1872– 1896 (electronic), 2005. [12] F. Brezzi, K. Lipnikov, M. Shashkov, and V. Simoncini. A new discretization methodology for diffusion problems on generalized polyhedral meshes. Comput. Methods Appl. Mech. Engrg., 196(37-40):3682–3692, 2007. [13] F. Brezzi, K. Lipnikov, and V. Simoncini. A family of mimetic finite difference methods on polygonal and polyhedral meshes. Math. Models Methods Appl. Sci., 15(10):1533–1551, 2005. [14] J. E. Castillo and R. D. Grone. A matrix analysis approach to higher-order approximations for divergence and gradients satisfying a global conservation law. SIAM J. Matrix Anal. Appl., 25(1):128–142, 2003. [15] J. E. Castillo and M. Yasuda. Linear systems arising for second-order mimetic divergence and gradient discretizations. J. Math. Model. Algorithms, 4(1):67–82, 2005. [16] W. Hundsdorfer and J. Verwer. Numerical solution of time-dependent advection-diffusionreaction equations, volume 33 of Springer Series in Computational Mathematics. SpringerVerlag, Berlin, 2003. [17] T. Koto. IMEX Runge-Kutta schemes for reaction-diffusion equations. J. Comput. Appl. Math., 215(1):182–195, 2008. [18] K. Lipnikov, G. Manzini, F. Brezzi, and A. Buffa. The mimetic finite difference method for the 3D magnetostatic field problems on polyhedral meshes. J. Comput. Phys., 230(2):305– 328, 2011. [19] K. Lipnikov, G. Manzini, and M. Shashkov. Mimetic finite difference method. J. Comput. Phys., 257:1163–1227, 2014. [20] L. Lopez. A method for the numerical solution of a class of nonlinear diffusion equations. Rocky Mountain J. Math., 21(3):1083–1097, 1991. 29

[21] A. Quarteroni and A. Valli. Numerical approximation of partial differential equations, volume 23 of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1994. [22] M. Shashkov. Conservative finite-difference methods on general grids. Symbolic and Numeric Computation Series. CRC Press, Boca Raton, FL, 1996. With 1 IBM-PC floppy disk (3.5 inch; HD). [23] M. Shashkov and S. Steinberg. Solving diffusion equations with rough coefficients in rough grids. J. Comput. Phys., 129(2):383–405, 1996. [24] B. P. Sommeijer, L. F. Shampine, and J. G. Verwer. RKC: an explicit solver for parabolic PDEs. J. Comput. Appl. Math., 88(2):315–326, 1998. [25] C. Talischi, G. H. Paulino, A. Pereira, and I. F.M . Menezes. Polymesher: a general-purpose mesh generator for polygonal elements written in matlab. Struct. Multidisc Optimiz., 45(3):309–328, 2012. [26] J. G. Verwer, B. P. Sommeijer, and W. Hundsdorfer. RKC time-stepping for advectiondiffusion-reaction problems. J. Comput. Phys., 201(1):61–79, 2004. [27] C. J. Zbinden. Partitioned Runge-Kutta-Chebyshev methods for diffusion-advectionreaction problems. SIAM J. Sci. Comput., 33(4):1707–1725, 2011.

30

E-COMPONENT CAPTIONS - ) fig1.fig: Sparsity of the matrix $\mathcal{R}_h$ for two polygonal decomposition of the domain $[0,1]\times [0,1]$: on the left the classical squares mesh decomposition, on the right Voronoi mesh decomposition. -) fig2.fig: Example of family of paths satisfying assumption (M1) for a Voronoi mesh decomposition. We denote with different colors the cells in different paths. -) fig3.fig: Asymptotic behaviour of the boundary eigenvalues as a functions of $h$ for the sequence of Voronoi meshes in Figure \ref{fig:7}. -) fig4.fig: Asymptotic behaviour of $\|\mathcal{R}_h\|_2$, $\|\mathcal{R}_h^{-1}\|_2$, $\mu(\mathcal{R}_h)_2$ as a functions of $h$ for the sequence of meshes in Figure \ref{fig:7}. -) fig5.fig: Behaviour of $\mu(\mathcal{R}_h)_2$ as a functions of $\sigma$ for the Voronoi mesh with $h=1/16$ in Figure \ref{fig:7}. -) fig6.fig: Asymptotic behaviour of $\lambda_1({\did}\,{\gd})$, $\lambda_{N_c}({\did}\, {\gd})$, $\mu({\did}\,{\gd})_2$ as a functions of $h$ for the sequence of meshes in Figure \ref{fig:7}. -) fig7.fig: Sequence of Voronoi meshes with $h=1/2, 1/4, 1/8, 1/16$.