A Fully Discrete Symmetric Finite Volume Element Approximation of Nonlocal Reactive Flows in Porous Media

Abstract

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 𝑢=𝑢(𝑡): 𝑢_𝑡+𝒜𝑢+^𝑡₀ℬ(𝑡,𝑠)𝑢(𝑠)𝑑𝑠=𝑓(𝑡),𝑢(0)=𝑢₀, (1) 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 [7–12].

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., [13–19]). 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 𝑢_𝑡−∇⋅(𝐴∇𝑢)−^𝑡₀((𝑠()(],𝑢()()(],𝑢()∇⋅𝐵∇𝑢))𝑑𝑠=𝑓,𝑥,𝑡∈Ω×0,𝑇𝑥,𝑡=0,𝑥,𝑡∈𝜕Ω×0,𝑇𝑥,0=𝑢₀(𝑥),𝑥∈Ω, (2) where Ω is a bounded polygonal domain in 𝑅² with boundary 𝜕Ω. 𝐴={𝑎_𝑖𝑗(𝑥)} is a 2×2 real-valued symmetric and uniformly positive definite matrix, 𝐵={𝑏_𝑖𝑗(𝑥,𝑡,𝑠)} is a 2×2 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 ℬ∶𝐻¹₀(Ω)→𝐻⁻¹(Ω): (𝒜𝑢,𝑣)=_Ω𝐴(𝑥)()=∇𝑢⋅∇𝑣𝑑𝑥,ℬ𝑢,𝑣_Ω𝐵(𝑥,𝑡,𝑠)∇𝑢⋅∇𝑣𝑑𝑥, (3) for all 𝑣∈𝐻¹₀(Ω) and for any 𝑡,𝑠∈(0,𝑇).

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: 𝑁_ℎ=𝑝∶𝑝isavertexofelement,𝐾∈𝒯_ℎand𝑝∈Ω. (4) 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).

Figure 1 

A control volume.

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 𝐶>0 such that 𝐶⁻¹ℎ²≤meas(𝑉_𝑖)≤𝐶ℎ², 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 𝐶ℎ²_𝐾≤meas(𝐾)≤ℎ²_𝐾,  ℎ_𝐾=diam(𝐾), 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 𝑆_ℎ: 𝑆_ℎ=𝑣∈𝐶(Ω)∶𝑣∣_𝐾islinearforall𝐾∈𝑇_ℎand𝑣∣_𝜕Ω=0, (5) and its dual volume element space 𝑆^∗_ℎ: 𝑆^∗_ℎ=𝑣∈𝐿²(Ω)∶𝑣∣_𝑉isconstantforall𝑉∈𝑇^∗_ℎand𝑣∣_𝜕Ω=0. (6) Obviously, 𝑆_ℎ=span{𝜙_𝑖(𝑥)∶𝑥_𝑖∈𝑁⁰_ℎ} and 𝑆^∗_ℎ=span{𝜒_𝑖(𝑥):𝑥_𝑖∈𝑁⁰_ℎ}, 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, 𝐼_ℎ𝑢=_{𝑥𝑖∈𝑁ℎ}𝑢𝑥_𝑖𝜙_𝑖(𝑥),𝐼^∗_ℎ𝑢=_{𝑥𝑖∈𝑁ℎ}𝑢𝑥_𝑖𝜒_𝑖(𝑥). (7)

The semidiscrete finite volume element approximation 𝑢_ℎ of (2) is a solution to the problem: find 𝑢_ℎ(𝑡)∈𝑆_ℎ, for 𝑡>0 such that 𝑢_ℎ,𝑡,𝑣_ℎ+𝐴𝑢_ℎ,𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢_ℎ(𝑠),𝑣_ℎ=𝑑𝑠𝑓,𝑣_ℎ,𝑣_ℎ∈𝑆^∗_ℎ,𝑢_ℎ(0)=𝑢_0,ℎ∈𝑆_ℎ. (8)

Here, the bilinear forms 𝐴(𝑢,𝑣) and 𝐵(𝑡,𝑠;𝑢,𝑣) are defined by 𝐴(𝑢,𝑣)=−_{𝑥𝑖∈𝑁ℎ}𝑣_{𝑖𝜕𝑉𝑖}𝐴(𝑥)()∇𝑢⋅𝐧𝑑𝑥,𝑢,𝑣∈𝑆_ℎ×𝑆^∗_ℎ,𝐵(𝑡,𝑠;𝑢,𝑣)=−_{𝑥𝑖∈𝑁ℎ}𝑣_{𝑖𝜕𝑉𝑖}𝐵()()𝑥,𝑡,𝑠∇𝑢⋅𝐧𝑑𝑥,𝑢,𝑣∈𝑆_ℎ×𝑆^∗_ℎ, (9) 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: 𝑢_ℎ,𝑡,𝐼^∗_ℎ𝑣_ℎ+𝐴𝑢_ℎ,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢_ℎ(𝑠),𝐼^∗_ℎ𝑣_ℎ=𝑑𝑠𝑓,𝐼^∗_ℎ𝑣_ℎ,𝑣_ℎ∈𝑆_ℎ. (10)

Next, we define the fully discrete time stepping schemes. Let Δ𝑡>0 be a time-step size and 𝑡_𝑛=𝑛Δ𝑡,  𝑢^𝑛=𝑢(𝑡_𝑛),  and  𝜕_𝑡𝑢^𝑛=(𝑢^𝑛−𝑢^𝑛−1)/Δ𝑡.𝑢_ℎ(0)∈𝑆_ℎ is an approximation of 𝑢₀, such that ‖𝑢_ℎ(0)−𝑢₀‖≤𝐶ℎ²‖𝑢₀‖₂.

The backward Euler scheme is defined to be the solution of 𝑢^𝑛_ℎ∈𝑆_ℎ such that 𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1Δ𝑡,𝐼^∗_ℎ𝑣_ℎ+𝐴𝑢^𝑛_ℎ,𝐼^∗_ℎ𝑣_ℎ+^𝑛−1_𝑘=0𝜔_𝑛,𝑘𝐵𝑡_𝑛,𝑡_𝑘;𝑢^𝑘_ℎ,𝐼^∗_ℎ𝑣_ℎ=𝑓^𝑛,𝐼^∗_ℎ𝑣_ℎ, (11) where 𝜔_𝑛,𝑗 are the weights, and the quadrature error is given by for any smooth functions 𝑔 and 𝑀 and its error 𝑞_𝑛𝑔=^𝑡𝑛₀𝑀𝑡_𝑛,𝑠𝑔(𝑠)𝑑𝑠−^𝑛−1_𝑘=0𝜔_𝑛,𝑗𝑀_𝑛,𝑗𝑔𝑡_𝑗 (12) satisfies 𝑞_𝑛𝑔≤𝐶Δ𝑡^𝑡𝑛₀𝑔+𝑔^′𝑑𝑡. (13)

Then we present fully discrete symmetric finite volume element scheme as the following: find 𝑢^𝑛_ℎ∈𝑆^ℎ, such that 𝜕_𝑡𝑢^𝑛_ℎ,𝐼^∗_ℎ𝑣_ℎ+𝐴𝑢^𝑛_ℎ,𝐼^∗_ℎ𝑣_ℎ+^𝑛−1_𝑘=0𝜔_𝑛,𝑘𝐵𝑡_𝑛,𝑡_𝑘;𝑢^𝑘_ℎ,𝐼^∗_ℎ𝑣_ℎ=𝑓^𝑛,𝐼^∗_ℎ𝑣_ℎ+𝐴𝑢_ℎ^𝑛−1,𝐼^∗_ℎ𝑣_ℎ−𝐴𝑢_ℎ^𝑛−1,𝐼^∗_ℎ𝑣_ℎ, (14) where 𝐴(𝑢,𝑣) is defined by 𝐴(𝑢,𝑣)=−_{𝑥𝑖∈𝑁ℎ}𝑣_{𝑖𝜕𝑉𝑖}𝐴(𝑥)()∇𝑢⋅𝐧𝑑𝑥,𝑢,𝑣∈𝑆_ℎ×𝑆^∗_ℎ,𝐴∣_𝐾=𝐴_𝐾,𝐴_𝐾=1meas(𝐾)_𝐾𝐴(𝑥)𝑑𝑥. (15)(16) We note that 𝐴(𝑢_ℎ^𝑛−1,𝐼^∗_ℎ𝑣_ℎ)−𝐴(𝑢_ℎ^𝑛−1,𝐼^∗_ℎ𝑣_ℎ) is a modified term.

3. Some Auxiliary Results

Lemma 1 (see [17]). If the matrix 𝐴=(𝑎_𝑖𝑗)_2×2 is constant over each element 𝐾∈𝒯_ℎ, then −_{𝑥𝑖∈𝑁ℎ}𝑣_{𝑖𝜕𝑉𝑖}𝐴∇𝑢⋅𝐧𝑑𝑠=_Ω𝐴∇𝑢⋅∇𝑣𝑑𝑥,∀𝑢,𝑣∈𝑆_ℎ. (17)

Proposition 2. Consider the following: −_{𝑥𝑖∈𝑁ℎ}𝑣_{𝑖𝜕𝑉𝑖}𝐴∇𝑢^𝑛_ℎ⋅𝐧𝑑𝑠=_Ω𝐴∇𝑢^𝑛_ℎ⋅∇𝑣𝑑𝑥,∀𝑢_ℎ,𝑣∈𝑆_ℎ. (18)

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 𝑃⁰_ℎ𝑣=(1/|𝐾|)∫_𝐾𝑣𝑑𝑥,  for all 𝐾∈𝒯_ℎ.

Lemma 4 (see [21]). Define 𝐴_ℎ=(𝑃⁰_ℎ𝑎_𝑖𝑗)_2×2, 𝑎_ℎ(𝑢_ℎ,𝑣)=∫_Ω𝐴_ℎ∇𝑢_ℎ⋅∇𝑣𝑑𝑥. Then there exist two positive constants 𝐶₁ and 𝛽, such that 𝑎_ℎ(𝑤,𝑣)≤𝐶₁‖𝑤‖₁‖𝑣‖₁,∀𝑤,𝑣∈𝑆_ℎ,𝑎_ℎ(𝑣,𝑣)≥𝛽‖𝑣‖²₁,∀𝑣∈𝑆_ℎ. (19)

Next we define some discrete norms on 𝑆_ℎ: 𝑢_ℎ²_0,ℎ=𝑢_ℎ,𝑢_ℎ0,ℎ=𝐼^∗_ℎ𝑢_ℎ,𝐼^∗_ℎ𝑢_ℎ,𝑢_ℎ²_1,ℎ=_{𝑥𝑖∈𝑁ℎ𝑥𝑗(𝑖)∈Π}meas𝑉_𝑖𝑢_𝑖−𝑢_𝑗𝑑_𝑖𝑗²,𝑢_ℎ²_1,ℎ=𝑢_ℎ²_0,ℎ+𝑢_ℎ²_1,ℎ,𝑢_ℎ²=𝑢_ℎ,𝐼^∗_ℎ𝑢_ℎ, (20) where 𝑑_𝑖𝑗=𝑑(𝑥_𝑖,𝑥_𝑗), the distance between 𝑥_𝑖 and 𝑥_𝑗. Obviously, these norms are well defined for 𝑢_ℎ∈𝑆^∗_ℎ as well and ‖𝑢_ℎ‖_0,ℎ=|||𝑢_ℎ|||.

Lemma 5 (see [20]). There exist two positive constants 𝐶₂,𝐶₃, independent of ℎ, such that 𝐶₂‖𝑣‖_0,ℎ≤𝑣_ℎ≤𝐶₃‖𝑣‖_0,ℎ,∀𝑣_ℎ∈𝑆_ℎ,𝐶₂𝑣_ℎ≤𝑣_ℎ≤𝐶₃𝑣_ℎ,∀𝑣_ℎ∈𝑆_ℎ,𝐶₂𝑣_ℎ1,ℎ≤𝑣_ℎ1≤𝐶₃‖𝑣‖_1,ℎ,∀𝑣_ℎ∈𝑆_ℎ. (21)

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 𝑊^1,∞(𝐾)-functions. Then,(a) there are positive constants ℎ₀ and 𝑐₀,𝑐₁, independent of ℎ and 𝑢, such that for all 0<ℎ≤ℎ₀, 𝐴𝑢_ℎ,𝐼^∗_ℎ𝑣_ℎ≤𝑐₁𝑢_ℎ1,ℎ𝑣_ℎ1,ℎ,∀𝑢_ℎ,𝑣_ℎ∈𝑆_ℎ,𝐴𝑣_ℎ,𝐼^∗_ℎ𝑣_ℎ≥𝑐₀𝑣_ℎ²_1,ℎ,∀𝑣_ℎ∈𝑆_ℎ; (22)(b) for 𝑇>0 is fixed there is a constant 𝐶=𝐶(𝑇)>0 independent of  ℎ and 𝑢, such that for 0<𝑡≤𝑇, 𝐵𝑡,𝑠;𝑢_ℎ,𝐼^∗_ℎ𝑣_ℎ≤𝐶𝑢_ℎ1,ℎ𝑣_ℎ1,ℎ,∀𝑢_ℎ,𝑣_ℎ∈𝑆_ℎ. (23)

4. 𝐿²-Norm Error Estimate

Now we introduce linear functions 𝑙_𝑖𝑗(𝑢), which are used in the error analysis of the finite volume element method: 𝑙_𝑖𝑗(𝑢)=−_𝛾𝑖𝑗𝐴∇𝐼_ℎ𝑢−𝑢⋅𝐧𝑑𝑥,𝛾_𝑖𝑗=𝑉_𝑖∩𝑉_𝑗. (24) The following estimate is a simple consequence of the Bramble-Hilbert lemma.

Lemma 7 (see [20]). If 𝑢∈𝐻²(Ω), then there is a positive constant 𝐶>0, independent of  ℎ, such that for 𝑒_𝑖𝑗=∪{𝐾∣𝐾∩𝛾_𝑖𝑗≠∅,𝐾∈𝒯_ℎ}, 𝑙_𝑖𝑗(𝑢)≤𝐶ℎ‖𝐴‖_0,∞|𝑢|_{2,𝑒𝑖𝑗}. (25)

Lemma 8 (see [20]). (a) If 𝑢∈𝐻²(Ω), then there exists a positive constant 𝐶>0, independent of ℎ and 𝑢, such that 𝐴𝑢−𝐼_ℎ𝑢,𝐼^∗_ℎ𝑣_ℎ≤𝐶ℎ‖𝑢‖₂𝑣_ℎ1,ℎ,𝑣_ℎ∈𝑆_ℎ. (26)
(b) If 𝑢(⋅)∈𝐿^∞(𝐻²), then, for 𝑇>0 fixed there is a constant 𝐶=𝐶(𝑇)>0, independent of ℎ and 𝑢, such that for 0<𝑡≤𝑇𝐵𝑡,𝑠;𝑢−𝐼_ℎ𝑢,𝐼^∗_ℎ𝑣_ℎ≤𝐶ℎ‖𝑢‖₂𝑣_ℎ1,ℎ,𝑣_ℎ∈𝑆_ℎ. (27) For any fixed 𝑡>0, one can define the Ritz projection function 𝑢(𝑥,𝑡) and the operator 𝑅_ℎ∶𝐻¹₀∩𝐻²→𝑆_ℎ, such that 𝐴𝑢−𝑅_ℎ𝑢,𝐼^∗_ℎ𝑣_ℎ=0,∀𝑣_ℎ∈𝑆_ℎ. (28)

Remark 9. If the partition 𝒯^∗_ℎ is regular (quasiuniform) and 𝑢 is 𝐻²-regular, then 𝑢−𝑅_ℎ𝑢₁≤𝐶ℎ‖𝑢‖₂. (29)

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 Ω×[0,𝑇]. The Ritz-Volterra projection 𝑉_ℎ∶𝐿^∞(𝐻¹₀∩𝐻²)→𝐿^∞(𝑆_ℎ) is defined for 𝑡≥0 by 𝐴𝑢−𝑉_ℎ𝑢,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢(𝑠)−𝑉_ℎ𝑢(𝑠),𝐼^∗_ℎ𝑣_ℎ𝑑𝑠=0,∀𝑣_ℎ∈𝑆_ℎ. (30)

Lemma 10 (see [20]). Assume that 𝑢∈𝐿^∞(𝐻¹₀∩𝐻²) for 𝑇>0 fixed, there is a constant 𝐶=𝐶(𝑇)>0 independent of ℎ and 𝑢, such that for all 0<𝑡≤𝑇, 𝑢−𝑉_ℎ𝑢(𝑡)₁≤𝐶ℎ‖𝑢(𝑡)‖₂+^𝑡₀‖𝑢(𝑠)‖₂𝑑𝑠. (31)

Now we consider an 𝐿²-estimate for the Ritz-Volterra projection. This estimate is optimal with respect to the order of convergence, but requires 𝑊^3,𝑝-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 𝑝>1, 𝑢∈𝐿^∞(𝑊^3,𝑝(Ω)). Then for 𝑇>0 fixed there exists a positive constant 𝐶=𝐶(𝑇)>0 independent of ℎ and 𝑢, such that, for all 0<𝑡≤𝑇, 𝑢−𝑉_ℎ𝑢(𝑡)≤𝐶ℎ²‖𝑢(𝑡)‖_3,𝑝+^𝑡₀‖𝑢(𝑠)‖_3,𝑝𝑑𝑠. (32)

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 𝐶>0 independent of ℎ, Δ𝑡, and 𝑢, such that 𝑢^𝑛−𝑢^𝑛_ℎ≤𝐶ℎ²𝑢_03,𝑝+^𝑡𝑛₀𝑢_𝑡3,𝑝𝑑𝑠+𝐶Δ𝑡^𝑡𝑛₀𝑢_𝑡𝑡+𝑢_𝑡1𝑑𝑠. (33)

Proof. Let 𝜃^𝑛=𝑢^𝑛_ℎ−𝑉_ℎ𝑢^𝑛,𝜌^𝑛=𝑢^𝑛−𝑉_ℎ𝑢^𝑛, from (2), (14), and (30), we have 𝜕_𝑡𝜃^𝑛,𝐼^∗_ℎ𝑣_ℎ+𝐴𝜃^𝑛,𝐼^∗_ℎ𝑣_ℎ+^𝑛−1_𝑘=0𝑤_𝑛,𝑘𝐵𝑡_𝑛,𝑡_𝑘;𝜃^𝑘,𝐼^∗_ℎ𝑣_ℎ=𝑢^𝑛_𝑡−𝜕_𝑡𝑢^𝑛,𝐼^∗_ℎ𝑣_ℎ+𝜕_𝑡𝜌^𝑛,𝐼^∗_ℎ𝑣_ℎ+𝐴𝜃^𝑛−1,𝐼^∗_ℎ𝑣_ℎ−𝐴𝜃^𝑛−1,𝐼^∗_ℎ𝑣_ℎ−Δ𝑡𝐴𝜕_𝑡𝑉_ℎ𝑢^𝑛,𝐼^∗_ℎ𝑣_ℎ−𝐴𝜕_𝑡𝑉_ℎ𝑢^𝑛,𝐼^∗_ℎ𝑣_ℎ+^𝑡𝑛₀𝐵𝑡_𝑛,𝑠;𝑉_ℎ𝑢^𝑛(𝑠),𝐼^∗_ℎ𝑣_ℎ−𝑑𝑠_{𝑛−1𝑘=0}𝑤_𝑛,𝑘𝐵𝑡_𝑛,𝑡_𝑘;𝑉_ℎ𝑢^𝑘,𝐼^∗_ℎ𝑣_ℎ. (34) Let 𝜏^𝑛=𝑢^𝑛_𝑡−𝜕_𝑡𝑢^𝑛+𝜕_𝑡𝜌^𝑛, 𝐿^𝐴−𝐴𝜃^𝑛−1,𝐼^∗_ℎ𝑣_ℎ=𝐴𝜃^𝑛−1,𝐼^∗_ℎ𝑣_ℎ−𝐴𝜃^𝑛−1,𝐼^∗_ℎ𝑣_ℎ,𝐿^𝐴−𝐴𝜕_𝑡𝑉_ℎ𝑢^𝑛,𝐼^∗_ℎ𝑣_ℎ=𝐴𝜕_𝑡𝑉_ℎ𝑢^𝑛,𝐼^∗_ℎ𝑣_ℎ−𝐴𝜕_𝑡𝑉_ℎ𝑢^𝑛,𝐼^∗_ℎ𝑣_ℎ,𝑞^𝑛𝑉_ℎ𝑢,𝐼^∗_ℎ𝑣_ℎ=^𝑡𝑛₀𝐵𝑡_𝑛,𝑠;𝑉_ℎ𝑢^𝑛(𝑠),𝐼^∗_ℎ𝑣_ℎ−𝑑𝑠_{𝑛−1𝑘=0}𝑤_𝑛,𝑘𝐵𝑡_𝑛,𝑡_𝑘;𝑉_ℎ𝑢^𝑘,𝐼^∗_ℎ𝑣_ℎ. (35) Seting 𝑣_ℎ=𝜃^𝑛 in (34), we get 𝜕_𝑡𝜃^𝑛,𝐼^∗_ℎ𝜃^𝑛+𝐴𝜃^𝑛,𝐼^∗_ℎ𝜃^𝑛=𝜏^𝑛,𝐼^∗_ℎ𝜃^𝑛+𝐿^𝐴−𝐴𝜃^𝑛−1,𝐼^∗_ℎ𝜃^𝑛−Δ𝑡𝐿^𝐴−𝐴𝜕_𝑡𝑉_ℎ𝑢^𝑛,𝐼^∗_ℎ𝜃^𝑛+𝑞^𝑛𝑉_ℎ𝑢,𝐼^∗_ℎ𝜃^𝑛−^𝑛−1_𝑘=0𝑤_𝑛,𝑘𝐵𝑡_𝑛,𝑡_𝑘;𝜃^𝑘,𝐼^∗_ℎ𝜃^𝑛. (36) We estimate the terms on the right-hand side: 𝐿^𝐴−𝐴𝜃^𝑛−1,𝐼^∗_ℎ𝜃^𝑛≤𝐶ℎ‖𝐴‖_1,∞Ω∇𝜃^𝑛−1∇𝜃^𝑛‖𝐴‖𝑑𝑥≤𝐶ℎ_1,∞∇𝜃^𝑛−1∇𝜃^𝑛≤𝐶𝜃^𝑛−12+𝜀𝜃^𝑛2₁,Δ𝑡𝐿^𝐴−𝐴𝜕_𝑡𝑉_ℎ𝑢^𝑛,𝐼^∗_ℎ𝜃^𝑛≤𝐶ℎΔ𝑡∇𝜕_𝑡𝑉_ℎ𝑢^𝑛∇𝜃^𝑛≤𝐶Δ𝑡∇𝜕_𝑡𝑉_ℎ𝑢^𝑛𝜃^𝑛≤𝐶Δ𝑡^𝑡𝑛_𝑡𝑛−1𝑢_𝑡²₁𝑑𝑠+𝐶𝜃^𝑛2,𝑞^𝑛𝑉_ℎ𝑢,𝐼^∗_ℎ𝜃^𝑛≤𝜀𝜃^𝑛2₁+𝐶Δ𝑡^𝑡𝑛₀𝐷_𝑡𝑉_ℎ𝑢(𝑠)₁𝑑𝑠²≤𝜀𝜃^𝑛2₁+𝐶Δ𝑡^𝑡𝑛₀𝑢_𝑡(𝑠)₁𝑑𝑠²,_{𝑛−1𝑘=0}𝑤_𝑛,𝑘𝐵𝑡_𝑛,𝑡_𝑘;𝜃^𝑘,𝐼^∗_ℎ𝜃^𝑛≤𝜀𝜃^𝑛2₁+𝐶^𝑛−1_𝑘=0Δ𝑡𝜃^𝑘2₁,𝜏_𝑛,𝐼^∗_ℎ𝜃^𝑛≤𝐶𝜃^𝑛2+𝐶Δ𝑡^𝑡𝑛_𝑡𝑛−1𝑢_𝑡𝑡(𝑠)²𝑑𝑠+𝐶ℎ⁴1Δ𝑡^𝑡𝑛_𝑡𝑛−1𝑢_𝑡(𝑠)²_3,𝑝𝑑𝑠. (37) Let 𝜀=𝐶₀/6, and use numerical quadrature error estimates to get 𝜃^𝑛2−𝜃^𝑛−122Δ𝑡+𝐶₀𝜃^𝑛2₁≤𝐶₀2𝜃^𝑛2₁+𝐶^𝑛−1_𝑘=0Δ𝑡𝜃^𝑘2₁+𝐶Δ𝑡^𝑡𝑛₀𝑢_𝑡(𝑠)₁𝑑𝑠²+𝐶𝜃^𝑛−12+𝐶𝜃^𝑛2+𝐶Δ𝑡^𝑡𝑛_𝑡𝑛−1𝑢_𝑡²₁𝑑𝑡+𝐶Δ𝑡^𝑡𝑛_𝑡𝑛−1𝑢_𝑡𝑡(𝑠)²𝑑𝑠+𝐶ℎ⁴1Δ𝑡^𝑡𝑛_𝑡𝑛−1𝑢_𝑡(𝑠)²_3,𝑝𝑑𝑠. (38) Thus, multiplying each term on both sides of (38) by Δ𝑡, summing over 𝑛, and employing Gronwall’s lemma, we obtain 𝜃^𝑛≤𝐶𝜃⁰+𝐶ℎ^2𝑡𝑛₀𝑢_𝑡3,𝑝𝑑𝑠+𝐶Δ𝑡^𝑡𝑛₀𝑢_𝑡𝑡+𝑢_𝑡1𝑑𝑠. (39) And ‖𝜃⁰‖≤𝐶ℎ²‖𝑢₀‖_3,𝑝; 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 Ω=(0,1)×(0,1), 𝑇=1, 𝑇_ℎ be Delaunay triangulation generated by EasyMesh [22] over Ω with mesh size ℎ as shown in Figure 2 and time step be Δ𝑡.

Figure 2 

Triangulation and circumcenter dual mesh.

We consider the case of 𝑎₁₁=𝑥₁+𝑥₂+3, 𝑎₁₂=𝑥₁+𝑥₂+4, 𝑎₁₂=𝑎₂₁=−(𝑥₁+𝑥₂), 𝑏₁₁=1, 𝑏₁₂=1/2, 𝑏₁₂=𝑏₂₁=0, the exact solution 𝑢(𝑥,𝑡)=(𝑡²+sin(𝜋𝑡))sin(𝜋𝑥₁)sin(𝜋𝑥₂), and the initial function 𝑢(𝑥,0)=0. Denote the numerical solution of 𝑢(𝑥,𝑇) by 𝑢_ℎ and set 𝛾_ℎ=‖𝑢(𝑥,𝑇)−𝑢_2ℎ‖/‖𝑢(𝑥,𝑇)−𝑢_ℎ‖ with fixed ratio Δ𝑡/ℎ². The numerical results are presented in Tables 1 and 2. It is observed that the results support our theory.

Table 1 

Error behavior for Δ𝑡/ℎ²=3.90625𝑒−4.

Table 2 

Error behavior for Δ𝑡/ℎ²=1.953125𝑒−4.

Acknowledgments

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