Mimetic finite difference operators and higher order quadratures

Mimetic finite difference operators D\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{D}$$\end{document}, G\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\textbf{G}$$\end{document} are discrete analogs of the continuous divergence (div) and gradient (grad) operators. In the discrete sense, these discrete operators satisfy the same properties as those of their continuum counterparts. In particular, they satisfy a discrete extended Gauss’ divergence theorem. This paper investigates the higher-order quadratures associated with the fourth- and sixth- order mimetic finite difference operators, and show that they are indeed numerical quadratures and satisfy the divergence theorem. In addition, extensions to curvilinear coordinates are treated. Examples in one and two dimensions to illustrate numerical results are presented that confirm the validity of the theoretical findings.


Introduction
The use of the finite difference (FD) method is a common approach in numerical solutions of partial differential equations (PDEs).The FD method discretizes the partial derivatives of PDEs into a set of algebraic equations, which is then solved.Although FD methods have some known drawbacks, one advantage of FD methods is its simplicity (especially when one can fit the domain into a box-shaped geometry), a straightforward implementation when compared to other methods such as the finite element or finite volume methods.
One of the drawbacks of the FD method is the sensitivity of the solution to boundary conditions (LeVeque 2007).FD methods derive stencils for the derivative operators using a Taylor's series approach.This approach has the advantage of being straight-forward, since one can easily implement matrices for the numerical derivatives.However, the underlying physics (or other mathematical characteristics) of the problem may not be adequately represented in this discretization process.Mimetic difference methods construct difference operators divergence, D, and gradient G, which discretely satisfy the extended Gauss' divergence theorem.These methods are called mimetic because the discrete difference operators mimic the properties of their continuum counterparts.Hence, numerical schemes obtained using the mimetic operators are more faithful to the physics of the problem under investigation (Castillo and Miranda 2013).
The classical divergence theorem states that the flux of a vector field v across the sectionally smooth boundary of a compact domain in two or three dimensions equals the surface or volume integral of div ( v), and this integral can be regarded as a functional.When dealing with quadrature approximations for integrals, one would then have a functional estimate for the domain integral.This functional estimate will be said to mimic the divergence theorem when this estimate is accurate and also turns out to be equal to the quadrature over the boundary.It is known that general high-order finite differences for the divergence over a uniform grid may be accurate, but they fail to mimic the divergence theorem.Hicken & Zingg (2013) have proven accuracy of the quadratures associated with summation-by-parts methods (SBP) (Kreiss and Scherer 1974;Strand 1994).SBP methods (derived on nodal grids) have been extended to staggered grids via interpolation, O'Reilly et al. (2017), but their order of approximation at the boundary is lower than the one in the interior of the domain.The construction of these operators requires a generalized inner product, referred to as the weight matrix, whose coefficients can be used for numerical integration of functions.
It is important to notice that the mimetic difference operators are built from a staggered grid (not from a nodal one) and as a result, their processes of constructions differ from those of the SBP approach.Our scalar functions are defined at the boundary and the center of cells, while our vector functions are defined at the edges or faces of cells.Moreover, the mimetic quadratures induce diagonal norms and our mimetic operators retain the same order of accuracy over the whole domain, including at the boundary.The SBP diagonal norm does not have the same order of accuracy at the boundary, i.e., their 4 th -order operator is only 3 rd -order accurate at the boundary.Navarro (2015) and Srinivasan & Castillo (2016) investigated the use of the coefficients obtained from the diagonal weight matrices Q and P, associated with mimetic difference operators divergence D and gradient G on staggered grids (Castillo and Grone 2003) as a tool for numerical integration.This preliminary study was inspired by the resemblance of the coefficients of the mimetic quadratures with those of Newton-Cotes'.In fact, the second order is exactly the 3/8, 9/8 Newton-Cotes quadrature, as previously noticed by Castillo et al. (2001).The results of Navarro and Srinivasan & Castillo demonstrated that mimetic quadrature coefficients are a viable alternative for numerical integration.
In this paper, a theoretical framework is provided for the fourth and sixth order mimetic quadratures (Castillo et al. 2001) associated with the corresponding D and G discrete mimetic operators.The mimetic coefficients considered in this paper are the ones from (Castillo and Grone 2003), which guarantee even order of accuracy at the boundaries and interior nodes for the derivative operators.The approach used for SBP quadratures reported by Hicken & Zingg (2013) is followed.
The novelty of this research lies in the demonstration that the weights obtained from the mimetic discretization method are a valid numerical quadrature formulation, and that these quadratures satisfy the divergence theorem.In addition all weights (coefficients) are positive-valued and result in diagonal matrices.Finally, high-order mimetic quadrature formulations retain accuracy when used in curvilinear coordinate systems.
This paper is organized as follows: after stating in Sect. 2 some identities that one would like to mimic in the discrete sense, Sect. 3 introduces mimetic operators, and the resulting high-order quadratures.Section 4 demonstrates that the mimetic quadratures are bonafide numerical quadratures so they can be used for numerical integration of functions.Section 5 presents the extension to curvilinear coordinates.Section 6 provides numerical implementations of the quadratures, along with their calculated accuracy orders.This section, also includes examples of how to solve PDEs using the mimetic operators.Finally, the "Appendix" shows details of the derivation of the fourth-order divergence and gradient operators.

The generalized inner product of functions
where U = (U i ), V = (V i ), i ∈ I , and W is a symmetric and positive definite matrix.It will turn out that one can restrict W to be diagonal with positive entries.(3) 4. The first derivative operator ∂ ∂ x over a smooth function V : [a, b] → R has a discrete analog D 1 which acts on discrete vector V = (V i ), i ∈ I and approximates

Mimetic finite difference
The higher dimension equivalent to the IBP is the extended Gauss' divergence theorem, In Eq. ( 5), ∇• is the divergence operator div, and ∇ the gradient operator grad.The integral on the right hand side of Eq. ( 5) represents the boundary integral operator.
The aim of mimetic discretizations is to seek discrete equivalents for the div and grad operators.
The mimetic discretization method utilizes a staggered grid.In two or three dimensions, the divergence differential operator acts on vector fields, and the gradient differential operator acts on scalar fields.PDE flux boundary conditions require that the flux is given in terms of a gradient.So, physically meaningful PDE discretizations are constrained to compute the result of a gradient on the boundary of a voxel or element.Similarly, the definition of the divergence (and that of a curl) of a vector field (as the limit of the average flux across the boundary of a region whose volume goes to zero) imposes the condition that when discretized, the result of a divergence should be computed on the interior of a discrete voxel.
On [a, b], we define 2N , i = 0, 1, . . ., 2N and consider: • the set of cell nodes is X = {x 0 , x 1 , . . ., x N }, and • the set of centers extended (cell centers and interval boundaries a and b) is The mimetic divergence and gradient operators are defined as linear maps +2) .A derivation of the mimetic fourth-order div and grad operators can be found in the "Appendix".The image of operator D is at the cell centers X .We extend the divergence operator to include the interval boundaries or X by setting the value of this extension to zero.We denote this new operator by D : X → X .Notice that D ∈ R (N +2)×(N +1) .As a linear map, D is extended to D by appending a first row and last row of zeros.
For simplicity, we identify [a, b] with [0, 1].In one dimension, Eq. ( 5) becomes where ∂v ∂ x corresponds to the div operator, and ∂ f ∂ x to the grad operator.Hence, v is seen as a discrete vector field and f is seen as a discrete scalar field.
Equation ( 6) can be rewritten as Using mimetic operators D and G, the discrete equivalent of Eq. ( 7), for general diagonal positive-definite weight matrices Q and P is (see (1)) Rather than satisfy Eq. ( 8) exactly, our constructed operators up to an error term of order h satisfies Eq. ( 8).Specifically, we construct operators that satisfy the following equality where E is defined in Eq. ( 14).One property of these mimetic operators of k th order is the zero row sum where 1 is a vector of ones of the appropriate size.Notice that 1 T DT = 0.
Similarly, if one replaces F = 1 in (8) then where q = ( qi ) = (Q ii ) is a vector that contains the diagonal of Q and vector b m+1 = (−1, 0, . . ., 0, 1) T .Matrices P and Q are defined as the solutions of the programming problems min +1) , the mimetic boundary operator, as Notice that if e 1 = (1, 0, . . ., 0) then the first row of B is given by where the last identity follows from ( 11) and that the first row of D (by construction since there are no divergence of boundary points, see "Appendix") is zero.Similarly, the last row of B is (0, . . ., 0, 1).In addition, for 2 ≤ i ≤ N + 1, the corresponding row sums of B are where δ i j is the Kronecker delta.The last identity follows from ( 10) and ( 8). , the first and last rows of E are zero.All other row sums of E are zero.From Eq. ( 12) we see that E is order h so it goes to zero as h goes to zero.
Table 1 shows the quadrature weights for orders k = 2 s = 2, 4, 6, obtained from the mimetic method in Castillo and Grone (2017).The λ i 's are the coefficients of the quadratures.The accuracy conditions that satisfy the truncation error for polynomials up to order k + 1 are 4 Quadrature properties of the mimetic coefficients The Euler-Maclaurin summation formula (Apostol 1999) relates the integral and the numerical derivative of a function.Consider f : A direct consequence of this formula can be stated as follows: 123 Proposition 1 Consider an interval x ∈ [0, 1] and a function g = f ∈ C 2 m+2 , with 2 m + 2 ≥ 2 s.Let W be a positive definite symmetric diagonal matrix of the form W = diag(W L , I , W R ), where W L = diag(λ 0 , λ 1 , . . ., λ 2s−1 ), W R = diag(λ 2 s−1 , . . ., λ 1 , λ 0 ) and I is the identity matrix.A quadrature for the function g of the form 1, hW g is a 2 s-order accurate approximation of f 1 − f 0 if and only if the 2s-quadrature weights satisfy the (2s − 1) Bernoulli conditions, with r = 2 s, and β j is the sequence of Bernoulli numbers, Hicken and Zingg (2013)).
Proposition 3 The algebraic system of 2s-equations is completed for the 2s-weights P 1 , P 2 , . . ., P 2s with the coefficients d i j 's (as noted in the "Appendix").These d i j 's are the coefficients that arise from the Toeplitz-structure for the centered-difference stencil on a staggered grid that is of 2s-order of accuracy.
Proof It is enough to take as the remaining complementary equation any appropriate simple column sum from h PG equated to zero.However, involving the conveniently adapted entries from the upper left corner of the stencil yields the desired 2 s-order of accuracy at the boundary nodes.This corresponds to the matrix I l of Eq. ( 27) in Castillo and Grone (2003).
Proposition 4 The quadrature obtained using the diagonal weight matrix P is (2s+1)order accurate.
) .So, the quadrature obtained by using the diagonal P matrix that is associated to G is (2s + 1)-order accurate.
The mimetic divergence operator and Q hold a similar property.It is well known that general higher-order FD schemes can accurately represent the divergence on a uniform grid.However, they fail to mimic the divergence theorem in the sense that the discrete quadrature of the volume integral does not produce a discrete quadrature for the surface integral.As noted in Hicken and Zingg (2013), the use of the corresponding diagonal matrix combined with the SBP operators produces a functional that mimics the divergence theorem.As shown in the previous section, P and Q, when combined with the mimetic G and D produce a functional that mimics the extended Gauss' divergence theorem on a staggered grid.
Mimetic operator D can be generalized to higher dimensions by means of Kronecker products (Castillo and Grone 2017).By construction, 1, h QDv = v n − v 0 in [0, 1].In a two dimensional grid, the exactness condition becomes I 2 , hD xy v Q xy = B xy v, where Q xy , D xy , B xy and I 2 are obtained from the Kronecker products using the one dimensional operators.The implementation of the Kronecker products is demonstrated in example 4 of the numerical results section.Therefore, one can conclude that the divergence theorem is satisfied because the discrete quadrature of Dv over the volume produces a discrete quadrature of the boundary flux v.

Mimetic quadratures and the divergence theorem
On a two dimensional grid with m and n equal size sub-intervals, the area integral on the left hand side of Eq. ( 5) can be discretized as where The significance of Eq. ( 16) is that the 2 s-order of accuracy is retained in higher dimensions for the extended Gauss' divergence theorem, since the higher dimension operators are formed, via Kronecker products, using the one-dimensional D and G matrices.Moreover, the discretization invokes the boundary operator, and therefore depends only on the boundary values for f and v.It can therefore be concluded that the mimetic quadratures accurately mimic the extended Gauss' divergence theorem.
A special case with f a constant function in Eq. ( 16) results in the divergence theorem.The mimetic quadratures therefore satisfy the discrete version of the divergence theorem.

Higher order mimetic quadratures on curvilinear coordinates
Evaluation of integrals by a coordinate transformation require the calculation of its Jacobian.In one dimension, the Jacobian of a coordinate transformation to a uniform staggered grid X is simply d f dx = f .If the original integrand is some function z, with the corresponding discrete vector Z , one can consider an expression of the form Z , h P f as an induced quadrature on the staggered grid.One can use proposition 4 for the quadrature of the form 1, h P f .If z and f are polynomials of degree i and j, with (i + j) ≤ 2 s, then the discrete equivalent of ).This is an accurate quadrature when the integrand is of the form z(x) f (x) [18].In higher dimensions, the Jacobian is a linear combination of products of first partial derivatives.Integrals in multiple coordinates on Cartesian domains can be expressed as iterated sums starting with one variable of integration.The mimetic operators in higher dimensions can be obtained from the one dimensional operators using the Kronecker products.Thus, the accuracy of the mimetic quadratures with curvilinear coordinates is retained in higher dimensions as well (Srinivasan et al. 2022).

Numerical examples
In this section, the implementation of higher order mimetic quadratures is illustrated in the first four numerical examples.The last two examples use the mimetic D and G operators for solving PDEs.

Example 1
The goal of this example is to evaluate the integral given in Schwartz (1969), and compare the numerical approximation proposed in this paper and its analytical solution.
The numerical solution was obtained for a = 0.5.A grid refinement study was performed to calculate the convergence rate for each of the methods considered.The errors are proportional to the p th power of the grid spacing h, where p is the order of convergence, i.e., E(h) = Ch p .
C and p for each of the methods, have been fitted by least squares.The results are summarized in Table 2. Figure 1 shows the log-log errors as a function of grid spacing.The numerical results achieve the desired order of accuracy for second, fourth and sixth orders.

Example 2
Consider the integral e2 from Bailey & Borwein (2006).Figure 2 shows the log-log errors.Since the sixth order numerical solution exhibits round-off errors as the grid size is decreased, the rate of convergence was computed at every step by halving the step size (as opposed to a least-square curve fit that was performed for the previous example).Table 3 shows the computed rates of convergence for each step size.

Example 3
Consider the double integral from Hicken & Zingg (2013).The computational domain (ξ, η) is mapped to the physical coordinates (x, y) via the transformation where ξ ∈ [0, 1] and η ∈ [0, 1].For each grid size, the physical coordinates were evaluated using a nonlinear solver.The physical coordinates were then used to compute the integrand function z in I 3 .The discrete equivalent of the double integral can be represented as where the Jacobian J is given by and where • denotes the element-wise product and I is the identity matrix (Table 4).
Figure 3 shows the log-log errors as a function of step size for example 3. The convergence for the sixth order scheme shows round-off errors since it quickly attains machine numeric precision for the calculations.16).Convenient functions v and f were chosen as Figure 4 shows the calculated log-log errors as a function of step size.Both B and B 1 achieve the desired order of accuracy as expected (Table 5).

Example 5 -Brusselator
In this example the solution of a one dimension system of diffusion Eq. (Gerhard and Hairer 1996) defined by the set of PDEs is approximated by mimetic operators.The initial and boundary conditions were set to u(0, t The boundary conditions for u and v were imposed at each step of numerical integration. The ∇ • ∇ operator was discretized using the fourth order mimetic Laplacian L = DG.The Brusselator system has been solved as a system of equations (both u and v simultaneously).No special treatment was done for the nonlinear u 2 v term here.As in, both u and v were defined simultaneously on staggered grids since this is a dissipative system.
The system of equations was solved using the adaptive fourth-order Runge Kutta time discretization after a mimetic fourth-order spatial discretization of the Laplacian was utilized.The numerical solution is shown in Fig. 5.

Example 6 -2D inviscid Burgers' equation
The 2D inviscid Burgers' equation in quasilinear form (Jameson 2008) is given by , σ x = σ y = 0.5 was used as initial condition in the domain x × y ∈ [−3, 3] 2 , with periodic boundary conditions.The ∇• operator was discretized using the sixth order D, along with the corresponding sixth order interpolant to have the quantities at the cell centers.The non-linear term was treated in quasilinear form as referenced in citejameson.Interpolation is neces-sary at each time step of integration, since the div operator maps u from nodes to cell centers.With I D as the interpolant, the discretized mimetic divergence was written as D I D at each time step.The sixth-order interpolant used is from (Srinivasan et al. 2022).The procedure to discretize ∇ • u using interpolant matrices is also mentioned in that report.
Figure 6 shows the solution obtained from a fourth-order Runge Kutta scheme with sixth-order mimetic spatial discretization.

Conclusion
This paper presents an introduction and a theoretical framework for higher-order quadratures arising from the mimetic methods.The mimetic quadrature weights are positive definite diagonal matrices that satisfy the divergence theorem.The weights satisfy also the Bernoulli conditions of the Euler-Maclaurin series, and can therefore be used as end corrections on a grid for numerical integration of functions.Mimetic methods satisfy vector calculus identities as well as a discrete version of the extended Gauss' divergence theorem.The mimetic D and G operators also exhibit uniform order of accuracy at all (interior and the boundary) grid points.In addition, the accuracy of the higher-order mimetic quadratures is preserved under curvilinear coordinates.
Funding No funding was received for conducting this study.Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material.If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

A Appendix: Mimetic fourth-order D and G
Consider a staggered grid on [0, 1] of N cells and grid spacing h = 1 N .Define

cell centers and boundary).
The mimetic divergence operator D : X → X uses data on X to compute derivatives at X .A staggered fourth-order finite difference scheme is given by [ 1 24 , − 9 8 , 9 8 , − 1 24 ].This scheme suffices for points x l ∈ X such x l−3/2 , x l−1/2 , x l+1/2 , x l+3/2 ∈ X .That is not the case for x 1/2 , x N −1/2 , where a special stencil is needed for the boundaries.Symmetry suggests that it is necessary to find a fourth-order scheme for x 1/2 .However, the approach taken in Castillo and Grone (2017) is to consider special schemes for Xb = {x 1/2 , x 3/2 , x 5/2 , x 7/2 } utilizing data from X b = {x 0 , x 1 , x 2 , x 3 , x 4 , x 5 } (similarly at the other boundary).Fourth-order Taylor expansions of the points in Xb to compute the first-order derivative of the points in X b generates a Vandermonde system of six unknowns (scheme weights) and five equations (degree of accuracy), with right hand side vector c = (0, 1, 0, 0, 0) T (since one wants to approximate the first derivative).Since the Vandermonde matrix V is full rank, then nullity(V ) = 1.
The resulting fourth-order divergence matrix D ∈ R (n,n+1) takes the form and it can be shown that the structure of V is such that its null space generator is (−1, 5, 10, 10, −5, 1).Therefore, for each grid point in Xb , there is one degree of freedom, and hence D has four degrees of freedom.
Let D be the top-left corner sub-matrix of D, then the general solution can be represented as If one takes α 2 = α 3 = α 4 = 0, one can extend the standard fourth-order scheme as much as possible.On the other hand, α 1 is chosen conveniently.
If one extends the image of D from X to X by padding with zeros the first and last rows then D ∈ R (n+2,n+1) is given by  The derivation for the coefficients of the fourth-order gradient follows a similar procedure as the one outlined above to obtain G ∈ R (m+1,m+2) as .
The matrix E is predominantly zeros with non-zero elements corresponding to the boundaries.E is O(h) and thus E → 0 as h → 0. This is illustrated in numerical example 4 of Sect.7.
3. The integration by parts (IBP) formula for U, V : [a, b] → R is given by
The authors have no relevant financial interests to disclose.Author José E. Castillo is an editor for Springer Nature collection Mathematical Modelling, Mimetic Methods and Simulations for Reactive Flows and Mechanics in Porous Media.The other authors have no relevant non-financial interests to disclose.

.
Symmetry of D ensures that the bottom right corner can be determined from the coefficients of the top left portion of the matrix.The weight matrix Q associated to D is given by

Table 1
Quadrature weights P and Q arising from G and D, respectively

Table 3
Calculated order of accuracy for example 2