We study symmetric finite volume element approximations for two-dimensional parabolic integrodifferential equations, arising in modeling of nonlocal reactive flows in porous media. It is proved that symmetric finite volume element approximations are convergent with optimal order in -norm. Numerical example is presented to illustrate the accuracy of our method.

1. Introduction

In this paper, we consider symmetric finite volume element discretizations of the following initial value problem for the operator equation for : where is a strongly elliptic differential operator and is a second-order elliptic differential operator in space. The operators and incorporate Dirichlet and Neumann boundary conditions. The problem (1) is an abstract form of an initial boundary value problem for a parabolic integrodifferential equation.

This model is very important in the transport of reactive and passive contaminats in aquifers, an area of active interdisciplinary research of mathematicians, engineers, and life scientists. From a mathematical point of view, the evolution of either a passive or reactive chemical within a velocity field exhibiting strong variation on many scales defies representation using classical Fickian theory. The evolution of a chemical in such a velocity field when modeled by Fickian-type theories leads to a dispersion tensor whose magnitude depends upon the time scales of observation. In order to avoid such difficulty, a new class of nonlocal models of transport have been derived. In this case, the constitutive relations involve either integrals or higher order derivatives, which take multiscales into consideration. We refer to [1, 2] for the derivation of the mathematical models and for the precise hypotheses and analysis.

Mathematical formulations of this kind also arise naturally in various engineering models, such as nonlocal reactive transport in underground water flows in porous media [3], heat conduction, radioactive nuclear decay in fluid flows [4], non-Newtonian fluid flows, or viscoelastic deformations of materials with memory (in particular polymers) [5], semiconductor modeling [6], and biotechnology. One very important characteristic of all these models is that they all express a conservation of a certain quantity (mass, momentum, heat, etc.) in any moment for any subdomain. This in many applications is the most desirable feature of the approximation method when it comes to the numerical solution of the corresponding initial boundary value problem.

This type of equations has been extensively treated by finite element, finite difference, and collocation methods in the last years [712].

The finite element method conserves the flux approximately; therefore, in the asymptotic limit (i.e., when the grid step-size tends to zero) it produces adequate results. However, this could be a disadvantage when relatively coarse grids are used. Perhaps the most important property of the finite volume method is that it exactly conserves the approximate flux (heat, mass, etc.) over each computational cell. This important property combined with adequate accuracy and ease of implementation has contributed to the recent renewed interest in the method. The theoretical framework and the basic tools for the analysis of the finite volume element methods have been developed in the last decade (see, e.g., [1319]). Reference [20] has given the finite volume element approximations of the problem (1). But in a general case the coefficient matrix of the linear system obtained from the finite volume element method is not symmetric. In this paper, we study symmetric modified finite volume element approximations for two-dimensional parabolic integrodifferential equations, arising in modeling of nonlocal reactive flows in porous media. It is proved that symmetric modified finite volume element approximations are convergent with optimal order in -norm.

Throughout this paper we use (without or with subscript or superscript) to denote a generic constant independent of the discretization parameter.

2. Fully Discrete Symmetric Finite Volume Element Scheme

Consider the following initial boundary value problem: find such that where is a bounded polygonal domain in with boundary . is a real-valued symmetric and uniformly positive definite matrix, is a matrix. and are known functions, which are assumed to be smooth so that problem (2) has a unique solution in a certain Sobolev space.

The problem (2) can be written in the form (1) by introducing the operators and : for all and for any .

We assume that is a convex polygonal domain. The domain is split into triangular finite elements . The elements are considered to be a closed set, and the triangulation is denoted by . Then and denotes all nodes or vertices: In order to accommodate Dirichlet boundary conditions, we will also need the set of vertices that are internal to , denoted by , that is, . For a given vertex , we define by the index set of all neighbors of in .

For a given triangulation , we construct a dual mesh based upon whose elements are called control volumes. In the finite volume methods, there are various ways to introduce the control volumes. Almost all approaches can be described in the following general scheme: in each triangle a point is selected; similarly, on each of the three edges of a point is selected; then is connected with the points by straight lines (see Figure 1).

Thus, around each vertex , we associate a control volume , which consists of the union of the subelements , which have as a vertex. Also, let denote the interface of two control volumes and . We call the partition regular or quasiuniform, if there exists a positive constant such that , for all . Here, is the maximal diameter of all elements . In the first (and most popular) control volume partition, the point is chosen to be the medicenter (the center of gravity or centroid) of the finite element and the points are chosen to be the midpoints of the edges of . This type of control volume can be introduced for any finite element partition and leads to relatively simple calculations. Besides, if the finite element partition is locally regular, that is, there is a constant such that ,  , for all elements , then the finite volume partition is also locally regular.

In this paper, we will also use the construction of the control volumes in which the point is the circumcenter of the element . Then obviously, are the perpendicular bisectors of the three edges of . This construction requires that all finite elements are triangles of acute type, which we will assume whenever such a triangulation is used.

We define the linear finite element space : and its dual volume element space : Obviously, and , where are the standard nodal linear basis functions associated with the node and are the characteristic functions of the volume . Let be the piecewise linear interpolation operator and let be the piecewise constant interpolation operator. That is,

The semidiscrete finite volume element approximation of (2) is a solution to the problem: find , for such that

Here, the bilinear forms and are defined by where denotes the outer-normal direction to the domain under consideration.

It is more convenient to rewrite (8) in the following form, which we use in the exposition:

Next, we define the fully discrete time stepping schemes. Let be a time-step size and ,  ,  and   is an approximation of , such that .

The backward Euler scheme is defined to be the solution of such that where are the weights, and the quadrature error is given by for any smooth functions and and its error satisfies

Then we present fully discrete symmetric finite volume element scheme as the following: find , such that where is defined by We note that is a modified term.

3. Some Auxiliary Results

Lemma 1 (see [17]). If the matrix is constant over each element , then

Proposition 2. Consider the following:

Proof. By the definition of (16), is piecewise constant on . Apply Lemma 1, and note that are constants on , we can get (18).

It can be deduced from Proposition 2 that the coefficient matrix of (14) is symmetric.

Remark 3. For each , if we select a point and replace (16) by ,  for all  ,   from the analysis given in the following we can get the same error estimates as in Theorem 12.
We define a space constant, for all , then introduce : by or ,  for all .

Lemma 4 (see [21]). Define , . Then there exist two positive constants and , such that

Next we define some discrete norms on : where , the distance between and . Obviously, these norms are well defined for as well and .

Lemma 5 (see [20]). There exist two positive constants , independent of , such that

Lemma 6 (see [20]). Assume that the jumps (if any) of coefficient matrices and are aligned with the finite element partition , and over each element their entries are -functions. Then,(a) there are positive constants and , independent of and , such that for all , (b) for is fixed there is a constant independent of  and , such that for ,

4. -Norm Error Estimate

Now we introduce linear functions , which are used in the error analysis of the finite volume element method: The following estimate is a simple consequence of the Bramble-Hilbert lemma.

Lemma 7 (see [20]). If , then there is a positive constant , independent of  , such that for ,

Lemma 8 (see [20]). (a) If , then there exists a positive constant , independent of and , such that
(b) If , then, for fixed there is a constant , independent of and , such that for For any fixed , one can define the Ritz projection function and the operator , such that

Remark 9. If the partition is regular (quasiuniform) and is -regular, then

However, these estimates for the Ritz projection lead to suboptimal error estimates for the finite volume element solution of the integrodifferential equation. In order to obtain optimal order estimates, we need a projection that also takes into account the integral term. This type of projection was called by Cannon and Lin [7] the Ritz-Volterra projection and has been used in the analysis of the finite element method for integrodifferential equations.

We define the Ritz-Volterra projection of a function defined on the cylinder . The Ritz-Volterra projection is defined for by

Lemma 10 (see [20]). Assume that for fixed, there is a constant independent of and , such that for all ,

Now we consider an -estimate for the Ritz-Volterra projection. This estimate is optimal with respect to the order of convergence, but requires -regularity of the solution. Therefore, it is suboptimal with respect to the regularity of the solution and can be useful for close to 1. Namely, we have the following result.

Lemma 11 (see [20]). Assume that, for some , . Then for fixed there exists a positive constant independent of and , such that, for all ,

Then we prove the following -norm error estimate.

Theorem 12. Let and be the solution of (2) and (14), respectively. Then we have a constant independent of , , and , such that

Proof. Let , from (2), (14), and (30), we have Let , Seting in (34), we get We estimate the terms on the right-hand side: Let , and use numerical quadrature error estimates to get Thus, multiplying each term on both sides of (38) by , summing over , and employing Gronwall’s lemma, we obtain And ; hence (33) follows from the above analysis and Lemma 11.

5. Numerical Experiments

In this section, we present a numerical example for solving the problem (2) by using the symmetric modified finite volume element scheme presented in Section 2. Let , , be Delaunay triangulation generated by EasyMesh [22] over with mesh size as shown in Figure 2 and time step be .

We consider the case of , , , , , , the exact solution , and the initial function . Denote the numerical solution of by and set with fixed ratio . The numerical results are presented in Tables 1 and 2. It is observed that the results support our theory.


This work is supported by the Excellent Young and Middle-aged Scientists Research Fund of Shandong Province (no. 2008BS09026); National Natural Science Foundation of China (no. 11171193); Natural Science Foundation of Shandong Province (no. ZR2011AM016).