A New Characteristic Nonconforming Mixed Finite Element Scheme for Convection-Dominated Diffusion Problem

Abstract

A characteristic nonconforming mixed finite element method (MFEM) is proposed for the convection-dominated diffusion problem based on a new mixed variational formulation. The optimal order error estimates for both the original variable 𝑢 and the auxiliary variable 𝜎 with respect to the space are obtained by employing some typical characters of the interpolation operator instead of the mixed (or expanded mixed) elliptic projection which is an indispensable tool in the traditional MFEM analysis. At last, we give some numerical results to confirm the theoretical analysis.

1. Introduction

Consider the following convection-dominated diffusion problem: 𝑢_𝑡𝑏,(),𝑢(),𝑢+𝐚𝑥,𝑦⋅∇𝑢−∇⋅𝑥,𝑦∇𝑢=𝑓𝑥,𝑦,𝑡𝑥,𝑦,𝑡∈Ω×0,𝑇𝑥,𝑦,𝑡=0,𝑥,𝑦,𝑡∈𝜕Ω×0,𝑇𝑥,𝑦,0=𝑢₀𝑥,𝑦,𝑥,𝑦∈Ω, (1) where Ω is a bounded polygonal domain in ℝ² with Lipschitz continuous boundary 𝜕Ω, 𝐽=(0,𝑇],0<𝑇<+∞. ∇ and ∇⋅ denote the gradient and the divergence operators, respectively.

Model (1) has been widely used to describe the conduction of heat in fluid, the diffusion of soluble minerals or pollutants in ground water, the incompressible miscible displacement in porous media, and so on. The parameters appearing in (1) satisfy the following assumptions [1, 2]:(A1)𝑢 denotes, for example, the concentration or saturation of soluble substances;(A2)𝐚(𝑥,𝑦)=(𝑎₁(𝑥,𝑦),𝑎₂(𝑥,𝑦)) represents Darcy velocity of mixed fluid, and 𝑓 a source term;(A3)𝑏(𝑥,𝑦) is sufficiently smooth and there exist constants 𝑏₁ and 𝑏₂, such that 0<𝑏₁≤𝑏𝑥,𝑦≤𝑏₂<+∞,∀𝑥,𝑦∈Ω. (2)

It is well known that convection dominated-diffusion problem (1) often presents serious numerical difficulties. The standard numerical methods, such as finite difference method (FDM), FEM and MFEM, usually produce numerical diffusion along sharp fronts. In order to overcome this fatal defect, Douglas et al. [3] combined the method of characteristics with FE procedures so as to reduce the truncation error, and it allows us to use large time steps without lose of accuracy. Moreover, there have appeared many effective discretization schemes concentrating on the hyperbolic nature of the equation, for example, characteristic FD streamline diffusion method [4, 5], Eulerian-Lagrangian method [6, 7], characteristic-finite volume element method [2, 8, 9], characteristics-mixed covolume method [10, 11], the modified method of characteristic-Galerkin FE procedure [12], characteristic nonconforming FEM [13–15], characteristic MFEM [16–19] and expanded characteristic MFEM [1, 20], and so forth.

As for the characteristic MFEM or expanded characteristic MFEM, the convergence rates of 𝑢 and 𝜎 in existing literature were suboptimal [11, 18, 21, 22] and the convergence analysis was valid only to the case of the lowest order MFE approximation [10, 17]. So far, to our best knowledge there are few studies on the optimal order error estimates except for [23], in which a family of characteristic MFEM with arbitrary degree of Raviart-Thomas-Nédélec space in [24, 25] for transient convection diffusion equations was studied.

Recently, based on the low regularity requirement of the flux variable in practical problems, a new mixed variational form for second elliptic problem was proposed in [26]. It has two typical advantages: the flux space belongs to the square integrable space instead of the traditional 𝐻(div;Ω), which makes the choices of MFE spaces sufficiently simple and easy; the LBB condition is automatically satisfied when the gradient of approximation space for the original variable is included in approximation space for the flux variable. Motivated by this idea, this paper will construct a characteristic nonconforming MFE scheme for (1) with a new mixed variational formulation. Similar to the expanded characteristic MFEM, the coefficient 𝑏 of (1) in this proposed scheme does not need to be inverted; therefore, it is also suitable for the case when 𝑏 is small. By employing some distinct characters of the interpolation operators on the element instead of the mixed or expanded mixed elliptic projection used in [1, 17, 20] which is an indispensable tool in the traditional characteristic MFEM analysis, the 𝑂(ℎ²) order error estimate in 𝐿²-norm for original variable 𝑢, which is one order higher than [1, 20] and half order higher than [18], is derived, and the optimal error estimates with order 𝑂(ℎ) for auxiliary variable 𝜎 in 𝐿²-norm and for 𝑢 in broken 𝐻¹-norm are obtained, respectively. It seems that the result for 𝑢 in broken 𝐻¹-norm has never been seen in the existing literature by making full use of the high-accuracy estimates of the lowest order Raviart-Thomas element proved by the technique of integral identities in [27] and the special properties of nonconforming 𝐸𝑄₁^rot element (see Lemma 1 below).

The paper is organized as follows. Section 2 is devoted to the introduction of the nonconforming FE approximation spaces and their corresponding interpolation operators. In Section 3, we will give the construction of the new characteristic nonconforming MFE scheme and two important lemmas, and the existence and uniqueness of the discrete scheme solution will be proved. In Section 4, the convergence analysis and optimal error estimates for both the original variable 𝑢 and the flux variable 𝜎 are obtained. In Section 5, some numerical results are provided to illustrate the effectiveness of our proposed method.

Throughout this paper, 𝐶 denotes a generic positive constant independent of the mesh parameters ℎ and Δ𝑡 with respect to domain Ω and time 𝑡.

2. Construction of Nonconforming MFEs

As in [28], we frequently employ the space 𝐿²(Ω) of square integrable functions with scalar product and norm (𝑢,𝑣)=(𝑢,𝑣)_𝐿²_(Ω)=_Ω𝑢𝑣𝑑𝑥𝑑𝑦^1/2,‖𝑣‖=‖𝑣‖_𝐿²_(Ω)=_Ω𝑣²𝑑𝑥𝑑𝑦^1/2. (3) We also employ the Sobolev space 𝐻^𝑚(Ω),𝑚≥1, of functions 𝑣 such that 𝐷^𝛽𝑣∈𝐿²(Ω) for all |𝛽|≤𝑚, equipped with the norm and seminorm ‖𝑣‖_𝑚,Ω=‖𝑣‖_𝐻^𝑚_(Ω)=_𝛽≤𝑚𝐷^𝛽𝑣^21/2,|𝑣|_𝑚,Ω=|𝑣|_𝐻^𝑚_(Ω)=_𝛽=𝑚𝐷^𝛽𝑣^21/2. (4) The space 𝐻¹₀(Ω) denotes the closure of the set of infinitely differentiable functions with compact supports in Ω. For any Sobolev space 𝑌, 𝐿^𝑝(0,𝑇;𝑌) is the space of measurable 𝑌-valued functions Φ of 𝑡∈(0,𝑇), such that ∫^𝑇₀‖Φ(⋅,𝑡)‖^𝑝_𝑌𝑑𝑡<∞ if 1≤𝑝<∞, or such that esssup_0<𝑡<𝑇‖Φ(⋅,𝑡)‖_𝑌<∞ if 𝑝=∞.

We now introduce the nonconforming MFE space described in [29] for and summarize it as follows.

Let Ω⊂ℝ² be a polygon domain with edges parallel to the coordinate axes on 𝑥𝑦 plane, and let 𝑇_ℎ be a rectangular subdivision of Ω satisfying the regular condition [30]. For a given element 𝑒∈𝑇_ℎ, denote the barycenter of element 𝑒 by (𝑥_𝑒,𝑦_𝑒), denote the length of edges parallel to 𝑥-axis and 𝑦-axis by 2ℎ_𝑥𝑒 and 2ℎ_𝑦𝑒, respectively, ℎ_𝑒=max_{𝑒∈𝑇ℎ}{ℎ_𝑥𝑒,ℎ_𝑦𝑒},ℎ=max_{𝑒∈𝑇ℎ}ℎ_𝑒.

Let 𝑒̂=[−1,1]×[−1,1] be the reference element on 𝑥𝑦 plane and four vertices 𝑑₁=(−1,−1), 𝑑₂=(1,−1),  𝑑₃=(1,1), and 𝑑₄=(−1,1), the four edges 𝑙₁̂=𝑑₁𝑑₂, 𝑙₂̂=𝑑₂𝑑₃, 𝑙₃̂=𝑑₃𝑑₄, and 𝑙₄̂=𝑑₄𝑑₁. Then there exists an affine mapping 𝐹_𝑒∶𝑒̂→𝑒 as 𝑥=𝑥_𝑒+ℎ_𝑥𝑒𝑥,𝑦=𝑦_𝑒+ℎ_𝑦𝑒𝑦. (5) Define the FE spaces (𝑒,̂𝑃^𝑖,∑^𝑖),  (𝑖=1,2,3) by ∑¹=𝑣₁,𝑣₂,𝑣₃,𝑣₄,𝑣₅,𝑃¹(𝑥)𝑦,∑=span1,𝑥,𝑦,𝜙,𝜙₂=𝑝₁,𝑝₂,𝑃²{𝑥},∑=span1,₃=𝑞₁,𝑞₂,𝑃³=span1,𝑦, (6) where 𝑣_𝑖=(1/|𝑙_𝑖̂̂𝑙|)∫_𝑖𝑣𝑑𝑠,̂(𝑖=1,2,3,4), 𝑣₅=(1/|𝑒|̂)∫𝑒̂𝑣𝑑𝑥𝑑𝑦, 𝜙(𝑡)=(1/2)(3𝑡²−1), 𝑝_𝑖=(1/|𝑙̂_2𝑖|) ∫̂𝑙_2𝑖𝑝𝑑𝑠̂, 𝑞_𝑖=(1/|𝑙̂_2𝑖−1̂𝑙|)∫_2𝑖−1𝑞𝑑𝑠,̂(𝑖=1,2).

The interpolation operators on 𝑒̂ are defined as follows: Π¹∶𝑣∈𝐻¹(𝑒)̂⟶Π¹𝑣∈𝑃¹,𝑙_𝑖̂Π¹𝑣𝑑̂(),𝑒̂Π𝑣−𝑠=0,𝑖=1,2,3,4₁𝑣𝑑Π𝑣−𝑥𝑑𝑦=0,₂∶𝑝∈𝐻¹(𝑒)̂⟶Π²𝑝∈𝑃²,𝑙̂_2𝑖Π²𝑝𝑑̂(),Π𝑝−𝑠=0,𝑖=1,2₃∶𝑞∈𝐻¹(𝑒)̂⟶Π³𝑞∈𝑃³,𝑙̂_2𝑖−1Π³𝑞−𝑞𝑑𝑠̂=0,(𝑖=1,2). (7) Then the associated nonconforming 𝐸𝑄₁^rot element space 𝑀_ℎ [29] and lowest order Raviart-Thomas element space 𝐕_ℎ [25, 27] are defined as 𝑀_ℎ=𝑣_ℎ∶𝑣_ℎ∣_𝑒=𝑣∘𝐹_𝑒⁻¹,𝑣∈𝑃¹,_𝐹𝑣_ℎ,𝐕𝑑𝑠=0,𝐹⊂𝜕𝑒_ℎ=𝐰_ℎ=𝑤_ℎ1,𝑤_ℎ2∶𝐰_ℎ∣_𝑒=𝑤₁∘𝐹_𝑒⁻¹,𝑤₂∘𝐹_𝑒⁻¹,𝐰=𝑤₁,𝑤₂∈𝑃²×𝑃³, (8) respectively, where [𝜑] represents the jump value of 𝜑 across the boundary 𝐹, and [𝜑]=𝜑 if 𝐹⊂𝜕Ω.

Similarly, the interpolation operators 𝜋¹_ℎ and 𝜋²_ℎ are defined as 𝜋¹_ℎ∶𝐻¹(Ω)⟶𝑀_ℎ,𝜋¹_ℎ𝑒=𝜋¹_𝑒,𝜋¹_𝑒𝑣=Π¹𝑣∘𝐹_𝑒⁻¹,∀𝑣∈𝐻¹(Ω),𝜋²_ℎ∶𝐻¹(Ω)²⟶𝐕_ℎ,𝜋²_ℎ|_𝑒=𝜋²_𝑒,𝜋²_𝑒𝐰=Π²𝑤₁∘𝐹_𝑒⁻¹,Π³𝑤₂∘𝐹_𝑒⁻¹,∀𝐰=𝑤₁,𝑤₂∈𝐻¹(Ω)². (9)

3. New Characteristic Nonconforming MFE Scheme and Two Lemmas

Let 𝜓(𝑥,𝑦)=(1+|𝐚(𝑥,𝑦)|²)^1/2 and 𝜏=𝜏(𝑥,𝑦) be the characteristic direction associated with 𝑢_𝑡+ 𝐚(𝑥,𝑦)⋅∇𝑢, such that 𝜕=1𝜕𝜏𝜓𝜕𝑥,𝑦+𝐚𝜕𝑡𝑥,𝑦𝜓𝑥,𝑦⋅∇. (10)

Then (1) can be put in the following system: 𝜓𝑥,𝑦𝜕𝑢𝑏,∀(],𝑢(],𝑢𝜕𝜏−∇⋅𝑥,𝑦∇𝑢=𝑓𝑥,𝑦,𝑡𝑥,𝑦,𝑡∈Ω×0,𝑇𝑥,𝑦,𝑡=0,𝑥,𝑦,𝑡∈𝜕Ω×0,𝑇𝑥,𝑦,0=𝑢₀𝑥,𝑦,𝑥,𝑦∈Ω. (11)

By introducing 𝜎=−𝑏(𝑥,𝑦)∇𝑢 and using Green's formula, we obtain the new characteristic mixed form of (11). Find (𝑢,𝜎)∶(0,𝑇]→𝐻¹₀(Ω)×(𝐿²(Ω))², such that 𝜓𝑥,𝑦𝜕𝑢−()=𝑓𝜕𝜏,𝑣𝜎,∇𝑣𝑥,𝑦,𝑡,𝑣∀𝑣∈𝐻¹₀(Ω),(𝜎,𝐰)+𝑏𝑥,𝑦∇𝑢,𝐰=0,∀𝐰∈𝐿²(Ω)². (12)

Let Δ𝑡>0,𝑁=𝑇/Δ𝑡∈ℤ, 𝑡^𝑛=𝑛Δ𝑡, and 𝜙^𝑛=𝜙(𝑥,𝑦,𝑡^𝑛). When solving 𝑢_ℎ^𝑛+1, we would like to make the scheme as implicit as possible by using of the characteristic vector 𝜏. Denote 𝑋=(𝑥,𝑦)∈Ω and 𝑋=𝑋−𝐚𝑥,𝑦Δ𝑡, (13) similar to [1, 3], and then we have the following approximation: 𝜓𝑥,𝑦𝜕𝑢𝜕𝜏_𝑡^𝑛𝑢≈𝜓𝑥,𝑦𝑋,𝑡^𝑛−𝑢𝑋,𝑡^𝑛−1𝑋−𝑋²+(Δ𝑡)²=𝑢𝑋,𝑡^𝑛−𝑢𝑋,𝑡^𝑛−1=𝑢Δ𝑡^𝑛−𝑢^𝑛−1.Δ𝑡 (14)

This leads to the following characteristic nonconforming MFE scheme. Find (𝑢_ℎ,𝜎_ℎ)∶{𝑡⁰,𝑡¹,…, 𝑡^𝑁}→𝑀_ℎ×𝐕_ℎ, such that 𝑢^𝑛_ℎ−𝑢_ℎ^𝑛−1Δ𝑡,𝑣_ℎ−𝜎^𝑛_ℎ,∇𝑣_ℎℎ=𝑓^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑀_ℎ,𝜎^𝑛_ℎ,𝐰_ℎ+𝑏∇𝑢^𝑛_ℎ,𝐰_ℎℎ=0,∀𝐰_ℎ∈𝐕_ℎ,𝑢⁰_ℎ=𝜋¹_ℎ𝑢₀𝑥,𝑦,𝜎⁰_ℎ=𝜋²_ℎ𝑏∇𝑢₀𝑥,𝑦,∀𝑥,𝑦∈Ω, (15a)(15b)(15c)where 𝑢^𝑛_ℎ=𝑢_ℎ(𝑋,𝑡^𝑛), (𝑢,𝑣)_ℎ=∑_{𝑒∈𝑇ℎ}∫_𝑒𝑢𝑣𝑑𝑥𝑑𝑦. Generally speaking, 𝑢_ℎ^𝑛−1(𝑛=2,…,𝑁) are not node values and should be derived by interpolation formulas on 𝑢_ℎ^𝑛−1.

Remark 1. In [1], the expanded characteristic MFE scheme was presented by introducing two new auxiliary variables which avoided the inversion of the coefficient 𝑏 when 𝑏 is small. The new mixed schemes (15a), (15b), and (15c) not only keep the advantage of expanded characteristic MFE scheme, but also donot need to solve three variables.

Now, we prove the existence and uniqueness of the solution of (15a), (15b), and (15c).

Theorem 1. Under assumption (A3), there exists a unique solution (𝑢_ℎ,𝜎_ℎ)∈𝑀_ℎ×𝐕_ℎ to the schemes (15a), (15b), and (15c).

Proof. The linear system generated by (15a), (15b), and (15c) is square, so the existence of the solution is implied by its uniqueness. From (15a), (15b), and (15c), we have 𝑢^𝑛_ℎΔ𝑡,𝑣_ℎ−𝜎^𝑛_ℎ,∇𝑣_ℎℎ=𝑢_ℎ^𝑛−1Δ𝑡,𝑣_ℎ+𝑓^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑀_ℎ,𝜎^𝑛_ℎ,𝐰_ℎ+𝑏∇𝑢^𝑛_ℎ,𝐰_ℎℎ=0,∀𝐰_ℎ∈𝐕_ℎ. (16)
Let 𝑢^𝑛_ℎ and 𝑓 be zero, and thus 𝑢^𝑛_ℎ is zero too; taking 𝑣_ℎ=𝑢^𝑛_ℎ,𝐰_ℎ=(1/𝑏)𝜎^𝑛_ℎ in (16) and adding them together, we have 1𝑢Δ𝑡^𝑛ℎ2+1𝑏𝜎^𝑛_ℎ,𝜎^𝑛_ℎ=0. (17) Thus assumption (A3) implies that 𝑢^𝑛_ℎ=𝜎^𝑛_ℎ=0. The proof is complete.

To get error estimates, we state the following two important lemmas.

Lemma 1 (see [27, 29, 31]). Assume that 𝑢∈𝐻¹(Ω),𝐩∈(𝐻²(Ω))², for all 𝑣_ℎ∈𝑀_ℎ,𝐰_ℎ∈𝐕_ℎ, and then there hold ∇𝑢−𝜋¹_ℎ𝑢,∇𝑣_ℎℎ=0,∇𝑢−𝜋¹_ℎ𝑢,𝐰_ℎℎ=0,𝐩−𝜋²_ℎ𝐩,𝐰_ℎ≤𝐶ℎ²|𝐩|_2,Ω𝐰_ℎ,_{𝑒∈𝑇ℎ𝜕𝑒}𝐩𝑣_ℎ⋅𝐧𝑑𝑠≤𝐶ℎ²|𝐩|_2,Ω𝑣_ℎ1,ℎ, (18)(19)(20) where ‖⋅‖_1,ℎ=(∑_{𝑒∈𝑇ℎ}|⋅|_1,𝑒)^1/2 is a norm on 𝑀_ℎ, and 𝐧 denotes the outward unit normal vector on 𝜕𝑒.

Lemma 2 (see [1, 3]). Let 𝜑∈𝐿²(Ω), and 𝜑=𝜑(𝑋−𝑔(𝑋)Δ𝑡), where function 𝑔 and its gradient ∇𝑔 are bounded, then 𝜑−𝜑₋₁≤𝐶𝜑Δ𝑡, (21) where ‖𝜑‖₋₁=sup_{𝜙∈𝐻1(Ω)}((𝜑,𝜙)/‖𝜙‖_1,Ω).

4. Convergence Analysis and Optimal Order Error Estimates

In this section, we aim to analyze the convergence analysis and error estimates of characteristic nonconforming MFEM. In order to do this, let 𝑢_ℎ−𝑢=𝑢_ℎ−𝜋¹_ℎ𝑢+𝜋¹_ℎ𝜎𝑢−𝑢=𝑒+𝜌,_ℎ−𝜎=𝜎_ℎ−𝜋²_ℎ𝜎+𝜋²_ℎ𝜎−𝜎=𝜉+𝜂. (22) Taking 𝑡=𝑡^𝑛 in (12) yields𝜓𝜕𝑢^𝑛𝜕𝜏,𝑣_ℎ−𝜎^𝑛,∇𝑣_ℎℎ+_{𝑒∈𝑇ℎ𝜕𝑒}𝜎^𝑛𝑣_ℎ⋅𝐧𝑑𝑠=𝑓^𝑛,𝑣_ℎ,∀𝑣_ℎ∈𝑀_ℎ,𝜎^𝑛,𝐰_ℎ+𝑏∇𝑢^𝑛,𝐰_ℎℎ=0,∀𝐰_ℎ∈𝐕_ℎ. (23a)(23b)From (23a), (23b), (15a), (15b), and (15c) we get𝑒^𝑛−𝑒^𝑛−1Δ𝑡,𝑣_ℎ−𝜉^𝑛,∇𝑣_ℎℎ=𝜓𝜕𝑢^𝑛−𝑢𝜕𝜏^𝑛−𝑢^𝑛−1Δ𝑡,𝑣_ℎ−𝜌^𝑛−𝜌^𝑛−1Δ𝑡,𝑣_ℎ+𝜂^𝑛,∇𝑣_ℎℎ+_{𝑒∈𝑇ℎ𝜕𝑒}𝜎^𝑛𝑣_ℎ⋅𝐧𝑑𝑠,∀𝑣_ℎ∈𝑀_ℎ,𝜉^𝑛,𝐰_ℎ+𝑏∇𝑒^𝑛,𝐰_ℎℎ=−𝜂^𝑛,𝐰_ℎ−𝑏∇𝜌^𝑛,𝐰_ℎℎ,∀𝐰_ℎ∈𝐕_ℎ. (24a)(24b)We are now in a position to prove the optimal order error estimates.

Theorem 2. Let (𝑢,𝜎) and (𝑢^𝑛_ℎ,𝜎^𝑛_ℎ) be the solutions of (12), (15a), (15b), and (15c), respectively, (𝜕²𝑢/𝜕𝜏²)∈𝐿²(0,𝑇;𝐿²(Ω)),𝑢_𝑡∈𝐿²(0,𝑇;𝐻²(Ω)),𝑢∈𝐿^∞(0,𝑇;𝐻²(Ω)),𝜎∈𝐿^∞(0,𝑇;𝐻²(Ω)) and assume that Δ𝑡=𝑂(ℎ²). Then under assumption (A3), we have max_{0≤𝑛≤𝑁}𝑢_ℎ−𝑢𝑡^𝑛_1,ℎ≤𝐶(Δ𝑡+ℎ),max_{0≤𝑛≤𝑁}𝑢_ℎ−𝑢𝑡^𝑛≤𝐶Δ𝑡+ℎ²,max_{0≤𝑛≤𝑁}𝜎_ℎ−𝜎𝑡^𝑛≤𝐶(Δ𝑡+ℎ). (25)(26)(27)

Proof. Taking 𝑣_ℎ=𝑒^𝑛 in (24a) and 𝐰_ℎ=∇𝑒^𝑛 in (24b), and adding them, we have 𝑒^𝑛−𝑒^𝑛−1Δ𝑡,𝑒^𝑛+𝑏∇𝑒^𝑛,∇𝑒^𝑛_ℎ=𝜓𝜕𝑢^𝑛−𝑢𝜕𝜏^𝑛−𝑢^𝑛−1Δ𝑡,𝑒^𝑛−𝜌^𝑛−𝜌^𝑛−1Δ𝑡,𝑒^𝑛−𝜌^𝑛−1−𝜌^𝑛−1Δ𝑡,𝑒^𝑛+_{𝑒∈𝑇ℎ𝜕𝑒}𝜎^𝑛𝑒^𝑛⋅𝐧𝑑𝑠−𝑏∇𝜌^𝑛,∇𝑒^𝑛_ℎ=⁵_𝑖=1(Err)_𝑖. (28) On the one hand, we consider the right hand of (28).
Using the method similar to [3], we have (Err)₁≤𝐶𝜓𝜕𝑢^𝑛−𝑢𝜕𝜏^𝑛−𝑢^𝑛−1Δ𝑡²+𝜀₁2𝑒^𝑛2𝜕≤𝐶Δ𝑡²𝑢𝜕𝜏²²_{𝐿2𝑡𝑛−1,𝑡𝑛;𝐿2(Ω)}+𝜀₁2𝑒^𝑛2. (29)(Err)₂ can be estimated as (Err)₂≤1Δ𝑡_Ω^𝑡𝑛_𝑡^𝑛−1𝜌_𝑡𝑑𝑠²𝑑𝑥𝑑𝑦^1/2𝑒^𝑛≤1√Δ𝑡^{Ω𝑡𝑛𝑡𝑛−1}𝜌²_𝑡𝑑𝑠𝑑𝑥𝑑𝑦^1/2𝑒^𝑛≤𝐶Δ𝑡^𝑡𝑛_𝑡^𝑛−1𝜌_𝑡²𝜀𝑑𝑠+¹2𝑒^𝑛2≤𝐶ℎ⁴Δ𝑡^𝑡𝑛_𝑡^𝑛−1𝑢_𝑡²_2,Ω𝜀𝑑𝑠+¹2𝑒^𝑛2. (30) By Lemma 2, we obtain (Err)₃≤1𝜌Δ𝑡^𝑛−1−𝜌^𝑛−1₋₁𝑒^𝑛_1,ℎ≤𝐶𝜌^𝑛−12+𝑏₁6𝑒^𝑛2_1,ℎ≤𝐶ℎ⁴𝑢^𝑛−12_2,Ω+𝑏₁6𝑒^𝑛2_1,ℎ. (31) It follows from Lemma 1 that (Err)₄≤𝐶ℎ⁴𝜎^𝑛2_2,Ω+𝑏₁6𝑒^𝑛2_1,ℎ. (32) Let 𝑏=(1/|𝑒|)∫_𝑒𝑏(𝑥,𝑦)𝑑𝑥𝑑𝑦. By Lemma 1, we have (Err)₅=−𝑏−𝑏∇𝜌^𝑛,∇𝑒^𝑛_ℎ≤𝐶ℎ|𝑏|_𝑊^1,∞_(Ω)𝜌^𝑛_1,ℎ𝑒^𝑛_1,ℎ≤𝐶ℎ⁴𝑢^𝑛2_2,Ω+𝑏₁6𝑒^𝑛2_1,ℎ. (33) On the other hand, the left hand of (28) can be bounded by 𝑒^𝑛−𝑒^𝑛−1Δ𝑡,𝑒^𝑛+𝑏∇𝑒^𝑛,∇𝑒^𝑛_ℎ≥1𝑒2Δ𝑡^𝑛,𝑒^𝑛−𝑒^𝑛−1,𝑒^𝑛−1+𝑏₁𝑒^𝑛2_1,ℎ≥1𝑒2Δ𝑡^𝑛2−(1+𝐶Δ𝑡)𝑒^𝑛−12+𝑏₁𝑒^𝑛2_1,ℎ, (34) where the inequality ‖𝑒^𝑛−1‖²≤(1+𝐶Δ𝑡)‖𝑒^𝑛−1‖² proved in [3] is used in the last step.
Combining (29)–(34) with (28) gives 1𝑒2Δ𝑡^𝑛2−𝑒^𝑛−12+𝑏₁𝑒^𝑛2_1,ℎ𝜕≤𝐶Δ𝑡²𝑢𝜕𝜏²²_{𝐿2(𝑡𝑛−1,𝑡𝑛;𝐿2(Ω))}+ℎ⁴Δ𝑡^𝑡𝑛_𝑡^𝑛−1𝑢_𝑡²_2,Ω𝑑𝑠+ℎ⁴𝑢^𝑛−12_2,Ω+𝑢^𝑛2_2,Ω+𝜎^𝑛2_2,Ω+𝜀₁𝑒^𝑛2+𝐶𝑒^𝑛−12+𝑏₁2𝑒^𝑛_1,ℎ. (35) Taking 1−2Δ𝑡𝜀₁>0, multiplying (35) by 2Δ𝑡, summing over from 𝑖=1 to 𝑖=𝑛, and noticing that 𝑒⁰=0, we obtain 𝑒^𝑛2+Δ𝑡^𝑛_𝑖=1𝑒^𝑖2_1,ℎ≤𝐶(Δ𝑡)²𝜕²𝑢𝜕𝜏²²_{𝐿20,𝑡𝑛;𝐿2(Ω)}+ℎ^4𝑡𝑛₀𝑢_𝑡²_2,Ω𝑑𝑠+Δ𝑡ℎ^4𝑛_𝑖=1𝑢^𝑖2_2,Ω+𝜎^𝑖2_2,Ω+𝐶^𝑛−1_𝑖=1𝑒^𝑖2. (36) It follows from discrete Gronwall’s lemma that 𝑒^𝑛2+Δ𝑡^𝑛_𝑖=1𝑒^𝑖2_1,ℎ≤𝐶(Δ𝑡)²𝜕²𝑢𝜕𝜏²²_{𝐿20,𝑡𝑛;𝐿2(Ω)}+ℎ⁴𝑢_𝑡²_{𝐿20,𝑡𝑛;𝐻2(Ω)}+‖𝑢‖²_{𝐿∞0,𝑡𝑛;𝐻2(Ω)}+‖𝜎‖²_{𝐿∞0,𝑡𝑛;𝐻2(Ω)2}. (37) From (37) we get the optimal order error estimate of ‖𝑒^𝑛‖ rather than ‖𝑒^𝑛‖_1,ℎ. So we start to reestimate ‖𝑒^𝑛‖_1,ℎ in the following manner and derive the estimation of ‖𝜉^𝑛‖ simultaneously.
Firstly, choosing 𝑣_ℎ=((𝑒^𝑛−𝑒^𝑛−1)/Δ𝑡) in (24a) and 𝐰_ℎ=∇((𝑒^𝑛−𝑒^𝑛−1)/Δ𝑡) in (24b), and adding them, we have 𝑒^𝑛−𝑒^𝑛−1,𝑒Δ𝑡^𝑛−𝑒^𝑛−1+Δ𝑡𝑏∇𝑒^𝑛𝑒,∇^𝑛−𝑒^𝑛−1Δ𝑡_ℎ=𝜓𝜕𝑢^𝑛−𝑢𝜕𝜏^𝑛−𝑢^𝑛−1,𝑒Δ𝑡^𝑛−𝑒^𝑛−1−𝜌Δ𝑡_𝑛−𝜌^𝑛−1,𝑒Δ𝑡^𝑛−𝑒^𝑛−1−𝜌Δ𝑡_𝑛−1−𝜌^𝑛−1,𝑒Δ𝑡^𝑛−𝑒^𝑛−1Δ𝑡+_{𝑒∈𝑇ℎ𝜕𝑒}𝜎^𝑛𝑒^𝑛−𝑒^𝑛−1Δ𝑡⋅𝐧𝑑𝑠−𝑏∇𝜌^𝑛𝑒,∇^𝑛−𝑒^𝑛−1Δ𝑡_ℎ=⁵_𝑖=1(Err)^′_𝑖. (38) The left hand can be estimated as 𝑒^𝑛−𝑒^𝑛−1,𝑒Δ𝑡^𝑛−𝑒^𝑛−1+Δ𝑡𝑏∇𝑒^𝑛𝑒,∇^𝑛−𝑒^𝑛−1Δ𝑡_ℎ≥𝑒^𝑛−𝑒^𝑛−1Δ𝑡²+12Δ𝑡𝑏∇𝑒^𝑛,∇𝑒^𝑛−𝑏∇𝑒^𝑛−1,∇𝑒^𝑛−1+𝑒^𝑛−1−𝑒^𝑛−1,𝑒Δ𝑡^𝑛−𝑒^𝑛−1,Δ𝑡 (39) and (Err)^′_𝑖,(𝑖=1,2,3,4,5) can be bounded by (Err)^′₁𝜕≤𝐶Δ𝑡²𝑢𝜕𝜏²²_{𝐿2𝑡𝑛−1,𝑡𝑛;𝐿2(Ω)}+14𝑒^𝑛−𝑒^𝑛−1Δ𝑡²,(Err)^′₂≤𝐶ℎ⁴Δ𝑡^𝑡𝑛_𝑡^𝑛−1𝑢_𝑡²_2,Ω1𝑑𝑠+4𝑒^𝑛−𝑒^𝑛−1Δ𝑡²,(Err)^′₃≤𝐶ℎ⁴𝑢Δ𝑡^{𝑛−122,Ω}+𝜀3𝑒Δ𝑡^𝑛−𝑒^𝑛−1Δ𝑡²_1,ℎ,(Err)^′₄≤𝐶ℎ⁴𝜎Δ𝑡^𝑛22,Ω+𝜀3𝑒Δ𝑡^𝑛−𝑒^𝑛−1Δ𝑡²_1,ℎ,(Err)^′₅≤𝐶ℎ⁴𝑢Δ𝑡^𝑛22,Ω+𝜀3𝑒Δ𝑡^𝑛−𝑒^𝑛−1Δ𝑡²_1,ℎ. (40) From (38)–(40), we get 12𝑒^𝑛−𝑒^𝑛−1Δ𝑡²+12Δ𝑡𝑏∇𝑒^𝑛,∇𝑒^𝑛_ℎ−𝑏∇𝑒^𝑛−1,∇𝑒^𝑛−1_ℎ𝜕≤𝐶Δ𝑡²𝑢𝜕𝜏²²_{𝐿2𝑡𝑛−1,𝑡𝑛;𝐿2(Ω)}+ℎ⁴Δ𝑡^𝑡𝑛_𝑡^𝑛−1𝑢_𝑡²_2,Ω𝑑𝑠+𝑢^𝑛2_2,Ω+𝑢^𝑛−12_2,Ω+𝜎^𝑛2_2,Ω𝑒+𝜀Δ𝑡^𝑛−𝑒^𝑛−1Δ𝑡²_1,ℎ+𝑒^𝑛−1−𝑒^𝑛−1,𝑒Δ𝑡^𝑛−𝑒^𝑛−1.Δ𝑡 (41) Multiplying (41) by 2Δ𝑡 and summing over in time from 𝑖=1 to 𝑖=𝑛 yield 𝑒Δ𝑡^𝑛−𝑒^𝑛−1Δ𝑡²+𝑏₁𝑒^𝑛2_1,ℎ≤𝐶(Δ𝑡)²𝜕²𝑢𝜕𝜏²²_{𝐿20,𝑡𝑛;𝐿2(Ω)}+ℎ⁴𝑢_𝑡²_{𝐿2(0,𝑡𝑛;𝐻2(Ω))}+ℎ^4𝑛_𝑖=1𝑢^𝑖2_2,Ω+𝜎^𝑖2_2,Ω+𝜀(Δ𝑡)^2𝑛_𝑖=1𝑒^𝑖−𝑒^𝑖−1Δ𝑡²_1,ℎ+^𝑛_𝑖=1𝑒^𝑖−1−𝑒^𝑖−1Δ𝑡,𝑒^𝑖−𝑒^𝑖−1. (42) Secondly, we take Δ𝑡→0 and Δ𝑡 must approach zero in such a way that Δ𝑡 and ℎ satisfy Δ𝑡=𝑂ℎ², (43) and by inverse inequality, we have (Δ𝑡)^2𝑛_𝑖=1𝑒^𝑖−𝑒^𝑖−1Δ𝑡²_1,ℎ≤𝐶Δ𝑡^𝑛_𝑖=1𝑒^𝑖−𝑒^𝑖−1Δ𝑡². (44) At the same time, using Lemma 2, we obtain ^𝑛_𝑖=1𝑒^𝑖−1−𝑒^𝑖−1Δ𝑡,𝑒^𝑖−𝑒^𝑖−1=𝑒^𝑛−1−𝑒^𝑛−1Δ𝑡,𝑒^𝑛+^𝑛−1_𝑖=1𝑒^𝑖−1−𝑒^𝑖−𝑒^𝑖−1−𝑒^𝑖Δ𝑡,𝑒^𝑖≤𝐶𝑒^𝑛−1𝑒^𝑛_1,ℎ+^𝑛−1_𝑖=1𝑒^𝑖−𝑒^𝑖−1𝑒^𝑖_1,ℎ≤𝐶𝑒^𝑛−12+𝑏₁2𝑒^𝑛2_1,ℎ+Δ𝑡^𝑛−1_𝑖=1𝑒^𝑖−𝑒^𝑖−1Δ𝑡²+𝐶Δ𝑡^𝑛−1_𝑖=1𝑒^𝑖2_1,ℎ. (45) From (42)–(45), taking suitable small 𝜀 such that 1−𝜀𝐶>0, we have 𝑒Δ𝑡^𝑛−𝑒^𝑛−1Δ𝑡²+𝑒^𝑛2_1,ℎ≤𝐶(Δ𝑡)²𝜕²𝑢𝜕𝜏²²_{𝐿20,𝑡𝑛;𝐿2(Ω)}+ℎ⁴𝑢_𝑡²_{𝐿20,𝑡𝑛;𝐻2(Ω)}+ℎ^4𝑛_𝑖=1𝑢^𝑖2_2,Ω+𝜎^𝑖2_2,Ω+𝑒^𝑛−12+𝐶Δ𝑡^𝑛−1_𝑖=1𝑒^𝑖−𝑒^𝑖−1Δ𝑡²+𝐶Δ𝑡^𝑛−1_𝑖=1𝑒^𝑖2_1,ℎ. (46) Finally, applying discrete Gronwall’s lemma yields 𝑒^𝑛2_1,ℎ≤𝐶(Δ𝑡)²𝜕²𝑢𝜕𝜏²²_{𝐿20,𝑡𝑛;𝐿2(Ω)}+ℎ⁴𝑢_𝑡²_{𝐿20,𝑡𝑛;𝐻2(Ω)}+ℎ²‖𝑢‖²_{𝐿∞0,𝑡𝑛;𝐻2(Ω)}+‖𝜎‖²_{𝐿∞0,𝑡𝑛;𝐻2(Ω)}. (47) In order to derive (27), set 𝐰_ℎ=𝜉^𝑛 in (24b) and employ Lemma 1 and assumption (A3) to give 𝜉^𝑛2=−𝑏∇𝑒^𝑛,𝜉^𝑛_ℎ−𝜂^𝑛,𝜉^𝑛−𝑏∇𝜌^𝑛,𝜉^𝑛_ℎ≤𝐶𝑒^𝑛2_1,ℎ+ℎ⁴𝜎^𝑛_2,Ω−𝑏−𝑏∇𝜌^𝑛,𝜉^𝑛_ℎ+14𝜉^𝑛2≤𝐶𝑒^𝑛2_1,ℎ+ℎ⁴𝜎^𝑛_2,Ω+𝑢^𝑛_2,Ω+12𝜉^𝑛2. (48) Combining (47) with (48) yields 𝜉^𝑛2≤𝐶(Δ𝑡)²𝜕²𝑢𝜕𝜏²²_{𝐿20,𝑡𝑛;𝐿2(Ω)}+ℎ⁴𝑢_𝑡²_{𝐿20,𝑡𝑛;𝐻2(Ω)}+ℎ²|𝑢|²_{𝐿∞0,𝑡𝑛;𝐻2(Ω)}+‖𝜎‖²_{𝐿∞0,𝑡𝑛;𝐻2(Ω)}. (49) By using of interpolation theory and the triangle inequality, (37), (47), and (49) lead to (25), (26), and (27), respectively, which are the desired results.

Remark 2. From (37), we have Δ𝑡^𝑛_𝑖=1𝑒^𝑖2_1,ℎ=Δ𝑡^𝑛_𝑖=1𝜋¹_ℎ𝑢−𝑢_ℎ^𝑖2_1,ℎ≤𝐶(Δ𝑡)²𝜕²𝑢𝜕𝜏²²_{𝐿20,𝑡𝑛;𝐿2(Ω)}+ℎ⁴𝑢_𝑡²_{𝐿20,𝑡𝑛;𝐻2(Ω)}+‖𝑢‖²_{𝐿∞0,𝑡𝑛;𝐻2(Ω)}+‖𝜎‖²_{𝐿∞(0,𝑡𝑛;𝐻2(Ω))2}. (50) This byproduct can be regarded as the superclose result between 𝜋¹_ℎ𝑢 and 𝑢_ℎ in mean broken 𝐻¹-norm. It seems that both (25) and (50) have never been seen in the existing studies. At the same time, by employing the new characteristic nonconforming MFE scheme, we can also obtain the same error estimate of (27) as traditional characteristic MFEM [10].

Remark 3. From the analysis of Theorem 2 in this paper, we may see that Lemma 1 is the key result leading to the successful optimal order error estimations. If we want to get higher order accuracy, similar to Lemma 1, the nonconforming finite elements for approximating 𝑢 should also possess a very special property, that is, the consistency error estimates with 𝑂(ℎ²) order, and satisfy (18). For the famous nonconforming Wilson element [32] whose shape function is span{1,𝑥,𝑦,𝑥²,𝑦²}, by a counter-example, it has been proven in [32] that its consistency error estimate is of 𝑂(ℎ) order and cannot be improved any more. For the rotated bilinear 𝑄₁ element [33] whose shape function is span{1,𝑥,𝑦,𝑥²−𝑦²}, although its consistency error with 𝑂(ℎ²) order and (∇(𝑢−𝜋¹_ℎ𝑢),∇𝑣_ℎ)_ℎ=0 on square meshes is satisfied, the second term of (18) is not valid. Thus when they are applied to (1) on new characteristic mixed finite element scheme, up to now, the optimal order error estimates of (25), (26), and (27) cannot be obtained directly.

5. Numerical Example

In order to verify our theoretical analysis in previous sections, we consider the convection-dominated diffusion problem (1) as follows: 𝑢_𝑡+𝑢_𝑥+𝑢_𝑦−10⁻⁴𝑢_𝑥𝑥+𝑢_𝑦𝑦,(),𝑢(),𝑢=𝑓𝑥,𝑦,𝑡𝑥,𝑦,𝑡∈Ω×0,𝑇𝑥,𝑦,𝑡=0,𝑥,𝑦,𝑡∈𝜕Ω×0,𝑇𝑥,𝑦,0=𝑢₀𝑥,𝑦,𝑥,𝑦∈Ω (51) with Ω=[0,1]×[0,1],𝐚(𝑥,𝑦)=(1,1), and 𝑏(𝑥,𝑦)=10⁻⁴.

The right hand term 𝑓(𝑥,𝑦,𝑡) is taken such that 𝑢=𝑒^−𝑡sin(𝜋𝑥)sin(2𝜋𝑦),𝜎=−10⁻⁴𝑒^−𝑡(𝜋cos(𝜋𝑥) sin(2𝜋𝑦),2𝜋sin(𝜋𝑥)cos(2𝜋𝑦)) are the exact solutions.

We first divide the domain Ω into 𝑚 and 𝑛 equal intervals along 𝑥-axis and 𝑦-axis and the numerical results at different times are listed in Tables 1, 2, and 3 and pictured in Figures 1, 2, 3, and 4, respectively. (𝑢_ℎ,𝐩_ℎ) denotes the characteristic nonconforming MFE solution of the problem (15a), (15b), and (15c). Δ𝑡 represents the time step and the experiment is done with Δ𝑡=ℎ². 𝛼 stands for the convergence order.

Table 1 

Numerical results of ‖𝑢−𝑢_ℎ‖_1,ℎ.

Table 2 

Numerical results of ‖𝑢−𝑢_ℎ‖.

Table 3 

Numerical results of ‖𝜎−𝜎_ℎ‖.

Figure 1 

Errors at 𝑡=0.4.

Figure 2 

Errors at 𝑡=0.5.

Figure 3 

Errors at 𝑡=0.8.

Figure 4 

Errors at 𝑡=0.9.

It can be seen from the above Tables 1, 2, and 3 that ‖𝑢−𝑢_ℎ‖_1,ℎ and ‖𝜎−𝜎_ℎ‖ are convergent at optimal rate of 𝑂(ℎ) and ‖𝑢−𝑢_ℎ‖ is convergent at optimal rate of 𝑂(ℎ²), respectively, which coincide with our theoretical investigation in Section 4.

Acknowledgments

The research was supported by the National Natural Science Foundation of China (Grant nos. 10971203, 11101384, and 11271340) and the Specialized Research Fund for the Doctoral Program of Higher Education (Grant no. 20094101110006). The author would like to thank the referees for their helpful suggestions.