On the Finite Volume Element Method for Self-Adjoint Parabolic Integrodifferential Equations

Abstract

Finite volume element schemes for non-self-adjoint parabolic integrodifferential equations are derived and stated. For the spatially discrete scheme, optimal-order error estimates in 𝐿², 𝐻¹, and 𝐿^𝑝, 𝑊^1,𝑝 norms for 2≤𝑝<∞, are obtained. In this paper, we also study the lumped mass modification. Based on the Crank-Nicolson method, a time discretization scheme is discussed and related error estimates are derived.

1. Introduction

The main purpose of this paper is to study semidiscrete and full discrete finite volume element method (FVE) for parabolic integrodifferential equation of the form 𝑢_𝑡−∇⋅(𝐴(𝑥,𝑡)∇𝑢)−^𝑡₀(𝐵()(𝑠()(],(],𝑢()∇⋅𝑥,𝑡,𝑠∇𝑢))𝑑𝑠=𝑓𝑥,𝑡,inΩ×0,𝑇𝑢=0,on𝜕Ω×0,𝑇⋅,0=𝑢₀,inΩ, (1) where Ω is a bounded domain in ℝ^𝑑,  𝑑=2,3, with smooth boundary 𝜕Ω, and 𝑇<∞. Here 𝐴(𝑡), a non-self-adjoint second-order strongly elliptic, and 𝐵(𝑡,𝑠), an arbitrary second-order linear partial differential operator, both with coefficients depending smoothly on 𝑥 and 𝑡,  𝑓=𝑓(𝑥,𝑡) and 𝑢₀(𝑥) are known functions, which are assumed to be smooth and satisfy certain compatibility conditions for 𝑥∈Ω and 𝑡=0, so that (1) has a unique solution in certain Sobolev space. Problem (1) occurs in nonlocal reactive flows in porous media, viscoelasticity, and heat conduction through materials with memory.

Finite volume method is an important numerical tool for solving partial differential equations. It has been widely used in several engineering fields, such as fluid mechanics, heat and mass transfer, and petroleum engineering. The method can be formulated in the finite difference framework or in the Petrov-Galerkin framework. Usually, the former one is called finite volume method [1], marker and cell (MAC) method [2], or cell-centered method [3], and the latter one is called finite volume element method (FVE) [4–9], covolume method [10], or vertex-centered method [11, 12]. We refer to the monographs [13, 14] for general presentation of these methods. The most important property of FVE method is that it can preserve the conservation laws (mass, momentum, and heat flux) on each control volume. This important property, combined with adequate accuracy and ease of implementation, has attracted more people to do research in this field.

Recently, the authors in [8, 15] studied FVE method for general self-adjoint elliptic problems. The authors in [16] presented and analyzed the semidiscrete and full discrete symmetric finite volume schemes for a class of parabolic problems. In [6, 7] the authors have studied FVE for one- and two-dimensional parabolic integrodifferential equations and have obtained an optimal-order estimate in the 𝐿²-norm. The regularity required on the exact solution 𝑢 is 𝑊^3,𝑝 for 𝑝>1 which is higher when compared to that for finite element methods.

The aim of this paper is to study the convergence of FVE discretization for a nonself-adjoint parabolic integrodifferential problem (1). Both spatially discrete scheme and discrete-in-time scheme are analyzed, and optimal error estimates in 𝐿² and 𝐻¹ norms are proved using only energy method. We also explore and generalize that idea to develop the lumped mass modification and 𝐿^𝑝 estimates, 2≤𝑝<∞. Our analysis avoids the use of semigroup theory, and the regularity requirement on the solution is the same of that of finite element method. Furthermore, based on the Crank-Nicolson method the fully discrete scheme is analyzed and the related optimal error estimates are established.

This paper is organized as follows. In Section 2, we introduce some notations and present some preliminary materials to be used later. The Ritz-Volterra projection to finite volume element spaces is introduced and related estimates are carried out in Section 3. In Section 4, we estimate the error of the finite volume element approximations derived in the previous section. In Section 5, the lumped mass is presented and optimal estimates in 𝐿² and 𝐻¹ norms are obtained Finally, the Crank-Nicolson scheme is studied in Section 6.

2. Finite Volume Element Scheme

In this section, we introduce some material which will be used repeatedly hareafter. Throughout this paper, 𝐶 (with or without index) denotes a generic positive constant which does not depend on the spatial and time discretization parameters ℎ and 𝑘,  respectively.

2.1. Notations

We will use ‖⋅‖_𝑚 and |⋅|_𝑚 (resp., ‖⋅‖_𝑚,𝑝 and |⋅|_𝑚,𝑝) to denote the norm and seminorm of the Sobolev space 𝐻^𝑚(Ω) (resp., 𝑊^𝑚,𝑝(Ω)). The scalar product and norm in 𝐿²(Ω) are denoted by (⋅,⋅) and ‖⋅‖, respectively. Let 𝐻⁰₁(Ω) be the standard Sobolev subspace of 𝐻¹(Ω) of functions vanishing on 𝜕Ω.

The weak form of (1) is used to find 𝑢(⋅,𝑡)∶[0,𝑇]→𝐻¹₀(Ω),  such that 𝑢_𝑡,𝑣+𝐴(𝑡;𝑢,𝑣)+^𝑡₀𝐵((𝑠))=𝑡,𝑠;𝑢,𝑣𝑑𝑠𝑓,𝑣,∀𝑣∈𝐻¹₀(Ω),𝑢(0)=𝑢₀, (2) where 𝐴(𝑡;𝑢,𝑣)=_Ω𝐴()𝐵((𝑠))=𝑥,𝑡∇𝑢⋅∇𝑣,𝑡,𝑠;𝑢,𝑣_Ω𝐵(𝑥,𝑡,𝑠)∇𝑢(𝑠)⋅∇𝑣. (3)

Let 𝒯_ℎ be a decomposition of Ω into triangles (for the 2D case) or tetrahedral (for the 3D case) with ℎ=maxℎ_𝐾, where ℎ_𝐾 is the diameter of the element 𝐾∈𝒯_ℎ.

In order to describe the FVE method for solving problem (1), we will introduce a dual partition 𝒯^∗_ℎ based upon the original partition 𝒯_ℎ whose elements are called control volumes. We construct the control volumes in the same way as in [7, 17]. Let 𝑧_𝐾 be a point of 𝐾∈𝒯_ℎ. In the 2D case, on each edge 𝑒 of 𝐾 a point 𝑞_𝑒 is selected; then we connect 𝑧_𝐾 with line segments to 𝑞_𝑒; thus, partitioning 𝐾 into three quadrilaterals 𝐾_𝑧,  𝑧∈𝑍_ℎ(𝐾), where 𝑍_ℎ(𝐾) are the vertices of 𝐾. Then with each vertex 𝑧∈𝑍_ℎ=∪_{𝐾∈𝒯ℎ}𝑍_ℎ(𝐾) we associate a control volume 𝑉_𝑧, which consists of the union of the subregions 𝐾_𝑧, sharing the vertex  𝑧  (see Figure 1).

(a)

(b)

Figure 1 

(a) A sample region with blue lines indicating the corresponding control volume 𝑉𝑧. (b) A triangle 𝐾 partitioned into three subregions 𝐾𝑧.

Similarly, in the 3D case, on each of the four faces 𝑆_𝑖,𝑖=1,…,4, a point 𝑞_𝑆𝑖,𝑖=1,…,4, is selected, and on each of the six edges 𝑒 a point 𝑞_𝑒 is selected. On each of the two faces 𝑆₁ and 𝑆₂ of 𝐾 sharing an edge 𝑒, we connect 𝑞_𝑆𝑖,𝑖=1,2,  with 𝑞_𝑒 and with 𝑧_𝐾 by line segments, thus, partitioning 𝐾 into twelve tetrahedron 𝐾_𝑧,  𝑧∈𝑍_ℎ(𝐾) (see Figure 2). Then for 𝑧∈𝑍_ℎ the control volume 𝑉_𝑧 consists of the union of the subregions 𝐾_𝑧 sharing the vertex 𝑧. Thus, we finally obtain a group of control volumes covering the domain Ω, which is called the dual partition 𝒯^∗_ℎ of the triangulation 𝒯_ℎ. We denote by 𝑍⁰_ℎ the set of interior vertices and 𝑁_ℎ=#𝑍⁰_ℎ. For a vertex 𝑧_𝑖∈𝑍⁰_ℎ,  let Π(𝑖) be the index set of those vertices that along with 𝑧_𝑖 are in some element of 𝑇_ℎ (Figure 2).

Figure 2 

A tetrahedron 𝐾 partitioned into twelve subregions 𝐾𝑧.

There are various ways to introduce a regular dual partition 𝒯^∗_ℎ. In this paper, we will also use the construction of the control volumes in which we let 𝑧_𝐾 be the barycenter of 𝐾∈𝒯_ℎ. In the 2D case, we choose 𝑞_𝑒 to be the midpoint of the edge 𝑒 (see Figure 3).

(a)

(b)

Figure 3 

𝑧𝐾 is the barycenter of 𝐾, and 𝑞𝑒 is to be the midpoint of the edge 𝑒.

In the 3D case, we choose 𝑞_𝑒 to be the midpoint of the edge 𝑒 and 𝑞_𝑆𝑖 to be the barycenter of the face 𝑆_𝑖  (Figure 4).

Figure 4 

𝑞𝑒 is the midpoint of the edge 𝑒, and  𝑞𝑆𝑖 is the barycenter of the face 𝑆𝑖.

We call the partition 𝒯^∗_ℎ regular or quasiuniform, if there exists a positive 𝐶>0 such that 𝐶⁻¹ℎ²≤meas𝑉_𝑧≤𝐶ℎ²,∀𝑉_𝑧∈𝒯^∗_ℎ. (4)

If the finite element triangulation 𝒯_ℎ is quasiuniform, that is, there exists a positive 𝐶>0 such that 𝐶⁻¹ℎ²≤meas(𝐾)≤𝐶ℎ²,∀𝐾∈𝒯_ℎ, (5) then the dual partition𝒯^∗_ℎis also quasiuniform.

Based on the triangulation𝑇_ℎ, let 𝑆_ℎ be the standard conforming finite element space of piecewise linear functions, defined on the triangulation 𝑇_ℎ as follows: 𝑆_ℎ=𝑣∈𝒞(Ω)∶𝑣|_𝐾islinear∀𝐾∈𝑇_ℎ,and𝑣|_Γ=0. (6)

Let 𝐼_ℎ∶𝒞(Ω)→𝑆_ℎ be the standard interpolation operators, such that 𝐼_ℎ𝑢=_{𝑧∈𝑍0ℎ}𝑣_𝑧(𝑡)𝜑_𝑧(𝑥),∀𝑣∈𝑆_ℎ, (7) where {𝜑_𝑧}_{𝑧∈𝑍0ℎ} are the standard basis functions of 𝑆_ℎ and 𝑣_𝑧(𝑡)=𝑣(𝑡;𝑧).

2.2. Construction of the FVE Scheme

We formulate the FVE method for the problem (1) as follows: Given a 𝑧∈𝑍⁰_ℎ, integrating (1)¹ over the associated control volume 𝑉_𝑧 and applying Green’s formula, we obtain an integral conservation as follows form: _𝑉𝑧𝑢_𝑡−_𝜕𝑉𝑧𝐴(𝑥,𝑡)∇𝑢⋅𝑛𝑑𝑠−_𝜕𝑉𝑧𝐵()=𝑥,𝑡,𝑠∇𝑢⋅𝑛𝑑𝑠_𝑉𝑧𝑓(𝑥,𝑡), (8) where𝑛denotes the unit outer normal vector to 𝜕𝑉_𝑧.

Let 𝐼^∗_ℎ∶𝒞(Ω)→𝑆^∗_ℎ be the transfer operator defined by 𝐼^∗_ℎ𝑣=_{𝑧∈𝑍0ℎ}𝑣(𝑧)𝜒_𝑧,∀𝑣∈𝑆_ℎ, (9) where 𝑆^∗_ℎ=𝑣∈𝐿²(Ω)∶𝑣_𝑖𝑉𝑧isconstant,∀𝑧∈𝑍⁰_ℎ, (10) and 𝜒_𝑧 is the characteristic function of the control volume 𝑉_𝑧.

Now for 𝑡>0 and for an arbitrary 𝐼^∗_ℎ𝑣, we multiply (8) by 𝑣(𝑧), and sum over all 𝑧∈𝑍⁰_ℎ. Then the semidiscrete FVE approximation 𝑢_ℎ of (1) is a solution to the following problem: find 𝑢_ℎ(𝑡)∈𝑆_ℎ for 𝑡>0 such that 𝑢_ℎ𝑡,𝑣_ℎ+𝐴𝑡;𝑢_ℎ,𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢_ℎ(𝑠),𝑣_ℎ=𝑑𝑠𝑓,𝑣_ℎ𝑣_ℎ∈𝑆^∗_ℎ,𝑢_ℎ(0)=𝑢_0ℎ∈𝑆_ℎ. (11)

Here the bilinear forms 𝐴(𝑡;𝑢,𝑣) and 𝐵(𝑡,𝑠;𝑢,𝑣) are defined by 𝐴()=𝑡;𝑢,𝑣−∑_{𝑧∈𝑍0ℎ}𝑣_𝑖∫_𝜕𝑉𝑧𝐴(𝑥,𝑡)∇𝑢⋅𝑛𝑑𝑠,(𝑢,𝑣)∈𝐻¹₀∩𝐻²∪𝑆_ℎ×𝑆^∗_ℎ,∫_Ω𝐴(𝑥)∇𝑢⋅∇𝑣𝑑𝑥,(𝑢,𝑣)∈𝐻¹₀×𝐻¹₀,𝐵()=𝑡,𝑠;𝑢,𝑣−∑_{𝑧∈𝑍0ℎ}𝑣_𝑧∫_𝜕𝑉𝑧𝐵(𝑥,𝑡,𝑠)∇𝑢⋅𝑛𝑑𝑠,(𝑢,𝑣)∈𝐻¹₀∩𝐻²∪𝑆_ℎ×𝑆^∗_ℎ,∫_Ω𝐵(𝑥,𝑡,𝑠)∇𝑢⋅∇𝑣𝑑𝑥,(𝑢,𝑣)∈𝐻¹₀×𝐻¹₀. (12)

Let 𝑢_ℎ=^𝑁ℎ_𝑗=1𝛼_𝑧(𝑡)𝜑_𝑧(𝑥),𝛼(𝑡)=𝛼₁(𝑡),𝛼₂(𝑡),…,𝛼_𝑁ℎ(𝑡)^𝑇. (13) Then, we can rewrite scheme (11)¹ as systems of ordinary differential equations as follows: 𝑀_ℎ𝛼^′(𝑡)+𝐴_ℎ(𝑡)𝛼(𝑡)+^𝑡₀𝐵_ℎ(𝑡,𝑠)𝛼(𝑠)𝑑𝑠=𝐹_ℎ(𝑡). (14) Here 𝐹_ℎ(𝑡)=(𝑓₁(𝑡),𝑓₂(𝑡),…,𝑓_𝑁ℎ(𝑡))^𝑇,  the mass matrix 𝑀_ℎ={𝑀_ℎ𝑖𝑗}={(𝜑_𝑖,𝜒_𝑗)} is tridiagonal, and both 𝐴_ℎ(𝑡)={𝐴(𝑡;𝜑_𝑖,𝜒_𝑗)} and 𝐵_ℎ(𝑡,𝑠)={𝐵(𝑡,𝑠;𝜑_𝑖,𝜒_𝑗)} are positive definites.

In order to describe features of the bilinear forms defined in (11), we introduce some discrete norms on 𝑆_ℎ in the same way as in [7]: 𝑣_ℎ²_0,ℎ=𝑣_ℎ,𝑣_ℎ0,ℎ=𝐼^∗_ℎ𝑣_ℎ,𝐼^∗_ℎ𝑣_ℎ,𝑣_ℎ²_1,ℎ=_{𝑥𝑖∈𝑍0ℎ𝑥𝑗(𝑖)∈Π}meas𝑉_𝑖𝑣_𝑖−𝑣_𝑗𝑑_𝑖𝑗²,𝑣_ℎ²_1,ℎ=𝑣_ℎ²_0,ℎ+𝑣_ℎ²_1,ℎ,𝑣_ℎ²=𝑣_ℎ,𝐼^∗_ℎ𝑣_ℎ, (15) where 𝑑_𝑖𝑗=𝑑(𝑥_𝑖,𝑥_𝑗),  the distance between 𝑥_𝑖 and 𝑥_𝑗. Obviously, these norms are well defined for 𝑣_ℎ∈𝑆^∗_ℎ as well and ‖𝑣_ℎ‖_0,ℎ=|||𝑣_ℎ|||.

Hereafter we state the equivalence of discrete norms ‖⋅‖_0,ℎ and ‖⋅‖_1,ℎ with usual norms ‖⋅‖ and ‖⋅‖₁ on 𝑆_ℎ, respectively.

Lemma 1 (see [7]). There exist two positive constants 𝐶₀ and 𝐶₁ such that for all 𝑣_ℎ∈𝑆_ℎ, we have 𝐶₀𝑣_ℎ0,ℎ≤𝑣_ℎ≤𝐶₁𝑣_ℎ0,ℎ,∀𝑣_ℎ∈𝑆_ℎ,𝐶₀𝑣_ℎ≤𝑣_ℎ≤𝐶₁𝑣_ℎ,∀𝑣_ℎ∈𝑆_ℎ,𝐶₀𝑣_ℎ1,ℎ≤𝑣_ℎ1≤𝐶₁𝑣_ℎ1,ℎ,∀𝑣_ℎ∈𝑆_ℎ. (16)

Next we recall some properties of the bilinear forms (see [7, 18])

Lemma 2 (see [7]). There exist two positive constants 𝐶 and 𝐶₀ such that for all 𝑢_ℎ,𝑣_ℎ∈𝑆_ℎ,  we have 𝐴𝑢_ℎ,𝐼^∗_ℎ𝑣_ℎ≤𝐶𝑢_ℎ1𝑣_ℎ1,∀𝑢_ℎ,𝑣_ℎ∈𝑆_ℎ,𝐴𝑣_ℎ,𝐼^∗_ℎ𝑣_ℎ≥𝐶₀𝑣_ℎ²₁,∀𝑣_ℎ∈𝑆_ℎ. (17)

The following lemmas are proved in [3, 7], which give the key feature of the bilinear forms in the FVE method.

Lemma 3 (see [3]). Assume that 𝜑∈𝑊₀^1,𝑝. Then, one has 𝐴𝑡;𝜑,𝑣_ℎ−𝐴𝑡;𝜑,𝐼^∗_ℎ𝑣_ℎ=_{𝐾∈𝜏ℎ𝜕𝐾}𝐴(𝑡)∇𝜑⋅𝐧𝑣_ℎ−𝐼^∗_ℎ𝑣_ℎ𝑑𝑠−_{𝐾∈𝜏ℎ𝐾}∇⋅𝐴(𝑡)∇𝜑𝑣_ℎ−𝐼^∗_ℎ𝑣_ℎ𝑑𝑠,∀𝑣_ℎ∈𝑆_ℎ. (18) The aforementioned identity holds true when 𝐴(⋅,⋅) is replaced by 𝐵(𝑡,𝑠;⋅,⋅).

Lemma 4 (see [3]). Assume that 𝜑∈𝑆_ℎ. Then, one has 𝐴𝑡;𝜑,𝜒−𝐴𝑡;𝜑,𝐼^∗_ℎ𝜒≤𝐶ℎ𝜑_1,𝑝𝜒_1,𝑞. (19) Furthermore, for 𝜑∈𝑊₀^1,𝑝∩𝑊^2,𝑝,  we have 𝐴𝑡;𝜑,𝜒−𝐴𝑡;𝜑,𝐼^∗_ℎ𝜒≤𝐶ℎ𝜑_2,𝑝𝜒_1,𝑞. (20)

3. Ritz-Volterra Projection and Related Estimates

Following [7, 19, 20], we define the Ritz-Volterra projection 𝑉_ℎ(𝑡)∶𝐻¹₀→𝑆_ℎ as follows: 𝐴𝑡;𝑢−𝑉_ℎ𝑢,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢(𝑠)−𝑉_ℎ𝑢(𝑠),𝐼^∗_ℎ𝑣_ℎ𝑑𝑠=0,𝑡>0,∀𝑣_ℎ∈𝑆_ℎ. (21)

This 𝑉_ℎ(𝑡) is an elliptic projection with memory of 𝑢 into 𝑆^∗_ℎ. It is easy to see that (21) is actually a system of integral equations of Volterra type. In fact, if 𝑉_ℎ(𝑡)𝑢=∑^𝑁ℎ_𝑗=1𝛼_𝑗(𝑡)𝜑_𝑗(𝑥), then (21) can be rewritten as 𝐴_ℎ(𝑡)𝛼(𝑡)+^𝑡₀𝐵_ℎ(𝑡,𝑠)𝛼(𝑠)𝑑𝑠=𝐹_ℎ(𝑡), (22) where 𝐴_ℎ(𝑡),  𝐵_ℎ(𝑡,𝑠) are matrices and 𝛼(𝑡),  𝐹_ℎ(𝑡) are vectors, defined via 𝛼(𝑡)=𝛼₁(𝑡),𝛼₂(𝑡),…,𝛼_𝑁ℎ(𝑡)^𝑇,𝐹_ℎ𝑘(𝑡)=𝐴𝑡;𝑢,𝜒_𝑘+^𝑡₀𝐵𝑡,𝑠;𝑢(𝑠),𝜒_𝑘𝑑𝑠,𝑘=1,2,…,𝑁_ℎ,𝐴_ℎ(𝑡)=𝐴𝑡,𝜑_𝑘(𝑥),𝜒_𝑙,𝐵_ℎ(𝑡,𝑠)=𝐵𝑡,𝑠;𝜑_𝑘(𝑥),𝜒_𝑙. (23)

From the positivity of 𝐴 (Lemma 2) and the linearity of (22), we see that the system (22) possesses a unique solution 𝛼(𝑡). Consequently, 𝑉_ℎ(𝑡)𝑢 in (21) is well defined.

Set 𝜌=𝑢−𝑉_ℎ(𝑡)𝑢. The following lemma was proved in [7], which shows the 𝐻¹ error estimate for 𝜌 and its temporal derivative.

Lemma 5 (see [7]). Assume that 𝐷^𝑛_𝑡𝑢∈𝐿^∞(𝐻¹₀∩𝐻²) for all 0≤𝑛≤𝑘, for some integer 𝑘≥0.  Then, for 𝑇>0 fixed there is a constant 𝐶=𝐶(𝑇;𝑘)>0, independent of ℎ and 𝑢, such that for all 0≤𝑛≤𝑘 and 0<𝑡<𝑇, 𝜌(𝑡)₁≤𝐶ℎ‖𝑢‖₂+^𝑡₀‖𝑢‖₂,𝐷𝑑𝑠_𝑛𝑡𝜌(𝑡)₁≤𝐶ℎ^𝑛_𝑖=0𝐷^𝑖_𝑡𝑢₂+^𝑡₀‖𝑢‖₂𝑑𝑠. (24)

Now we establish 𝐿² error estimate for 𝜌 and its temporal derivative which improves Theorem  2.2 in [7]. This estimate is optimal with respect to the order.

Lemma 6. Assume that, for some integer 𝑘≥0,  𝐷^𝑛_𝑡𝑢∈𝐿^∞(𝐻¹₀∩𝐻²)  for all 0≤𝑛≤𝑘. Then, for 𝑇>0 fixed there is a constant 𝐶=𝐶(𝑇;𝑘)>0, independent of ℎ and 𝑢, such that for all 0≤𝑛≤𝑘 and 0<𝑡<𝑇, 𝜌(𝑡)≤𝐶ℎ²‖𝑢‖₂+^𝑡₀‖𝑢‖₂,𝐷𝑑𝑠_𝑛𝑡𝜌(𝑡)≤𝐶ℎ^2𝑛_𝑖=0𝐷^𝑖_𝑡𝑢₂+^𝑡₀‖𝑢‖₂𝑑𝑠. (25)

Proof. The proof will proceed by duality argument. Let 𝜓∈𝐻²(Ω)∩𝐻¹₀(Ω) be the solution of 𝐴^∗(𝑡)𝜓=𝜌,inΩ𝜓=0,in𝜕Ω. (26) The solution 𝜓∈𝐻²(Ω)∩𝐻¹₀(Ω) satisfies the following regularity estimate: 𝜓₂≤𝐶𝜌. (27) Multiplying this equation by 𝜌 and then taking 𝐿² innerproduct over Ω, we obtain the following: 𝜌²=𝐴𝑡;𝜌,𝜓=𝐴𝑡;𝜌,𝜓−𝑅_ℎ𝜓+𝐴𝑡;𝜌,𝑅_ℎ𝜓−𝐼^∗_ℎ𝑅_ℎ𝜓−^𝑡₀𝐵𝑡,𝑠;𝜌(𝑠),𝐼^∗_ℎ𝑅_ℎ𝜓−𝑅_ℎ𝜓−𝑑𝑠_𝑡0𝐵𝑡,𝑠;𝜌(𝑠),𝑅_ℎ−𝜓−𝜓𝑑𝑠_𝑡0𝐵𝑡,𝑠;𝜌(𝑠),𝜓𝑑𝑠=𝐼₁+𝐼₂+𝐼₃+𝐼₄+𝐼₅. (28) We have 𝐼₁+𝐼₄≤𝐶ℎ²‖𝑢‖₂+^𝑡₀‖𝑢‖₂𝑑𝑠𝜓₂. (29) Applying Lemma 4, we obtain 𝐼₂+𝐼₃≤𝐶ℎ²‖𝑢‖₂+^𝑡₀‖𝑢‖₂𝑑𝑠𝜓₂. (30) Finally, we have 𝐼₅≤^𝑡₀𝜌(𝑠),𝐵^∗(𝑡,𝑠)𝜓𝑑𝑠≤𝐶^𝑡₀𝜌𝑑𝑠𝜓₂; (31) then, we have 𝜌≤𝐶ℎ²‖𝑢‖₂+^𝑡₀‖𝑢‖₂𝑑𝑠+𝐶^𝑡₀𝜌𝑑𝑠. (32)
Finally, an application of Gronwall’s lemma yields the first estimate.
The second inequality follows in a similar fashion.

Lemma 7. There exists a constant 𝐶 independent of ℎ such that 𝜌_0,𝑝+ℎ𝜌_1,𝑝≤𝐶ℎ²‖𝑢‖_2,𝑝+^𝑡₀‖𝑢‖_2,𝑝𝑑𝑠. (33)

Proof. Let 𝜌_𝑥 be an arbitrary component of ∇𝜌, with 𝑝 and 𝑞 conjugate indices; we have ‖𝜌_𝑥‖_𝑝=sup{(𝜌_𝑥,𝜑);𝜑∈𝒞^∞₀(Ω),‖𝜑‖_𝑞=1}.
For any such 𝜑, let 𝜓 be the solution of 𝐴^∗𝑡;𝜓,𝑣=−𝜑_𝑥,𝑣∀𝑣∈𝐻¹₀(Ω),𝜓=0,on𝜕Ω. (34) It follows from the regularity theory for the elliptic problem that 𝜓_1,𝑞≤𝐶_𝑝𝜑_𝑞=𝐶_𝑝. (35) We then have by application of (21) that 𝜌_𝑥,𝜑=𝐴𝑡;𝜌,𝜓=𝐴𝑡;𝜌,𝜓−𝑅_ℎ𝜓+𝐴𝑡;𝜌,𝑅_ℎ𝜓−𝐼^∗_ℎ𝑅_ℎ𝜓+^𝑡₀𝐵𝑡,𝑠;𝜌(𝑠),𝐼^∗_ℎ𝑅_ℎ𝜓𝑑𝑠=𝐼₁+𝐼₂+𝐼₃,𝐴𝑡;𝜌,𝜓−𝑅_ℎ𝜓=𝐴𝑡;𝑅_ℎ𝑅𝑢−𝑢,𝜓=−_ℎ𝑢−𝑢_𝑥,𝜑≤𝐶ℎ‖𝑢‖_2,𝑝. (36) Applying Lemma 4, we have 𝐼₂=𝐴𝑡;𝑢,𝑅_ℎ𝜓−𝐼^∗_ℎ𝑅_ℎ𝜓−𝐴𝑡;𝑉_ℎ𝑢,𝑅_ℎ𝜓−𝐼^∗_ℎ𝑅_ℎ𝜓≤𝐶ℎ‖𝑢‖_2,𝑝. (37) Finally, 𝐼₃ is estimated as follows: 𝐼₃=^𝑡₀𝐵𝑡,𝑠;𝜌(𝑠),𝐼^∗_ℎ𝑅_ℎ𝜓𝑑𝑠≤𝐶_𝑝^𝑡₀𝜌_1,𝑝𝑑𝑠. (38) Combining these estimates, we get 𝜌_1,𝑝≤𝐶ℎ‖𝑢‖_2,𝑝+𝐶_𝑝^𝑡₀𝜌_1,𝑝𝑑𝑠, (39) hence by Gronwall’s lemma 𝜌_1,𝑝≤𝐶ℎ‖𝑢‖_2,𝑝+^𝑡₀‖𝑢‖_2,𝑝𝑑𝑠. (40) The derivation of the error estimate in 𝐿^𝑝 is similar to the case when 𝑝=2.

4. Error Estimates for Semidiscrete Approximations

We split the error 𝑒(𝑡)=𝑢(𝑡)−𝑢_ℎ(𝑡) as follows: 𝑒(𝑡)=𝑢(𝑡)−𝑉_ℎ𝑢(𝑡)+𝑉_ℎ𝑢(𝑡)−𝑢_ℎ(𝑡)=𝜌+𝜃. (41)

It is easy to see that 𝜃=𝑉_ℎ𝑢(𝑡)−𝑢_ℎ(𝑡)∈𝑆_ℎ satisfies an error equation of the form 𝜃_𝑡,𝐼^∗_ℎ𝑣_ℎ+𝐴𝑡;𝜃,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝑣_ℎ𝜌𝑑𝑠=−_𝑡,𝐼^∗_ℎ𝑣_ℎ,𝑣_ℎ∈𝑆_ℎ. (42)

Since the estimates of 𝜌 are already known, it is enough to have estimates for 𝜃.

We will prove a sequence of lemmas which lead to the following result.

Lemma 8. There is a positive constant 𝐶 independent of ℎ such that |||𝜃(𝑡)|||≤𝐶|||𝜃(0)|||²+^𝑡₀𝜌_𝑡𝑑𝑠. (43)

Proof. Since 𝜃∈𝑆_ℎ we may take 𝑣_ℎ=𝜃 in (42) to obtain 12𝑑𝜃(𝑡𝑑𝑡|||)|||²+𝑐‖𝜃‖²₁≤𝜌_𝑡‖𝜃‖+𝐶^𝑡₀‖𝜃‖₁𝑑𝑠‖𝜃‖₁≤𝜌_𝑡𝜃‖+1‖2𝑐‖𝜃‖²₁+𝐶^𝑡₀‖𝜃‖²₁𝑑𝑠, (44) and hence by integration and Lemma 1, we have ||𝜃(𝑡)||²+^𝑡₀‖𝜃‖²₁𝜃(0𝑑𝑠≤𝐶|||)|||²+^𝑡₀𝜌_𝑡‖𝜃‖𝑑𝑠+^𝑡₀^𝑠₀‖𝜃(𝜏)‖²₁𝑑𝜏𝑑𝑠. (45) Gronwall’s lemma now implies the following: |||𝜃(𝑡)|||²+^𝑡₀‖𝜃‖²₁𝑑𝑠≤𝐶|||𝜃(0)|||²+^𝑡₀𝜌_𝑡𝜃‖𝜃(0‖𝑑𝑠≤𝐶|||)|||²+12sup_𝑠≤𝑡‖𝜃(𝑠)‖²+^𝑡₀𝜌_𝑡𝑑𝑠². (46) Since this holds for all ∈𝐽, we may conclude that ||𝜃(𝑡)||≤𝐶|||𝜃(0)|||+^𝑡₀𝜌_𝑡𝑑𝑠. (47)

Remark 9. If the initial value was chosen so that ‖𝑢_0ℎ−𝑢₀‖≤𝐶ℎ²‖𝑢₀‖₂,  then ‖𝜃(0)‖≤‖𝑢_0ℎ−𝑢₀‖+‖𝑉_ℎ𝑢₀−𝑢₀‖≤𝐶ℎ²‖𝑢₀‖₂. One can derive |||𝜃(𝑡)|||≤𝐶ℎ²𝑢₀₂+^𝑡₀𝑢_𝑡2𝑑𝑠. (48)

Lemma 10. There is a positive constant 𝐶 independent of ℎ such that ^𝑡₀𝜃_𝑡²𝑑𝑠+‖𝜃‖²₁≤𝐶‖𝜃(0)‖²₁+^𝑡₀𝜌_𝑡²𝑑𝑠. (49)

Proof. Set 𝑣_ℎ=𝜃_𝑡 in (42) to get 𝜃_𝑡²+12𝑑𝐴𝑑𝑡𝑡;𝜃,𝐼^∗_ℎ𝜃=−𝜌_𝑡,𝐼^∗_ℎ𝜃_𝑡−^𝑡₀𝐵𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝜃_𝑡(𝑡)+1𝑑𝑠2𝐴_𝑡𝑡;𝜃,𝐼^∗_ℎ𝜃+12𝐴𝑡;𝜃_𝑡,𝐼^∗_ℎ𝜃−𝐴𝑡;𝜃,𝐼^∗_ℎ𝜃_𝑡≤12𝜌_𝑡²+12𝜃_𝑡²+𝐴_𝑡𝑡;𝜃,𝐼^∗_ℎ𝜃+12𝐴𝑡;𝜃_𝑡,𝐼^∗_ℎ𝜃−𝐴𝑡;𝜃,𝐼^∗_ℎ𝜃_𝑡−𝑑𝑑𝑡^𝑡₀𝐵𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝜃(𝑡)(𝑡)𝑑𝑠+𝐵𝑡,𝑡;𝜃,𝐼^∗_ℎ𝜃(𝑡)+^𝑡₀𝐵_𝑡𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝜃(𝑡)𝑑𝑠. (50)
Then 𝜃_𝑡²+𝑑𝐴𝑑𝑡𝑡;𝜃,𝐼^∗_ℎ𝜃≤𝜌_𝑡²−𝑑𝑑𝑡^𝑡₀𝐵𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝜃𝜃‖𝑑𝑠+𝐶‖²₁+^𝑡₀‖𝜃(𝑠)‖²₁+1𝑑𝑠2𝐴𝑡;𝜃_𝑡,𝐼^∗_ℎ𝜃−𝐴𝑡;𝜃,𝐼^∗_ℎ𝜃_𝑡. (51)
In addition, recall that 𝐴𝑡;𝑢_ℎ,𝐼^∗_ℎ𝑣_ℎ−𝐴𝑡;𝑣_ℎ,𝐼^∗_ℎ𝑢_ℎ≤𝐶ℎ𝑢_ℎ1𝑣_ℎ1,∀𝑢_ℎ,𝑣_ℎ∈𝑆_ℎ; (52) then applying an inverse inequality and using kickback argument, we obtain 𝐴𝑡;𝜃_𝑡,𝐼^∗_ℎ𝜃−𝐴𝑡;𝜃,𝐼^∗_ℎ𝜃_𝑡≤𝐶ℎ𝜃_𝑡1‖𝜃‖₁≤𝐶𝜃_𝑡‖𝜃‖₁≤𝜀𝜃_𝑡²+𝐶‖𝜃‖²₁. (53) Combining these estimates, we derive 𝜃_𝑡²+𝑑𝐴𝑑𝑡𝑡;𝜃,𝐼^∗_ℎ𝜃≤𝜌_𝑡²−𝑑𝑑𝑡^𝑡₀𝐵𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝜃𝜃‖𝑑𝑠+𝐶‖²₁+^𝑡₀‖𝜃(𝑠)‖²₁𝑑𝑠. (54) So after integration in time and using the weak coercivity of 𝐴(𝑡;𝜃,𝐼^∗_ℎ𝜃), we get ^𝑡₀𝜃_𝑡²𝑑𝑠+𝑐₀‖𝜃‖²₁≤𝑐₀‖𝜃(0)‖²₁+^𝑡₀𝜌_𝑡²+𝑑𝑠_𝑡0𝐵𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝜃𝑑𝑠+𝐶^𝑡₀‖𝜃(𝑠)‖²₁𝑑𝑠≤𝑐₀‖𝜃(0)‖²₁+𝑐2‖𝜃‖²₁+𝐶^𝑡₀𝜌_𝑡²+‖𝜃(𝑠)‖²₁𝑑𝑠, (55) and by Gronwall’s lemma, ^𝑡₀𝜃_𝑡²𝑑𝑠+𝑐‖𝜃‖²₁≤𝐶‖𝜃(0)‖²₁+^𝑡₀𝜌_𝑡²𝑑𝑠. (56)

Remark 11. If 𝜃(0)=0, then ^𝑡₀𝜃_𝑡²𝑑𝑠+𝑐‖𝜃‖²₁≤𝐶ℎ^2𝑡₀𝑢_𝑡²₂𝑑𝑠. (57)

Theorem 12 (error estimates in 𝐿² and 𝐻¹-norms). Let 𝑢,  𝑢_ℎ be the solutions of (2) and (11), respectively. Assume that 𝑢,𝑢_𝑡∈𝐿^∞(𝐻¹₀∩𝐻²).(a)Let 𝑢_0ℎ be chosen so that ‖𝑢_0ℎ−𝑢₀‖≤𝐶ℎ²‖𝑢₀‖₂. Then for 𝑇>0 fixed there is a constant 𝐶=𝐶(𝑇) independent of ℎ, such that for all 0<𝑡<𝑇, 𝑢_ℎ(𝑡)−𝑢(𝑡)≤𝐶ℎ²𝑢₀₂+^𝑡₀𝑢_𝑡2𝑑𝑠. (58)(b)Let 𝑢_0ℎ be chosen so that ‖𝑢_0ℎ−𝑢₀‖₁≤𝐶ℎ‖𝑢₀‖₂. Then for 𝑇>0 fixed there is a constant 𝐶=𝐶(𝑇) independent of ℎ, such that for all 0<𝑡<𝑇, 𝑢_ℎ(𝑡)−𝑢(𝑡))₁≤𝐶ℎ𝑢₀₂+^𝑡₀𝑢_𝑡2𝑑𝑠. (59)

We now prove error estimates for FVE approximations in 𝐿^𝑝 and 𝑊^1,𝑝-norms.

Theorem 13 (error estimates in𝐿^𝑝and𝑊^1,𝑝-norms). Let 𝑢,  𝑢_ℎ be the solutions of (2) and (11), respectively and 𝑢_0ℎ=𝑉_ℎ𝑢₀. Assume that 𝑢,𝑢_𝑡∈𝐿^∞(𝐻¹₀∩𝑊^2,𝑝). For ℎ sufficiently small, we have 𝑢−𝑢_ℎ0,𝑝≤𝐶ℎ²𝑢₀₂+^𝑡₀𝑢_𝑡2𝑑𝑠,𝑢−𝑢_ℎ1,𝑝≤𝐶ℎ𝑢₀₂+‖𝑢‖_2,𝑝+^𝑡₀𝑢_𝑡2𝑑𝑠. (60)

Proof. If 2≤𝑝<∞, by the following Sobolev embedding inequality ‖𝜃‖_0,𝑝≤𝐶‖𝜃‖₁, (61) then the first desired estimate follows from Lemmas 7 and 10.
Given 𝜑∈𝒞^∞₀(Ω), find 𝜓∈𝐻¹₀(Ω) such that 𝐴(𝑡)^∗𝜓=−𝜑_𝑥𝜓,inΩ,𝜓=0,on𝜕Ω,_1,𝑞≤𝜑_0,𝑞. (62) We have 𝑢−𝑢_ℎ𝑥,𝜑=𝐴𝑡;𝑢−𝑢_ℎ,𝜓=𝐴𝑡;𝑢−𝑢_ℎ,𝜓−𝑅_ℎ𝜓+𝐴𝑡;𝑢−𝑢_ℎ,𝑅_ℎ𝜓−𝐼^∗_ℎ𝑅_ℎ𝜓−^𝑡₀𝐵𝑡,𝑠;𝑢−𝑢_ℎ(𝑠),𝐼^∗_ℎ𝑅_ℎ𝜓−𝑑𝑠𝑢−𝑢_ℎ𝑡,𝐼^∗_ℎ𝑅_ℎ𝜓=𝐼₁+𝐼₂+𝐼₃+𝐼₄,𝐼₁≤𝐴𝑡;𝑢−𝑅_ℎ𝑢,𝜓≤𝐶𝑢−𝑅_ℎ𝑢_1,𝑝𝜓_1,𝑞≤𝐶ℎ‖𝑢‖_2,𝑝𝜓_1,𝑞. (63) By Lemma 4𝐼₂≤𝐴𝑡;𝑢−𝑢_ℎ,𝑅_ℎ𝜓−𝐼^∗_ℎ𝑅_ℎ𝜓≤𝐶ℎ𝑢−𝑢_ℎ1,𝑝+|𝑢|_2,𝑝𝜓_1,𝑞,𝐼₃≤^𝑡₀𝑢−𝑢_ℎ1,𝑝𝑑𝑠𝜓_1,𝑞,𝐼₄≤𝑢−𝑢_ℎ𝜓≤𝐶ℎ²𝑢₀₂+^𝑡₀𝑢_𝑡2𝑑𝑠𝜓_1,𝑞, (64) where we have used the fact ‖𝜓‖≤‖𝜓‖_1,𝑟,  𝑟>1. Combining these estimates, we get 𝑢−𝑢_ℎ𝑥,𝜑≤𝐶ℎ𝑢₀₂+‖𝑢‖_2,𝑝+^𝑡₀𝑢_𝑡2𝑑𝑠𝜓_1,𝑞,𝑢−𝑢_{ℎ𝑥0,𝑝}=sup𝑢−𝑢_ℎ𝑥,𝜑𝜑_0,𝑞≤𝐶ℎ𝑢−𝑢_ℎ1,𝑝+𝐶ℎ𝑢₀₂+‖𝑢‖_2,𝑝+^𝑡₀𝑢_𝑡2𝑑𝑠. (65) Hence using the Poincaré inequality, we have forℎsufficiently small 𝑢−𝑢_ℎ1,𝑝≤𝐶ℎ𝑢₀₂+‖𝑢‖_2,𝑝+^𝑡₀𝑢_𝑡2𝑑𝑠. (66)

We compare the relationship between covolume solution and the Galerkin finite element solution.

Corollary 14. Let𝑢_ℎbe the finite element solution to (2); that is, 𝑢_ℎ𝑡,𝑣_ℎ+𝐴𝑡;𝑢_ℎ,𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢_ℎ(𝑠),𝑣_ℎ𝑑𝑠=𝑓,𝑣_ℎ,𝑣_ℎ∈𝑆_ℎ,𝑢_ℎ(0)=𝑅_ℎ𝑢₀. (67) For ℎ sufficiently small, we have 𝑢_ℎ−𝑢_ℎ1,𝑝ℎ≤𝐶𝑢−𝑢_ℎ1,𝑝+𝑢−𝑢_ℎ𝑡+𝑢_ℎ−𝑢_𝑡+∫^𝑡₀𝑢−𝑢_ℎ(𝑠)_1,𝑝+𝑢−𝑢_ℎ(𝑠)_1,𝑝(𝑢)𝑑𝑠≤𝐶ℎ. (68)

Proof. By (2) and (67), 𝑢_ℎ−𝑢_𝑡,𝑣_ℎ+𝐴𝑡;𝑢_ℎ−𝑢,𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢_ℎ−𝑢(𝑠),𝑣_ℎ𝑑𝑠=0,𝑣_ℎ∈𝑆_ℎ. (69) Consider the following auxiliary problem. For any such 𝜑, let 𝜓 be the solution of the following: 𝐴(𝑡)^∗𝜓=−𝜑_𝑥,inΩ,𝜓=0,on𝜕Ω, (70) with 𝜓_1,𝑞≤𝜑_0,𝑞,𝑢_ℎ−𝑢_ℎ𝑥𝑢,𝜑=𝐴𝑡;_ℎ−𝑢_ℎ𝑢,𝜓=𝐴𝑡;_ℎ−𝑢_ℎ,𝜓−𝑅_ℎ𝜓+𝐴𝑡;𝑢−𝑢_ℎ,𝑅_ℎ𝜓−𝐴𝑡;𝑢−𝑢_ℎ,𝐼^∗_ℎ𝑅_ℎ𝜓−𝑢−𝑢_ℎ𝑡,𝐼^∗_ℎ𝑅_ℎ𝜓−^𝑡₀𝐵𝑡,𝑠;𝑢−𝑢_ℎ(𝑠),𝐼^∗_ℎ𝑅_ℎ𝜓𝑢𝑑𝑠+𝐴𝑡;_ℎ−𝑢,𝑅_ℎ𝜓=𝐴𝑡;𝑢−𝑢_ℎ,𝑅_ℎ𝜓−𝐴𝑡;𝑢−𝑢_ℎ,𝐼^∗_ℎ𝑅_ℎ𝜓−𝑢−𝑢_ℎ𝑡,𝐼^∗_ℎ𝑅_ℎ𝜓−𝑢_ℎ−𝑢_𝑡,𝑅_ℎ𝜓−^𝑡₀𝐵𝑡,𝑠;𝑢−𝑢_ℎ(𝑠),𝐼^∗_ℎ𝑅_ℎ𝜓−𝑑𝑠_𝑡0𝐵𝑡,𝑠;𝑢_ℎ−𝑢(𝑠),𝑅_ℎ𝜓𝑑𝑠=𝐼₁+𝐼₂+𝐼₃. (71) On the other hand, 𝐼₁≤𝐶ℎ𝑢−𝑢_ℎ1,𝑝𝜓_1,𝑞,𝐼₂≤𝐶𝑢−𝑢_ℎ𝑡+𝑢_ℎ−𝑢_𝑡𝜓≤𝐶𝑢−𝑢_ℎ𝑡+𝑢_ℎ−𝑢_𝑡𝜓_1,𝑞, (72) where we have used the fact ‖𝜓‖≤‖𝜓‖_1,𝑟,  𝑟>1𝐼₃≤^𝑡₀𝑢−𝑢_ℎ(𝑠)_1,𝑝+𝑢−𝑢_ℎ(𝑠)_1,𝑝𝑑𝑠𝜓_1,𝑞𝑢_ℎ−𝑢_{ℎ𝑥0,𝑝}=sup_{𝜑∈𝒞∞0}𝑢_ℎ−𝑢_ℎ𝑥,𝜑𝜑_0,𝑞ℎ≤𝐶𝑢−𝑢_ℎ1,𝑝+𝑢−𝑢_ℎ𝑡+𝑢_ℎ−𝑢_𝑡+∫^𝑡₀𝑢−𝑢_ℎ(𝑠)_1,𝑝+𝑢−𝑢_ℎ(𝑠)_1,𝑝.𝑑𝑠 (73) We deduce the result from the known finite element estimates.

Remark 15. In order to estimate ‖(𝑢−𝑢_ℎ)_𝑡‖, by differentiating (42) with respect to 𝑡, we obtain 𝜃_𝑡𝑡,𝐼^∗_ℎ𝑣_ℎ+𝐴𝑡;𝜃_𝑡,𝐼^∗_ℎ𝑣_ℎ+𝐴_𝑡𝑡;𝜃_𝑡,𝐼^∗_ℎ𝑣_ℎ+𝐵𝑡,𝑡;𝜃,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵_𝑡𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝑣_ℎ𝜌𝑑𝑠=−_𝑡𝑡,𝐼^∗_ℎ𝑣_ℎ. (74) Setting 𝑣_ℎ=𝜃_𝑡, we obtain 12𝑑𝜃𝑑𝑡^𝑡2+𝑐𝜃_𝑡²₁≤𝜌_𝑡𝑡𝜃_𝑡+12𝑐𝜃_𝑡²₁+𝐶‖𝜃‖²₁+^𝑡₀‖𝜃‖²₁≤𝜌𝑑𝑠_𝑡𝑡𝜃_𝑡+12𝑐𝜃_𝑡²₁+𝐶^𝑡₀𝜃_𝑡²₁𝑑𝑠. (75) Using kickback argument, integrating and applying Gronwall’s lemma, we deduce 𝜃_𝑡≤𝐶𝜃_𝑡(0)+^𝑡₀𝜌_𝑡𝑡1𝑑𝑠≤𝐶ℎ²𝑢₀₂+𝑢₁₂+^𝑡₀𝑢_𝑡𝑡2𝑑𝑠. (76)

5. The Lumped Mass Finite Volume Element Method

In this section, we restrict our study to the 2D case. A simple way to define the lumped mass scheme [21] is to replace the mass matrix 𝑀_ℎ in (14) by the diagonal matrix 𝑀_ℎ obtained by taking for its diagonal elements the numbers 𝑀_ℎ𝑖𝑖=∑^𝑁ℎ_𝑗=1𝑀_ℎ𝑖𝑗 or by lumping all masses in one row into the diagonal entry. This makes the inversion of the matrix in front of 𝛼^′(𝑡) a triviality. We will therefore study the matrix problem 𝑀_ℎ𝛼^′(𝑡)+𝐴_ℎ(𝑡)𝛼(𝑡)+^𝑡₀𝐵_ℎ(𝑡,𝑠)𝛼(𝑠)𝑑𝑠=𝐹_ℎ(𝑡). (77)

We know that the lumped mass method defined by (77) above is equivalent to 𝐼^∗_ℎ𝑢_ℎ𝑡,𝐼^∗_ℎ𝑣_ℎ+𝐴𝑡;𝑢_ℎ,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢_ℎ(𝑠),𝐼^∗_ℎ𝑣_ℎ=𝑑𝑠𝑓,𝐼^∗_ℎ𝑣_ℎ,𝑣_ℎ∈𝑆_ℎ. (78)

Our alternative interpretation of this procedure will be to think of (77) as being obtained by evaluating the first term in (78) by numerical quadrature. Let 𝐾 be a triangle of the triangulation 𝑇_ℎ, let 𝑥_𝑗,𝑗=1,2,3,  be its vertices, and consider the quadrature formula 𝑄_𝐾,ℎ𝑓=13area𝐾³_𝑗=1𝑓𝑥_𝑗≃_𝐾𝑓𝑑𝑥. (79)

We may then define the associated bilinear form in𝑆_ℎ×𝑆^∗_ℎ, using the quadrature scheme, by the following: 𝑣_ℎ,𝜂_ℎℎ=_{𝐾∈𝑇ℎ}𝑄_𝐾,ℎ𝑣_ℎ𝜂_ℎ=_{𝑥𝑖∈𝑁̀𝑎ℎ}𝑣_ℎ𝑥_𝑖𝜂_ℎ𝑥_𝑖𝑉_𝑥𝑖,∀𝑣_ℎ∈𝑆_ℎ,𝜂_ℎ∈𝑆^∗_ℎ. (80)

We note that ‖𝑣_ℎ‖²_ℎ=(𝑣_ℎ,𝐼^∗_ℎ𝑣_ℎ)_ℎ is a norm in 𝑆_ℎ which is equivalent to the 𝐿²-norm uniformly in ℎ; that is, there exist two positive constants 𝐶₁ and 𝐶₂ such that for all 𝑣_ℎ∈𝑆_ℎ, we have 𝐶₀𝑣_ℎ≤𝑣_ℎℎ≤𝐶₁𝑣_ℎ,∀𝑣_ℎ∈𝑆_ℎ. (81)

We note that the aforementioned definition (𝑣_ℎ,𝜂_ℎ)_ℎ may be used also for 𝜂_ℎ∈𝑆_ℎ  and that (𝑣_ℎ,𝑤_ℎ)_ℎ=(𝑣_ℎ,𝐼^∗_ℎ𝑤_ℎ)_ℎ for 𝑣_ℎ,𝑤_ℎ∈𝑆_ℎ.

The lumped mass method defined by (78) is equivalent to 𝑢_ℎ𝑡,𝐼^∗_ℎ𝑣_ℎℎ+𝐴𝑡;𝑢_ℎ,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢_ℎ(𝑠),𝐼^∗_ℎ𝑣_ℎ=𝑑𝑠𝑓,𝐼^∗_ℎ𝑣_ℎ,𝑣_ℎ∈𝑆_ℎ. (82)

We introduce the quadrature error 𝜀_ℎ𝑣_ℎ,𝑤_ℎ=𝑣_ℎ,𝑤_ℎℎ−𝑣_ℎ,𝑤_ℎ, (83)

Lemma 16 (see [21]). Let 𝑣_ℎ,𝑤_ℎ∈𝑆_ℎ. Then 𝜀_ℎ𝑣_ℎ,𝑤_ℎ≤𝐶ℎ²∇𝑣_ℎ∇𝑤_ℎ. (84)

Theorem 17. Let 𝑢_ℎ and 𝑢 be the solutions of (82) and (2), respectively, and assume 𝑢_ℎ(0)=𝑅_ℎ𝑢₀. Then we have for the error in the lumped mass semidiscrete method, for 𝑡≥0, the following: 𝑢_ℎ(𝑡)−𝑢(𝑡)≤𝐶ℎ²𝑢₀₂+‖𝑢‖₂+^𝑡₀𝑢_𝑡2𝑑𝑠. (85)

Proof. In order to estimate ‖𝜃‖, we write 𝜃_𝑡,𝐼^∗_ℎ𝑣_ℎℎ+𝐴𝑡;𝜃,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝑣_ℎ=𝑢𝑑𝑠_ℎ𝑡,𝐼^∗_ℎ𝑣_ℎℎ+𝐴𝑡;𝑢_ℎ,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝑢_ℎ(𝑠),𝐼^∗_ℎ𝑣_ℎ−𝑉𝑑𝑠_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎℎ−𝐴𝑡;𝑉_ℎ𝑢,𝐼^∗_ℎ𝑣_ℎ−^𝑡₀𝐵𝑡,𝑠;𝑉_ℎ𝑢(𝑠),𝐼^∗_ℎ𝑣_ℎ=𝑑𝑠𝑓,𝐼^∗_ℎ𝑣_ℎ−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎℎ−𝐴𝑡;𝑢,𝐼^∗_ℎ𝑣_ℎ−^𝑡₀𝐵𝑡,𝑠;𝑢(𝑠),𝐼^∗_ℎ𝑣_ℎ=𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎ−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎℎ=−𝜌_𝑡,𝐼^∗_ℎ𝑣_ℎ−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎℎ+𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎ. (86) We rewrite 𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎℎ−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎ=𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎℎ−𝑉_ℎ𝑢_𝑡,𝑣_ℎ+𝑉_ℎ𝑢_𝑡,𝑣_ℎ−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎ=𝜀_ℎ𝑉_ℎ𝑢_𝑡,𝑣_ℎ+𝑉_ℎ𝑢_𝑡,𝑣_ℎ−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎ.𝜃_𝑡,𝐼^∗_ℎ𝑣_ℎℎ+𝐴𝑡;𝜃,𝐼^∗_ℎ𝑣_ℎ+^𝑡₀𝐵𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝑣_ℎ𝜌𝑑𝑠=−_𝑡,𝐼^∗_ℎ𝑣_ℎ+𝜀_ℎ𝑉_ℎ𝑢_𝑡,𝑣_ℎ+𝑉_ℎ𝑢_𝑡,𝑣_ℎ−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝑣_ℎ. (87) Setting 𝑣_ℎ=𝜃 in (87), we obtain 12𝑑‖𝜃‖𝑑𝑡^2ℎ+𝑐₀‖𝜃‖²₁≤𝜌_𝑡𝜃‖+1‖2𝑐₀‖𝜃‖²₁+𝐶^𝑡₀‖𝜃‖²₁𝑑𝑠+𝜀_ℎ𝑉_ℎ𝑢_𝑡,𝜃+𝑉_ℎ𝑢_𝑡,𝜃−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝜃. (88) Using Lemma 16 and the inverse estimate, we get 𝜀_ℎ𝑉_ℎ𝑢_𝑡,𝜃≤𝐶ℎ²∇𝑉_ℎ𝑢_𝑡‖∇𝜃‖≤𝐶ℎ²∇𝑢_𝑡‖∇𝜃‖≤𝐶ℎ∇𝑢_𝑡‖𝜃‖, (89) we have 𝑉_ℎ𝑢_𝑡,𝜃−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝜃≤𝐶ℎ∇𝑢_𝑡‖𝜃‖. (90) Using Young’s inequality and Gronwall’s lemma to eliminate ‖𝜃‖₁ on the right-hand side and using integration in 𝑡, we get the result. 12𝑑‖𝜃‖𝑑𝑡^2ℎ+𝑐₀‖𝜃‖≤𝜌_𝑡‖𝜃‖+𝐶ℎ∇𝑢_𝑡‖𝜃‖. (91) Using Young’s inequality to eliminate ‖𝜃‖ on the right hand side it becomes.
Using integration in 𝑡, we get the result.

We will now show that the 𝐻¹-norm error bound of theorem remains valid for the lumped mass method (82).

Theorem 18. Let 𝑢_ℎ and 𝑢 be the solutions of (82) and (2), respectively, and assume 𝑢_ℎ(0)=𝑅_ℎ𝑢₀,𝑢_1ℎ(0)−𝑢₁≤𝐶ℎ²𝑢₁₂. (92) Then, we have for the error in the lumped mass semidiscrete method, for 𝑡≥0, the following: 𝑢_ℎ(𝑡)−𝑢(𝑡)₁≤𝐶ℎ²𝑢₀₂+𝑢₁₂+^𝑡₀𝑢_𝑡𝑡2𝑑𝑠. (93)

Proof. Setting 𝑣_ℎ=𝜃_𝑡 in (87), we obtain 𝜃_𝑡²_ℎ+12𝑑𝐴𝑑𝑡𝑡;𝜃,𝐼^∗_ℎ𝜃=12𝐴_𝑡𝑡;𝜃,𝐼^∗_ℎ𝜃+12𝐴𝑡;𝜃_𝑡,𝐼^∗_ℎ𝜃−𝐴𝑡;𝜃,𝐼^∗_ℎ𝜃_𝑡𝐵𝑡,𝑡;𝜃(𝑡),𝐼^∗_ℎ𝜃(𝑡)+^𝑡₀𝐵_𝑡𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝜃(𝑡)−𝑑𝑑𝑠𝑑𝑡^𝑡₀𝐵𝑡,𝑠;𝜃(𝑠),𝐼^∗_ℎ𝜃(𝑡)𝑑𝑠−𝜌_𝑡,𝐼^∗_ℎ𝜃_𝑡−𝜀_ℎ𝑉_ℎ𝑢_𝑡,𝜃_𝑡+𝑉_ℎ𝑢_𝑡,𝜃−𝑉_ℎ𝑢_𝑡,𝐼^∗_ℎ𝜃. (94) It follows, thus, that using integration in 𝑡 and Gronwall’s lemma, we have ^𝑡₀𝜃_𝑡²_ℎ+‖𝜃‖²₁≤𝐶‖∇𝜃(0)‖²+𝐶^𝑡₀𝜌_𝑡𝜃_𝑡𝑑𝑠+𝐶ℎ^2𝑡₀𝑢_𝑡²₁𝑑𝑠. (95)

6. Full Discretization

Let 𝜕𝑈^𝑛=(𝑈^𝑛−𝑈^𝑛−1)/𝑘 be the backward difference quotient of 𝑈^𝑛; assume that 𝐴_ℎ=𝑃_ℎ𝐴 is a discrete analogue of 𝐴 (similarly 𝐵_ℎ=𝑃_ℎ𝐵), where 𝑃_ℎ∶𝐿²(Ω)→𝑆^∗_ℎ the 𝐿² projection is defined by 𝑃_ℎ𝑣,𝐼^∗_ℎ𝑣_ℎ=𝑣,𝐼^∗_ℎ𝑣_ℎ,𝑣∈𝐿²(Ω),𝑣_ℎ∈𝑆_ℎ. (96)

In order to define fully discrete approximation of (11), we discretize the time by taking 𝑡_𝑛=𝑛𝑘,  𝑘>0,  𝑛=1,2,…   and use the numerical quadrature ^{𝑡𝑛−1/2}₀𝑔(𝑠)𝑑𝑠≈^𝑛_𝑘=1𝜔_𝑛,𝑘𝑔𝑡_𝑘−1/2,𝑡_𝑛−1/2=1𝑛−2𝑘. (97)

Here {𝜔_𝑛,𝑘} are the integration weights and we assume that the following error estimate is valid: 𝑞^𝑛𝑔=^{𝑡𝑛−1/2}₀𝑔(𝑠)𝑑𝑠−^𝑛_𝑘=1𝜔_𝑛,𝑘𝑔𝑡_𝑘−1/2≤𝐶𝑘^2𝑡𝑛₀𝑔^′+𝑔^′′𝑑𝑠. (98)

Now, define our complete discrete FVE approximation of (11) by the following: find 𝑈^𝑛∈𝑆_ℎ for 𝑛=1,2,…, such that for all 𝑣_ℎ∈𝑆_ℎ𝜕𝑈^𝑛,𝐼^∗_ℎ𝑣_ℎ+𝐴𝑡_𝑛−1/2;𝑈^𝑛−1/2,𝐼^∗_ℎ𝑣_ℎ+^𝑛_𝑘=1𝜔_𝑛,𝑘𝐵𝑡_𝑛−1/2,𝑡_𝑘−1/2,𝑈^𝑘−1/2,𝐼^∗_ℎ𝑣_ℎ=𝑓^𝑛−1/2,𝐼^∗_ℎ𝑣_ℎ,𝑈⁰in𝑆_ℎ, (99) where 𝑈^𝑛−1/2=(𝑈^𝑛+𝑈^𝑛−1)/2.

Theorem 19. Let 𝑢(𝑡) and 𝑈^𝑛 be the solutions of problem (2) and its complete discrete scheme (99), respectively. Then for any 𝑇>0 there exists a positive constant 𝐶=𝐶(𝑇)>0, independent of ℎ, such that for 0<𝑡_𝑛≤𝑇𝑢𝑡_𝑛−𝑈^𝑛≤𝐶ℎ²𝑢₀₂+^𝑡𝑛₀𝑢_𝑡2𝑑𝑠+𝐶𝑘^2𝑡𝑛₀‖𝑢‖₂+𝑢_𝑡2+𝑢_𝑡𝑡2+𝑢_𝑡𝑡𝑡𝑑𝑠. (100)

Proof. Let us split the error into two parts: 𝑢(𝑡_𝑛)−𝑈^𝑛=𝜌^𝑛+𝜃^𝑛,  where  𝜌^𝑛=𝑢(𝑡_𝑛)−𝑉_ℎ𝑢(𝑡_𝑛) and 𝜃^𝑛=𝑉_ℎ𝑢(𝑡_𝑛)−𝑈^𝑛, and let 𝑊=𝑉_ℎ𝑢(𝑡)∈𝑆_ℎ be the Ritz-Volterra projection of 𝑢. Then from (2) and (99) we have for all 𝑣_ℎ∈𝑆_ℎ the following: 𝜕𝜃^𝑛,𝐼^∗_ℎ𝑣_ℎ+𝐴𝑡_𝑛−1/2;𝜃^𝑛−1/2,𝐼^∗_ℎ𝑣_ℎ+^𝑛_𝑘=1𝜔_𝑛,𝑘𝐵𝑡_𝑛−1/2,𝑡_𝑘−1/2,𝜃^𝑘−1/2,𝐼^∗_ℎ𝑣_ℎ=−𝑟_𝑛,𝐼^∗_ℎ𝑣_ℎ,∀𝑣_ℎ∈𝑆_ℎ, (101) where 𝑟_𝑛=𝑟¹_𝑛+𝑟²_𝑛+𝑟³_𝑛+𝑟⁴_𝑛,𝑟¹_𝑛=𝜕𝜌^𝑛,𝑟²_𝑛=𝜕𝑢𝑡_𝑛−𝑢_𝑡𝑡_𝑛−1/2,𝑟³_𝑛=𝐴𝑡_𝑛−1/2;𝑢𝑡_𝑛+𝑢𝑡_𝑛−12−𝑢𝑡_𝑛−1/2,𝑟⁴_𝑛=𝑞^𝑛𝐵_ℎ𝑊=^𝑛_𝑘=1𝜔_𝑛,𝑘𝐵_ℎ𝑡_𝑛−1/2,𝑡_𝑘−1/2,𝑊^𝑘−1/2,𝐼^∗_ℎ𝑣_ℎ−^{𝑡𝑛−1/2}₀𝐵𝑡_𝑛,𝑠,𝑊(𝑠),𝐼^∗_ℎ𝑣_ℎ𝑑𝑠. (102)
In fact, by Taylor expansion, 𝑢^𝑛+1=𝑢^𝑛+𝑘𝑢^′𝑡_𝑛+^𝑡𝑛+1_𝑡𝑛𝑢^′′(𝑠)𝑡_𝑛+1−𝑠𝑑𝑠=𝑢^𝑛+𝑘𝑢^′𝑡_𝑛+𝑘²2𝑢^′′𝑡_𝑛+𝑘³6𝑢⁽³⁾𝑡_𝑛+16^𝑡𝑛+1_𝑡𝑛𝑢⁽⁴⁾(𝑠)𝑡_𝑛+1−𝑠³𝑑𝑠, (103) we have 𝑟¹_𝑛=𝜕𝜌^𝑛≤1𝑘^𝑡𝑛_𝑡𝑛−1𝜌_𝑡ℎ𝑑𝑠≤𝐶²𝑘^𝑡𝑛_𝑡𝑛−1𝑢_𝑡𝑡2𝑟𝑑𝑠,_2𝑛=𝜕𝑢𝑡_𝑛−𝑢_𝑡𝑡_𝑛−1/2=1𝑘^𝑡𝑛_𝑡𝑛−1𝑢_𝑡(𝑠)−𝑢_𝑡𝑡_𝑛−1/2𝑑𝑠≤𝐶𝑘^𝑡𝑛_𝑡𝑛−1𝑢⁽³⁾(𝑠𝑟)𝑑𝑠,_3𝑛=𝐴𝑡_𝑛−1/2;𝑢𝑡_𝑛+𝑢𝑡_𝑛−12−𝑢𝑡_𝑛−1/2,𝐼^∗_ℎ𝑣_ℎ≤𝐶𝑘^𝑡𝑛_𝑡𝑛−1𝐴𝑢_𝑡𝑡(𝑠)𝑑𝑠≤𝐶𝑘^𝑡𝑛_𝑡𝑛−1𝑢_𝑡𝑡2𝑑𝑠. (104) In addition, the quadrature error satisfies 𝑟⁴_𝑛=𝑞^𝑛−1/2𝐵_ℎ𝑊=^𝑛_𝑘=1𝜔_𝑛,𝑘𝐵𝑡_𝑛−1/2,𝑡_𝑘−1/2,𝑊^𝑘−1/2,𝐼^∗_ℎ𝑣_ℎ−^{𝑡𝑛−1/2}₀𝐵𝑡_𝑛,𝑠,𝑊(𝑠),𝐼^∗_ℎ𝑣_ℎ𝑑𝑠≤𝐶𝑘^2𝑡𝑛₀𝐵_ℎ𝑊_𝑠𝑠𝑑𝑠≤𝐶𝑘^2𝑡𝑛₀‖𝑢‖₂+𝑢_𝑡2+𝑢_𝑡𝑡2𝑘𝑑𝑠._𝑁𝑛=1𝑟_𝑛≤𝐶ℎ^2𝑡𝑛₀𝑢_𝑡𝑡2𝑑𝑠+𝐶𝑘^2𝑡𝑛₀‖𝑢‖₂+𝑢_𝑡2+𝑢_𝑡𝑡2+𝑢⁽³⁾𝑑𝑠. (105) Taking 𝑣_ℎ=𝜃^𝑛−1/2 in (101) and noting that (𝜕𝜃^𝑛,𝐼^∗_ℎ𝜃^𝑛−1/2)=(1/2)𝜕|||𝜃^𝑛|||², there is 𝜃^𝑛2−𝜃^𝑛−12+2𝑘𝑐𝜃^𝑛−1/22₁≤𝐶𝑘^2𝑛_𝑘=1𝜃^𝑘−1/2₁𝜃^𝑛−1/2₁+𝐶𝑘𝑟_𝑛𝜃^𝑛−1/2≤𝑘𝑐𝜃^𝑛−1/22₁+𝐶𝑘^2𝑛_𝑘=1𝜃^𝑘−1/22₁+𝐶𝑘𝑟_𝑛𝜃^𝑛−1/2. (106) Summing from 𝑛=1 to 𝑁 and then, after cancelling the common factor and using Gronwall’s lemma, we obtain 𝜃^𝑁2≤𝐶𝜃⁰²+𝐶𝑘^𝑁_𝑘=1𝑟_𝑛𝜃^𝑘+𝜃^𝑘−1/2, (107) and then 𝜃^𝑁≤𝐶𝜃⁰+𝐶𝑘^𝑁_𝑛=1𝑟_𝑛, (108) the theorem follows from the estimates of 𝜌^𝑛 and 𝑟^𝑛.