A New Analysis of an Implicit Mimetic Scheme for the Heat Equation

Mariela Castillo Nava^1,2, Juan Manuel Guevara-Jordan²

¹Universidad Metropolitana, Caracas, Venezuela.
²Universidad Central de Venezuela, Caracas, Venezuela.

DOI: 10.4236/jamp.2023.113056                

Abstract

A mimetic finite difference scheme for the transient heat equation under Robin’s conditions is presented. The scheme uses second order gradient and divergence mimetic operators, on a staggered grid, to approximate the space derivatives. The temporal derivative is replaced by a first order backward difference approximation to obtain an implicit formulation. The resulting scheme contains nonstandard finite difference stencils. An original convergence analysis by the matrix’s method shows that the proposed scheme is unconditionally stable. A comparative study against standard finite difference schemes, based on central difference or first order one side approximations, reveals the advantages of our scheme without being its implementation more expensive or difficult to achieve.

Keywords

Mimetic, Finite Difference, Robin’s Condition, Implicit, Stable, Transient, Heat

1. Introduction

The mimetic methods for partial differential Equations (PDEs) have their origin in the Soviet Union atomic program with the seminal work of academicians Tikhonov and Samarskii [1] , which describes discretizations of differential operators satisfying discrete versions of integral identities. Their ideas were extended and tested in the eighties, see for example [2] [3] , when it was realized the importance of these discretizations for solving PDEs with discontinuous coefficients. In the nineties, after the end of the Soviet Union, research on mimetic methods moved to the west, in particular to the national labs in the United States. Most of the initial research on mimetic methods in the United States was summarized in [4] , where the mimetic methods are called or named support operators methods. Starting this century research on mimetic methods has been extremely intensive and several versions of them are available at this time. Each one of them has its advantages and disadvantages, so there is not the best formulation yet. Good summaries and description of them can be found in the research books [5] [6] [7] . The main idea behind all the mimetic methods is to provide discretizations of gradient, divergence, and curl operators and generalized inner products in such a way that they satisfy discrete version of main integral identities from vector analysis. These mimetic operators are easily combined to obtain mimetic schemes associated to specific PDEs. This approach is different from traditional methods whose main concern is the numerical schemes for specific equations rather than the operators.

In this article, a second order mimetic operators described in [8] [9] will be analyzed in the context of the transient diffusivity equation. This version of the mimetic operators has been successfully applied to elliptic equations [9] [10] , transient diffusivity equations [11] - [19] , reservoir flow problems [20] [21] , the acoustic wave equation [22] [23] , the biharmonic equation [24] and the biharmonic wave equation [25] . From all these references, the more relevant to this article are [11] - [18] [20] [21] , which presents mimetic schemes for the heat or diffusivity equations, and a brief review of their content will be described. The presentation of the unpublished results of reference [19] will be developed in this paper.

The series of articles [11] [12] [13] [14] presents a complete analysis of an original mimetic Crank Nicolson scheme for the heat equations. Their convergence analysis shows the stability and consistence of the scheme, and it is similar to our analysis but ours differ in detail and content. Specifically, our analysis of convergence studies the fully implicit scheme and makes complete use of truncation errors at the border, which are not used in the analysis presented for the Crank Nicolson scheme.

The thesis [15] presents for first time the algorithm for the implicit mimetic scheme analyzed in this paper. Results from numerical implementation in the thesis gave evidence of the quadratic convergence in space of the scheme and its advantage over standard finite difference method based on ghost point. No formal convergence analysis of the scheme is presented in the thesis.

Reference [16] applies the method of the lines by using second order mimetic discretization in space and solving the system of ordinary differential equations by computing the exponential of the matrix associated to the system, and no analytic analysis of convergence is provided. Numerical results show the advantage of the mimetic approximations against standard finite difference scheme. All the theoretical and numerical work is presented for Dirichlet conditions.

Recent articles [17] [18] apply an implicit mimetic finite difference scheme for the heat equation to adaptative grid refinement techniques and imaging processing. Their scheme is similar to the one presented in this article but they differ in details and it does not have an analytical analysis of convergence. More comments about their scheme are presented in the conclusions section of this paper.

In [20] [21] , they develop mimetic schemes for the convection diffusion equations; these schemes are explicit and they reduce to numerical schemes for the heat equation when the velocity is null; no formal analysis of convergence are provided and the stabilities are computed numerically or they are taken from standard finite difference analysis.

Our work fully analyzes and extends the numerical scheme presented in [15] . In particular, a complete convergence analysis is given by showing, analytically, the consistence and stability of the scheme based on the matrix approach, which is the best adapted for the mimetic scheme analyzed in this work. Moreover, numerical implementation of the scheme is tested and evaluated in terms L₂ errors and numerical convergence rate in space and time. To the best of our knowledge, no complete analysis of the second order implicit scheme for the heat o diffusivity equations has been previously reported in the technical literature based on the mimetic operators described in [8] [9] .

The content of this article is divided in seven sections. First section is this introduction, second section is describes the second order mimetic operators used to deduct the scheme, third section deducts and states the implicit mimetic for the diffusivity or heat equations, fourth section provides the convergence analysis, and fifth section present a complete numerical comparative study of the standard finite differences schemes based on lateral derivatives and the ghost point central difference formulation against the implicit scheme analyzed in the article. Sixth section gives the conclusions and recommendations based on the results of this work. Finally, seventh section is an Appendix. It contains an extension of the implicit mimetic scheme for the heat equation with variable coefficients.

2. Mimetic Operators

In this article all the mimetic discretizations are developed for uniform staggered grids. One dimensional second order mimetic discretizations for the gradient and divergence operators have been well documented in [6] [8] [9] . Their description will be briefly presented, for sake of completeness, in reference to the one dimensional uniform staggered grid displayed in Figure 1.

In this figure crosses represent block’s centers and black circles are the block’s edges. One dimensional mimetic discretization of divergence operator is located at the block’s center and it is defined by standard central difference

$D{u}_{i+1/2}\equiv \left(u_{i+1}-u_i\right)/h$ (1)

Figure 1. One dimensional staggered grid.

with $h\equiv x_{i+1}-x_i$ . Its matrix’s formulation will be denoted by D, which contains two additional rows, the first and last, filled with zeros in order to obtain matrix’s dimensional consistency. The gradient mimetic discretizations are defined at the block’s edges by a central finite difference

$G{u}_i\equiv \left(u_{i+1/2}-u_{i-1/2}\right)/h$ , (2)

$h\equiv x_{i+1/2}-x_{i-1/2}$ , and one side second order approximations at boundary edges

$G{u}_0\equiv \left(-\left(8/3\right)u_0+3u_{1/2}-\left(1/3\right)u_{3/2}\right)/h$ , (3)

$G{u}_n\equiv \left(\left(8/3\right)u_n-3u_{n-1/2}+\left(1/3\right)u_{n-3/2}\right)/h$ . (4)

Its associated matrix’s will be represented by G. It can be easily deducted from this description that D and G have second order truncation error at all grid points.

The above gradient and divergence mimetic discretizations satisfy a discrete versions of the Green-Gauss-Stokes’ theorem, which takes the following form

${\langle \text{D}u\mathrm{,}v\rangle }_{hQ}+{\langle \text{G}u\mathrm{,}v\rangle }_{hP}={\langle \text{B}u\mathrm{,}v\rangle }_I\mathrm{.}$ (5)

In this equation B is the boundary operator and the bracket ${\langle \cdot \mathrm{,}\cdot \rangle }_A$ is a generalized inner product of the form ${\langle u\mathrm{,}v\rangle }_A\equiv u^{t}Av$ with weight A. Matrix Q is diagonal with the midpoint quadrature weights as its entries. P is a diagonal

matrix with the 3/8’s Newton Cotes quadrature weights, I is the identity matrix,

and h the grid’s block size. The Green-Gauss-Stokes’ theorem (5) provides an explicit expression for the boundary operator B,

$\text{B}=hQ\text{D}+h{\left(\text{G}\right)}^{t}P$ (6)

which is an immediate consequence of the generalized inner products properties.

Mimetic operators D, G, and B are the only ones needed for discretization of the heat equation in a one dimensional context. Their description and deduction of mimetic operators in a general multidimensional case can be easily achieved by systematic use of Kronecker’s tensor product as described in [6] [7] [10] [22] . However, our analysis will be restricted to the one dimensional context because it allows a simpler, elegant, and unique presentation. Analysis of multidimensional mimetic scheme for the heat equation can be performed in many ways, but all them will be based in the one dimensional case reported in this paper.

3. Mimetic Scheme

In this section an implicit mimetic scheme for the heat equation is presented. The model of the heat equation used in this article takes the form

$\frac{\partial u}{\partial t}-\Delta u=f$ (7)

where u is the temperature, f is a source term, $\partial$ is a partial derivative, t is time and Δ the Laplace’s operator or laplacian. For simplicity all the variables in (7) will be assumed as dimensionless. The laplacian can be expressed as a composition of the divergence, $\nabla \cdot$ , and the gradient, $\nabla$ , therefore Equation (7) can be written as

$\frac{\partial u}{\partial t}-\nabla \cdot \nabla u=f\mathrm{.}$ (8)

This representation is the most convenient for describing the mimetic scheme. Let us denote by T a square matrix of order $n+2$ , whose entries are all null with the exception of the diagonal entries ${\text{T}}_{i\mathrm{,}i}=1/dt$ for $i=2,\cdots ,n+1$ and dt denotes the size of the time step of the mimetic scheme that will be developed. If U_m represents the mimetic approximations to u at the time level m, F_m is the projections of f on the staggered grid at the time m, the continuous gradient and divergence operators in (8) are replaced by their mimetic approximations, and the partial time derivative in (8) is replaced by the vectorial expression for the backward finite difference approximation $\text{T}\left(U_{m+1}-U_m\right)$ then

$\text{T}\left(U_{m+1}-U_m\right)-\text{DG}{U}_{m+1}=F_{m+1}$ (9)

represents an implicit mimetic approximation to the heat Equation (8). The above expression requires boundary and initial conditions to complete the formulation of the new mimetic scheme.

Initial condition give the temperature u at $t=0$ . Its discretization is the standard one, projection on the staggered combined with a backward finite difference in time, which provides approximations to U₀ needed to start the numerical scheme.

Standard boundary conditions for the heat equation are Dirichlet’s or Neumann’s condition, which can be easily combined in a general Robin’s condition

$\beta \frac{\partial u}{\partial n}+\alpha u=h$ (10)

where $\beta$ , $\alpha$ , and h are known boundary’s functions or constant. Partial derivative respect n is the directional derivative respect to the exterior normal unitary vector on the boundary. If $\beta$ or $\alpha$ is null then the Dirichlet or Neumann condition are obtained respectively. Equation (7) under Robin’s boundary conditions produces a well posed initial boundary value problem [26] [27] . Mimetic discretizations of (10) is

$B\text{BG}{U}_{m+1}+A{U}_{m+1}=H_{m+1}$ (11)

where B, A, and $H_{m+1}$ are the projections of $\beta$ , $\alpha$ , and $h\left(\cdot \mathrm{,}\left(m+1\right)dt\right)$ on the staggered grid. B and A are, for simplification, time independent diagonal matrices and $H_{m+1}$ is a time dependent vector. By the superposition principle the mimetic discretizations of Robin’s condition (11) and the implicit mimetic discretization of the heat Equation (9) can be combined in a single vectorial or matricial expression

$\left(T-\text{DG}+B\text{BG}+A\right)U_{m+1}=\text{T}{U}_m+F_{m+1}+H_{m+1}\mathrm{.}$ (12)

which represents the complete implicit mimetic scheme for the heat equation. The mathematical and numerical analysis of this scheme are the main objectives of this paper.

In order to appreciate the new features of the proposed schemes, equations defined as entries of (12) will be describe on the grid in Figure 1. First equation represents the left side boundary condition

$\left(\frac{{\beta }_{o}8}{3h}+{\alpha }_o\right)U_0^{m+1}-\frac{3{\beta }_o}{h}{U}_{\frac{1}{2}}^{m+1}+\frac{{\beta }_o}{3h}{U}_{\frac{3}{2}}^{m+1}=h_0^{m+1}$ (13)

which is a one side approximation with three terms to ensure second order accuracy. Next equation at the first inner node is given by

$-\left(\frac{1}{3h}+\frac{8}{3h^2}\right)U_o^{m+1}+\left(\frac{4}{h^2}+\frac{1}{dt}+\frac{1}{2h}\right)U_{\frac{1}{2}}^{m+1}-\left(\frac{4}{3h^2}+\frac{1}{6h}\right)U_{\frac{3}{2}}^{m+1}=\frac{U_{\frac{1}{2}}^m}{dt}+F_{\frac{1}{2}}^{m+1}$ (14)

which is first order in time and space. This reduction of accuracy in space was expected due to the non uniform distribution of staggered grid nodes near the boundary. Notice the unusual linear terms for h in the denominators, they come from the superposition of the boundary conditions. The equation associated to the second inner node takes the following expression

$\frac{1}{3h}{U}_o^{m+1}-\left(\frac{1}{2h}+\frac{1}{h^2}\right)U_{\frac{1}{2}}^{m+1}+\left(\frac{1}{dt}+\frac{1}{6h}+\frac{2}{h^2}\right)U_{\frac{3}{2}}^{m+1}-\frac{1}{h^2}{U}_{\frac{5}{2}}^{m+1}=\frac{U_{\frac{3}{2}}^m}{dt}+F_{\frac{3}{2}}^{m+1}$ (15)

This equation is also first order in time and space. This first order of accuracy in space is due the four points non symmetric laplacian approximation around $x_{3/2}$ . After the second inner node discretizations take the standard finite difference heat equation approximation

$-\frac{1}{h^2}{U}_{i-\frac{1}{2}}^{m+1}+\left(\frac{2}{h^2}+\frac{1}{dt}\right)U_{i+\frac{1}{2}}^{m+1}-\frac{1}{h^2}{U}_{i+\frac{3}{2}}^{m+1}=\frac{U_{i+\frac{1}{2}}^m}{dt}+F_{i+\frac{1}{2}}^{m+1}$ (16)

They have second and first order of accuracy in space and time, respectively. The equations associated with nodes $x_{n-1/2}$ , $x_{n-3/2}$ , and $x_n$ are reflection of those associated with the first and second inner nodes, (14) and (15), and the left boundary condition, (13), respectively. All these equations generates an almost tridiagonal linear system represented in symbols by (12), which contains two submatrices of third order at the upper left and lower right corners of the system, and it is structurally symmetric.

Sometimes the heat Equation (7) is too simple to model many situations in applications. An alternative heat equation is developed for those cases and it takes the form $c\cdot {\partial }_{t}u-\nabla \cdot \left(k\nabla u\right)=f$ , where c and k are positive scalar functions which represents specific physical parameters in each application [27] . Note that in our work k is a scalar function, but in a multidimensional context it represents a tensor. The mimetic scheme (12) can be modified to produce an implicit mimetic scheme for the alternative heat equation, it takes the expression

$\left(\text{CT}-\text{DKG}+B\text{BG}+A\right)U_{m+1}=\text{CT}{U}_m+F_{m+1}+H_{m+1}\mathrm{.}$ (17)

where C and K are square diagonal matrices and they represent the projections of c and k on the staggered grid. The matrix C has order n + 2, it contains in its diagonal the projection of c on the block centers in the staggered grid and its first and last rows are null in order to obtain dimensional consistency. Matrix K has order n + 1 and it contains the projections of k at the block edges. Description of the equations generated by (17) will be presented in Appendix. They won’t be used in the sequel.

4. Convergence Analysis

The convergence analysis presented in this section is based on the well-known Lax-Richtmyer or Lax’s equivalence theorem [28] [29] , which states that consistency and stability of time dependent numerical schemes of a well-posed initial boundary value problem imply its convergence. Therefore, the consistency and stability conditions associated with the implicit mimetic scheme described in this paper will be established in the next two subsections.

4.1. Consistency

The consistency analysis will be based on the systematic application of Taylor’s expansions to each one of the finite difference equations of the implicit mimetic scheme developed in the previous section. Let us define $u_k^m\equiv u\left(x_k\mathrm{,}m\cdot dt\right)$ , where u is the exact solution of the continuous one-dimensional heat equation under Robin’s condition, $m\cdot dt$ denote the m time step of size dt and $x_k$ represents a generic node in the staggered grid. In order to establish the consistence of the mimetic scheme the Taylor expansion of u at the block center or edge, in the case of boundary conditions, will be evaluated in the mimetic finite difference equations to obtain the truncation error. In that sense the expansion for the boundary condition (13) provides the following truncation error

$\left(\frac{{\beta }_{o}8}{3h}+{\alpha }_o\right)u_0^{m+1}-\frac{3{\beta }_o}{h}{u}_{\frac{1}{2}}^{m+1}+\frac{{\beta }_o}{3h}{u}_{\frac{3}{2}}^{m+1}-\left({\beta }_o{\partial }_n+\alpha \right)u_0^{m+1}=O\left(h^2\right)$ (18)

where

$O\left(h^2\right)=\left({\beta }_o{h}^2/8\right){\partial }_x^{3}u\left(\mathrm{0,}\left(m+1\right)\cdot dt\right)+\left({\beta }_o{h}^3/16\right){\partial }_x^{4}u\left(\mathrm{0,}\left(m+1\right)\cdot dt\right)+\cdots$ (19)

is a second order consistent approximation in space. The same analysis for Equation (14) produces the truncation error

$\begin{array}{l}\frac{u_{\frac{1}{2}}^{m+1}-u_{\frac{1}{2}}^m}{dt}-\left(\frac{1}{3h}+\frac{8}{3h^2}\right)u_o^{m+1}+\left(\frac{4}{h^2}+\frac{1}{2h}\right)u_{\frac{1}{2}}^{m+1}-\left(\frac{4}{3h^2}+\frac{1}{6h}\right)u_{\frac{3}{2}}^{m+1}\\ -\left({\partial }_t-{\partial }_x^2\right)u_{\frac{1}{2}}^{m+1}=O\left(dt\right)+O\left(h\right)\end{array}$ (20)

in which

$O\left(dt\right)+O\left(h\right)=-\left(h/8\right){\partial }_x^{2}u\left(h/2\mathrm{,}\left(m+1\right)\cdot dt\right)+\left(dt/2\right){\partial }_t^{2}u\left(h/2\mathrm{,}\left(m+1\right)\cdot dt\right)+\cdots$ (21)

and it represents a first order consistent approximation in space and time. Next equation in the mimetic scheme is given by (15) and its truncation error is

$\begin{array}{l}\frac{u_{\frac{3}{2}}^{m+1}-u_{\frac{3}{2}}^m}{dt}+\frac{1}{3h}{u}_o^{m+1}-\left(\frac{1}{2h}+\frac{1}{h^2}\right)u_{\frac{1}{2}}^{m+1}+\left(\frac{1}{6h}+\frac{2}{h^2}\right)u_{\frac{3}{2}}^{m+1}-\frac{1}{h^2}{u}_{\frac{5}{2}}^{m+1}\\ -\left({\partial }_t-{\partial }_x^2\right)u_{\frac{3}{2}}^{m+1}=O\left(dt\right)+O\left(h\right)\end{array}$ (22)

in this equation

$O\left(dt\right)+O\left(h\right)=\left(h/8\right){\partial }_x^{2}u\left(3h/2\mathrm{,}\left(m+1\right)\cdot dt\right)+\left(dt/2\right){\partial }_t^{2}u\left(3h/2\mathrm{,}\left(m+1\right)\cdot dt\right)+\cdots$ (23)

and it also gives a first order consistent approximation in space and time. At the internal blocks in the staggered which are those with at least two blocks away from the boundaries the standard implicit finite difference approximation for the heat equation is given by Equation (16) and its truncation error takes the usual form

$\frac{u_{i+\frac{1}{2}}^{m+1}-u_{i+\frac{1}{2}}^m}{dt}-\frac{1}{h^2}{u}_{i-\frac{1}{2}}^{m+1}+\frac{2}{h^2}{u}_{i+\frac{1}{2}}^{m+1}-\frac{1}{h^2}{u}_{i+\frac{3}{2}}^{m+1}-\left({\partial }_t-{\partial }_x^2\right)u_{i+\frac{1}{2}}^{m+1}=O\left(dt\right)+O\left(h^2\right)$ (24)

where

$\begin{array}{c}O\left(dt\right)+O\left(h^2\right)=-\left(h^2/12\right){\partial }_x^{4}u\left(\left(i+3/2\right)h\mathrm{,}\left(m+1\right)\cdot dt\right)\\ +\left(dt/2\right){\partial }_t^{2}u\left(\left(i+3/2\right)h\mathrm{,}\left(m+1\right)\cdot dt\right)+\cdots \end{array}$ (25)

which represents a standard second order $O\left(h^2\right)$ in space and first order in time $O\left(dt\right)$ consistent approximation. The consistency analysis associated with equations at nodes $x_n$ , $x_{n-1/2}$ , and $x_{n-3/2}$ are identically to those presented above to the left boundary condition along with the equations associated to the first and second inner nodes, so they are omitted.

4.2. Stability

The stability analysis will prove that our implicit mimetic scheme satisfies the Lax-Richtmyer condition. This is achieved by subtracting the mimetic finite difference equations from the consistence equations in order to evaluate the error $e_k^m\equiv u_k^m-U_k^m$ between the analytic, $u_k^m$ , and the mimetic approximations, $U_k^m$ , obtained by the proposed implicit scheme. The subtraction generates a linear system

$\left(\text{T}-\text{DG}+B\text{BG}+A\right)e_{m+1}=\text{T}{e}_m+T{R}_{m+1}$ (26)

where $e_{m+1}$ , $e_m$ and $T{R}_{m+1}$ , $T{R}_m$ are the column vectors for the errors and the truncation errors at the m + 1 and m time steps, respectively. Some transformations and assumptions are required to simplify the stability analysis. The first simplification is that boundary conditions coefficients are set equal to one. This assumption does not affect the generality of the analysis because it is always possible to normalize the coefficients of directional derivatives and the Dirichlet coefficient in the Robin’s conditions does not play an important role in the analysis. The second assumption or transformation is that all the equations in the system (26) will be multiplied by dt and it takes the form

$\left(\text{I}-dt\text{DG}+dt\left(B\text{BG}+A\right)\right)e_{m+1}=\text{I}{e}_m+dt\cdot T{R}_{m+1}$ (27)

where I represents an identity matrix whose first and last rows are null because of the definition of T. Taking into account such transformations and solving for $e_0^{m+1}$ and $e_n^{m+1}$ in the first and last equations respectively results

$e_0^{m+1}=\frac{3e_{\frac{1}{2}}^{m+1}-\frac{1}{3}{e}_{\frac{3}{2}}^{m+1}}{h+\frac{8}{3}}+\frac{O\left(h^3\right)}{h+\frac{8}{3}}\approx e_0^{m+1}=\frac{9e_{\frac{1}{2}}^{m+1}-e_{\frac{3}{2}}^{m+1}}{8}+\frac{O\left(h^3\right)}{h+\frac{8}{3}}$ (28)

and a similar expression for $e_n^{m+1}$ . Notice that dt does not appear in (28) because the first and last equations in system (27) represent the boundary conditions. Since $e_0^{m+1}$ is only presented in the equations associated to the first and second inner nodes in system (27) then substitution of (28) in those equations results in the following relations.

$\left(\frac{dt}{8h}+\frac{dt}{h^2}+1\right)e_{\frac{1}{2}}^{m+1}-\left(\frac{dt}{8h}+\frac{dt}{h^2}\right)e_{\frac{3}{2}}^{m+1}=e_{\frac{1}{2}}^m+O\left(d{t}^2\right)+O\left(dt\cdot h\right)$ (29)

$-\left(\frac{dt}{8h}+\frac{dt}{h^2}\right)e_{\frac{1}{2}}^{m+1}+\left(1+\frac{dt}{8h}+\frac{2dt}{h^2}\right)e_{\frac{3}{2}}^{m+1}-\frac{dt}{h^2}{e}_{\frac{5}{2}}^{m+1}=e_{\frac{3}{2}}^m+O\left(d{t}^2\right)+O\left(dt\cdot h\right)$ (30)

Similar expressions are obtained for equations associated to nodes $x_{n-1/2}$ , and $x_{n-3/2}$ after substitution in them of the analog Equation (20) for $e_n^{m+1}$ . The other equations in system (27) take the well known finite difference expression.

$-\frac{dt}{h^2}{e}_{i-\frac{1}{2}}^{m+1}+\left(1+\frac{2dt}{h^2}\right)e_{i+\frac{1}{2}}^{m+1}-\frac{dt}{h^2}{e}_{i+\frac{3}{2}}^{m+1}=e_{i+\frac{1}{2}}^m+O\left(d{t}^2\right)+O\left(dt\cdot h^2\right)$ (31)

The reduced linear system obtained by substitution of (28) and its similar expression for $e_n^{m+1}$ in the equations associated to inner nodes in the staggered grid are represented by Equations (29), (30), and (31). The eigenvalues $\lambda$ of this reduced linear system can be easily estimated by Greschgorin theorem [28] [29] , which allows us to state that they satisfy the estimate $1<\left|\lambda \right|$ . This estimate implies the stability of the proposed mimetic scheme. The main contribution of this stability analysis is that it takes into account the truncation errors from the equations for the boundary condition, which are systematically dropped in some published analysis [11] [12] [13] [14] in the context of Crank-Nicolson schemes.

Our analysis of convergence includes systematically the truncation errors from the equations associated to the boundary conditions and avoids their use as simple constrains [19] . Although dropping truncation errors from the boundaries conditions is acceptable for Dirichlet problems that are not the case for general Robin’s conditions considered in this work.

5. Numerical Study

This section will present and analyze a comparative study of three numerical solutions for an initial boundary value (IBVP) problem for the heat equation in a one dimensional context, whose analytic solutions are known. Three different discretizations techniques will be used to obtain the numerical solutions: the standard implicit finite difference method with first order lateral or one side finite difference approximation for the boundary conditions, the standard implicit finite difference method with second order central difference approximation based on ghost points for the boundary conditions, and the implicit mimetic scheme developed in this article. All these methods will be discretized on the staggered grid described in Figure 1.

Some comments about the standard finite difference methods, [27] [28] [29] , are needed because they are being used on a staggered. The standard finite difference discretizations of the one dimensional heat equation are identical to those obtained for the mimetic scheme but dropping the terms (k/h) in the coefficients, where k is a constant. This observation applies specifically to Equations (14) and (15) and their analogous expression at the grid points $x_{n-1/2}$ and $x_{n-3/2}$ . In the case of the implicit finite difference scheme based on ghost points two additional finite difference equations for the heat equation are needed at the points $x_o$ and $x_n$ respectively. Those equations are identical to (16) but their coefficients are multiply by four. The first order lateral or one side approximation for Robin’s condition (10) takes, on the left boundary of the grid, the form

$\left(\frac{{\beta }_{o}2}{h}+{\alpha }_o\right)U_o^{m+1}-\frac{{\beta }_{o}2}{h}{U}_{\frac{1}{2}}^{m+1}=h_0^{m+1}$ (32)

and similar expression is obtained as function of $U_n^{m+1}$ and $U_{n-1/2}^{m+1}$ for the

boundary condition at the right end of the grid. The second order central finite difference approximation with ghost point for (10) at left boundary of the grid is given by

$\frac{{\beta }_o}{h}{U}_{-\frac{1}{2}}^{m+1}+{\alpha }_o{U}_o^{m+1}-\frac{{\beta }_o}{h}{U}_{\frac{1}{2}}^{m+1}=h_0^{m+1}$ (33)

where $U_{-1/2}^{m+1}$ is the variable associated with the ghost point, which is the reflection point of $x_{1/2}$ against $x_o$ . An analogous equations results as a function of $U_n^{m+1}$ , $U_{n-1/2}^{m+1}$ and the ghost variable $U_{n+1/2}^{m+1}$ associated to the ghost point obtained by reflection of $x_{n-1/2}$ against $x_n$ . Truncation errors for Equations (32) and (33) are $O\left(h\right)$ and $O\left(h^2\right)$ respectively. Notice that truncation errors for the heat equation discretizations in all the numerical methods described in this article are first order in time and second or first order in space. In particular, at the nodes or grid points near the boundary they are first order in space and time, $O\left(dt\right)+O\left(h\right)$ , while they are first order in time and second order in space, $O\left(dt\right)+O\left(h^2\right)$ , for inner grid points two blocks away from the boundary.

The IBVP considered as test problem in this study will be evaluated as a function of accuracy and convergence rates in time and space. The linear or algebraic systems generated by their discretizations are sparse and banded, they will be solved by LU decomposition [28] [29] . The test analyzes the numerical solution of the nonhomogeneous heat equation

${\partial }_{t}u-{\partial }_x^{2}u=-\left[{\text{e}}^{-t-5x}-{\text{e}}^{-t}+25{\text{e}}^{-t-5x}\right]/\left({\text{e}}^{-5}-1\right)$ (34)

on the plane region $\left(\mathrm{0,1}\right)\times \left(\mathrm{0,}\infty \right)$ , under Robin’s conditions

$\left(\left(-{\text{e}}^{-5}\right)u-\left[\left({\text{e}}^{-5}-1\right)/\left(-5\right)\right]{\partial }_{x}u\right)\left(\mathrm{0,}t\right)=-{\text{e}}^{-t}$ (35)

and

$\left(\left(-{\text{e}}^{-5}\right)u-\left[\left({\text{e}}^{-5}-1\right)/\left(-5\right)\right]{\partial }_{x}u\right)\left(\mathrm{1,}t\right)=0$ , (36)

initial condition

$u\left(x\mathrm{,0}\right)=\left({\text{e}}^{-5x}-1\right)/\left({\text{e}}^{-5}-1\right)$ , (37)

and analytical solution

$u\left(x\mathrm{,}t\right)={\text{e}}^{-t}\left({\text{e}}^{-5x}-1\right)/\left({\text{e}}^{-5}-1\right)$ . (38)

This IBVP problem is discretized for the three methods described in this section for grid sizes ranging from six to seventy grid blocks in space, where the block size $h=1/n$ with n denoting the number of grid blocks, and a time step size $dt=h^2=1/n^2$ . All the runs were stopped at $t=1$ . The numerical results of all these runs are displayed for some specific representative values or parameters at the Table 1 and Figure 2 which will be analyzed and described in a sequential order.

Table 1 contains two types of information separated in two regions by a horizontal line. In the upper region there are four columns, they contain information on the spatial grid size while second, third and fourth columns show the associated L² errors for the three numerical methods or schemes of discretization considered in this paper. Each column for the errors shows that their magnitude decreases monotonically as the number of grid blocks increases, which implies that the computer implementation of the three methods is convergent undergrid refinement. It can be observed that the mimetic scheme, proposed and analyzed in this paper, is one order of magnitude more accurate than the central finite difference based on ghost points and two orders of magnitude better than the standard finite difference based on one-side first order approximation at the

Table 1. Grid sizes, errors in L₂ norm, convergence rates in space and time.

Figure 2. Errors in L² norm vs. spatial grid sizes (left) and time grid sizes (right).

borders without ghost points. In the lower region of the table there are also four columns with the convergence rates associated to each one of the methods under analysis, they are presented in the second to the fourth columns. The convergence rate as a function of space grid size shows the expected rates for the standard lateral and central finite difference approximations, namely first and second order rates in space [28] [29] . However, the second order convergence rate for the mimetic method cannot be easily deducted from the analysis presented in previous sections, because it was proved that truncations errors for the mimetic discretization of the heat equation is just first order in space for grid blocks near or close to the border. This second order convergence rate for the mimetic scheme is the result of the second order approximation in the Robin’s condition which enforces the modulus maximum principle [27] [28] [30] , in the L² errors for the mimetic scheme. Moreover, it is known in the numerical analysis literature [28] [30] [31] , that numerical schemes associated to the heat equation inherit the converge rate of the same numerical scheme for the Laplacian or the elliptic operator in space present in the heat or parabolic equations. Since the mimetic scheme for the static diffusion equations has an optimum second order convergence rate [9] , then that rate in space is observed for the mimetic scheme in Table 1. In the last row of the table we find the convergence rates in time for the three schemes. The mimetic and the central difference scheme with ghost point present first order convergence rates as it was expected from the truncation error analysis. The one half or sub-linear convergence rate observed in the table for the lateral finite difference scheme is a natural result of the time step size. In effect, since the overall convergence rate of the method is theoretically $O\left(dt\right)+O\left(h\right)$ and $dt=h^2$ then $O\left(dt\right)+O\left(\sqrt{dt}\right)=O\left(\sqrt{dt}\right)$ which is the sub-linear rate observed in the table and implies the first order convergence rate in time.

Figure 2 contains two graphs. The left graph has the log of L² errors in the vertical axis and the log of spatial grid block size in the horizontal axis. The upper, middle, and lower curves represent the errors associated to finite different scheme with one-side approximation at the border, the second order central finite difference scheme with ghost point, and the mimetic scheme, respectively. That order reveals that the best approximation to the analytic solution is achieved by the mimetic scheme. The three curves are straight lines whose slopes provide the convergence rate of each method or scheme in space. It is quite evident that the slopes of the lower two lines are graphically identical and they correspond to the second order convergence rate in space for the central finite difference approximation with ghost points and the mimetic scheme. The slope of the upper line is half of the slope of the lower lines and it represents the first order convergence rate in space for finite difference scheme with lateral approximation at the border. The right graph has the same information, as the left graphs, in the vertical axis but the horizontal axis represents the log of the time grid block size. Since the time step size $dt=h^2$ then the horizontal axis is larger and it is almost twice the length of the one in the left graph. Consequently, the slopes of the lines in the right graph are reduced by a half of those found in the left graph. The curves in the right graphs represent the same L² errors of the left graph but their reduced slopes are the convergence rates in time of their associated numerical methods. Specifically, the lower lines have slopes approximately equal to 1 and the upper line has slope equal to 0.5, all in agreement with the information displayed in Table 1.

The main conclusion of this numerical study is that the implicit mimetic scheme analyzed in this paper belongs to the set or class of numerical schemes whose truncations errors do not provide an optimum convergence rate estimate. The test does not allow us to conclude that the mimetic scheme is better than standard finite difference. In fact, in our experience, the mimetic scheme always achieves comparable or similar accuracy as the standard finite difference method based on ghost points for easy problems. However, the implicit mimetic scheme allows a simpler computational implementation or program development and generates linear system with a smaller number of unknowns.

6. Conclusions

This paper has presented an original mathematical proof of convergence for an implicit mimetic finite difference scheme for the heat equation [19] . The scheme was originally proposed, without analytical proof of convergence, in [15] . The numerical study in this paper extends the results presented in [15] by analyzing the convergence rates in time, explaining the second order convergence rate for the mimetic scheme in space, and adding one more numerical method in the comparative study, namely the standard finite difference method with lateral or one-side approximation at the boundaries. All the theoretical and numerical results presented in this paper are completely justified and they provide a solid foundation for the proposed mimetic scheme. The theoretical results show that the proposed mimetic scheme converges without given optimum upper bounds for its convergence rate. Numerical results give strong evidence that the convergence rate obtained by the truncation error analysis is too conservative and an optimum second order convergence rate in space is achieved numerically and justified on the modulus maximum principle for the heat equation. To the best of our knowledge, there is no analytical proof of optimum convergence rate for the implicit mimetic scheme presented in this paper and it is an open research problem.

In recent years, several applications of a similar, not equal, implicit mimetic scheme for the heat equation based on the mimetic operators described in this paper have been published in academic journals [17] [18] . In our notation, the implicit mimetic scheme proposed in those articles, in its simplest form, takes the expression

$\text{I}-\left({\text{T}}^{-1}\text{DG}+B\text{BG}+A\right)U_{m+1}=U_m+{\text{T}}^{-1}{F}_{m+1}+H_{m+1}$ (39)

which is certainly different from our scheme in Equation (12). It would be interesting to perform the convergence analysis for (39) and a comparative study of the mimetic schemes (12) and (39) on uniform staggered grids under Robin’s conditions. However, many authors of mimetic papers make use of simplified versions of the mimetic operator B in their numerical implementations. The use of a simplified B in (12) or (39) produces an identical conservative finite difference scheme but it is not mimetic. The conservative schemes associated with a simplified B were fully studied for the heat equation in [32] .

Acknowledgements

We thank God.

Appendix

This appendix presents the finite difference equations associated to the implicit mimetic scheme in the case of a general transient diffusion equation $c\cdot {\partial }_{t}u-\nabla \cdot \left(k\nabla u\right)=f$ described by (17). Following the style of the paper the equations associated to the first inner node, the second inner node, and the generic inner node will be presented. Analogous expression could be easily found for the others nodes. Equation for the Robin’s boundary condition is identical as (13) so it won’t be presented here. As usual all the sub-index are referred to the staggered grid in Figure 1.

Equation for the first inner node:

$\begin{array}{l}-\left(\frac{1}{3h}+\frac{8k_o}{3h^2}\right)U_o^{m+1}+\left(\frac{3k_o+k_1}{h^2}+\frac{c_{\frac{1}{2}}}{dt}+\frac{1}{2h}\right)U_{\frac{1}{2}}^{m+1}-\left(\frac{k_o}{3h^2}+\frac{k_1}{h^2}+\frac{1}{6h}\right)U_{\frac{3}{2}}^{m+1}\\ =\frac{c_{\frac{1}{2}}{U}_{\frac{1}{2}}^m}{dt}+F_{\frac{1}{2}}^{m+1}\end{array}$ (40)

Equation for the second inner node:

$\begin{array}{l}\frac{1}{3h}{U}_o^{m+1}-\left(\frac{1}{2h}+\frac{k_1}{h^2}\right)U_{\frac{1}{2}}^{m+1}+\left(\frac{c_{\frac{3}{2}}}{dt}+\frac{1}{6h}+\frac{k_1+k_2}{h^2}\right)U_{\frac{3}{2}}^{m+1}-\frac{k_2}{h^2}{U}_{\frac{5}{2}}^{m+1}\\ =\frac{c_{\frac{3}{2}}{U}_{\frac{3}{2}}^m}{dt}+F_{\frac{3}{2}}^{m+1}\end{array}$ (41)

Equation for the general inner node away from the boundaries:

$-\frac{k_i}{h^2}{U}_{i-\frac{1}{2}}^{m+1}+\left(\frac{k_i+k_{i+1}}{h^2}+\frac{c_{i+\frac{1}{2}}}{dt}\right)U_{i+\frac{1}{2}}^{m+1}-\frac{k_{i+1}}{h^2}{U}_{i+\frac{3}{2}}^{m+1}=\frac{c_{i+\frac{1}{2}}{U}_{i+\frac{1}{2}}^m}{dt}+F_{i+\frac{1}{2}}^{m+1}$ (42)

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Tikhonov, A.N. and Samarskii, A.A. (1962) Homogeneous Difference Schemes. URSS Computational Mathematics and Mathematical Physics, 1, 5-67. https://doi.org/10.1016/0041-5553(62)90005-8
[2]	Samarskii, A.A., Tishkin, V.F., Favorskii, A.P. and Shashkov, M. (1981) Operational Finite Difference Schemes. Differential Equations, 17, 854-862.
[3]	Samarskii, A.A., Tishkin, V.F., Favorskii, A.P. and Shashkov, M. (1983) Employment of the Reference Operators in the Construction of Finite Difference Analogs of Tensor Operations. Differential Equations, 18, 881-885.
[4]	Shashkov, M. (1996) Conservative Finite Difference Methods. CRC Press, Boca Raton, FL.
[5]	Beirao da Veiga, L., Lipnikov, K. and Manzini, G. (2014) The Mimetic Finite Difference Method for Elliptic Problems. Springer, New York. https://doi.org/10.1007/978-3-319-02663-3
[6]	Castillo, J.E. and Miranda, G. (2013) Mimetic Discretization Method. CRC Press, Boca Raton, FL. https://doi.org/10.1201/b14575
[7]	Justo, D. (2009) Higher Order Mimetic Methods and Absorbing Boundary Conditions. VDM Verlag, Riga.
[8]	Castillo, J.E. and Yasuda, M. (2005) Linear Systems Arising in Second Order Mimetic Divergence and Gradient Discretizations. Journal of Mathematical Modelling and Algorithms, 4, 67-82. https://doi.org/10.1007/s10852-004-3523-1
[9]	Guevara-Jordan, J.M., Rojas, S., Freites-Villegas, M. and Castillo, J.E. (2007) Convergence of a Mimetic Finite Difference Method for Static Diffusion Equations. Advances in Difference Equations, 2007, Article No. 012303. https://doi.org/10.1155/2007/12303
[10]	Blanco, J., Rojas, O., Chacon-Suescun, C., Guevara-Jordan, J.M. and Castillo, J.E. (2017) Tensor Formulation for 3D Mimetic Operators Applied to Elliptic Equations. Electronic Transactions in Numerical Analysis (ETNA), 45, 457-475.
[11]	Mannarino, I. (2009) A Mimetic Finite Difference Method Using Crank Nicolson Scheme for Unsteady Diffusion Equation. Revista de Matematicas: Teoria y Aplicaciones, 16, 221-230. https://doi.org/10.15517/rmta.v16i2.302
[12]	Mannarino, I. (2007) Un Método Mimético de Diferencias Finitas para la Ecuación No Estática de Difusión. Master’s Thesis in Mathematics, Universidad de Central de Venezuela, Caracas.
[13]	Mannarino, I., Rojas, O. and Zenaida, C. (2011) A Multidimensional Mimetic Scheme for the Heat Equation. Memorias del VIII Congreso Colombiano de Métodos Numéricos, Editorial de la Universidad EAFIT, Medellin.
[14]	Mannarino, I., Guevara-Jordan, J.M. and Quintana, Y. (2008) A Numerical Study of a Mimetic Scheme for Unsteady Heat Equation. Revista Faraute Ciencia Y Tecnología, 3, 64-70.
[15]	Fernandez, F. and Rocha, F. (2006) Evaluación del Desempeño de un Método Mimético Aplicado a la Ecuación de Conducción Unidimensional de Calor en Régimen Transitorio. Thesis in Mechanical Engineering, Universidad de Central de Venezuela, Caracas.
[16]	Hernandez-Walls, R., Castillo, J.E., Rojas-Mayoral, E.M. and Giron Nava, J.A. (2014) Semidiscrete Numeric Solution for the Non Stationary Heat Equation Using Mimetic Techniques. European Journal of Physics, 35, Article 065013. https://doi.org/10.1088/0143-0807/35/6/065013
[17]	Villamizar, J., Calderon, G., Carrillo, J., Bautista, L., Carrillo, J., Rueda, J. and Castillo, J. (2022) Mimetic Finite Difference Methods for Restauration of Fundus Images for Automatic Detection of Glaucoma Suspects. Computer Methods in Biomechanics and Biomedical Engineering: Imaging and Visualization, 10, 492-499. https://doi.org/10.1080/21681163.2021.1914733
[18]	Calderon, G., Carrillo, J.C., Villamizar, J., Torres, C. and Castillo, J. (2019) Estimation of Errors and Adaptability in Mimetic Differences for a Parabolic Type Problem. Revista UIS Ingenierias, 18, 141-149. https://doi.org/10.18273/revuin.v18n3-2019014
[19]	Castillo-Nava, M.C. (2019) Un Esquema Implícito de Diferencias Finitas Mimético para la Ecuación de Difusividad Transitoria. Inner Report, Universidad de Central de Venezuela, Caracas.
[20]	Corbino-Delgado, J. (2018) Subflow: Simulating Geological Storage of CO2 Using Mimetic Operators. Ph.D. Thesis, San Diego State University, San Diego.
[21]	Guevara-Jordan, J.M. and Arteaga-Arispe, J. (2007) A Second Order Mimetic Approach for Tracer Flow in Oil Reservoirs. Proceedings of the LAPEC Conference, Buenos Aires, April 2007, SPE-107366-MS. https://doi.org/10.2118/107366-MS
[22]	Solano-Feo, F., Guevara-Jordan, J.M., Rojas, O., Otero, B. and Rodriguez, R. (2016) A New Mimetic Scheme for the Acoustic Wave Equation. Journal of Computational and Applied Mathematics, 295, 2-12. https://doi.org/10.1016/j.cam.2015.09.037
[23]	Solano-Feo, F., Guevara-Jordan, J.M., Gonzalez, C., Rojas, O. and Otero, B. (2017) Modeling Seismic Wave Propagation Using Staggered Grid Mimetic Finite Difference. Bulletin of Computational and Applied Mathematics, 5, 9-28. https://doi.org/10.1016/j.cam.2015.09.037
[24]	Gomez-Polanco, A., Guevara-Jordan, J.M. and Molina, B. (2013) A Mimetic Iterative Scheme for Biharmonic Equations. Mathematical and Computer Modelling, 57, 2132-2139. https://doi.org/10.1016/j.mcm.2011.03.015
[25]	Solano-Feo, F., Guevara-Jordan, J.M. and Rojas, O. (2012) An Explicit Mimetic Method for Transient Beam Equation. In: Montero, B. and Poltronieri, J., Eds., Avances en Simulación Computacional y Modelaje Numérico, Memorias de Cimenics 2012, Davila, E., Uzcategui, G. and Cerrolaza, Y.M., Venezuela, 1-8.
[26]	Mikhailov, V.P. (1978) Partial Differential Equations. Mir Publishers, Moscow.
[27]	Tikhonov, A.N. and Samarskii, A.A. (2011) Equations of Mathematical Physics. Dover Publishers, New York.
[28]	Morton, K.W. and Mayers, D.F. (2001) Numerical Solution of Partial Differential Equations. Cambridge University Press, Cambridge.
[29]	Smith, G.D. (1978) Numerical Solution of Partial Differential Equations: Finite Difference Methods. Oxford University Press, Oxford.
[30]	Larsson, S. and Thomee, V. (2009) Partial Differential Equations with Numerical Methods. Springer Verlag, New York.
[31]	Varga, R.S. (1991) Functional Analysis and Approximation Theory in Numerical Analysis. SIAM, Philadelphia.
[32]	Arteaga-Arispe, J. and Guevara-Jordan, J.M. (2013) A New Conservative Method for the Unsteady Heat Equation. Mathematical and Computer Modelling, 57, 2302-2315. https://doi.org/10.1016/j.mcm.2011.08.024