A Novel Characteristic Expanded Mixed Method for Reaction-Convection-Diffusion Problems

Abstract

A novel characteristic expanded mixed finite element method is proposed and analyzed for reaction-convection-diffusion problems. The diffusion term ∇⋅(𝑎(𝐱,𝑡)∇𝑢) is discretized by the novel expanded mixed method, whose gradient belongs to the square integrable space instead of the classical 𝐇(div;Ω) space and the hyperbolic part 𝑑(𝐱)(𝜕𝑢/𝜕𝑡)+𝐜(𝐱,𝑡)⋅∇𝑢 is handled by the characteristic method. For a priori error estimates, some important lemmas based on the novel expanded mixed projection are introduced. The fully discrete error estimates based on backward Euler scheme are obtained. Moreover, the optimal a priori error estimates in 𝐿²- and 𝐻¹-norms for the scalar unknown 𝑢 and a priori error estimates in (𝐿²)²-norm for its gradient 𝜆 and its flux 𝜎 (the coefficients times the negative gradient) are derived. Finally, a numerical example is provided to verify our theoretical results.

1. Introduction

In this paper, we consider the following reaction-convection-diffusion problems: 𝑑(𝐱)𝜕𝑢()(𝑎())()𝑢(),()𝑢()()𝜕𝑡+𝐜𝐱,𝑡⋅∇𝑢−∇⋅𝐱,𝑡∇𝑢+𝑅𝐱,𝑡=𝑓𝐱,𝑡𝐱,𝑡∈Ω×𝐽,𝐱,𝑡=0,𝐱,𝑡∈𝜕Ω×𝑢()𝐽,𝐱,0=𝑢₀(𝐱),𝐱∈Ω, (1) where Ω is a bounded convex polygonal domain in ℝ² with Lipschitz continuous boundary 𝜕Ω and 𝐽=(0,𝑇] is the time interval with 0<𝑇<∞. The initial value 𝑢₀(𝐱) and the source term 𝑓(𝐱,𝑡) are given functions. Throughout this paper, we assume that the accumulation, diffusion, and reaction coefficients, 𝑑=𝑑(𝐱), 𝑎=𝑎(𝐱,𝑡), and 𝑅=𝑅(𝐱,𝑡), satisfy 𝐇₁∶0<𝑑₀≤𝑑(𝐱)≤𝑑₁𝐇<+∞,₂∶0<𝑎₀≤𝑎(𝐱,𝑡)≤𝑎₁<+∞,𝑎_𝑡(𝐱,𝑡)≤𝑎₁,𝐇₃∶0≤𝑅₀≤𝑅(𝐱,𝑡)≤𝑅₁<+∞,𝑅_𝑡(𝐱,𝑡)≤𝑅₁, (2) and the bounded vector 𝐜(𝐱,𝑡)=(𝑐₁(𝐱,𝑡),𝑐₂(𝐱,𝑡)) satisfies 𝐇₄∶0<𝑏₀≤²_𝑖=1𝑐²_𝑖(𝐱,𝑡)^1/2≤𝑏₁<+∞. (3)

Reaction-convection-diffusion equations are a class of important evolution partial differential equations and have a lot of applications in many physical problems, such as the infiltration of liquid, the proliferation of gas, the conduction of heat, and the spread of impurities in semiconductor materials. In recent years, a lot of numerical methods, such as mixed finite element methods [1–3], Pod methods [4, 5], characteristic-mixed covolume methods [6], space-time discontinuous Galerkin methods [7, 8], least-squares finite element methods [9, 10], nonconforming finite element method [11], conforming rectangular element method [12], and two-grid expanded mixed methods [13–16], have been studied for reaction-convection-diffusion equations.

In 1994, Chen [17, 18] proposed an expanded mixed finite element method for second-order linear elliptic equation. Compared to standard mixed element methods, the expanded mixed method can approximate three variables simultaneously, namely, the scalar unknown, its gradient, and its flux (the tensor coefficient times the negative gradient). From then on, the expanded mixed method was applied to many evolution equations [2, 19, 20], and some new numerical methods based on the Chen’s expanded mixed method were proposed, for example, expanded mixed covolume method [21], expanded mixed hybrid methods [22], positive definite expanded mixed method [23], and so on. From the above literature on the study of expanded mixed method, we can find that all papers were studied based on Chen’s expanded mixed method [17, 18].

In 2011, we developed a new expanded mixed finite element method [24] based on the mixed weak formulations [25–27]. In this paper, we develop and analyze a novel characteristic expanded mixed finite element method, which combines the novel expanded mixed method [24] applied to approximating the diffusion term and the characteristic method that handled the hyperbolic part, for reaction-convection-diffusion equations. The gradient for our method belongs to the square integrable space instead of the classical 𝐇(div;Ω) space of Chen’s expanded mixed method. We derive a priori error estimates based on backward Euler method. Moreover, we prove the optimal a priori error estimates in 𝐿²- and 𝐻¹-norms for the scalar unknown 𝑢 and a priori error estimates in (𝐿²)²-norm for its gradient 𝜆 and its flux 𝜎 (the coefficients times the negative gradient).

Throughout this paper, 𝐶 will denote a generic positive constant which does not depend on the spatial mesh parameter ℎ and time discretization parameter Δ𝑡. At the same time, we denote the natural inner product in 𝐿²(Ω) or (𝐿²(Ω))² by (⋅,⋅) with the corresponding norm ‖⋅‖. The other notations and definitions of Sobolev spaces as in [28] are used.

2. Novel Characteristic Expanded Mixed System

Introducing the two auxiliary variables 𝜆=∇𝑢, 𝜎=−(𝑎(𝐱,𝑡)∇𝑢)=−(𝑎(𝐱,𝑡)𝜆), we obtain the following first-order system for (1): 𝑑𝜕𝑢(),𝜕𝑡+𝐜⋅∇𝑢+∇⋅𝜎+𝑅𝑢=𝑓𝐱,𝑡𝜆−∇𝑢=0,𝜎+𝑎𝜆=0. (4) Let the characteristic direction corresponding to the hyperbolic part 𝑑𝜕𝑢/𝜕𝑡+𝐜⋅∇𝑢 be denoted by 𝜏=𝜏(𝐱), which is subject to 𝜕(𝐱)=𝑑(𝐱)𝜕𝜏𝜓()𝜕𝐱,𝑡+𝐜(𝐱)𝜕𝑡𝜓(𝐱,𝑡)⋅∇, (5) where 𝜓(𝐱,𝑡)=(𝑑²(𝐱)+|𝐜(𝐱,𝑡)|²)^1/2.

Then the first-order system (4) can be rewritten as 𝜓𝜕𝑢(),𝜕𝜏+∇⋅𝜎+𝑅𝑢=𝑓𝐱,𝑡𝜆−∇𝑢=0,𝜎+𝑎𝜆=0. (6) Using the similar method to the one in [29], the expanded mixed weak formulation for problem (1) is to find {𝑢,𝜆,𝜎}∶[0,𝑇]↦𝑊×𝐕×𝐗 such that𝜓𝜕𝑢+()+()=()+()()+()𝜕𝜏,𝑣∇⋅𝜎,𝑣𝑅𝑢,𝑣𝑓,𝑣,∀𝑣∈𝑊,𝜆,𝐰𝑢,∇⋅𝐰=0,∀𝐰∈𝐕,𝜎,𝐳𝑎𝜆,𝐳=0,∀𝐳∈𝐗, (7a)(7b)(7c) where 𝑊=𝐿²(Ω) or 𝑊={𝑤∈𝐿²(Ω)|𝑤|_𝜕Ω=0}, 𝐕=𝐇(div;Ω)={𝐯∈(𝐿²(Ω))²∣∇⋅𝐯∈𝐿²(Ω)}, and 𝐗=(𝐿²(Ω))².

In this paper, we propose and discuss a novel characteristic expanded mixed method. The new characteristic expanded mixed weak formulation is to find {𝑢,𝜆,𝜎}∶[0,𝑇]↦𝐻¹₀×(𝐿²(Ω))²×(𝐿²(Ω))² such that𝜓𝜕𝑢−()+()=𝜕𝜏,𝑣𝜎,∇𝑣𝑅𝑢,𝑣𝑓,𝑣,∀𝑣∈𝐻¹₀()−()𝐿,𝜆,𝐰∇𝑢,𝐰=0,∀𝐰∈²(Ω)²,(𝜎,𝐳)+(𝑎𝜆,𝐳)=0,∀𝐳∈𝐿²(Ω)². (8a)(8b)(8c)For approximating the solution at time 𝑡_𝑛=𝑛Δ𝑡, the characteristic derivative will be approximated by 𝐱_𝑛=𝐱−Δ𝑡𝑑(𝐱)𝑐𝐱,𝑡_𝑛, (9) and then we have the following approximation: 𝜓𝐱,𝑡_𝑛𝜕𝑢𝜕𝜏_𝑡𝑛≈𝜓𝐱,𝑡_𝑛𝑢𝐱,𝑡_𝑛−𝑢𝐱_𝑛,𝑡_𝑛−1𝐱−𝐱_𝑛²+(Δ𝑡)²(𝐱)𝑢=𝑑𝐱,𝑡_𝑛−𝑢𝐱_𝑛,𝑡_𝑛−1.Δ𝑡 (10)

Let (𝑉_ℎ,𝐖_ℎ)⊂(𝐻¹₀,(𝐿²(Ω))²) be defined by the following finite element pair 𝑃₁−𝑃²₀ [25, 26]: 𝑉_ℎ=𝑣_ℎ∈𝐶⁰(Ω)∩𝐻¹₀∣𝑣_ℎ∈𝑃₁(𝐾),∀𝐾∈𝒦_ℎ,𝐖_ℎ=𝐰_ℎ=𝑤_1ℎ,𝑤_2ℎ∈𝐿²(Ω)²∣𝑤_1ℎ,𝑤_2ℎ∈𝑃₀(𝐾),∀𝐾∈𝒦_ℎ. (11)

Now the novel characteristic expanded mixed finite element procedure for (8a), (8b), and (8c) is to find {𝑢^𝑛_ℎ,𝜆^𝑛_ℎ,𝜎^𝑛_ℎ}∶[0,𝑇]↦𝑉_ℎ×𝐖_ℎ×𝐖_ℎ satisfying𝑑𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1Δ𝑡,𝑣_ℎ−𝜎^𝑛_ℎ,∇𝑣_ℎ+𝑅𝑢^𝑛_ℎ,𝑣_ℎ=𝑓^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝜆^𝑛_ℎ,𝐰_ℎ−∇𝑢^𝑛_ℎ,𝐰_ℎ=0,∀𝐰_ℎ∈𝐖_ℎ,𝜎^𝑛_ℎ,𝐳_ℎ+𝑎^𝑛𝜆^𝑛_ℎ,𝐳_ℎ=0,∀𝐳_ℎ∈𝐖_ℎ, (12a)(12b)(12c) where 𝑢^𝑛_ℎ=𝑢_ℎ(𝑡_𝑛), 𝑢_ℎ^𝑛−1=𝑢_ℎ(𝐱_𝑛,𝑡_𝑛−1)=𝑢_ℎ(𝐱−(Δ𝑡/𝑑(𝐱))𝑐(𝐱,𝑡_𝑛),𝑡_𝑛−1).

Remark 1. From [25, 26], we find that (𝑉_ℎ,𝐖_ℎ) satisfies the so-called discrete Ladyzhenskaya-Babuska-Brezzi condition.

Remark 2. Compared to the scheme (7a), (7b), and (7c) based on Chen’s expanded mixed element method, the gradient in the scheme (8a), (8b), and (8c) belongs to the simple square integrable space instead of the classical 𝐇(div;Ω) space. For 𝐇(div;Ω)⊂(𝐿²(Ω))², we easily find that our method reduces the regularity requirement on the gradient solution 𝜆=∇𝑢.

Remark 3. Based on finite element space (𝑉_ℎ,𝐖_ℎ) in (11), the number of total degrees of freedom for our scheme ((12a), (12b), and (12c)) is less than that for the scheme in [29]. By the same discussion as Remark  1 in [30], we can obtain the detailed analysis for degrees of freedom.

3. Some Lemmas and Error Estimates

3.1. Novel Expanded Mixed Projection and Lemmas

We first introduce the novel expanded mixed elliptic projection [24] associated with our equations to derive a priori error estimates for the proposed method.

Let (𝑢_ℎ,𝜆_ℎ,𝜎_ℎ)∶[0,𝑇]→𝑉_ℎ×𝐖_ℎ×𝐖_ℎ be given by the following mixed relations:𝜎−𝜎_ℎ,∇𝑣_ℎ=0,∀𝑣_ℎ∈𝑉_ℎ,𝜆−𝜆_ℎ,𝐰_ℎ−∇𝑢−𝑢_ℎ,𝐰_ℎ=0,∀𝐰_ℎ∈𝐖_ℎ𝜎,𝜎−_ℎ,𝐳_ℎ+𝑎𝜆−𝜆_ℎ,𝐳_ℎ=0,∀𝐳_ℎ∈𝐖_ℎ. (13a)(13b)(13c) In the following discussion, we will give some important lemmas based on the novel expanded mixed projection.

Lemma 4 (see [24]). There is a constant 𝐶 independent of ℎ such that ‖𝛿‖≤𝐶ℎ‖𝜆‖_(𝐻1(Ω))²+‖𝑢‖_𝐻²,𝜌≤𝐶ℎ‖𝑢‖_𝐻²+‖𝜎‖_(𝐻1(Ω))²+‖𝜆‖_(𝐻1(Ω))²,∇𝜂≤𝐶ℎ‖𝑢‖_𝐻²+‖𝜆‖_(𝐻1(Ω))², (14)(15)(16) where 𝜂=𝑢−𝑢_ℎ, 𝛿=𝜆−𝜆_ℎ, and 𝜌=𝜎−𝜎_ℎ.

Lemma 5 (see [24]). There is a constant 𝐶 independent of ℎ such that 𝜂≤𝐶ℎ²‖𝑢‖_𝐻²,𝜂_𝑡≤𝐶ℎ²𝑢_𝑡𝐻2,𝜂₁≤𝐶ℎ‖𝑢‖_𝐻²+‖𝜆‖_(𝐻1(Ω))². (17)

Lemma 6 (see [6]). For each 𝑛 one has (𝑑𝑣,𝑣)−(𝑑𝑣,𝑣)≤𝐶Δ𝑡(𝑑𝑣,𝑣),∀𝑣∈𝐿²(Ω), (18) where 𝑣=𝑣(𝐱−𝐜(𝐱,𝑡_𝑛)Δ𝑡/𝑑(𝐱)).

Remark 7. For proving Lemmas 4 and 5, we introduce two linear operators [25, 26] Π_ℎ∶(𝐿²(Ω))²→𝐖_ℎ and 𝑃_ℎ∶𝐻¹₀(Ω)→𝑉_ℎ and consider the following auxiliary elliptic problem: −∇⋅𝑎∇𝜒=𝜂,inΩ,𝜒=0,on𝜕Ω. (19) The detailed proofs for Lemmas 4 and 5 are shown in [24].

3.2. A Priori Error Estimates

In the following discussion, we will derive a priori error estimates based on fully discrete backward Euler method. Let 0=𝑡₀<𝑡₁<𝑡₂<⋯<𝑡_𝑀=𝑇 be a given partition of the time interval [0,𝑇] with step length Δ𝑡=𝑇/𝑀 and nodes 𝑡_𝑛=𝑛Δ𝑡, for some positive integer 𝑀. For a smooth function 𝜙 on [0,𝑇], define 𝜙^𝑛=𝜙(𝑡_𝑛).

In order to derive a priori error estimates, we now write 𝑢𝑡_𝑛−𝑢^𝑛_ℎ=𝑢𝑡_𝑛−𝑢^𝑛_ℎ+𝑢^𝑛_ℎ−𝑢^𝑛_ℎ=𝜂^𝑛+𝜍^𝑛,𝜆𝑡_𝑛−𝜆^𝑛_ℎ=𝜆𝑡_𝑛−𝜆^𝑛_ℎ+𝜆^𝑛_ℎ−𝜆^𝑛_ℎ=𝛿^𝑛+𝜃^𝑛,𝜎𝑡_𝑛−𝜎^𝑛_ℎ=𝜎𝑡_𝑛−𝜎^𝑛_ℎ+𝜎^𝑛_ℎ−𝜎^𝑛_ℎ=𝜌^𝑛+𝜉^𝑛. (20) Combining (8a), (8b), and (8c), (12a), (12b), and (12c) and (13a), (13b), and (13c) at 𝑡=𝑡_𝑛, we get the error equations𝑑𝜍^𝑛−𝜍̂^𝑛−1Δ𝑡,𝑣_ℎ−𝜉^𝑛,∇𝑣_ℎ+𝑅𝜍^𝑛,𝑣_ℎ=𝜓^𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1Δ𝑡,𝑣_ℎ−𝑑𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝑣_ℎ−𝑅𝜂^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑉_ℎ,𝜃^𝑛,𝐰_ℎ−∇𝜍^𝑛,𝐰_ℎ=0,∀𝐰_ℎ∈𝐖_ℎ,𝜉^𝑛,𝐳_ℎ+𝑎^𝑛𝜃^𝑛,𝐳_ℎ=0,∀𝐳_ℎ∈𝐖_ℎ. (21a)(21b)(21c) Based on the error equations (21a), (21b), and (21c), we obtain the following theorem for the fully discrete error estimates.

Theorem 8. Assume that 𝑢⁰_ℎ=𝑢_ℎ(0); then there exist some positive constants independent of ℎ=𝑂(Δ𝑡) such that 𝑢𝑡_𝐽−𝑢^𝐽_ℎ≤𝐶𝑎₀,𝑑₁𝑂ℎ²𝑎,Δ𝑡+𝐶₀,𝑑₁,𝑅₁ℎ²‖𝑢‖_𝐿^∞_(𝐻²₎,𝑢𝑡_𝐽−𝑢^𝐽_ℎ1≤𝐶𝑎₀,𝑑₁𝑂ℎ¹𝐶𝑎,Δ𝑡+ℎ₀,𝑑₁,𝑅₁‖𝑢‖_𝐿^∞_(𝐻²₎+𝐶‖𝜆‖_{𝐿∞((𝐻1)}²₎,𝜆𝑡_𝐽−𝜆^𝐽_ℎ≤𝐶𝑎₀,𝑑₁,𝑅₁𝑂ℎ¹𝐶𝑎,Δ𝑡+ℎ₀,𝑑₁,𝑅₁,𝑅₁‖𝑢‖_𝐿^∞_(𝐻²₎+𝐶‖𝜆‖_{𝐿∞((𝐻1)}²₎,𝜎𝑡_𝐽−𝜎^𝐽_ℎ≤𝐶𝑎₀,𝑑₁,𝑅₁𝑂ℎ¹𝐶‖𝜆‖,Δ𝑡+ℎ_{𝐿∞((𝐻1)}²₎+‖𝜎‖_{𝐿∞((𝐻1)}²₎+𝐶𝑎₀,𝑑₁,𝑅₁,𝑅₁‖𝑢‖_𝐿^∞_(𝐻²₎, (22) where |‖𝑂(ℎ^2−𝑗,Δ𝑡)‖|=ℎ^2−𝑗(||𝑢||_𝐿^∞_(𝐻²₎+||𝑢_𝑡||_𝐿²_(𝐻²₎)+Δ𝑡||𝑢_𝑡𝑡||_𝐿²_(𝐿²₎,𝑗=0,1.

Proof. Set 𝑣_ℎ=𝜍^𝑛 in (21a), 𝐰_ℎ=𝜉^𝑛 in (21b), and 𝐳_ℎ=𝜃^𝑛 in (21c) to obtain𝑑𝜍^𝑛−𝜍̂^𝑛−1Δ𝑡,𝜍^𝑛−𝜉^𝑛,∇𝜍^𝑛+𝑅^𝑛𝜍^𝑛,𝜍^𝑛=𝜓^𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1Δ𝑡,𝜍^𝑛−𝑑𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝜍^𝑛−𝑅^𝑛𝜂^𝑛,𝜍^𝑛,∇𝜍^𝑛,𝜉^𝑛−𝜃^𝑛,𝜉^𝑛𝜉=0,_𝑛,𝜃^𝑛+𝑎^𝑛1/2𝜃^𝑛2=0. (23a)(23b)(23c)Adding the above three equations, we obtain 𝑑𝜍^𝑛−𝜍̂^𝑛−1Δ𝑡,𝜍^𝑛+𝑎^𝑛1/2𝜃^𝑛2+𝑅^𝑛1/2𝜍^𝑛2=𝜓^𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1Δ𝑡,𝜍^𝑛−𝑑𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝜍^𝑛−𝑅^𝑛𝜂^𝑛,𝜍^𝑛. (24) Noting that −𝑑𝜍̂^𝑛−1𝜍^𝑛/Δ𝑡=(||𝑑^1/2(𝜍^𝑛−𝜍̂^𝑛−1)||²−(||𝑑^1/2𝜍̂^𝑛−1||²+||𝑑^1/2𝜍^𝑛||²))/2Δ𝑡, (24) may be written as 𝑑^1/2𝜍^𝑛2−𝑑^1/2𝜍̂^𝑛−122Δ𝑡+𝑎₀𝜃^𝑛2+𝑅^𝑛1/2𝜍^𝑛2≤𝑑^1/2𝜍^𝑛2−𝑑^1/2𝜍̂^𝑛−12+𝑑^1/2𝜍^𝑛−𝜍̂^𝑛−12+𝑎2Δ𝑡_𝑛1/2𝜃^𝑛2+𝑅^𝑛1/2𝜍^𝑛2=𝜓^𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1Δ𝑡,𝜍^𝑛−𝑑𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝜍^𝑛−𝑅^𝑛𝜂^𝑛,𝜍^𝑛. (25) Using the same analysis as that in [31], the right-hand side of (25) can be rewritten as 𝜓^𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1Δ𝑡,𝜍^𝑛−𝑑𝜂^𝑛−𝜂^𝑛−1Δ𝑡,𝜍^𝑛−𝑅^𝑛𝜂^𝑛,𝜍^𝑛=𝜓^𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1Δ𝑡,𝜍^𝑛−𝑑𝜂^𝑛−𝜂^𝑛−1+𝜂^𝑛−1−𝜂^𝑛−1Δ𝑡,𝜍^𝑛−𝑅^𝑛𝜂^𝑛,𝜍^𝑛≤𝜓^𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1+𝑑𝜂Δ𝑡^𝑛−𝜂^𝑛−1+𝑑𝜂Δ𝑡_𝑛−1−𝜂^𝑛−1+𝑅Δ𝑡^𝑛𝜂^𝑛𝜍^𝑛≤𝐶𝑑₁Δ𝑡^𝑡𝑛_𝑡𝑛−1𝜕²𝑢𝜕𝜏²²1𝑑𝑠+Δ𝑡^𝑡𝑛_𝑡𝑛−1𝜕𝜂𝜕𝜏²𝑑𝑠+𝜍^𝑛2+𝑅^𝑛𝜂^𝑛2+𝑑𝜂^𝑛−1−𝜂^𝑛−1𝜍Δ𝑡^𝑛. (26) Using the assumption ℎ=𝑂(Δ𝑡), we obtain 𝑑𝜂^𝑛−1−𝜂^𝑛−1𝜍Δ𝑡^𝑛≤𝐶ℎ+Δ𝑡𝜂Δ𝑡^𝑛−1𝜃^𝑛≤𝐶𝑎₀𝜂²_𝐿∞𝐿2+𝑎₀2‖𝜃‖². (27) By Lemma 6, we have 𝑑^1/2𝜍^𝑛−12−𝑑^1/2𝜍̂^𝑛−12≥−𝐶Δ𝑡𝑑^1/2𝜍^𝑛−12. (28) Using (28), we obtain 𝑑^1/2𝜍^𝑛2−𝑑^1/2𝜍̂^𝑛−122Δ𝑡+𝑎₀𝜃^𝑛2+𝑅^𝑛1/2𝜍^𝑛2≥𝑑^1/2𝜍^𝑛2−(1+𝐶Δ𝑡)𝑑^1/2𝜍^𝑛−122Δ𝑡+𝑎₀𝜃^𝑛2+𝑅^𝑛1/2𝜍^𝑛2=𝑑^1/2𝜍^𝑛2−𝑑^1/2𝜍^𝑛−12−𝐶Δ𝑡𝑑^1/2𝜍^𝑛−122Δ𝑡+𝑎₀𝜃^𝑛2+𝑅^𝑛1/2𝜍^𝑛2. (29) Multiplying by 2Δ𝑡, summing (29) from 𝑛=1 to 𝐽, and using (25)–(29), the resulting inequality becomes 𝑑₀𝜍^𝐽2+Δ𝑡^𝐽_𝑛=1𝑎₀𝜃^𝑛2+𝑅^𝑛1/2𝜍^𝑛2≤𝑑₁𝜍⁰²+𝐶𝑎₀,𝑑₁×Δ𝑡^𝐽_𝑛=1Δ𝑡^𝑡𝑛_𝑡𝑛−1𝜕²𝑢𝜕𝜏²²1𝑑𝑠+Δ𝑡^𝑡𝑛_𝑡𝑛−1𝜕𝜂𝜕𝜏²𝑎𝑑𝑠+𝐶₀,𝑑₁,𝑅₁𝜂²_𝐿∞𝐿2+𝐶𝑑₁Δ𝑡^𝐽_𝑛=1𝜍^𝑛−12=𝑑₁𝜍⁰²+𝐶𝑎₀,𝑑₁×(Δ𝑡)^2𝑡𝐽_𝑡0𝜕²𝑢𝜕𝜏²²𝑑𝑠+^𝑡𝐽_𝑡0𝜕𝜂𝜕𝜏²𝑎𝑑𝑠+𝐶₀,𝑑₁,𝑅₁𝜂²_𝐿∞𝐿2+𝐶𝑑₁Δ𝑡^𝐽_𝑛=1𝜍^𝑛−12. (30) Using Gronwall lemma, we have 𝑑₀𝜍^𝐽2+Δ𝑡^𝐽_𝑛=1𝑎₀𝜃^𝑛2+𝑅^𝑛1/2𝜍^𝑛2≤𝐶𝑎₀,𝑑₁(Δ𝑡)^2𝑡𝐽_𝑡0𝜕²𝑢𝜕𝜏²²𝑑𝑠+^𝑡𝐽_𝑡0𝜕𝜂𝜕𝜏²𝑎𝑑𝑠+𝐶₀,𝑑₁,𝑅₁𝜂²_𝐿∞𝐿2. (31) From (21b), we get 𝜃^𝑛−𝜃^𝑛−1Δ𝑡,𝐰_ℎ−∇𝜍^𝑛−∇𝜍^𝑛−1Δ𝑡,𝐰_ℎ=0. (32) Choose 𝐰_ℎ=𝜉^𝑛 in (32), 𝑣_ℎ=(𝜍^𝑛−𝜍^𝑛−1)/Δ𝑡 in (21a), and 𝐳_ℎ=(𝜃^𝑛−𝜃^𝑛−1)/Δ𝑡 in (21c) to obtain𝑑𝜍^𝑛−𝜍̂^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1−𝜉Δ𝑡^𝑛,∇𝜍^𝑛−∇𝜍^𝑛−1+𝑅Δ𝑡_𝑛𝜍^𝑛,𝜍^𝑛−𝜍^𝑛−1=𝜓Δ𝑡_𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1−𝑑𝜂Δ𝑡_𝑛−𝜂^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1−𝑅Δ𝑡^𝑛𝜂^𝑛,𝜍^𝑛−𝜍^𝑛−1,Δ𝑡∇𝜍^𝑛−∇𝜍^𝑛−1Δ𝑡,𝜉^𝑛−𝜃^𝑛−𝜃^𝑛−1Δ𝑡,𝜉^𝑛𝜉=0,_𝑛,𝜃^𝑛−𝜃^𝑛−1+𝑎Δ𝑡^𝑛𝜃^𝑛,𝜃^𝑛−𝜃^𝑛−1Δ𝑡=0. (33a)(33b)(33c)Adding the three equations and using Cauchy-Schwarz inequality and Young inequality, we obtain 𝑑𝜍^𝑛−𝜍̂^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1+𝑎Δ𝑡^𝑛𝜃^𝑛,𝜃^𝑛−𝜃^𝑛−1+𝑅Δ𝑡_𝑛𝜍^𝑛,𝜍^𝑛−𝜍^𝑛−1=𝜓Δ𝑡_𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1−𝑑𝜂Δ𝑡_𝑛−𝜂^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1−𝑅Δ𝑡^𝑛𝜂^𝑛,𝜍^𝑛−𝜍^𝑛−1≤𝜓Δ𝑡_𝑛𝜕𝑢^𝑛𝑢𝜕𝜏−𝑑^𝑛−𝑢^𝑛−1+𝑑𝜂Δ𝑡^𝑛−𝜂^𝑛−1+𝑑𝜂Δ𝑡_𝑛−1−𝜂^𝑛−1+𝑅Δ𝑡^𝑛𝜂^𝑛𝜍^𝑛−𝜍^𝑛−1𝑑Δ𝑡≤𝐶₁Δ𝑡^𝑡𝑛_𝑡𝑛−1𝜕²𝑢𝜕𝜏²²𝐶𝑑𝑑𝑠+¹Δ𝑡^𝑡𝑛_𝑡𝑛−1𝜕𝜂𝜕𝜏²𝑅𝑑𝑠+𝐶₁𝜂²_𝐿∞𝐿2+12𝜍^𝑛−𝜍^𝑛−1Δ𝑡². (34) The left-hand side of (34) can be written as 𝑑𝜍^𝑛−𝜍̂^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1+𝑎Δ𝑡^𝑛𝜃^𝑛,𝜃^𝑛−𝜃^𝑛−1+𝑅Δ𝑡_𝑛𝜍^𝑛,𝜍^𝑛−𝜍^𝑛−1=𝑑Δ𝑡_1/2𝜍^𝑛−𝜍^𝑛−1Δ𝑡²+𝑑𝜍^𝑛−1−𝜍̂^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1+𝑎Δ𝑡_𝑛1/2𝜃^𝑛−𝜃^𝑛−12+𝑎2Δ𝑡_𝑛1/2𝜃^𝑛2−𝑎^𝑛−11/2𝜃^𝑛−12−𝑎2Δ𝑡_𝑛−𝑎^𝑛−1𝜃2Δ𝑡^𝑛−1,𝜃^𝑛−1+𝑅^𝑛1/2𝜍^𝑛−𝜍^𝑛−12+𝑅2Δ𝑡_𝑛1/2𝜍^𝑛2−𝑅^𝑛−11/2𝜍^𝑛−12−𝑅2Δ𝑡_𝑛−𝑅^𝑛−1𝜍2Δ𝑡^𝑛−1,𝜍^𝑛−1. (35) Substitute (35) into (34) and multiply by 2Δ𝑡 to get Δ𝑡𝑑^1/2𝜍^𝑛−𝜍^𝑛−1Δ𝑡²+𝑎^𝑛1/2𝜃^𝑛−𝜃^𝑛−12+𝑅^𝑛1/2𝜍^𝑛−𝜍^𝑛−12+𝑎^𝑛1/2𝜃^𝑛2−𝑎^𝑛−11/2𝜃^𝑛−12+𝑅^𝑛1/2𝜍^𝑛2−𝑅^𝑛−11/2𝜍^𝑛−12≤𝐶𝑑₁(Δ𝑡)^2𝑡𝑛_𝑡𝑛−1𝜕²𝑢𝜕𝜏²²𝑑𝑑𝑠+𝐶₁^𝑡𝑛_𝑡𝑛−1𝜕𝜂𝜕𝜏²𝑑𝑠+𝐶𝑅₁Δ𝑡𝜂²_𝐿∞𝐿2𝑑𝜍−Δ𝑡^𝑛−1−𝜍̂^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1∫Δ𝑡+Δ𝑡^{𝑡𝑛𝑡𝑛−1}𝑎_𝑡𝑑𝑠𝜃2Δ𝑡^𝑛−1,𝜃^𝑛−1∫−Δ𝑡^{𝑡𝑛𝑡𝑛−1}𝑅_𝑡𝑑𝑠𝜍2Δ𝑡^𝑛−1,𝜍^𝑛−1. (36) Summing (36) from 𝑛=1 to 𝐽, the resulting inequality becomes Δ𝑡^𝐽_𝑛=1𝑑^1/2𝜍^𝑛−𝜍^𝑛−1Δ𝑡²+^𝐽_𝑛=1𝑎^𝑛1/2𝜃^𝑛−𝜃^𝑛−12+𝑅^𝑛1/2𝜍^𝑛−𝜍^𝑛−12+𝑎^𝐽1/2𝜃^𝐽2+𝑅^𝐽1/2𝜍^𝐽2≤𝑎^01/2𝜃⁰²+𝐶𝑑₁(Δ𝑡)^2𝑡𝐽_𝑡0𝜕²𝑢𝜕𝜏²²𝑑𝑠+^𝑡𝐽_𝑡0𝜕𝜂𝜕𝜏²𝑅𝑑𝑠+𝐶₁𝜂²_𝐿∞𝐿2+𝐶𝑎₁Δ𝑡^𝐽−1_𝑛=0𝜃^𝑛2+𝐶𝑅₁Δ𝑡^𝐽−1_𝑛=0𝜍^𝑛2−Δ𝑡^𝐽_𝑛=1𝑑𝜍^𝑛−1−𝜍̂^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1.Δ𝑡 (37) Using the same analysis as that in [29], we obtain Δ𝑡^𝐽_𝑛=1𝑑𝜍^𝑛−1−𝜍̂^𝑛−1,𝜍Δ𝑡^𝑛−𝜍^𝑛−1Δ𝑡=Δ𝑡^𝐽_𝑛=1𝑑𝜍^𝑛−1−𝜍̂^𝑛−1,𝜍Δ𝑡Δ^𝑛𝑡−Δ𝑡^𝐽_𝑛=1𝑑𝜍^𝑛−1−𝜍̂^𝑛−1,𝜍Δ𝑡Δ𝑡^𝑛−1=Δ𝑡^𝐽_𝑛=1𝑑𝜍^𝑛−1−𝜍^𝑛−𝜍̂^𝑛−1+𝜍^𝑛̂,𝜍Δ𝑡Δ^𝑛𝑡+Δ𝑡^𝐽_𝑛=1𝑑𝜍^𝑛−𝜍^𝑛̂,𝜍Δ𝑡Δ^𝑛𝑡−Δ𝑡^𝐽_𝑛=1𝑑𝜍^𝑛−1−𝜍̂^𝑛−1,𝜍Δ𝑡Δ𝑡^𝑛−1=Δ𝑡^𝐽_𝑛=1𝑑𝜍^𝑛−1−𝜍^𝑛−𝜍̂^𝑛−1+𝜍^𝑛̂,𝜍Δ𝑡^𝑛+𝑑𝜍Δ𝑡^𝐽−𝜍^𝐽̂Δ𝑡,𝜍^𝐽≤Δ𝑡^𝐽_𝑛=1𝜍^𝑛−1−𝜍^𝑛−𝜍̂^𝑛−1+𝜍^𝑛̂Δ𝑡𝑑𝜍^𝑛+𝑑𝜍Δ𝑡^𝐽−𝜍^𝐽̂Δ𝑡,𝜍^𝐽≤𝜀Δ𝑡^𝐽_𝑛=1𝜍^𝑛−1−𝜍^𝑛Δ𝑡²+𝐶(𝜀)Δ𝑡^𝐽_𝑛=1𝜍^𝑛2+𝐶𝑑₁𝜍^𝐽2. (38) Substitute (31) and (38) into (37) to obtain 𝑑₀−𝜀Δ𝑡^𝐽_𝑛=1𝜍^𝑛−𝜍^𝑛−1Δ𝑡²+^𝐽_𝑛=1𝑎^𝑛1/2𝜃^𝑛−𝜃^𝑛−12+𝑅^𝑛1/2𝜍^𝑛−𝜍^𝑛−12+𝑎^𝐽1/2𝜃^𝐽2+𝑅^𝐽1/2𝜍^𝐽2≤𝐶𝑎₀,𝑑₁,𝑅₁(Δ𝑡)^2𝑡𝐽_𝑡0𝜕²𝑢𝜕𝜏²²𝑑𝑠+^𝑡𝐽_𝑡0𝜕𝜂𝜕𝜏²𝑎𝑑𝑠+𝐶₀,𝑑₁,𝑅₁,𝑅₁𝜂²_𝐿∞𝐿2+𝐶𝑎₁Δ𝑡^𝐽−1_𝑛=0𝜃^𝑛2. (39) Using Gronwall lemma, we have Δ𝑡^𝐽_𝑛=1𝜍^𝑛−𝜍^𝑛−1Δ𝑡²+𝜃^𝐽2+𝑅^𝐽1/2𝜍^𝐽2≤𝐶𝑎₀,𝑑₁,𝑅₁(Δ𝑡)^2𝑡𝐽_𝑡0𝜕²𝑢𝜕𝜏²²𝑑𝑠+^𝑡𝐽_𝑡0𝜕𝜂𝜕𝜏²𝑎𝑑𝑠+𝐶₀,𝑑₁,𝑅₁,𝑅₁𝜂²_𝐿∞𝐿2. (40)
Taking 𝐳_ℎ=𝜉^𝑛 in (21a), (21b), and (21c), we get 𝜉^𝑛2≤𝐶𝜃^𝑛2. (41) Substitute (40) into (41) to get 𝜉^𝐽2≤𝐶𝑎₀,𝑑₁,𝑅₁(Δ𝑡)^2𝑡𝐽_𝑡0𝜕²𝑢𝜕𝜏²²𝑑𝑠+^𝑡𝐽_𝑡0𝜕𝜂𝜕𝜏²𝑎𝑑𝑠+𝐶₀,𝑑₁,𝑅₁,𝑅₁𝜂²_𝐿∞𝐿2. (42)
Taking 𝐰_ℎ=∇𝜍^𝑛 in (21b), we have ∇𝜍^𝑛2≤𝜃^𝑛2≤𝐶𝑎₀,𝑑₁,𝑅₁(Δ𝑡)^2𝑡𝐽_𝑡0𝜕²𝑢𝜕𝜏²²𝑑𝑠+^𝑡𝐽_𝑡0𝜕𝜂𝜕𝜏²𝑎𝑑𝑠+𝐶₀,𝑑₁,𝑅₁,𝑅₁𝜂²_{𝐿∞(𝐿2)}. (43) Combining (14)–(16), (31), (40), (43), (42), and the triangle inequality, we complete the proof.

Remark 9. Compared to the results in [29], we can obtain the optimal a priori error estimates in 𝐻¹-norm for the scalar unknown 𝑢 in Theorem 8.

4. Numerical Experiment

In this section, in order to confirm our theoretical results for the novel characteristic expanded mixed finite element method, we consider the following test problem: 𝑑(𝐱)𝜕𝑢()(𝑎())()𝑢(),()𝜕𝑡+𝐜𝐱,𝑡⋅∇𝑢−∇⋅𝐱,𝑡∇𝑢+𝑅𝐱,𝑡=𝑓𝐱,𝑡𝐱,𝑡∈Ω×𝐽, (44) with boundary condition 𝑢(𝐱,𝑡)=0,(𝐱,𝑡)∈𝜕Ω×𝐽, (45) and initial condition 𝑢(𝐱,0)=𝑥₁𝑥₂𝑥₁−1𝑥₂−12𝑥₂−1,𝐱∈Ω, (46) where Ω=[0,1]×[0,1], 𝐱=(𝑥₁,𝑥₂), 𝐽=(0,1], 𝑎(𝐱,𝑡)=1+2𝑥²₁+𝑥²₂, 𝐜(𝐱,𝑡)=(1+𝑥²₁+𝑥²₂+𝑡²,1+𝑥²₁+𝑥²₂+𝑡²), 𝑅(𝐱,𝑡)=1+𝑥²₁+2𝑥²₂, and 𝑑(𝐱)=1 and 𝑓(𝐱,𝑡) is chosen so that the exact solution for the scalar unknown function is 𝑢(𝐱,𝑡)=𝑒^−𝑡𝑥₁𝑥₂𝑥₁−1𝑥₂−12𝑥₂−1. (47) The corresponding exact gradient function is 𝜆=∇𝑢=𝑒^−𝑡𝑥₂2𝑥₁−1𝑥₂−12𝑥₂𝑥−1₁𝑥₁−16𝑥²₂−6𝑥₂+1, (48) and its exact flux function is 𝜎=−𝑎∇𝑢=−𝑒^−𝑡1+2𝑥²₁+𝑥²₂𝑥₂2𝑥₁−1𝑥₂−12𝑥₂𝑥−1₁𝑥₁−16𝑥²₂−6𝑥₂.+1 (49)

We divide the domain Ω=[0,1]×[0,1] into the triangulations of spatial mesh parameter ℎ uniformly and use the backward Euler procedure with uniform time discretization parameter Δ𝑡. We consider the piecewise linear space 𝑉_ℎ with index 𝑘=1 and the corresponding piecewise constant space 𝐖_ℎ with index 𝑘=0.

In Table 1, we get the optimal a priori error estimate in 𝐿²-norm for the scalar unknown 𝑢 with ℎ=2√2Δ𝑡=√2/16,√2/32,√2/64. At the same time, we also obtain the optimal a priori error estimate in 𝐻¹-norm for the scalar unknown 𝑢.

Table 1 

The errors and convergence rate for 𝑢.

In Table 2, we obtain some convergence results in (𝐿²(Ω))²-norm for the gradient 𝜆 with ℎ=2√2Δ𝑡=√2/16, √2/32, √2/64 in Table 2. The similar results are obtained for the flux 𝜎 in Table 2. From the data obtained in Table 2, we can find that the numerical results confirm the theoretical results of Theorem 8.

Table 2 

The errors and convergence rate for 𝜆 and 𝜎.

Figure 1 shows the surface for the exact solution 𝑢 at 𝑡=1, and Figure 2 describes the corresponding surface for the numerical solution 𝑢_ℎ with ℎ=2√2Δ𝑡=√2/16 at 𝑡=1. In Figures 1 and 2, we can find easily that the exact solution 𝑢 is approximated very well by the numerical solution 𝑢_ℎ.

Figure 1 

The surface for the exact solution 𝑢.

Figure 2 

The surface for the numerical solution 𝑢ℎ.

Figures 3 and 4 show the surface of the exact gradient function 𝜆 and the surface of the numerical gradient function 𝜆_ℎ with ℎ=2√2Δ𝑡=√2/16 at 𝑡, respectively. A similar comparison between the surface of the exact flux function 𝜎 and the surface of the numerical flux function 𝜎_ℎ in Figures 5 and 6 also is made.

(a)

(b)

Figure 3 

The surface for the exact solution 𝜆=(𝜆1,𝜆2).

(a)

(b)

Figure 4 

The surface for the numerical solution 𝜆ℎ=(𝜆1ℎ,𝜆2ℎ).

(a)

(b)

Figure 5 

The surface for the exact solution 𝜎=(𝜎1,𝜎2).

(a)

(b)

Figure 6 

The surface for the numerical solution 𝜎ℎ=(𝜎1ℎ,𝜎2ℎ).

From the convergence results for ‖𝑢−𝑢_ℎ‖_𝐿^∞_(𝐿²_(Ω)), ‖𝑢−𝑢_ℎ‖_𝐿^∞_(𝐻¹_(Ω)), ‖𝜆−𝜆_ℎ‖_{𝐿∞((𝐿2(Ω))}²₎, and ‖𝜎−𝜎_ℎ‖_{𝐿∞((𝐿2(Ω))}²₎ in Tables 1 and 2 and Figures 1 and 2, we find that the numerical results confirm our theoretical analysis.

5. Concluding Remarks

In this paper, we propose and study a novel characteristic expanded mixed finite element method, which combines the novel expanded mixed method [24] applied to approximating the diffusion term and the characteristic method that handled the hyperbolic part, for reaction-convection-diffusion equations. Compared to Chen’s expanded mixed method, the gradient for our method belongs to the square integrable space instead of the classical 𝐇(div;Ω) space. We derive a priori error estimates based on backward Euler method. Moreover, we prove the optimal a priori error estimates in 𝐿²- and 𝐻¹-norms for the scalar unknown 𝑢 and a priori error estimates in (𝐿²)²-norm for its gradient 𝜆 and its flux 𝜎. Finally, we choose a test problem to confirm our theoretical results. In the near future, the proposed characteristic expanded mixed scheme will be applied to other linear/nonlinear evolution equations, such as nonlinear reaction-diffusion equations, linear/nonlinear convection-dominated Sobolev equations, and time-dependent convection-diffusion optimal control problems.

Acknowledgments

The authors thank the anonymous referees and editors for their helpful comments and suggestions, which greatly improve the paper. This work is supported by the National Natural Science Fund (11061021), Natural Science Fund of Inner Mongolia Autonomous Region (2012MS0108, 2012MS0106, and 2011BS0102), Scientific Research Projection of Higher Schools of Inner Mongolia (NJZZ12011, NJ10006, NJ10016, and NJZY13199), Key Project of Chinese Ministry of Education (12024), and the Program of Higher-Level Talents of Inner Mongolia University (125119, Z200901004, and 30105-125132).