A New Positive Definite Expanded Mixed Finite Element Method for Parabolic Integrodifferential Equations

Abstract

A new positive definite expanded mixed finite element method is proposed for parabolic partial integrodifferential equations. Compared to expanded mixed scheme, the new expanded mixed element system is symmetric positive definite and both the gradient equation and the flux equation are separated from its scalar unknown equation. The existence and uniqueness for semidiscrete scheme are proved and error estimates are derived for both semidiscrete and fully discrete schemes. Finally, some numerical results are provided to confirm our theoretical analysis.

1. Introduction

In this paper, we consider the following initial-boundary value problem of parabolic partial integrodifferential equations: 𝑢_𝑡−∇⋅𝑎(𝑥,𝑡)∇𝑢+𝑏(𝑥,𝑡)^𝑡₀(),()𝑢()()∇𝑢𝑑𝑠=𝑓𝑥,𝑡𝑥,𝑡∈Ω×𝐽,𝑥,𝑡=0,𝑥,𝑡∈𝜕Ω×𝑢()𝐽,𝑥,0=𝑢₀(𝑥),𝑥∈Ω, (1.1) where Ω is a bounded convex polygonal domain in ℝ^𝑑,(𝑑=1,2,3) with a smooth boundary 𝜕Ω, 𝐽=(0,𝑇] is the time interval with 0<𝑇<∞. The coefficients 𝑎=𝑎(𝑥,𝑡), 𝑏=𝑏(𝑥,𝑡) are two functions, which satisfy the property that there exist some positive constants 𝑎_min, 𝑎_max, 𝑏_min, and 𝑏_max such that 0<𝑎_min≤𝑎(𝑥,𝑡)≤𝑎_max and 0<𝑏_min≤𝑏(𝑥,𝑡)≤𝑏_max.

Parabolic integrodifferential equations are a class of very important evolution equations which describe many physical phenomena such as heat conduction in material with memory, compression of viscoelastic media, and nuclear reactor dynamics. In recent years, a lot of researchers have studied the numerical methods for parabolic integrodifferential equations, such as finite element methods [1–5], mixed finite element methods [6–9], and finite volume element method [10] and so forth.

In 1994, Chen [11, 12] proposed a new mixed method, which is called a expanded mixed finite element method and proved some mathematical theories for second-order linear elliptic equation. Compared to standard mixed element methods, the expanded mixed method is expanded in the sense that three variables are explicitly approximated, namely, the scalar unknown, its gradient, and its flux (the tensor coefficient times the gradient). In 1997, Arbogast et al. [13] derived and exploited a connection between the expanded mixed method and a certain cell-centered finite difference method. And Chen proved some mathematical theories for second-order quasilinear elliptic equation [14] and fourth-order elliptic problems [15]. With the development of the expanded mixed finite element method, the method was applied to many evolution equations. In [16], some error estimates of the expanded mixed element for a kind of parabolic equation were given. Woodward and Dawson [17] studied the expanded mixed finite element method for nonlinear parabolic equation. Wu and Chen et al. [18–22] studied the two-grid methods for expanded mixed finite-element solution of semilinear reaction-diffusion equations. Song and Yuan [23] proposed the expanded upwind-mixed multistep method for the miscible displacement problem in three dimensions. Guo and Chen [24] developed and analysed an expanded characteristic-mixed finite element method for a convection-dominated transport problem. In 2010, Chen and Wang [25] proposed an 𝐻¹-Galerkin expanded mixed method for a nonlinear parabolic equation in porous medium flow, and Liu and Li [26] studied the 𝐻¹-Galerkin expanded mixed method for pseudo-hyperbolic equation. Liu [27], studied the 𝐻¹-Galerkin expanded mixed method for RLW-Burgers equation and proved semidiscrete and fully discrete optimal error estimates. Jiang and Li [28] studied an expanded mixed semidiscrete scheme for the problem of purely longitudinal motion of a homogeneous bar. In [29, 30], the expanded mixed covolume method was studied for the linear integrodifferential equation of parabolic type and elliptic problems, respectively. In [31], a posteriori error estimator for expanded mixed hybrid methods was studied and analysed.

In 2001, Yang [32] proposed a new mixed finite element method called the splitting positive definite mixed finite element procedure to treat the pressure equation of parabolic type in a nonlinear parabolic system describing a model for compressible flow displacement in a porous medium. Compared to standard mixed methods whose numerical solutions have been quite difficult because of losing positive definite properties, the proposed one does not lead to some saddle point problems. From then on, the method was applied to the hyperbolic equations [33] and pseudo-hyperbolic equations [34].

In this paper, our purpose is to propose and analyse a new expanded mixed method based on the positive definite system [32–34] for parabolic integrodifferential equations. Compared to expanded mixed methods, the proposed mixed element system is symmetric positive definite and avoids some saddle point problems. What is more, both the gradient equation and the flux equation are separated from its scalar unknown equation. The existence and uniqueness for semidiscrete scheme are proved and error estimates are derived for both semidiscrete and fully discrete schemes.

The layout of the paper is as follows. In Section 2, the positive definite expanded mixed weak formulation and semidiscrete mixed scheme are formulated, and the proof of the existence and uniqueness of the discrete solutions is given. Error estimates are derived for both semidiscrete and fully discrete schemes for problems, respectively, in Sections 3 and 4. In Section 5, some numerical results are provided to illustrate the efficiency of our method. Finally in Section 6, we will give some concluding remarks about the positive definite expanded mixed finite element method.

Throughout this paper, 𝐶 will denote a generic positive constant which does not depend on the spatial mesh parameters ℎ_𝑢,ℎ_𝜎 and time-discretization parameter 𝛿 and may be different at their occurrences. Usual definitions, notations, and norms of the Sobolev spaces as in [35–37] are used. We denote the natural inner product in 𝐿²(Ω) or [𝐿²(Ω)]^𝑑 by (⋅,⋅) with norm ‖⋅‖_𝐿2(Ω) or ‖⋅‖_𝐋2(Ω) and introduce the function space 𝐖=𝐻(div;Ω)={𝜔∈[𝐿²(Ω)]^𝑑;∇⋅𝜔∈𝐿²(Ω)}.

2. A New Expanded Mixed Variational Formulation

Introducing the auxiliary variables: 𝜆=∇𝑢,𝜎=𝑎(𝑥,𝑡)∇𝑢+𝑏(𝑥,𝑡)^𝑡₀∇𝑢𝑑𝑠=𝑎𝜆+𝑏^𝑡₀𝜆𝑑𝑠, (2.1) Then we obtain the equivalent system of parabolic partial integrodifferential equations for the problem (1.1): (a)𝑢_𝑡(),()()−∇⋅𝜎=𝑓𝑥,𝑡𝑥,𝑡∈Ω×𝐽,(b)𝜆−∇𝑢=0,𝑥,𝑡∈Ω×𝐽,(c)𝜎−𝑎𝜆−𝑏^𝑡₀𝜆𝑑𝑠=0,(𝑥,𝑡)∈Ω×𝐽, (2.2) with the initial values 𝜆(𝑥,0)=∇𝑢₀(𝑥),𝜎(𝑥,0)=𝑎∇𝑢₀(𝑥), and 𝑢(𝑥,0)=𝑢₀(𝑥).

Then, the following expanded mixed weak formulation of (2.2) can be given by: (a)𝑢_𝑡,𝑣−(∇⋅𝜎,𝑣)=𝑓(𝑥,𝑡),𝑣,∀𝑣∈𝐿²(Ω),(b)+()(c)−()−𝑏)(𝜆,𝐰𝑢,∇⋅𝐰=0,∀𝐰∈𝐖,)(𝜎,𝐳𝑎𝜆,𝐳^𝑡₀𝜆𝑑𝑠,𝐳=0,∀𝐳∈𝐖. (2.3)

From (2.3)(b) we derive 𝜆_𝑡,𝐰+𝑢_𝑡,∇⋅𝐰=0. (2.4) Taking 𝑣=∇⋅𝐰 in (2.3)(a) for 𝐰∈𝐖 and then substituting it into (2.4), we derive a new equivalent expanded mixed weak formulation of the system (2.3): (a)𝜆_𝑡+()𝑓()(b)−()−𝑏,𝐰∇⋅𝜎,∇⋅𝐰=−𝑥,𝑡,∇⋅𝐰,∀𝐰∈𝐖,)(𝜎,𝐳𝑎𝜆,𝐳^𝑡₀(c)𝑢𝜆𝑑𝑠,𝐳=0,∀𝐳∈𝐖,_𝑡,𝑣−(∇⋅𝜎,𝑣)=𝑓(𝑥,𝑡),𝑣,∀𝑣∈𝐿²(Ω). (2.5)

Let 𝒯_ℎ𝑢 and 𝒯_ℎ𝜎 be two families of quasi-regular partitions of the domain Ω, which may be the same one or not, such that the elements in the partitions have the diameters bounded by ℎ_𝑢 and ℎ_𝜎, respectively. Let 𝑋_ℎ𝑢⊂𝐿²(Ω) and 𝐕_ℎ𝜎⊂𝐖 be finite element spaces defined on the partitions 𝒯_ℎ𝑢 and 𝒯_ℎ𝜎.

Now the semidiscrete positive definite expanded mixed finite element method for (2.5) consists in determining (𝑢_ℎ,𝜆_ℎ,𝜎_ℎ)∈𝑋_ℎ𝑢×𝐕_ℎ𝜎×𝐕_ℎ𝜎 such that (a)𝜆_ℎ𝑡,𝐰_ℎ+∇⋅𝜎_ℎ,∇⋅𝐰_ℎ=−𝑓(𝑥,𝑡),∇⋅𝐰_ℎ,∀𝐰_ℎ∈𝐕_ℎ𝜎,(b)𝜎_ℎ,𝐳_ℎ−𝑎𝜆_ℎ,𝐳_ℎ−𝑏^𝑡₀𝜆_ℎ𝑑𝑠,𝐳_ℎ=0,∀𝐳_ℎ∈𝐕_ℎ𝜎,(c)𝑢_ℎ𝑡,𝑣_ℎ−∇⋅𝜎_ℎ,𝑣_ℎ=𝑓(𝑥,𝑡),𝑣_ℎ,∀𝑣_ℎ∈𝑋_ℎ𝑢, (2.6) with given an initial approximation (𝑢⁰_ℎ,𝜆⁰_ℎ,𝜎⁰_ℎ)∈𝑋_ℎ𝑢×𝐕_ℎ𝜎×𝐕_ℎ𝜎.

Remark 2.1. Compared to expanded mixed weak formulation (2.3), the new expanded mixed element system (2.6) is symmetric positive definite, that is to say the gradient function and the flux function system (2.6)(a,b) is symmetric positive definite. And both the gradient equation and the flux equation are separated from its scalar unknown equation (2.6)(c).

Theorem 2.2. There exists a unique discrete solution to the system (2.6).

Proof. Let {𝜓_𝑖(𝑥)}^𝑁1_𝑖=1 and {𝜑_𝑗(𝑥)}^𝑁2_𝑗=1 be bases of 𝑋_ℎ𝑢 and 𝐕_ℎ𝜎, respectively. Then, we have the following expressions: 𝑢_ℎ=^𝑁1_𝑖=1𝑢_𝑖(𝑡)𝜓_𝑖(𝑥),𝜆_ℎ=^𝑁2_𝑗=1𝜆_𝑗(𝑡)𝜑_𝑗(𝑥),𝜎_ℎ=^𝑁2_𝑗=1𝜎_𝑗(𝑡)𝜑_𝑗(𝑥). (2.7) Substituting these expressions into (2.6) and choosing 𝑣_ℎ=𝜓_𝑚,𝐰_ℎ=𝐳_ℎ=𝜑_𝑙, then the problems (2.6) can be written in vector matrix form as: find {𝑈(𝑡),Λ(𝑡),Σ(𝑡)} such that, for all 𝑡∈(0,𝑇](a)𝐴Λ^′(𝑡)(𝑡)(𝑡),(b)(𝑡)(𝑡)+𝐵Σ=𝐹𝐴Σ−𝐶Λ−𝐻^𝑡₀Λ(𝑠)(c)𝑑𝑠=𝟎,𝐷𝑈^′(𝑡)−𝐸Σ(𝑡)=𝐺(𝑡), (2.8) where 𝐴=𝜑_𝑗,𝜑_{𝑙𝑁2×𝑁2},𝐵=∇⋅𝜑_𝑗,∇⋅𝜑_{𝑙𝑁2×𝑁2},𝐶=𝑎𝜑_𝑗,𝜑_{𝑙𝑁2×𝑁2},𝐷=𝜓_𝑖,𝜓_{𝑚𝑁1×𝑁1},𝐻=𝑏𝜑_𝑗,𝜑_{𝑙𝑁2×𝑁2},𝐸=∇⋅𝜑_𝑗,𝜓_{𝑚𝑁1×𝑁2},𝑈(𝑡)=𝑢₁(𝑡),𝑢₂(𝑡),…,𝑢_𝑁1(𝑡)^𝑇,Λ(𝑡)=𝜆₁(𝑡),𝜆₂(𝑡),…,𝜆_𝑁2(𝑡)^𝑇,Σ(𝑡)=𝜎₁(𝑡),𝜎₂(𝑡),…,𝜎_𝑁2(𝑡)^𝑇,𝐹(𝑡)=−𝑓(𝑡),∇⋅𝜑_𝑙^𝑇_1×𝑁2,𝐺(𝑡)=𝑓(𝑡),𝜓_𝑚^𝑇_1×𝑁1. (2.9) It is easy to see that both 𝐴 and 𝐷 are symmetric positive definite. From (2.8), the problems can be written as follows: (a)Λ^′(𝑡)=−𝐴⁻¹𝐵𝐴⁻¹𝐶Λ(𝑡)+𝐻^𝑡₀Λ(𝑠)𝑑𝑠+𝐴⁻¹𝐹(𝑡),(b)Σ(𝑡)=𝐴⁻¹𝐶Λ(𝑡)+𝐴⁻¹𝐻^𝑡₀Λ(𝑠)(c)𝑈𝑑𝑠,_′(𝑡)=𝐷⁻¹𝐸Σ(𝑡)+𝐷⁻¹𝐺(𝑡). (2.10) Thus, by the theory of differential equations [38, 39], (2.10) has a unique solution, and equivalently (2.6) has a unique solution.

Remark 2.3. It is easy to see that the coefficient matrixes 𝐴,𝐵,𝐶, and 𝐻 of system (2.8) are symmetric positive definite. In view of this, the new expanded mixed element system (2.6) is symmetric positive definite.

3. Semidiscrete Error Estimates

Let 𝑋_ℎ𝑢 and 𝐕_ℎ𝜎 be finite dimensional subspaces of 𝐿²(Ω) and 𝐖, respectively, with the inverse property (see [36]) and the following approximation properties (see [40–44]): for 0≤𝑝≤+∞ and 𝑟,𝑟^∗,𝑘 positive integers inf_{𝐰ℎ∈𝐕ℎ𝜎}𝐰−𝐰_{ℎ𝐋𝑝(Ω)}≤𝐶ℎ_𝜎^𝑟+1‖𝐰‖_{𝐖𝑟+1,𝑝(Ω)},∀𝐰∈𝐻(div;Ω)∩𝑊^𝑟+1,𝑝(Ω)^𝑑,inf_{𝐰ℎ∈𝐕ℎ𝜎}∇⋅𝐰−𝐰_{ℎ𝐿𝑝(Ω)}≤𝐶ℎ^𝑟∗_𝜎‖∇⋅𝐰‖_𝑊𝑟^∗_,𝑝(Ω),∀𝐰∈𝐻(div;Ω)∩𝑊^𝑟+1,𝑝(Ω)^𝑑,inf_{𝑣ℎ∈𝑋ℎ𝑢}𝑣−𝑣_{ℎ𝐿𝑝(Ω)}≤𝐶ℎ_𝑢^𝑘+1‖𝑣‖_{𝑊𝑘+1,𝑝(Ω)},∀𝑣∈𝐿²(Ω)∩𝑊^𝑘+1,𝑝(Ω), (3.1) where 𝑟^∗=𝑟+1 for the Brezzi-Douglas-Fortin-Marini spaces [43] and the Raviart-Thomas spaces [42] and 𝑟^∗=𝑟 for the Brezzi-Douglas-Marini spaces [40, 43].

For our subsequent error analysis, we introduce two operators. It is well known that, in any one of the classical mixed finite element spaces, there exists an operator 𝑅_ℎ from 𝐻(div;Ω) onto 𝐕_ℎ𝜎, see [40–44], such that, for 1≤𝑝≤+∞, ∇⋅𝜎−𝑅_ℎ𝜎,𝜙_ℎ=0,∀𝜙_ℎ∈∇⋅𝐕_ℎ𝜎=𝜙_ℎ=∇⋅𝐰_ℎ,𝐰_ℎ∈𝐕_ℎ𝜎;𝜎−𝑅_ℎ𝜎_𝐋𝑝(Ω)≤𝐶ℎ_𝜎^𝑟+1‖𝜎‖_{𝐖𝑟+1,𝑝(Ω)};∇⋅𝜎−𝑅_ℎ𝜎_𝐿𝑝(Ω)≤𝐶ℎ^𝑟∗_𝜎‖∇⋅𝜎‖_𝑊𝑟^∗_,𝑝(Ω). (3.2) We also define the 𝐿²-project operator 𝑃_ℎ from 𝐿²(Ω) onto 𝑋_ℎ𝑢 such that 𝑣−𝑃_ℎ𝑣,𝑣_ℎ=0,∀𝑣∈𝐿²(Ω),𝑣_ℎ∈𝑋_ℎ𝑢;𝑣−𝑃_ℎ𝑣_𝐿𝑝(Ω)≤𝐶ℎ_𝑢^𝑘+1‖𝑣‖_{𝐻𝑘+1,𝑝(Ω)},∀𝑣∈𝐻^𝑘+1(Ω). (3.3) Using the definitions of the operators 𝑅_ℎ and 𝑃_ℎ, we can easily obtain the following lemma.

Lemma 3.1. Assume that the solution of system (2.5) has the regular properties that 𝑢_𝑡∈𝐿²(𝐻^𝑘+1(Ω)),    𝜆_𝑡,𝜆_𝑡𝑡,𝜎_𝑡∈𝐿²(𝐇^𝑟+1(Ω)), then one has the following estimates: 𝜆−𝑅_ℎ𝜆_{𝑡𝐋𝑝(Ω)}≤𝐶ℎ_𝜎^𝑟+1𝜆_{𝑡𝐖𝑟+1,𝑝(Ω)},𝜆−𝑅_ℎ𝜆_{𝑡𝑡𝐋𝑝(Ω)}≤𝐶ℎ_𝜎^𝑟+1𝜆_{𝑡𝑡𝐖𝑟+1,𝑝(Ω)},𝜎−𝑅_ℎ𝜎_{𝑡𝐋𝑝(Ω)}≤𝐶ℎ_𝜎^𝑟+1𝜎_{𝑡𝐖𝑟+1,𝑝(Ω)},𝑢−𝑃_ℎ𝑢_{𝑡𝐿𝑝(Ω)}≤𝐶ℎ_𝑢^𝑘+1𝑢_{𝑡𝐻𝑘+1,𝑝(Ω)}. (3.4)
Let 𝑢−𝑢_ℎ=𝑢−𝑃_ℎ𝑢+𝑃_ℎ𝑢−𝑢_ℎ=𝜂+𝜍,𝜆−𝜆_ℎ=𝜆−𝑅_ℎ𝜆+𝑅_ℎ𝜆−𝜆_ℎ=𝜌+𝜉,𝜎−𝜎_ℎ=𝜎−𝑅_ℎ𝜎+𝑅_ℎ𝜎−𝜎_ℎ=𝛾+𝜃. (3.5) Subtracting (2.6) from (2.5) and using projections (3.2) and (3.3), one obtains 𝜉_𝑡,𝐰_ℎ+∇⋅𝜃,∇⋅𝐰_ℎ=−𝜌_𝑡,𝐰_ℎ,∀𝐰_ℎ∈𝐕_ℎ𝜎,𝜃,𝐳_ℎ−𝑎𝜉,𝐳_ℎ−𝑏^𝑡₀𝜉𝑑𝑠,𝐳_ℎ=−𝛾,𝐳_ℎ+𝑎𝜌,𝐳_ℎ+𝑏^𝑡₀𝜌𝑑𝑠,𝐳_ℎ,∀𝐳_ℎ∈𝐕_ℎ𝜎,𝜍_𝑡,𝑣_ℎ−∇⋅𝜎−𝜎_ℎ,𝑣_ℎ=0,∀𝑣_ℎ∈𝑋_ℎ𝑢. (3.6)

Theorem 3.2. Assume that the approximate properties (3.1) hold, and the solution of system (2.5) has regular properties that 𝑢_𝑡,𝑢_𝑡𝑡∈𝐿²(𝐻^𝑘+1(Ω)),𝜆_𝑡,𝜆_𝑡𝑡,𝜎_𝑡,𝜎_𝑡𝑡∈𝐿²(𝐇^𝑟+1(Ω)). Then one has the error estimates 𝜆−𝜆_{ℎ𝐿∞(𝐋2(Ω))}+𝜎−𝜎_{ℎ𝐿∞(𝐋2(Ω))}≤𝐶ℎ_𝜎^𝑟+1,∇⋅𝜎−𝜎_{ℎ𝐿∞(𝐿2(Ω))}≤𝐶ℎ^𝑟∗_𝜎,𝑢−𝑢_{ℎ𝐿∞(𝐿2(Ω))}≤𝐶ℎ^𝑟∗_𝜎+ℎ_𝑢^𝑘+1. (3.7)

Proof . Choose 𝐰_ℎ=𝜃 in (3.6)(a) and 𝐳_ℎ=−𝜉_𝑡 in (3.6)(b), and add the two equations to obtain 12𝑑𝑎𝑑𝑡^1/2𝜉²_𝐋2(Ω)+‖∇⋅𝜃‖²_𝐿2(Ω)=12𝑎_𝑡𝜉,𝜉+𝛾,𝜉_𝑡−𝑎𝜌,𝜉_𝑡−𝜌_𝑡,𝜃−𝑏^𝑡₀𝜉𝑑𝑠,𝜉_𝑡−𝑏^𝑡₀𝜌𝑑𝑠,𝜉_𝑡=12𝑎_𝑡+𝑑𝜉,𝜉−𝛾𝑑𝑡𝛾,𝜉_𝑡−𝑑,𝜉+𝑑𝑡𝑎𝜌,𝜉𝑎𝜌_𝑡+𝑎_𝑡−𝜌𝜌,𝜉_𝑡−𝑑,𝜃𝑏𝑑𝑡^𝑡₀𝜉𝑑𝑠,𝜉+𝑏_𝑡^𝑡₀+()−𝑑𝜉𝑑𝑠,𝜉𝑏𝜉,𝜉𝑏𝑑𝑡^𝑡₀𝜌𝑑𝑠,𝜉+𝑏_𝑡^𝑡₀𝜌𝑑𝑠,𝜉+𝑏𝜌,𝜉. (3.8) Integrate with respect to time from 0 to 𝑡 and apply the Cauchy-Schwarz's inequality and the Young's inequality to obtain 𝑎^1/2𝜉²_𝐋2(Ω)+2^𝑡₀‖∇⋅𝜃‖²_𝐿2(Ω)≤𝑎𝑑𝑠_1/2𝜉(0)²_𝐋2(Ω)+𝜌²_𝐋2(Ω)+𝛾²_𝐋2(Ω)+𝜀𝜉²_𝐋2(Ω)+𝐶^𝑡₀𝜌²_𝐋2(Ω)+𝜌_𝑠²_𝐋2(Ω)+𝛾_𝑠²_𝐋2(Ω)+𝜉²_𝐋2(Ω)+‖𝜃‖²_𝐋2(Ω)𝑑𝑠. (3.9)
Choose 𝐳_ℎ=𝜃 in (3.6)(b) to get ‖𝜃‖²_𝐋2(Ω)≤𝐶𝑎^1/2𝜉²_𝐋2(Ω)+𝛾²_𝐋2(Ω)+𝜌²_𝐋2(Ω)+^𝑡₀𝜉²_𝐋2(Ω)+𝜌²_𝐋2(Ω)𝑑𝑠. (3.10) Combining (3.9) and (3.10), we obtain 𝜉²_𝐋2(Ω)+‖𝜃‖²_𝐋2(Ω)+2^𝑡₀‖∇⋅𝜃‖²_𝐿2(Ω)≤𝑎𝑑𝑠_1/2𝜉(0)²_𝐋2(Ω)+𝜌²_𝐋2(Ω)+𝛾²_𝐋2(Ω)+𝜀𝜉²_𝐋2(Ω)+𝐶^𝑡₀𝜌²_𝐋2(Ω)+𝜌_𝑠²_𝐋2(Ω)+𝛾_𝑠²_𝐋2(Ω)+𝜉²_𝐋2(Ω)+‖𝜃‖²_𝐋2(Ω)𝑑𝑠. (3.11) Using the fact that 𝜉(0)=0 and Gronwall's lemma, we obtain 𝜉²_𝐋2(Ω)+‖𝜃‖²_𝐋2(Ω)+2^𝑡₀‖∇⋅𝜃‖²_𝐿2(Ω)≤𝜌𝑑𝑠_2𝐋2(Ω)+𝛾²_𝐋2(Ω)+𝐶^𝑡₀𝜌²_𝐋2(Ω)+𝜌_𝑠²_𝐋2(Ω)+𝛾_𝑠²_𝐋2(Ω)𝑑𝑠. (3.12)
Differentiating (3.6)(b) and taking 𝐳_ℎ=𝜉_𝑡, we obtain 𝜃_𝑡,𝜉_𝑡=𝑎𝜉_𝑡²_𝐋2(Ω)+𝑎_𝑡𝜉,𝜉_𝑡−𝛾_𝑡,𝜉_𝑡+𝑎_𝑡𝜌+𝑎𝜌_𝑡,𝜉_𝑡+𝑏_𝑡^𝑡₀𝜉𝑑𝑠,𝜉_𝑡+𝑏𝜉,𝜉_𝑡+𝑏_𝑡^𝑡₀𝜌𝑑𝑠,𝜉_𝑡+𝑏𝜌,𝜉_𝑡. (3.13) Choosing 𝐰_ℎ=𝜃_𝑡 in (3.6)(a), we obtain 𝜉_𝑡,𝜃_𝑡+12𝑑‖‖𝑑𝑡∇⋅𝜃²_𝐿2(Ω)𝑑=−𝜌𝑑𝑡^𝑡,𝜃+𝜌_𝑡𝑡,𝜃. (3.14) Add (3.13) and (3.14) to obtain 𝑎𝜉_𝑡²_𝐋2(Ω)+12𝑑‖‖𝑑𝑡∇⋅𝜃²_𝐿2(Ω)𝑑=−𝜌𝑑𝑡^𝑡,𝜃+𝜌_𝑡𝑡,𝜃−𝑎_𝑡𝜉,𝜉_𝑡+𝛾_𝑡,𝜉_𝑡−𝑎_𝑡𝜌+𝑎𝜌_𝑡,𝜉_𝑡−𝑏_𝑡^𝑡₀𝜉𝑑𝑠,𝜉_𝑡−𝑏𝜉,𝜉_𝑡−𝑏_𝑡^𝑡₀𝜌𝑑𝑠,𝜉_𝑡−𝑏𝜌,𝜉_𝑡. (3.15) Integrate (3.15) with respect to time from 0 to 𝑡 to obtain ^𝑡₀𝜉_𝑠²_𝐋2(Ω)𝑑𝑠+‖∇⋅𝜃‖²_𝐿2(Ω)≤𝜌²_𝐋2(Ω)+𝛾²_𝐋2(Ω)+𝐶𝜉²_𝐋2(Ω)+^𝑡₀𝜉²_𝐋2(Ω)𝑑𝑠+𝐶^𝑡₀𝜌²_𝐋2(Ω)+𝜌_𝑠²_𝐋2(Ω)+𝜌_𝑠𝑠²_𝐋2(Ω)+𝛾_𝑠²_𝐋2(Ω)𝑑𝑠. (3.16) Substitute (3.12) into (3.16) to have ^𝑡₀𝜉_𝑠²_𝐋2(Ω)𝑑𝑠+‖∇⋅𝜃‖²_𝐿2(Ω)≤𝜌²_𝐋2(Ω)+𝛾²_𝐋2(Ω)+𝐶^𝑡₀𝜌²_𝐋2(Ω)+𝛾²_𝐋2(Ω)+𝜌_𝑠²_𝐋2(Ω)+𝜌_𝑠𝑠²_𝐋2(Ω)+𝛾_𝑠²_𝐋2(Ω)𝑑𝑠. (3.17) Choosing 𝑣_ℎ=𝜍 in (3.6)(b) and applying Cauchy-Schwarz's inequality, we obtain ‖𝜍‖_𝐿2(Ω)≤∇⋅𝜎−𝜎_ℎ𝐿2(Ω). (3.18) Using Lemma 3.1, (3.17), and Gronwall's lemma, we get ‖𝜍‖_{𝐿∞(𝐿2(Ω))}≤𝐶ℎ^𝑟∗_𝜎+ℎ_𝑢^𝑘+1. (3.19) Using (3.12), (3.17), (3.19), (3.2), (3.3), and Lemma 3.1, we apply the triangle inequality to complete the proof.

4. Fully Discrete Error Estimates

In this section, we get the error estimates of fully discrete schemes. For the backward Euler procedure, let 0=𝑡₀<𝑡₁<𝑡₂<⋯<𝑡_𝑀=𝑇 be a given partition of the time interval [0,𝑇] with step length 𝛿=𝑇/𝑀, for some positive integer 𝑀. For a smooth function 𝜙 on [0,𝑇], define 𝜙^𝑛=𝜙(𝑡_𝑛) and 𝜕_𝑡𝜙^𝑛=(𝜙^𝑛−𝜙^𝑛−1)/𝛿. For approximating the integrals, we use the composite left rectangle rule 𝛿^𝑛−1_𝑗=0𝜙^𝑗≈^𝑡𝑛₀𝜙(𝑠)𝑑𝑠. (4.1) Note that 𝜙∈𝐶¹[0,𝑇], the quadrature error satisfies 𝛿^𝑛−1_𝑗=0𝜙^𝑗−^𝑡𝑛₀𝜙(𝑠)𝑑𝑠≤𝐶𝛿^𝑡𝑛₀𝜙_𝑠(𝑠)𝑑𝑠. (4.2)

Equation (2.5) has the following equivalent formulation: (a)𝜕_𝑡𝜆^𝑛,𝐰+∇⋅𝜎^𝑛,∇⋅𝐰=−𝑓^𝑛,∇⋅𝐰+𝑅^𝑛₁(b𝜎,𝐰,∀𝐰∈𝐖,)^𝑛,𝐳−𝑎^𝑛𝜆^𝑛,𝐳−𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜆^𝑗,𝐳=−𝑅^𝑛₃(c)𝜕,𝐳,∀𝐳∈𝐖,_𝑡𝑢^𝑛,𝑣−∇⋅𝜎^𝑛,𝑣=𝑓^𝑛,𝑣+𝑅^𝑛₂,𝑣,∀𝑣∈𝐿²(Ω), (4.3) where 𝑅^𝑛₁=𝜕_𝑡𝜆^𝑛−𝜆_𝑡𝛿𝜕=𝑂²𝜆𝜕𝑡²,𝑅^𝑛₂=𝜕_𝑡𝑢^𝑛−𝑢_𝑡𝛿𝜕=𝑂²𝑢𝜕𝑡²,𝑅^𝑛₃=𝛿^𝑛−1_𝑗=0𝜆^𝑗−^𝑡𝑛₀𝜆(𝑠)𝑑𝑠=𝑂𝛿𝜕𝜆.𝜕𝑡 (4.4) Now we can formulate a fully discrete scheme based on (4.3).

Fully discrete scheme: find (𝑢^𝑛_ℎ,𝜆^𝑛_ℎ,𝜎^𝑛_ℎ)∈𝑋_ℎ𝑢×𝐕_ℎ𝜎×𝐕_ℎ𝜎,(𝑛=1,2,…,𝑀−1) such that (a)𝜕_𝑡𝜆^𝑛_ℎ,𝐰_ℎ+∇⋅𝜎^𝑛_ℎ,∇⋅𝐰_ℎ=−𝑓^𝑛,∇⋅𝐰_ℎ,∀𝐰_ℎ∈𝐕_ℎ𝜎,(b)𝜎^𝑛_ℎ,𝐳_ℎ−𝑎^𝑛𝜆^𝑛_ℎ,𝐳_ℎ−𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜆^𝑗_ℎ,𝐳_ℎ=0,∀𝐳_ℎ∈𝐕_ℎ𝜎,(c)𝜕_𝑡𝑢^𝑛_ℎ,𝑣_ℎ−∇⋅𝜎^𝑛_ℎ,𝑣_ℎ=𝑓^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑋_ℎ𝑢, (4.5) with given an initial approximation (𝑢⁰_ℎ,𝜆⁰_ℎ,𝜎⁰_ℎ)∈𝑋_ℎ𝑢×𝐕_ℎ𝜎×𝐕_ℎ𝜎.

For fully discrete error estimates, we now split the errors 𝑢^𝑛−𝑢^𝑛_ℎ=𝑢^𝑛−𝑃_ℎ𝑢^𝑛+𝑃_ℎ𝑢^𝑛−𝑢^𝑛_ℎ=𝜂^𝑛+𝜍^𝑛,𝜆^𝑛−𝜆^𝑛_ℎ=𝜆^𝑛−𝑅_ℎ𝜆^𝑛+𝑅_ℎ𝜆^𝑛−𝜆^𝑛_ℎ=𝜌^𝑛+𝜉^𝑛,𝜎^𝑛−𝜎^𝑛_ℎ=𝜎^𝑛−𝑅_ℎ𝜎^𝑛+𝑅_ℎ𝜎^𝑛−𝜎^𝑛_ℎ=𝛾^𝑛+𝜃^𝑛. (4.6) From (4.3) to (4.5), we then obtain (a)𝜕_𝑡𝜉^𝑛,𝐰_ℎ+∇⋅𝜃^𝑛,∇⋅𝐰_ℎ=−𝜕_𝑡𝜌^𝑛,𝐰_ℎ+𝑅^𝑛₁,𝐰_ℎ,∀𝐰_ℎ∈𝐕_ℎ𝜎,(b)𝜃^𝑛,𝐳_ℎ−𝑎^𝑛𝜉^𝑛,𝐳_ℎ−𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜉^𝑗,𝐳_ℎ=−𝛾^𝑛,𝐳_ℎ+𝑎^𝑛𝜌^𝑛,𝐳_ℎ+𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜌^𝑗,𝐳_ℎ−𝑅^𝑛₃,𝐳_ℎ,∀𝐳_ℎ∈𝐕_ℎ𝜎,(c)𝜕_𝑡𝜍^𝑛,𝑣_ℎ−∇⋅𝜃^𝑛+𝛾^𝑛,𝑣_ℎ=𝑅^𝑛₂,𝑣_ℎ,∀𝑣_ℎ∈𝑋_ℎ𝑢. (4.7)

Lemma 4.1. Assume that the solution of system (2.5) has regular properties that 𝑢_𝑡∈𝐿²(𝐻^𝑘+1(Ω)),𝜆_𝑡, 𝜎_𝑡∈𝐿²(𝐇^𝑟+1(Ω)). Then one has the estimates max_{0≤𝑛≤𝑀}𝜕_𝑡𝜆−𝑅_ℎ𝜆^𝑛_𝐋2(Ω)+max_{0≤𝑛≤𝑀}𝜕_𝑡𝜎−𝑅_ℎ𝜎^𝑛_𝐋2(Ω)≤𝐶ℎ_𝜎^𝑟+1,max_{0≤𝑛≤𝑀}𝜕_𝑡𝑢−𝑃_ℎ𝑢^𝑛_𝐿2(Ω)≤𝐶ℎ_𝑢^𝑘+1. (4.8)

Theorem 4.2. Assume that 𝜕²𝑢/𝜕𝑡²,𝜕𝑢/𝜕𝑡∈𝐿²(𝐻^𝑘+1(Ω)),𝜕𝜆/𝜕𝑡,𝜕²𝜆/𝜕𝑡²,𝜕𝜎/𝜕𝑡,𝜕²𝜎/𝜕𝑡²∈𝐿²(𝐇^𝑟+1(Ω)), 𝑢∈𝐿^∞(𝐻^𝑘+1(Ω)), and 𝜆,𝜎∈𝐿^∞(𝐇^𝑟+1(Ω)), then there exists a constant C such that max_{0≤𝑛≤𝑀}𝜆−𝜆_ℎ^𝑛_𝐋2(Ω)+max_{0≤𝑛≤𝑀}𝜎−𝜎_ℎ^𝑛_𝐋2(Ω)≤𝐶ℎ_𝜎^𝑟+1+𝛿,max_{0≤𝑛≤𝑀}∇⋅𝜎−𝜎_ℎ^𝑛_𝐿2(Ω)≤𝐶ℎ^𝑟∗_𝜎+𝛿,max_{0≤𝑛≤𝑀}𝑢−𝑢_ℎ^𝑛_𝐿2(Ω)≤𝐶ℎ_𝑢^𝑘+1+ℎ^𝑟∗_𝜎+𝛿. (4.9)

Proof . Set 𝐰_ℎ=𝜃^𝑛 in (4.7)(a) and 𝐳_ℎ=−𝜕_𝑡𝜉^𝑛 in (4.7)(b) and add the two equations to obtain ∇⋅𝜃^𝑛2_𝐿2(Ω)+𝑎^𝑛𝜕_𝑡𝜉^𝑛,𝜉^𝑛=−𝜕_𝑡𝜌^𝑛,𝜃^𝑛+𝑅^𝑛₁,𝜃^𝑛+𝛾^𝑛,𝜕_𝑡𝜉^𝑛−𝑎^𝑛𝜌^𝑛,𝜕_𝑡𝜉^𝑛−𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜉^𝑗,𝜕_𝑡𝜉^𝑛−𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜌^𝑗,𝜕_𝑡𝜉^𝑛+𝑅^𝑛₃,𝜕_𝑡𝜉^𝑛. (4.10) Note that 𝜕_𝑡𝑎^𝑛𝜉^𝑛2_𝐋2(Ω)=𝑎^𝑛𝜉^𝑛,𝜉^𝑛−𝑎^𝑛−1𝜉^𝑛−1,𝜉^𝑛−1𝛿=𝑎^𝑛𝜉^𝑛,𝜉^𝑛−𝑎^𝑛𝜉^𝑛,𝜉^𝑛−1+𝑎^𝑛𝜉^𝑛,𝜉^𝑛−1−𝑎^𝑛−1𝜉^𝑛−1,𝜉^𝑛−1𝛿=𝑎^𝑛𝜕_𝑡𝜉^𝑛,𝜉^𝑛+𝑎^𝑛𝜕_𝑡𝜉^𝑛,𝜉^𝑛−1+𝜕_𝑡𝑎^𝑛𝜉^𝑛−1,𝜉^𝑛−1=2𝑎^𝑛𝜕_𝑡𝜉^𝑛,𝜉^𝑛−𝑎^𝑛1/2𝜉^𝑛−𝜉^𝑛−12_𝐋2(Ω)𝛿+𝜕_𝑡𝑎^𝑛𝜉^𝑛−1,𝜉^𝑛−1. (4.11) So, we get 𝑎^𝑛𝜕_𝑡𝜉^𝑛,𝜉^𝑛=12𝜕_𝑡𝑎^𝑛𝜉^𝑛2_𝐋2(Ω)+𝑎^𝑛1/2𝜉^𝑛−𝜉^𝑛−12_𝐋2(Ω)−12𝛿2𝜕_𝑡𝑎^𝑛𝜉^𝑛−1,𝜉^𝑛−1. (4.12) Note that 𝛾^𝑛,𝜕_𝑡𝜉^𝑛=𝜉^𝑛,𝛾^𝑛−𝜉^𝑛−1,𝛾^𝑛−1𝛿−𝜕_𝑡𝛾^𝑛,𝜉^𝑛−1,𝑎^𝑛𝜌^𝑛,𝜕_𝑡𝜉^𝑛=𝜉^𝑛,𝑎^𝑛𝜌^𝑛−𝜉^𝑛−1,𝑎^𝑛−1𝜌^𝑛−1𝛿−𝜕_𝑡𝑎^𝑛𝜌^𝑛,𝜉^𝑛−1,𝑅^𝑛₃,𝜕_𝑡𝜉^𝑛=𝜉^𝑛,𝑅^𝑛₃−𝜉^𝑛−1,𝑅₃^𝑛−1𝛿−𝜕_𝑡𝑅^𝑛₃,𝜉^𝑛−1.𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜉^𝑗,𝜕_𝑡𝜉^𝑛=𝑏^𝑛𝛿∑^𝑛−1_𝑗=0𝜉^𝑗,𝜉^𝑛−𝑏^𝑛−1𝛿∑^𝑛−2_𝑗=0𝜉^𝑗,𝜉^𝑛−1𝛿−𝜕_𝑡𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜉^𝑗,𝜉^𝑛−1−𝑏^𝑛−1𝜉^𝑛−1,𝜉^𝑛−1,𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜌^𝑗,𝜕_𝑡𝜉^𝑛=𝑏^𝑛𝛿∑^𝑛−1_𝑗=0𝜌^𝑗,𝜉^𝑛−𝑏^𝑛−1𝛿∑^𝑛−2_𝑗=0𝜌^𝑗,𝜉^𝑛−1𝛿−𝜕_𝑡𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜌^𝑗,𝜉^𝑛−1−𝑏^𝑛−1𝜌^𝑛−1,𝜉^𝑛−1. (4.13) Substitute (4.12)-(4.13) into (4.10) to get ∇⋅𝜃^𝑛2_𝐿2(Ω)+12𝜕_𝑡𝑎^𝑛𝜉^𝑛2_𝐋2(Ω)+𝑎^𝑛1/2𝜉^𝑛−𝜉^𝑛−12_𝐋2(Ω)𝜕2𝛿=−_𝑡𝜌^𝑛,𝜃^𝑛+𝑅^𝑛₁,𝜃^𝑛+12𝜕_𝑡𝑎^𝑛𝜉^𝑛−1,𝜉^𝑛−1+𝜉^𝑛,𝛾^𝑛−𝜉^𝑛−1,𝛾^𝑛−1𝛿−𝜕_𝑡𝛾^𝑛,𝜉^𝑛−1−𝜉^𝑛,𝑎^𝑛𝜌^𝑛−𝜉^𝑛−1,𝑎^𝑛−1𝜌^𝑛−1𝛿+𝜕_𝑡𝑎^𝑛𝜌^𝑛,𝜉^𝑛−1−𝑏^𝑛𝛿∑^𝑛−1_𝑗=0𝜉^𝑗+𝜌^𝑗,𝜉^𝑛−𝑏^𝑛−1𝛿∑^𝑛−2_𝑗=0𝜉^𝑗+𝜌^𝑗,𝜉^𝑛−1𝛿+𝜕_𝑡𝑏^𝑛𝛿^𝑛−1_𝑗=0𝜉^𝑗+𝜌^𝑗,𝜉^𝑛−1+𝑏^𝑛−1𝜉^𝑛−1+𝜌^𝑛−1,𝜉^𝑛−1+𝜉^𝑛,𝑅^𝑛₃−𝜉^𝑛−1,𝑅₃^𝑛−1𝛿−𝜕_𝑡𝑅^𝑛₃,𝜉^𝑛−1. (4.14) Summing from 1 to 𝑛, we find that 𝑎^𝑛1/2𝜉^𝑛2_𝐋2(Ω)+2𝛿^𝑛_𝑗=1∇⋅𝜃^𝑗2_𝐿2(Ω)+^𝑛_𝑗=1𝑎^𝑛1/2𝜉^𝑛−𝜉^𝑛−12_𝐋2(Ω)≤𝑎^01/2𝜉⁰²_𝐋2(Ω)−𝜉⁰,𝛾⁰+𝜉⁰,𝑎⁰𝜌⁰+𝜀𝜉^𝑛2_𝐋2(Ω)+𝐶𝜌^𝑛2_𝐋2(Ω)+𝐶𝑅^𝑛₃²_𝐋2(Ω)+𝐶𝛾^𝑛2_𝐋2(Ω)+𝐶𝛿^𝑛−1_𝑗=1𝜉^𝑗2_𝐋2(Ω)+𝐶𝛿^𝑛_𝑗=1𝜃^𝑗2_𝐋2(Ω)+𝐶𝛿^𝑛_𝑗=1𝜕_𝑡𝜌^𝑗2_𝐋2(Ω)+𝜕_𝑡𝛾^𝑗2_𝐋2(Ω)+𝜕_𝑡𝑎^𝑗𝜌^𝑗2_𝐋2(Ω)+𝜌^𝑗2_𝐋2(Ω)+𝑅^𝑗₁²_𝐋2(Ω)+𝜕_𝑡𝑅^𝑗₃²_𝐋2(Ω). (4.15) Choose 𝐳_ℎ=𝜃^𝑛 in (3.6)(b) to get 𝜃^𝑛2_𝐋2(Ω)≤𝐶𝑎^𝑛1/2𝜉^𝑛2_𝐋2(Ω)+𝛾^𝑛2_𝐋2(Ω)+𝜌^𝑛2_𝐋2(Ω)+𝐶𝛿^𝑛−1_𝑗=1𝜉^𝑗2_𝐋2(Ω)+𝜌^𝑗2_𝐋2(Ω). (4.16) Substitute (4.16) into (4.15) and note that 𝜉⁰=0 to get 𝜉^𝑛2_𝐋2(Ω)+𝜃^𝑛2_𝐋2(Ω)+2𝛿^𝑛_𝑗=1∇⋅𝜃^𝑗2_𝐿2(Ω)≤𝐶𝜌^𝑛2_𝐋2(Ω)+𝐶𝛾^𝑛2_𝐋2(Ω)+𝐶𝑅^𝑛₃²_𝐋2(Ω)+𝐶𝛿^𝑛−1_𝑗=1𝜉^𝑗2_𝐋2(Ω)+𝐶𝛿^𝑛_𝑗=1𝜃^𝑗2_𝐋2(Ω)+𝐶𝛿^𝑛_𝑗=1𝜕_𝑡𝜌^𝑗2_𝐋2(Ω)+𝜕_𝑡𝛾^𝑗2_𝐋2(Ω)+𝜕_𝑡𝑎^𝑗𝜌^𝑗2_𝐋2(Ω)+𝜌^𝑗2_𝐋2(Ω)+𝑅^𝑗₁²_𝐋2(Ω)+𝜕_𝑡𝑅^𝑗₃²_𝐋2(Ω). (4.17) Using Gronwall's lemma, we obtain 𝜉^𝑛2_𝐋2(Ω)+(1−𝐶𝛿)𝜃^𝑛2_𝐋2(Ω)+2𝛿^𝑛_𝑗=1∇⋅𝜃^𝑗2_𝐋2(Ω)≤𝐶𝜌^𝑛2_𝐋2(Ω)+𝐶𝛾^𝑛2_𝐋2(Ω)+𝐶𝑅^𝑛₃²_𝐋2(Ω)+𝐶𝛿^𝑛_𝑗=1𝜕_𝑡𝜌^𝑗2_𝐋2(Ω)+𝜕_𝑡𝛾^𝑗2_𝐋2(Ω)+𝜕_𝑡𝑎^𝑗𝜌^𝑗2_𝐋2(Ω)+𝜌^𝑗2_𝐋2(Ω)+𝑅^𝑗₁²_𝐋2(Ω)+𝜕_𝑡𝑅^𝑗₃²_𝐋2(Ω). (4.18)
Note that (a)𝛿^𝑛_𝑗=1𝑅^𝑗₁²_𝐋2(Ω)≤𝐶𝛿^3𝑛_𝑗=1𝜕²𝜆𝜕𝑡²²_{𝐿2𝐋2(Ω)}≤𝐶𝛿²𝜕²𝜆𝜕𝑡²²_{𝐿∞𝐋2(Ω)},(b)𝛿^𝑛_𝑗=1𝜕_𝑡𝜌^𝑗2_𝐋2(Ω)≤𝐶𝛿^𝑛_𝑗=1𝜕𝜌^𝑗𝜕𝑡²_𝐋2(Ω)𝑇≤𝐶𝑀⋅𝑛ℎ_𝜎^2𝑟+2≤𝐶ℎ_𝜎^2𝑟+2,(c)𝛿^𝑛_𝑗=1𝜕_𝑡𝛾^𝑗2_𝐋2(Ω)≤𝐶𝛿^𝑛_𝑗=1𝜕𝛾^𝑗𝜕𝑡²_𝐋2(Ω)𝑇≤𝐶𝑀⋅𝑛ℎ_𝜎^2𝑟+2≤𝐶ℎ_𝜎^2𝑟+2,(d)𝛿^𝑛_𝑗=1𝜕_𝑡𝑅^𝑗₃²_𝐋2(Ω)=𝐶𝛿^𝑛_𝑗=1∫^𝑡𝑗_𝑡𝑗−1𝜆^𝑗−1−𝜆𝑑𝑠𝛿²_𝐋2(Ω)≤𝐶𝛿²𝜕𝜆𝜕𝑡²_𝐿^∞_𝐋²_(Ω). (4.19) Therefore, substituting the above estimates into (4.18) and choosing 𝛿₀ in such a way that for 0<𝛿≤𝛿₀,(1−𝐶𝛿)>0, we obtain 𝜃^𝑛2_𝐋2(Ω)+𝜉^𝑛2_𝐋2(Ω)+𝛿^𝑛_𝑗=1∇⋅𝜃^𝑗2_𝐿2(Ω)≤𝐶ℎ_𝜎^2𝑟+2+𝛿². (4.20)
By (4.7)(b), we obtain 𝜕_𝑡𝜃^𝑛,𝐳_ℎ=𝑎^𝑛𝜉^𝑛−𝑎^𝑛−1𝜉^𝑛−1𝛿,𝐳_ℎ−𝜕_𝑡𝛾^𝑛,𝐳_ℎ+𝑎^𝑛𝜌^𝑛−𝑎^𝑛−1𝜌^𝑛−1𝛿,𝐳_ℎ+𝑏^𝑛𝛿∑^𝑛−1_𝑗=0𝜉^𝑗+𝜌^𝑗,𝐳_ℎ−𝑏^𝑛−1𝛿∑^𝑛−2_𝑗=0𝜉^𝑗+𝜌^𝑗,𝐳_ℎ𝛿−𝜕_𝑡𝑅^𝑛₃,𝐳_ℎ=𝑎^𝑛𝜉^𝑛−𝑎^𝑛−1𝜉^𝑛−1𝛿,𝐳_ℎ−𝜕_𝑡𝛾^𝑛,𝐳_ℎ+𝑎^𝑛𝜌^𝑛−𝑎^𝑛−1𝜌^𝑛−1𝛿,𝐳_ℎ+𝑏^𝑛𝜉^𝑛−1+𝜌^𝑛−1,𝐳_ℎ+𝜕_𝑡𝑏^𝑛𝛿^𝑛−2_𝑗=0𝜉^𝑗+𝜌^𝑗,𝐳_ℎ−𝜕_𝑡𝑅^𝑛₃,𝐳_ℎ. (4.21) Set 𝐳_ℎ=𝜕_𝑡𝜉^𝑛 in (4.21) to obtain 𝜕_𝑡𝜃^𝑛,𝜕_𝑡𝜉^𝑛=𝑎^𝑛𝜉^𝑛−𝑎^𝑛−1𝜉^𝑛−1𝛿,𝜕_𝑡𝜉^𝑛−𝜕_𝑡𝛾^𝑛,𝜕_𝑡𝜉^𝑛+𝑎^𝑛𝜌^𝑛−𝑎^𝑛−1𝜌^𝑛−1𝛿,𝜕_𝑡𝜉^𝑛+𝑏^𝑛𝜉^𝑛−1+𝜌^𝑛−1,𝜕_𝑡𝜉^𝑛+𝜕_𝑡𝑏^𝑛𝛿^𝑛−2_𝑗=0𝜉^𝑗+𝜌^𝑗,𝜕_𝑡𝜉^𝑛−𝜕_𝑡𝑅^𝑛₃,𝜕_𝑡𝜉^𝑛=𝑎^𝑛1/2𝜕_𝑡𝜉^𝑛2_𝐋2(Ω)+𝜕_𝑡𝑎^𝑛𝜉^𝑛−1,𝜕_𝑡𝜉^𝑛+𝑎^𝑛𝜕_𝑡𝜌^𝑛,𝜕_𝑡𝜉^𝑛+𝜕_𝑡𝑎^𝑛𝜌^𝑛−1,𝜕_𝑡𝜉^𝑛−𝜕_𝑡𝛾^𝑛,𝜕_𝑡𝜉^𝑛+𝑏^𝑛𝜉^𝑛−1+𝜌^𝑛−1,𝜕_𝑡𝜉^𝑛+𝜕_𝑡𝑏^𝑛𝛿^𝑛−2_𝑗=0𝜉^𝑗+𝜌^𝑗,𝜕_𝑡𝜉^𝑛−𝜕_𝑡𝑅^𝑛₃,𝜕_𝑡𝜉^𝑛. (4.22) Set 𝐰_ℎ=𝜕_𝑡𝜃^𝑛 in (4.7)(a) to obtain 𝜕_𝑡𝜉^𝑛,𝜕_𝑡𝜃^𝑛+12𝜕_𝑡∇⋅𝜃^𝑛2_𝐿2(Ω)+∇⋅𝜃^𝑛−𝜃^𝑛−12_𝐿2(Ω)𝜕2𝛿=−_𝑡𝜌^𝑛,𝜕_𝑡𝜃^𝑛+𝑅^𝑛₁,𝜕_𝑡𝜃^𝑛. (4.23) Substitute (4.22) into (4.23) to get 𝑎^𝑛1/2𝜕_𝑡𝜉^𝑛2_𝐋2(Ω)+12𝜕_𝑡∇⋅𝜃^𝑛2_𝐿2(Ω)+∇⋅𝜃^𝑛−𝜃^𝑛−12_𝐿2(Ω)𝜕2𝛿=−_𝑡𝜌^𝑛,𝜕_𝑡𝜃^𝑛−𝜕_𝑡𝑎^𝑛𝜉^𝑛−1,𝜕_𝑡𝜉^𝑛−𝑎^𝑛𝜕_𝑡𝜌^𝑛,𝜕_𝑡𝜉^𝑛−𝜕_𝑡𝑎^𝑛𝜌^𝑛−1,𝜕_𝑡𝜉^𝑛+𝜕_𝑡𝛾^𝑛,𝜕_𝑡𝜉^𝑛−𝑏^𝑛𝜉^𝑛−1+𝜌^𝑛−1,𝜕_𝑡𝜉^𝑛+𝑅^𝑛₁,𝜕_𝑡𝜃^𝑛−𝜕_𝑡𝑏^𝑛𝛿^𝑛−2_𝑗=0𝜉^𝑗+𝜌^𝑗,𝜕_𝑡𝜉^𝑛+𝜕_𝑡𝑅^𝑛₃,𝜕_𝑡𝜉^𝑛. (4.24) Take 𝐳_ℎ=𝜕_𝑡𝜃^𝑛 in (4.21) to obtain 𝜕_𝑡𝜃^𝑛2_𝐋2(Ω)=𝑎^𝑛𝜉^𝑛−𝑎^𝑛−1𝜉^𝑛−1𝛿,𝜕_𝑡𝜃^𝑛−𝜕_𝑡𝛾^𝑛,𝜕_𝑡𝜃^𝑛+𝑎^𝑛𝜌^𝑛−𝑎^𝑛−1𝜌^𝑛−1𝛿,𝜕_𝑡𝜃^𝑛+𝑏^𝑛𝜉^𝑛−1+𝜌^𝑛−1,𝜕_𝑡𝜃^𝑛+𝜕_𝑡𝑏^𝑛𝛿^𝑛−2_𝑗=0𝜉^𝑗+𝜌^𝑗,𝜕_𝑡𝜃^𝑛−𝜕_𝑡𝑅^𝑛₃,𝜕_𝑡𝜃^𝑛=𝑎^𝑛𝜕_𝑡𝜉^𝑛,𝜕_𝑡𝜃^𝑛+𝜕_𝑡𝑎^𝑛𝜉^𝑛−1,𝜕_𝑡𝜃^𝑛−𝜕_𝑡𝛾^𝑛,𝜕_𝑡𝜃^𝑛+𝑎^𝑛𝜕_𝑡𝜌^𝑛,𝜕_𝑡𝜃^𝑛+𝜕_𝑡𝑎^𝑛𝜌^𝑛−1,𝜕_𝑡𝜃^𝑛+𝑏^𝑛𝜉^𝑛−1+𝜌^𝑛−1,𝜕_𝑡𝜃^𝑛+𝜕_𝑡𝑏^𝑛𝛿^𝑛−2_𝑗=0𝜉^𝑗+𝜌^𝑗,𝜕_𝑡𝜃^𝑛−𝜕_𝑡𝑅^𝑛₃,𝜕_𝑡𝜃^𝑛. (4.25) Add (4.24) and (4.25) to get 𝑎^𝑛1/2𝜕_𝑡𝜉^𝑛2_𝐋2(Ω)+𝜕_𝑡𝜃^𝑛2_𝐋2(Ω)+12𝜕_𝑡∇⋅𝜃^𝑛2_𝐿2(Ω)+∇⋅𝜃^𝑛−𝜃^𝑛−12_𝐿2(Ω)𝜕2𝛿=−_𝑡𝜌^𝑛,𝜕_𝑡𝜃^𝑛−𝜕_𝑡𝑎^𝑛𝜉^𝑛−1,𝜕_𝑡𝜉^𝑛−𝑎^𝑛𝜕_𝑡𝜌^𝑛,𝜕_𝑡𝜉^𝑛−𝜕_𝑡𝑎^𝑛𝜌^𝑛−1,𝜕_𝑡𝜉^𝑛+𝜕_𝑡𝛾^𝑛,𝜕_𝑡𝜉^𝑛+𝑎^𝑛𝜕_𝑡𝜉^𝑛,𝜕_𝑡𝜃^𝑛+𝜕_𝑡𝑎^𝑛𝜉^𝑛−1,𝜕_𝑡𝜃^𝑛−𝜕_𝑡𝛾^𝑛,𝜕_𝑡𝜃^𝑛+𝑎^𝑛𝜕_𝑡𝜌^𝑛,𝜕_𝑡𝜃^𝑛+𝜕_𝑡𝑎^𝑛𝜌^𝑛−1,𝜕_𝑡𝜃^𝑛+𝑅^𝑛₁,𝜕_𝑡𝜃^𝑛−𝜕_𝑡𝑅^𝑛₃,𝜕_𝑡𝜃^𝑛+𝜕_𝑡𝜉^𝑛+𝑏^𝑛𝜉^𝑛−1+𝜌^𝑛−1,𝜕_𝑡𝜃^𝑛+𝜕_𝑡𝜉^𝑛+𝜕_𝑡𝑏^𝑛𝛿^𝑛−2_𝑗=0𝜉^𝑗+𝜌^𝑗,𝜕_𝑡𝜃^𝑛+𝜕_𝑡𝜉^𝑛. (4.26) Apply Cauchy-Schwarz's inequality and Young's inequality to obtain 𝑎_min𝜕_𝑡𝜉^𝑛2_𝐋2(Ω)+𝜕_𝑡𝜃^𝑛2_𝐋2(Ω)+12𝜕_𝑡∇⋅𝜃^𝑛2_𝐿2(Ω)≤𝜀𝜕_𝑡𝜉^𝑛2_𝐋2(Ω)+𝜕_𝑡𝜃^𝑛2_𝐋2(Ω)+𝐶𝜉^𝑛−12_𝐋2(Ω)+𝜕_𝑡𝜌^𝑛2_𝐋2(Ω)+𝜕_𝑡𝛾^𝑛2_𝐋2(Ω)+𝜌^𝑛−12_𝐋2(Ω)+𝑅^𝑛₁²_𝐋2(Ω)+𝜕_𝑡𝑅^𝑛₃²_𝐋2(Ω)+𝐶𝛿^2𝑛−2_𝑗=0𝜉^𝑗2_𝐋2(Ω)+𝜌^𝑗2_𝐋2(Ω). (4.27) Using (4.19) and (4.20) and summing from 1 to 𝑛, we obtain ∇⋅𝜃^𝑛2_𝐿2(Ω)+𝛿^𝑛_𝑗=1𝜕_𝑡𝜉^𝑛2_𝐋2(Ω)+𝜕_𝑡𝜃^𝑛2_𝐋2(Ω)≤𝐶ℎ_𝜎^2𝑟+2+𝛿². (4.28) Choosing 𝑣_ℎ=𝜍^𝑛 in (4.7)(c) and applying Cauchy-Schwarz's inequality, Young's inequality, and (4.28), we have 𝜍^𝑛2_𝐿2(Ω)−𝜍^𝑛−12_𝐿2(Ω)=2𝛿∇⋅𝜎−𝜎_ℎ^𝑛,𝜍^𝑛+2𝛿𝑅^𝑛₂,𝜍^𝑛≤𝐶𝛿∇⋅𝜎−𝜎_ℎ^𝑛2_𝐿2(Ω)+𝑅^𝑛₂²_𝐿2(Ω)+𝜍^𝑛2_𝐿2(Ω). (4.29) Summing from 1 to 𝑛 and using the Gronwall lemma, we obtain 𝜍^𝑛2_𝐿2(Ω)≤𝜍⁰²_𝐿2(Ω)+𝐶𝛿^𝑛_𝑗=1∇⋅𝜎−𝜎_ℎ^𝑗2_𝐿2(Ω)+𝑅^𝑗₂²_𝐿2(Ω). (4.30) Note that 𝛿^𝑛_𝑗=1𝑅^𝑗₂²_𝐿2(Ω)≤𝐶𝛿^3𝑛_𝑗=1𝜕²𝑢𝜕𝑡²²_{𝐿2𝐿2(Ω)}≤𝐶𝛿²𝜕²𝑢𝜕𝑡²²_{𝐿∞𝐿2(Ω)}. (4.31) Substituting (4.31) into (4.30) and using (3.2), (4.20), and the triangle inequality, we get 𝜍^𝑛2_𝐿2(Ω)≤𝐶ℎ_𝑢^2𝑘+2+ℎ_𝜎^2𝑟∗+𝛿². (4.32) Combining (3.2), (3.3), (4.20), (4.28), (4.32), and Lemma 4.1, we apply the triangle inequality to complete the proof.

5. Numerical Example

In this section, we analyse some numerical results to illustrate the efficiency of the proposed method. We consider the following 2D parabolic partial integrodifferential equations with initial-boundary value condition: 𝑢_𝑡−∇⋅𝑎(𝑥,𝑡)∇𝑢+𝑏(𝑥,𝑡)^𝑡₀(),()𝑢()()∇𝑢𝑑𝑠=𝑓𝑥,𝑡𝑥,𝑡∈Ω×𝐽,𝑥,𝑡=0,𝑥,𝑡∈𝜕Ω×𝑢()𝐽,𝑥,0=sin𝜋𝑥₁sin𝜋𝑥₂,𝑥∈Ω, (5.1) where Ω=[0,1]×[0,1], 𝐽=(0,1], 𝑎(𝑥,𝑡)=1+𝑥²₁+2𝑥²₂, 𝑏(𝑥,𝑡)=1+2𝑥²₁+𝑥²₂, 𝑥=(𝑥₁,𝑥₂), and 𝑓(𝑥,𝑡) is chosen so that the exact solution for the scalar unknown function is 𝑢(𝑥,𝑡)=𝑒^−𝑡sin𝜋𝑥₁sin𝜋𝑥₂. (5.2) The corresponding exact gradient is 𝜆(𝑥,𝑡)=𝜆₁,𝜆₂=𝜋𝑒^−𝑡cos𝜋𝑥₁sin𝜋𝑥₂,𝜋𝑒^−𝑡sin𝜋𝑥₁cos𝜋𝑥₂, (5.3) and its exact flux is 𝜎(𝑥,𝑡)=𝜎₁,𝜎₂=𝜋1+2−𝑒^−𝑡𝑥²₁+1+𝑒^−𝑡𝑥²₂cos𝜋𝑥₁sin𝜋𝑥₂,𝜋1+2−𝑒^−𝑡𝑥²₁+1+𝑒^−𝑡𝑥²₂sin𝜋𝑥₁cos𝜋𝑥₂. (5.4)

Dividing the domain Ω into the triangulations of mesh size ℎ_𝑢=ℎ_𝜎=ℎ uniformly, considering the piecewise constant space 𝑋_ℎ𝑢 with index 𝑘=0 for the scalar unknown function 𝑢 and the lowest-order Raviart-Thomas triangular space 𝐕_ℎ𝜎 [42, 45] for the gradient 𝜆 and the flux 𝜎 and using the backward Euler procedure with uniform time step length 𝛿=1/𝑀, we obtain some convergence results for ‖𝑢−𝑢_ℎ‖_{𝐿∞(𝐿2(Ω))}, ‖𝜆−𝜆_ℎ‖_{𝐿∞(𝐋2(Ω))}, ‖𝜎−𝜎_ℎ‖_{𝐿∞(𝐋2(Ω))} and ‖𝜎−𝜎_ℎ‖_{𝐿∞(𝐻(div;Ω))} with ℎ=√2𝛿=√2/8,√2/16,√2/32 in Table 1. With time 𝑡=1, ℎ=√2𝛿=√2/32, the exact solution 𝑢, 𝜆, 𝜎 is shown in Figures 1, 3, and 5, respectively, and the corresponding numerical solution 𝑢_ℎ, 𝜆_ℎ, 𝜎_ℎ is shown in Figures 2, 4, and 6, respectively.

Table 1 

The errors and convergence order.

Figure 1 

The exact solution 𝑢.

Figure 2 

The numerical solution 𝑢ℎ.

(a)

(b)

Figure 3 

The exact gradient 𝜆1(a) and 𝜆2(b).

(a)

(b)

Figure 4 

The numerical gradient 𝜆1ℎ (a) and 𝜆2ℎ (b).

(a)

(b)

Figure 5 

The exact flux 𝜎1 (a) and 𝜎2 (b).

(a)

(b)

Figure 6 

The numerical flux 𝜎1ℎ(a) and 𝜎2ℎ(b).

We can see from Table 1 that the convergence rate is order 1 which confirms the theoretical results of Theorem 4.2 for the above chosen spaces 𝑋_ℎ𝑢 and 𝐕_ℎ𝜎. The numerical results in Table 1 and Figures 1–6 show that new positive definite expanded mixed scheme is efficient.

6. Concluding Remarks

In the paper, a new expanded mixed finite element method based on a positive definite system is proposed for parabolic partial integrodifferential equation. Compared to expanded mixed method and standard mixed methods, the new expanded mixed element system is symmetric positive definite and both the gradient equation and the flux equation are separated from its scalar unknown equation. The existence and uniqueness for semidiscrete scheme are proved and error estimates are derived for both semidiscrete and fully discrete schemes. Finally, some numerical results are provided to confirm our theoretical analysis. In the near further, we will study the others evolution equations such as hyperbolic wave equation, and miscible displacement of compressible flow in porous media.

Acknowledgments

The authors would like to thank the Editor and the referees for their valuable suggestions and comments. This work is supported by the National Natural Science Fund of China (no. 11061021), the Scientific Research Projection of Higher Schools of Inner Mongolia (no. NJ10006, no. NJ10016, no. NJZZ12011), the Program of Higher-level talents of Inner Mongolia University (SPH-IMU, no. Z200901004, no. 125119), and the YSF of Inner Mongolia University (no. ND0702).