1 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, \(\textbf{D}\), and gradient \(\textbf{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 \(\vec v\) across the sectionally smooth boundary of a compact domain in two or three dimensions equals the surface or volume integral of \(div(\vec 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 \(\textbf{D}\) and gradient \(\textbf{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 \(\textbf{D}\) and \(\textbf{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.

2 Identities

Let [a, b] be an interval. Let \(\{ x_0, x_1, \cdots , x_N \}\) a homogeneous partition of [a, b], i.e., \(x_i = a + ih, \, i \in I, \, h = \frac{b-a}{N-1}, I = \{ 0, 1, \cdots , N \}\).

  1. 1.

    A quadrature for the numerical integration of function \(\mathcal {U}:[a,b] \rightarrow {\mathbb {R}}\) is given by the closed Newton-Cotes formula,

    $$\begin{aligned} \int _{a}^{b} \mathcal {U} \ \textrm{d}x \approx \sum _{i=0}^{N} {{\tilde{w}}}_i \,U_i = \big<U, \mathbbm {1}\big>_{{{\tilde{W}}}} = h \, \big <U,\mathbbm {1}\big >_{{W}} = h \, {U}^T W \mathbbm {1} , \end{aligned}$$
    (1)

    where \({{\bar{w}}}_i, \, i \in I\), are its integration weights, and discrete vector \(U = (U_i) = (U(x_i)), i \in I\). In addition, \(w_i = \frac{{{\bar{w}}}_i}{h}, \, i \in I\), and \({\bar{W}} = \text {diag } \{ {{\bar{w}}}_i \} = h \, \text {diag } \{ w_i \} = h \, W\), and \(\mathbbm {1} = [1,1,...,1]^T\). Notice that \(\big <\cdot ,\cdot \big >_{V}\) for \(V = {\bar{W}}, W\) refers to the generalized inner product.

  2. 2.

    The generalized inner product of functions \(\mathcal {U}, \mathcal {V}: [a,b] \rightarrow {\mathbb {R}}\) is given by

    $$\begin{aligned} \int _{a}^{b} \mathcal {U} \mathcal {V}\ \textrm{d}x \approx {\big <{U},{V} \big >_{{{\mathbb {W}}}} = U^T {\mathbb {W}} V = V^T {\mathbb {W}} U}, \end{aligned}$$
    (2)

    where \(U = (U_i), \, V = (V_i), \, i \in I\), and \({\mathbb {W}}\) is a symmetric and positive definite matrix. It will turn out that one can restrict \({\mathbb {W}}\) to be diagonal with positive entries.

  3. 3.

    The integration by parts (IBP) formula for \(\mathcal {U}, \mathcal {V}: [a,b] \rightarrow {\mathbb {R}}\) is given by

    $$\begin{aligned} \int _{a}^{b} \mathcal {U} \mathcal {V}_x \ \textrm{d}x = \mathcal{U}\mathcal{V}\bigg |_{a}^{b} - \int _{a}^{b} \mathcal {V} \mathcal {U}_x \ \textrm{d}x. \end{aligned}$$
    (3)
  4. 4.

    The first derivative operator \(\frac{\partial }{\partial x}\) over a smooth function \(\mathcal {V}:[a,b] \rightarrow {\mathbb {R}}\) has a discrete analog \(\mathcal {D}_1\) which acts on discrete vector \(V = (V_i), \, i \in I\) and approximates \(\frac{\partial \mathcal {V}}{\partial x} = \mathcal {V}_x\) by

    $$\begin{aligned} \mathcal {V}_x = \frac{\partial \mathcal {V}}{\partial x} \approx \mathcal {D}_1 \, V. \end{aligned}$$
    (4)

3 Mimetic finite difference

The higher dimension equivalent to the IBP is the extended Gauss’ divergence theorem,

$$\begin{aligned} \int _{\Omega } f \ \nabla \cdot \vec {v}\ \ \textrm{d}V + \int _{\Omega } \vec {v} \cdot \left( \nabla f \right) \ \textrm{d}V = \int _{\partial {\Omega }} f \vec {v} \cdot \vec {n} \ \textrm{d}S \end{aligned}$$
(5)

In Eq. (5), \(\nabla \cdot \) is the divergence operator div, and \(\nabla \) 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 \(x_i = a + \frac{i(b-a)}{2N}, \, i = 0,1,\ldots ,2N\) and consider:

  • the set of cell centers is \({\tilde{X}} = \{ x_{1/2},x_{3/2},\ldots x_{N-1/2} \}\),

  • the set of cell nodes is \(X = \{ x_0,x_1,\ldots ,x_N \}\), and

  • the set of centers extended (cell centers and interval boundaries a and b) is \({\hat{X}} = \{ x_0,x_{1/2},x_{3/2},\ldots ,x_{N-1/2},x_N \}\).

The mimetic divergence and gradient operators are defined as linear maps

$$\begin{aligned} \textbf{D}:X \rightarrow {\tilde{X}}, \qquad \textbf{G}: {\hat{X}} \rightarrow X. \end{aligned}$$

Therefore, \(\textbf{D} \in {{\mathbb {R}}}^{N \times (N+1)}\) and \(\textbf{G} \in {{\mathbb {R}}}^{(N+1) \times (N+2)}\). A derivation of the mimetic fourth-order div and grad operators can be found in the “Appendix”.

The image of operator \(\textbf{D}\) is at the cell centers \(\tilde{X}\). We extend the divergence operator to include the interval boundaries or \({\hat{X}}\) by setting the value of this extension to zero. We denote this new operator by \(\hat{\textbf{D}}:X \rightarrow {\hat{X}}\). Notice that \(\hat{\textbf{D}} \in {{\mathbb {R}}}^{(N+2) \times (N+1)}\). As a linear map, \(\textbf{D}\) is extended to \(\hat{\textbf{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

$$\begin{aligned} \int _{0}^{1} f \frac{\partial v}{\partial x} \ \textrm{d}x + \int _{0}^{1} \frac{\partial f}{\partial x} v \ \textrm{d}x = v(1)f(1) - v(0)f(0), \end{aligned}$$
(6)

where \(\dfrac{\partial v}{\partial x} \) corresponds to the div operator, and \(\dfrac{\partial f}{\partial 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

$$\begin{aligned} \big<\dfrac{\partial v}{\partial x},f \big> + \big <v, \frac{\partial f}{\partial x} \big > = v(1)f(1) - v(0)f(0). \end{aligned}$$
(7)

Let \(V = v\vert _X = (V_i) = V(x_i), \, x_i \in X\) be the discrete version of the vector field v (V is v restricted to X), and let \(F = f\vert _{{\hat{X}}} = (F_i) = F(x_i), \, x_i \in {\hat{X}}\) be the discrete version of the scalar field f (F is f restricted to \({\hat{X}}\)).

Using mimetic operators \(\hat{\textbf{D}}\) and \(\textbf{G}\), the discrete equivalent of Eq. (7), for general diagonal positive-definite weight matrices Q and P is (see (1))

$$\begin{aligned} h \big<\hat{\textbf{D}}V,F \big>_Q + h \big <V, \textbf{G} F \big >_P = h V^T \hat{\textbf{D}}^T Q F + h V^T P \textbf{G} F = V_N F_N - V_0 F_0. \end{aligned}$$
(8)

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

$$\begin{aligned} h V^T \hat{\textbf{D}} Q \hat{F} + h V^T P \textbf{G} \hat{F} = h V^T \textbf{B}^T \hat{F} = V_N F_N - V_1 F_1 + h V^T E \hat{F}. \end{aligned}$$
(9)

where E is defined in Eq. (14).

One property of these mimetic operators of \(k^{th}\) order is the zero row sum

$$\begin{aligned} \hat{\textbf{D}} \, \mathbbm {1} = 0, \qquad \textbf{G} \, \mathbbm {1} =0, \end{aligned}$$
(10)

where \(\mathbbm {1}\) is a vector of ones of the appropriate size. Notice that \({\mathbbm {1}}^T \hat{\textbf{D}}^T = 0\).

For \(V = {\mathbbm {1}}\), (8) becomes \(h {\mathbbm {1}}^T P \textbf{G} F = F_N - F_0 = (-1,0,\ldots ,0,1)^T F\), for any F, or equivalently, \(h \textbf{G}^T P {\mathbbm {1}} = b_{m+2} = (-1,0,\ldots ,0,1)\). Therefore,

$$\begin{aligned} h \, {\textbf{G}}^T p = b_{m+2}, \end{aligned}$$
(11)

where \(p = (p_i) = (P_{ii})\) is a vector that contains the diagonal of P.

Similarly, if one replaces \(F = {\mathbbm {1}}\) in (8) then

$$\begin{aligned} h \, {\hat{\textbf{D}}}^T {\hat{q}} = b_{m+1}, \end{aligned}$$
(12)

where \({\hat{q}} = ({{\hat{q}}}_i) = (Q_{ii})\) is a vector that contains the diagonal of Q and vector \(b_{m+1} = (-1,0,\ldots ,0,1)^T\).

Matrices P and Q are defined as the solutions of the programming problems

$$\begin{aligned} \min \quad \sum _{i=1}^{N+1} p_i, \quad \text {subject to } \quad h \, {\textbf{G}}^T p = b_{m+2}, \quad p_i \ge 0, \end{aligned}$$

and

$$\begin{aligned} \min \quad \sum _{i=1}^N q_i, \quad \text {subject to } \quad h \, {\textbf{D}}^T q = b_{m+1}, \quad q_i \ge 0, \end{aligned}$$

respectively.

From (8), define \(\textbf{B} \in {{\mathbb {R}}}^{(N+2) \times (N+1)}\), the mimetic boundary operator, as

$$\begin{aligned} \textbf{B} = h \, (Q \hat{\textbf{D}} + \textbf{G}^T P). \end{aligned}$$
(13)

Notice that if \(e_1 =(1,0,\ldots ,0)\) then the first row of \(\textbf{B}\) is given by

$$\begin{aligned} e_1^T \textbf{B} = h \, e_1^T \textbf{G}^T P + h \, e_1^T Q \hat{\textbf{D}} = h \, {G}_{1\ldots }^T P + h \, Q_{1 \ldots } \hat{\textbf{D}} = (-1,0,\ldots ,0), \end{aligned}$$

where the last identity follows from (11) and that the first row of \(\hat{\textbf{D}}\) (by construction since there are no divergence of boundary points, see “Appendix”) is zero. Similarly, the last row of \(\textbf{B}\) is \((0,\ldots ,0,1)\).

In addition, for \(2 \le i \le N+1\), the corresponding row sums of \(\textbf{B}\) are

$$\begin{aligned} \sum _{j=1}^{N+1} {\textbf{B}}_{ij}= & {} h \, \sum _{j=1}^{N+1} (\textbf{G}^T P + Q \hat{\textbf{D}})_{ij} = h \, \sum _{j=1}^{N+1} \sum _{l=1}^{N+2} ({\textbf{G}}_{il}^T P_{lj} + Q_{il} {\hat{\textbf{D}}}_{lj}) \\= & {} h \, \sum _{j=1}^{N+1} \sum _{l=1}^{N+2} ({\textbf{G}}_{il}^T \, p_l \, \delta _{lj} + q_i \, \delta _{il} \, {\hat{\textbf{D}}}_{lj}) = h \, \sum _{j=1}^{N+1} {\textbf{G}}_{ij}^T p_j + h \, q_i \sum _{j=1}^{N+1}{\hat{\textbf{D}}}_{ij} = 0, \end{aligned}$$

where \( \delta _{ij}\) is the Kronecker delta. The last identity follows from (10) and (8).

If one defines \({\textbf{B}}_1 \in {{\mathbb {R}}}^{(N+2) \times (N+1)}\) by

$$\begin{aligned} {\textbf{B}}_1 = \left( \begin{array}{ccccc} -1 &{} 0 &{} \ldots &{} 0 &{} 0 \\ 0 &{} 0 &{} \ldots &{} 0 &{} 0 \\ \vdots &{} \vdots &{} \vdots &{} \vdots &{} \vdots \\ 0 &{} 0 &{} \ldots &{} 0 &{} 0 \\ 0 &{} 0 &{} \ldots &{} 0 &{} 1 \end{array} \right) , \end{aligned}$$

and

$$\begin{aligned} \textbf{E} = \textbf{B} - {\textbf{B}}_1, \end{aligned}$$
(14)

the first and last rows of \(\textbf{E}\) are zero. All other row sums of \(\textbf{E}\) are zero. From Eq. (12) we see that \(\textbf{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 \(\lambda _i\)’s are the coefficients of the quadratures. The accuracy conditions that satisfy the truncation error for polynomials up to order \(k+1\) are

$$\begin{aligned} \hat{\textbf{D}} v^{j} = j v^{j-1} \quad \text {and } \quad \textbf{G} \hat{f}^{j} = j \hat{f}^{j-1}, \quad j = 0, 1, \dots , k \end{aligned}$$
(15)
Table 1 Quadrature weights P and Q arising from \(\textbf{G}\) and \(\textbf{D}\), respectively
Full size table

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:[0,1] \rightarrow {\mathbb {R}}, \, f \in C^{2m+2}\). The Euler-Maclaurin summation formula is

$$\begin{aligned} \int _{0}^{1} f(x) \ \textrm{d}x = h \sum _{v=0}^{n} f_v + \sum _{k=1}^{2m} \frac{\beta _k}{k!} h^k \Big (f_0^{(k-1)} - (-1)^k f_n^{(k-1)} \Big ) + \mathcal {O} (h^{2m+2}) \end{aligned}$$

A direct consequence of this formula can be stated as follows:

Proposition 1

Consider an interval \(x \in [0,1]\) and a function \(g = f' \in C^{2\,m+2}\), with \(2\,m+2 \ge 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(\lambda _0, \lambda _1, \dots , \lambda _{2s-1})\), \( W_R = diag(\lambda _{2\,s-1}, \dots , \lambda _1, \lambda _0)\) and I is the identity matrix. A quadrature for the function g of the form \(\big < \mathbbm {1}, hWg \big>\) 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,

$$\begin{aligned} j \sum _{v=0}^{r-1} \lambda _v (r-v)^{j-1} = r^j - (-1)^j \beta _j, \quad j = 1,2,\dots ,(2s-1), \end{aligned}$$

with \(r = 2\,s\), and \(\beta _j\) is the sequence of Bernoulli numbers, \(\beta _1 = -\frac{1}{2}, \beta _2 = \frac{1}{6}, \beta _4 = -\frac{1}{30}, \beta _6 = \frac{1}{42}, \beta _3 = \beta _5 = \dots =0\) (refer to theorem 1 in Hicken and Zingg (2013)).

Proposition 2

The \(P_i\)’s satisfy the \((2s-1)\) Bernoulli conditions.

Proof

By construction, \(\big <\mathbbm {1}, h \textbf{G}f \big >_P = f_1 - f_0\) is satisfied. Hence, the weights of P satisfy the \((2s-1)\) Bernoulli conditions. The exactness condition corresponds to Eqs. (27) and (29) in Castillo and Grone (2003). \(\square \)

Proposition 3

The algebraic system of 2s-equations is completed for the 2s-weights \(P_1, P_2, \dots , P_{2s}\) with the coefficients \(d_{ij}\)’s (as noted in the “Appendix”). These \(d_{ij}\)’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 P\textbf{G}\) 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). \(\square \)

Proposition 4

The quadrature obtained using the diagonal weight matrix P is \((2s+1)\)-order accurate.

Proof

By construction \(\big < \mathbbm {1}, hP \textbf{G}f \big > = f_1 - f_0\) holds. Since \(\textbf{G}\) is 2s-order accurate, it follows that \(\textbf{G} f = f' + O(h^{2s})\). Thus, \(hP \textbf{G}f = hPf' + hP[O(h^{2\,s})]\). It follows that \(f_1 - f_0 = \big< \mathbbm {1}, hP \textbf{G}f \big> = \big< \mathbbm {1}, hPf' \big> + \big < \mathbbm {1}, O(h^{2\,s+1}) \big > \). So, the quadrature obtained by using the diagonal P matrix that is associated to \(\textbf{G}\) is \((2s+1)\)-order accurate. \(\square \)

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 \(\textbf{G}\) and \(\textbf{D}\) produce a functional that mimics the extended Gauss’ divergence theorem on a staggered grid.

Mimetic operator \(\textbf{D}\) can be generalized to higher dimensions by means of Kronecker products (Castillo and Grone 2017). By construction, \(\big < \mathbbm {1}, hQ \textbf{D} v \big > = v_n -v_0\) in [0, 1]. In a two dimensional grid, the exactness condition becomes \(\big < I_2, h\textbf{D}_{\textbf{xy}} v \big > _{Q_{xy}} = \mathbf {B_{xy}}v\), where \(Q_{xy}, \textbf{D}_{\textbf{xy}}, \textbf{B}_{\textbf{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 \(\textbf{D}v\) over the volume produces a discrete quadrature of the boundary flux v.

5 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

$$\begin{aligned} \big<h {{\hat{{\textbf{D}}}}_{xy}} v,\hat{f} \big>_{Q_{xy}} + \big <h\textbf{G}^T_{\textbf{xy}}v, \hat{f} \big >_{P_{xy}} = h \hat{f}^T [Q_{xy} {{\hat{\textbf{D}}}_{xy}} + \textbf{G}^T_{\textbf{xy}} P_{xy}] v = h\hat{f}^T {\textbf{B}_{xy}} v, \end{aligned}$$
(16)

where \({\hat{\textbf{D}}_{xy}} = \big [ {{\hat{\textbf{I}}_n}} \otimes {{\hat{\textbf{D}}_m}}, {{\hat{\textbf{D}}_n}} \otimes {\hat{\textbf{I}}}_{m} \big ]\), \(\textbf{G}_{\textbf{xy}} = \big [ {{\hat{\textbf{I}}_n}}^T \otimes \mathbf {G_m}; \mathbf {G_n} \otimes {{\hat{\textbf{I}}_m}}^T \big ]\), \(\mathbf {B_{xy}} = \big [ {{\hat{\textbf{I}}_n}} \otimes \mathbf {B_m}, \mathbf {B_n} \otimes {{\hat{\textbf{I}}_m}} \big ]\), \(P_{xy} = diag([\mathbf {I_n} \otimes P_{m+1}, P_{n+1} \otimes \mathbf {I_m}])\), \(Q_{xy} = diag([\mathbf {I_{n+2}} \otimes Q_{m+2}, Q_{n+2} \otimes \mathbf {I_{m+2}}])\), \({\hat{\textbf{I}}_n} = [0; \mathbf {{I}_n}; 0]\).

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 \(\textbf{D}\) and \(\textbf{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.

6 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 \({\hat{X}}\) is simply \(\dfrac{\textrm{d}f}{\textrm{d}x} = f^{'}\). If the original integrand is some function z, with the corresponding discrete vector Z, one can consider an expression of the form \(\big < Z, hPf' \big>\) as an induced quadrature on the staggered grid. One can use proposition 4 for the quadrature of the form \(\big < \mathbbm {1}, hPf' \big>\). If z and f are polynomials of degree i and j, with \((i+j) \le 2\,s\), then the discrete equivalent of \( \int _{0}^{1} z(x) f'(x) \ \textrm{d}x\) can be expressed as \(\big < Z, h\textbf{D}f \big >_Q = Z^T hQ \textbf{D}f\). It now follows from the construction of \(\textbf{D}\) that if \(z(x)f' \in C^{2\,s} [0,1]\), then \(Z^T hQ \textbf{D}f = \int _{0}^{1} z(x) f'(x) \ \textrm{d}x + O(h^{2\,s})\). 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).

7 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 \(\textbf{D}\) and \(\textbf{G}\) operators for solving PDEs.

7.1 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.

$$\begin{aligned} \mathcal {I}_1 = \int _{0}^{1} \frac{1}{(x-0.5)^2+a^2} \ \textrm{d}x = \frac{1}{a}\Bigg (atan\Big (\frac{0.5}{a}\Big )-atan\Big (\frac{-0.5}{a}\Big )\Bigg ) \end{aligned}$$
Table 2 Calculated least-square order of accuracy for example 1
Full size table
Fig. 1
figure 1

Log-log error versus step-size for example 1

Full size image

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.

7.2 Example 2

Consider the integral \({\textbf {e2}}\) from Bailey & Borwein (2006).

$$\begin{aligned} \mathcal {I}_2 = \int _{0}^{1} x^2 arctan(x) \ \textrm{d}x = \frac{\pi - 2 + 2 ln2}{12} \end{aligned}$$

This integral involves a coordinate transformation \(\xi = arctan(x), \xi \in [0,\frac{\pi }{4}]\). The numerical integration uses \(\mathcal {I}_{n2} = z' hP\textbf{G}x\), where z is the integrand in \(\mathcal {I}_2\). \(\textbf{G}x\) refers to the inverse of the Jacobian, which is the first derivative of the transformation \(\dfrac{\textrm{d}x}{\textrm{d}\xi }\).

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.

Fig. 2
figure 2

Log-log error versus step-size for example 2

Full size image
Table 3 Calculated order of accuracy for example 2
Full size table

7.3 Example 3

Consider the double integral

$$\begin{aligned} \mathcal {I}_3 = \int \int _{\Omega _x} (x^2 + y^2) e^{\frac{1}{3}(1-x^2+y^2)} sin \Big ( {\frac{xy-1}{2}} \Big ) \ \textrm{d}x \textrm{d}y = 3 (1-e^{-1})(1-cos(1)) \end{aligned}$$

from Hicken & Zingg (2013). The computational domain \((\xi , \eta )\) is mapped to the physical coordinates (x, y) via the transformation

$$\begin{aligned} \xi = \frac{x^2-y^2-1}{3}, \quad \text {and} \quad \eta = \frac{xy-1}{2}, \end{aligned}$$

where \(\xi \in [0,1]\) and \(\eta \in [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 \(\mathcal {I}_3\). The discrete equivalent of the double integral can be represented as

$$\begin{aligned} \mathcal {I}_{n3} = z^T (P \otimes P) J, \end{aligned}$$

where the Jacobian J is given by

$$\begin{aligned} J = [(\textbf{I} \otimes \textbf{G})x] \circ [(\textbf{G} \otimes \textbf{I})y] - [(\textbf{I} \otimes \textbf{G})y] \circ [(\textbf{G} \otimes \textbf{I})x], \end{aligned}$$

and where \(\circ \) denotes the element-wise product and \(\textbf{I}\) is the identity matrix (Table 4).

Table 4 Calculated order of accuracy for example 3
Full size table
Fig. 3
figure 3

Log-log error versus step-size for example 3

Full size image

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.

7.4 Example 4

This example demonstrates the numerical convergence of the boundary operator \(\textbf{B}\) and \({\textbf{B}}_1\) in 2D as shown in Eq. (16). Convenient functions v and f were chosen as

$$\begin{aligned} v(x,y) = sin(x)cos(y), \quad f(x,y) = 2e^x y^2, \quad x \in [0,1], y \in [0,1] \end{aligned}$$
(17)

Figure 4 shows the calculated log-log errors as a function of step size. Both \(\textbf{B}\) and \({\textbf{B}}_1\) achieve the desired order of accuracy as expected (Table 5).

Fig. 4
figure 4

Log-log error versus step-size for example 4

Full size image
Table 5 Calculated order of accuracy for example 4
Full size table
Fig. 5
figure 5

Example 5, solution u(x, t).

Full size image

7.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

$$\begin{aligned} u_t= & {} 1 + u^2v - 4u + 1/50 \ (\nabla \cdot \nabla u) \\ v_t= & {} 3u - u^2v + 1/50 \ (\nabla \cdot \nabla v) \end{aligned}$$

is approximated by mimetic operators.

The initial and boundary conditions were set to \(u(0,t) = u(1,t) = 1\), \(v(0,t) = v(1,t) = 3\), \(u(x,0) = 1+sin(2\pi x)\), \(v(x,0) = 3\) in the domain \(x \in [0,1]\) and \(t \in [0,10]\). The boundary conditions for u and v were imposed at each step of numerical integration.

The \(\nabla \cdot \nabla \) operator was discretized using the fourth order mimetic Laplacian \(\textbf{L} = \hat{\textbf{D}} \textbf{G}\). 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^2v\) term here. As in, both u and v were defined simultaneously on staggered grids since this is a dissipative system.

Fig. 6
figure 6

2D inviscid Burgers’ equation, solution u(x, y, t).

Full size image

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.

7.6 Example 6 - 2D inviscid Burgers’ equation

The 2D inviscid Burgers’ equation in quasilinear form (Jameson 2008) is given by

$$\begin{aligned} u_t + \frac{2}{3} \nabla \cdot \left( \frac{u^2}{2} \right) + \frac{1}{3} u \nabla \cdot u = 0 \end{aligned}$$

A 2D Gaussian input \(u(x,y,0) = exp \left( -\dfrac{x^2}{2 \sigma ^2_x} - \dfrac{y^2}{2 \sigma ^2_y } \right) \), \(\sigma _x = \sigma _y = 0.5\) was used as initial condition in the domain \(x \times y \in [-3,3]^2\), with periodic boundary conditions.

The \(\nabla \cdot \) operator was discretized using the sixth order \(\hat{\textbf{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 necessary 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 \(\nabla \cdot 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.

8 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 \(\hat{\textbf{D}}\) and \(\textbf{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.