The Upwind Finite Volume Element Method for Two-Dimensional Burgers Equation

Abstract

A finite volume element method for approximating the solution to two-dimensional Burgers equation is presented. Upwind technique is applied to handle the nonlinear convection term. We present the semi-discrete scheme and fully discrete scheme, respectively. We show that the schemes are convergent to order one in space in 𝐿²-norm. Numerical experiment is presented finally to validate the theoretical analysis.

1. Introduction

We consider the following two-dimensional Burgers equation [1–3]: (a)𝜕𝑢𝜕𝑡+𝑢𝜕𝑢𝜕𝑥₁+𝑣𝜕𝑢𝜕𝑥₂=𝜁Δ𝑢,𝑥=𝑥₁,𝑥₂(],(b)∈Ω,𝑡∈𝐽=0,𝑇𝜕𝑣𝜕𝑡+𝑢𝜕𝑣𝜕𝑥₁+𝑣𝜕𝑣𝜕𝑥₂()(c)𝑢()(𝑥)()(𝑥)(d)=𝜁Δ𝑣,𝑥,𝑡∈Ω×𝐽,𝑥,0=𝜙,𝑣𝑥,0=𝜓,𝑥∈Ω,𝑢=𝑔₁,𝑣=𝑔₂,(𝑥,𝑡)∈𝜕Ω×𝐽, (1) for the unknown functions 𝑢 and 𝑣 in a bounded spatial domain Ω⊂ℝ², over a time interval [0,𝑇]. The coefficient 𝜁 is a positive number.

Burgers equation is the simplest nonlinear convection-diffusion model [1]. It is often used in modeling such physical phenomena as turbulence, shocks, and so forth. The study of Burgers equation has been a very active area because of its importance.

It is well known that strictly parabolic discretization schemes applied to Burgers equation do not work well when it is advection dominated. Effective discretization schemes recognize to some extent the hyperbolic nature of the equation.

The finite volume element method (FVEM) [4–12] is an important discretization technique for partial differential equations, especially those that arise from physical conservation laws. FVEM has ability to be faithful to the physics in general and conservation in particular, to produce simple stencils, and to treat effectively Neumann boundary conditions and nonuniform grids, and so forth.

Liang [11, 12] combined the upwind technique and the FVEM to handle the linear convection-dominated problems. In this paper, we will consider upwind finite volume element method for the approximation of (1). Upwind approximation is applied to handle the nonlinear convection term. The semi-discrete and fully discrete schemes are defined, respectively. We prove that they are both convergent to order one in space. Numerical experiments are presented finally to validate the theoretical analysis.

In this paper, we use the following Sobolev spaces and the norms associated with these spaces: 𝐿²(Ω)=𝑓∶_Ω𝑓²𝑑𝑥<∞,𝑓=_Ω𝑓²𝑑𝑥^1/2,𝐿^∞(Ω)=𝑓∶esssup_Ω𝑓<∞,𝑓_∞=esssup_Ω𝑓,𝐻^𝑚(Ω)=𝜕𝑓∶^|𝛼|𝑓𝜕𝑥^𝛼∈𝐿²(Ω),|𝛼|,𝑓≤𝑚𝑚≥0,_𝑚=_{|𝛼|≤𝑚}𝜕^|𝛼|𝑓𝜕𝑥^𝛼21/2,𝑊^𝑚_∞(Ω)=𝜕𝑓∶^|𝛼|𝑓𝜕𝑥^𝛼∈𝐿^∞(Ω),|𝛼|,𝑓≤𝑚𝑚≥0._𝑚,∞=max_{|𝛼|≤𝑚}𝜕^|𝛼|𝑓𝜕𝑥^𝛼_𝐿∞, (2) In particular, 𝐻⁰(Ω)=𝐿²(Ω),  𝑊⁰_∞(Ω)=𝐿^∞(Ω). Let [𝑎,𝑏]⊂[0,𝑇] and let 𝑋 be any of the spaces just defined. If 𝑓(𝑥,𝑡) represents functions on Ω×[𝑎,𝑏], we set 𝐻^𝑚(𝑎,𝑏;𝑋)=𝑓∶^𝑏_𝑎𝜕^𝛼𝑓𝜕𝑡^𝛼(⋅,𝑡)²_𝑋,𝑓𝑑𝑡<∞,𝛼≤𝑚_{𝐻𝑚(𝑎,𝑏;𝑋)}=^𝑚_𝛼=0^𝑏_𝑎𝜕^𝛼𝑓𝜕𝑡^𝛼(⋅,𝑡)²_𝑋𝑑𝑡^1/2𝑊,𝑚≥0,_𝑚∞(𝑎,𝑏;𝑋)=𝑓∶esssup_[𝑎,𝑏]𝜕^𝛼𝑓𝜕𝑡^𝛼(⋅,𝑡)_𝑋,𝑓<∞,𝛼≤𝑚_{𝑊𝑚∞(𝑎,𝑏;𝑋)}=max_{0≤𝛼≤𝑚}esssup_[]𝑎,𝑏𝜕^𝛼𝑓𝜕𝑡^𝛼(⋅,𝑡)_𝑋𝐿,𝑚≥0,₂(𝑎,𝑏;𝑋)=𝐻⁰(),𝐿𝑎,𝑏;𝑋_∞(𝑎,𝑏;𝑋)=𝑊⁰_∞(𝑎,𝑏;𝑋). (3) If [𝑎,𝑏]=[0,𝑇], we drop it from the notation. We also drop Ω; thus, we write 𝐿^∞(𝑊¹_∞) for 𝐿^∞(0,𝑇;𝑊¹_∞(Ω)).

If 𝑤=(𝑤₁,𝑤₂) is a vector function, we say that 𝑤∈𝑋² if 𝑤₁∈𝑋 and 𝑤₂∈𝑋.

An outline of the paper follows. In the next section we define the upwind finite volume element schemes for (1). Some Lemmas are presented in Section 3. We derive the 𝐿²-norm error estimates for the semi-discrete scheme and the fully discrete scheme in Sections 4 and 5, respectively. Finally in Section 6, we give some numerical experiments.

Throughout the paper we will denote by 𝐶 and 𝐶_𝑖(𝑖=1,2,…) generic constants independent of the mesh parameters, which may take different values in different occurrences.

2. The Approximation Schemes

In order to rewrite (1) as the vector form we define some vector notations. The gradient of a vector function 𝑤=(𝑤₁,𝑤₂)∶ℝ²→ℝ² is a matrix, and the divergence of a matrix function 𝐴=(𝑎_𝑖𝑗)_{1≤𝑖,𝑗≤2}∶ℝ²→ℝ^2×2 is a vector ∇𝑤=𝜕𝑤_𝑖𝜕𝑥_{𝑗1≤𝑖,𝑗≤2},∇⋅𝐴=²_𝑗=1𝜕𝑎_1𝑗𝜕𝑥_𝑗,²_𝑗=1𝜕𝑎_2𝑗𝜕𝑥_𝑗. (4) Consequently, we have for a vector function 𝑤=(𝑤₁,𝑤₂)Δ𝑤=∇⋅∇𝑤=Δ𝑤₁,Δ𝑤₂. (5) Let 𝜃=(𝑢,𝑣), 𝜃₀(𝑥)=(𝑢₀(𝑥),𝑣₀(𝑥)), and let 𝑔=(𝑔₁,𝑔₂); then the system (1) can be written as the following vector form: (a)𝜕𝜃()(b)𝜃()𝜕𝑡+𝜃⋅∇𝜃−𝜁Δ𝜃=0,𝑥,𝑡∈Ω×𝐽,𝑥,0=𝜃₀(𝑥)(c)𝜃()(),𝑥∈Ω,𝑥,𝑡=𝑔,𝑥,𝑡∈𝜕Ω×𝐽, (6) where 𝜃⋅∇𝜃=𝑢𝜕𝜃𝜕𝑥₁+𝑣𝜕𝜃𝜕𝑥₂. (7)

Let 𝒯_ℎ={𝐾} be a triangulation of the domain Ω, and as usual, we assume the triangles 𝐾 to be shape regular. Denote by Ω_ℎ={𝑃_𝑖} the set of the vertices of all the triangles 𝐾, and let Ω_ℎ=Ω_ℎ⧵𝜕Ω. For a given triangulation 𝒯_ℎ, we construct a dual mesh 𝒯^∗_ℎ whose elements are called control volumes. Each triangle 𝐾∈𝒯_ℎ can be divided into three subdomains by connecting an inner point of the triangle to the midpoints of the three edges. Around each 𝑃_𝑖∈Ω_ℎ, we associate a control volume 𝐾^∗_𝑖=𝐾^∗_𝑃𝑖, which consists of the union of subregions having 𝑃_𝑖 as a vertex. For a vertex 𝑃_𝑖∈𝜕Ω, we can define its control volume in a similar way. Then we define the dual partition 𝒯^∗_ℎ={𝐾^∗_𝑃𝑖,𝑃_𝑖∈Ω_ℎ} to be the union of all the control volumes. Usually we can choose the inner points as the barycenters or the circum centers, and in the later case we assume that all the inner angles of each triangle are not larger than 𝜋/2. We will use the barycenters duel mesh in this paper, while, with some trivial changes, our analysis can be also applied to the case when the circum centers are used.

We now characterize the finite-dimensional spaces which will be employed in approximating (6). For the sake of simplicity, we will assume that 𝑔₁=𝑔₂=0. We define the following finite dimensional spaces: 𝑃_ℎ=𝜔_ℎ∈𝐻¹₀(Ω),𝜔_ℎ𝐾∈𝒫₁(𝐾),𝐾∈𝑇_ℎ,𝑌_ℎ=𝜑_ℎ∈𝐿²(Ω),𝜑_{ℎ𝐾∗𝑖}∈𝒫₀𝐾^∗,𝑃_𝑖∈Ω_ℎ;𝜑_{ℎ𝐾∗𝑖}=0,𝑃_𝑖,𝑈∈𝜕Ω_ℎ=𝑃²_ℎ,𝑉_ℎ=𝑌²_ℎ, (8) where 𝒫_𝑙(𝑉)(𝑙=0,1) denotes the set of polynomials on 𝑉 with a degree of not more than 𝑙.

Multiplying (6a) by test function 𝑧∈𝑉_ℎ and integrating by parts yield 𝜕𝜃𝜕𝑡,𝑧+𝐵₁(𝜃;𝜃,𝑧)+𝐵₂()()𝜃;𝜃,𝑧+𝐴𝜃,𝑧=0,∀𝑧∈𝑉_ℎ, (9) where 𝐵₁𝜑;𝑤,𝑧=−_{𝑃𝑖∈Ωℎ}𝑧𝑃_𝑖⋅_𝐾∗𝑖()𝐵∇⋅𝑤𝜑𝑑𝑥,₂𝜑;𝑤,𝑧=_{𝑃𝑖∈Ωℎ}𝑧𝑃_𝑖⋅_{𝜕𝐾∗𝑖}𝐴()𝜑⋅𝜈𝑤𝑑𝑠,𝑤,𝑧=_{𝑃𝑖∈Ωℎ}𝑧𝑃_𝑖⋅_{𝜕𝐾∗𝑖}𝜁(𝜈⋅∇𝑤)𝑑𝑠, (10) here 𝜈 is the unit outward normal vector of 𝜕𝐾^∗_𝑖.

Now we approximate 𝐵₂(𝜑;𝑤,𝑧) by using the upwind technique.

Let Λ_𝑖={𝑗∶𝑃_𝑗 is adjoint with 𝑃_𝑖}. Assuming that 𝑗∈Λ_𝑖, let Γ_𝑖𝑗=𝜕𝐾^∗_𝑖∩𝜕𝐾^∗_𝑗 and 𝛾_𝑖𝑗 is the length of Γ_𝑖𝑗. Denote by 𝜈_𝑖𝑗 the unit outward normal vector of Γ_𝑖𝑗 when Γ_𝑖𝑗 is regarded as the boundary of 𝐾^∗_𝑖. Define 𝛽_𝑖𝑗𝜑=_Γ𝑖𝑗𝜑⋅𝜈_𝑖𝑗𝑑𝑠. (11) Let 𝛽⁺_𝑖𝑗𝜑=max𝛽_𝑖𝑗𝜑,0,𝛽⁻_𝑖𝑗𝜑=max−𝛽_𝑖𝑗𝜑,0_{𝜕𝐾∗𝑖},𝜑⋅𝜈𝑤𝑑𝑠≈_{𝑗∈Λ𝑖}𝛽⁺_𝑖𝑗𝜑𝑤𝑃_𝑖−𝛽⁻_𝑖𝑗𝜑𝑤𝑃_𝑗. (12) The upwind discretization of the nonlinear term 𝐵₂(𝜑;𝑤,𝑧) is defined by the form 𝐵_2ℎ𝜑;𝑤,𝑧=_{𝑃𝑖∈Ωℎ𝑗∈Λ𝑖}𝛽⁺_𝑖𝑗𝜑𝑤𝑃_𝑖−𝛽⁻_𝑖𝑗𝜑𝑤𝑃_𝑗⋅𝑧𝑃_𝑖. (13) Using the heaviside function 𝐻(𝑟)=1,𝑟≥0,0,𝑟<0, (14) we can write 𝐵_2ℎ(𝜑;𝑤,𝑧) as 𝐵_2ℎ=𝜑;𝑤,𝑧_{𝑃𝑖∈Ωℎ𝑗∈Λ𝑖}𝛽_𝑖𝑗𝜑×𝐻𝛽_𝑖𝑗𝜑𝑤𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜑𝑤𝑃_𝑗⋅𝑧𝑃_𝑖. (15)

Introduce the interpolation operators Π_ℎ∶𝐻¹₀(Ω)→𝑃_ℎ and Π^∗_ℎ∶𝑃_ℎ→𝑌_ℎ, respectively. For 𝑤=(𝑤₁,𝑤₂), define Π_ℎ𝑤=(Π_ℎ𝑤₁,Π_ℎ𝑤₂) and Π^∗_ℎ𝑤=(Π^∗_ℎ𝑤₁,Π^∗_ℎ𝑤₂). Assuming that 𝑤∈𝐻²(Ω)², we can easily get the following interpolation estimates: 𝑤−Π_ℎ𝑤_𝑠≤ℎ^2−𝑠‖𝑤‖₂,𝑠=0,1. (16)

The semi-discrete upwind finite volume scheme of (6) is as follows: find 𝜃_ℎ∶[0,𝑇]→𝑈_ℎ such that 𝜕𝜃_ℎ𝜕𝑡,Π^∗_ℎ𝑧_ℎ+𝐵₁𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝑧_ℎ+𝐵_2ℎ𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝑧_ℎ+𝐴𝜃_ℎ,Π^∗_ℎ𝑧_ℎ=0,∀𝑧_ℎ∈𝑈_ℎ,𝜃_ℎ(𝑥,0)=𝜃_0ℎ(𝑥), (17) where 𝜃_0ℎ(𝑥) is the interpolation projection of 𝜃₀, that is, 𝜃_0ℎ(𝑥)=Π_ℎ𝜃₀.

Partition [0,𝑇] into 0=𝑡⁰<𝑡¹<⋯<𝑡^𝑁=𝑇, with 𝜏^𝑛=𝑡^𝑛−𝑡^𝑛−1. Our analysis is valid for variable time steps, but we drop the superscript from 𝜏 for convenience. For functions 𝑓 on Ω×𝐽, we write 𝑓^𝑛(𝑥) for 𝑓(𝑥,𝑡^𝑛). By approximating 𝜕𝜃_ℎ/𝜕𝑡 at the time 𝑡=𝑡_𝑛 with the backward difference 𝜕_𝑡𝜃^𝑛_ℎ=(𝜃^𝑛_ℎ−𝜃_ℎ^𝑛−1)/𝜏, we define the fully discrete upwind finite volume scheme for (6) as follows: find 𝜃^𝑛_ℎ∈𝑈_ℎ, such that 𝜕_𝑡𝜃^𝑛_ℎ,Π^∗_ℎ𝑧_ℎ+𝐵_1ℎ𝜃_ℎ^𝑛−1;𝜃^𝑛_ℎ,Π^∗_ℎ𝑧_ℎ+𝐵_2ℎ𝜃_ℎ^𝑛−1;𝜃^𝑛_ℎ,Π^∗_ℎ𝑧_ℎ+𝐴𝜃^𝑛_ℎ,Π^∗_ℎ𝑧_ℎ=0,𝑛≥1,∀𝑧_ℎ∈𝑈_ℎ,𝜃⁰_ℎ=𝜃_0ℎ. (18)

3. Some Lemmas

Now we present several lemmas. Let 𝑤_ℎ=(𝑤₁,𝑤₂)∈𝑈_ℎ,𝑤_ℎ=(𝑤₁,𝑤₂)∈𝑈_ℎ.

Lemma 1. (i)   Π^∗_ℎ is a self-adjoint operator, that is, 𝑤_ℎ,Π^∗_ℎ𝑤_ℎ=𝑤_ℎ,Π^∗_ℎ𝑤_ℎ,∀𝑤_ℎ,𝑤_ℎ∈𝑈_ℎ. (19)
(ii) Let |‖𝑤_ℎ‖|=(𝑤_ℎ,Π^∗_ℎ𝑤_ℎ)^1/2. Then, for some positive constants 𝐶₁ and 𝐶₂ that are independent of ℎ, 𝐶₁𝑤_ℎ≤𝑤_ℎ≤𝐶₂𝑤_ℎ,∀𝑤_ℎ∈𝑈_ℎ. (20)

Proof. It is easy to know that 𝑤_ℎ,Π^∗_ℎ𝑤_ℎ=𝑤₁,Π^∗_ℎ𝑤₁+𝑤₂,Π^∗_ℎ𝑤₂,𝑤_ℎ²=𝑤_ℎ,Π^∗_ℎ𝑤_ℎ=𝑤₁,Π^∗_ℎ𝑤₁+𝑤₂,Π^∗_ℎ𝑤₂. (21) From [4] we know that for 𝑤₁, 𝑤₂, 𝑤₁, 𝑤₂∈𝑃_ℎ, 𝑤₁,Π^∗_ℎ𝑤₁=𝑤₁,Π^∗_ℎ𝑤₁,𝑤₂,Π^∗_ℎ𝑤₂=𝑤₂,Π^∗_ℎ𝑤₂,𝐶²₁𝑤₁²≤𝑤₁,Π^∗_ℎ𝑤₁≤𝐶²₂𝑤₁²,𝐶²₁𝑤₂²≤𝑤₂,Π^∗_ℎ𝑤₂≤𝐶²₂𝑤₂², (22) where 𝐶₁ and 𝐶₂ are some positive constants that are independent of ℎ. Thus we obtain (19) and (20) immediately.

Lemma 2. For the bilinear form 𝐴(⋅,Π^∗_ℎ⋅), one has the following conclusions:(i)For 𝑤_ℎ,𝑤_ℎ∈𝑈_ℎ, one has 𝐴𝑤_ℎ,Π^∗_ℎ𝑤_ℎ=𝐴𝑤_ℎ,Π^∗_ℎ𝑤_ℎ. (23)(ii)There exists a positive constant 𝐶 such that 𝐴𝑤−Π_ℎ𝑤,Π^∗_ℎ𝑤_ℎ≤𝐶ℎ‖𝑤‖₂𝑤_ℎ1,∀𝑤∈𝐻²(Ω)²,𝑤_ℎ∈𝑈_ℎ. (24)(iii)There exists a positive constant 𝛼 such that 𝐴𝑤_ℎ,Π^∗_ℎ𝑤_ℎ≥𝛼𝑤_ℎ²₁,∀𝑤_ℎ∈𝑈_ℎ. (25)

Proof. For 𝜙,𝜓∈𝑃_ℎ, define the bilinear form 𝑎𝜙,Π^∗_ℎ𝜓=−_{𝑃𝑖∈Ωℎ}𝜓𝑃_{𝑖𝜕𝐾∗𝑖}𝜁∇𝜙⋅𝜈𝑑𝑠. (26) Then, we get 𝐴𝑤_ℎ,Π^∗_ℎ𝑤_ℎ=−_{𝑃𝑖∈Ωℎ}𝑤_ℎ𝑃_𝑖⋅_{𝜕𝐾∗𝑖}𝜁𝜈⋅∇𝑤_ℎ𝑑𝑠=−_{𝑃𝑖∈Ωℎ}𝑤₁𝑃_{𝑖𝜕𝐾∗}𝜁∇𝑤₁⋅𝜈𝑑𝑠−_{𝑃𝑖∈Ωℎ}𝑤₂𝑃_{𝑖𝜕𝐾∗𝑖}𝜁∇𝑤₂𝑤⋅𝜈𝑑𝑠=𝑎₁,Π^∗_ℎ𝑤₁+𝑎𝑤₂,Π^∗_ℎ𝑤₂,𝐴𝑤−Π_ℎ𝑤,Π^∗_ℎ𝑤_ℎ=𝑎𝑤₁−Π_ℎ𝑤₁,Π^∗_ℎ𝑤₁+𝑎𝑤₂−Π_ℎ𝑤₂,Π^∗_ℎ𝑤₂. (27) By combining the above results and the corresponding conclusions for 𝑎(⋅,Π^∗_ℎ⋅) in [4], we can obtain (23)–(25).

Lemma 3. For 𝜑∈(𝑊⁰_∞(Ω))²,𝜃∈(𝐻¹₀(Ω))²,𝜑_ℎ∈𝑈_ℎ, and 𝑧_ℎ∈𝑈_ℎ, one has 𝐵₂𝜑;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜑_ℎ;Π_ℎ𝜃,Π^∗_ℎ𝑧_ℎ≤𝑧_ℎ1ℎ𝜑_∞|𝜃|₁+‖𝜃‖_∞𝜑−𝜑_ℎ+ℎ𝜑−𝜑_ℎ1. (28)

Proof. First we have 𝐵₂𝜑;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜑_ℎ;Π_ℎ𝜃,Π^∗_ℎ𝑧_ℎ=𝐵₂𝜑;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜑_ℎ;𝜃,Π^∗_ℎ𝑧_ℎ+𝐵_2ℎ𝜑_ℎ;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜑_ℎ;Π_ℎ𝜃,Π^∗_ℎ𝑧_ℎ. (29) Noting that 𝜃(𝑃_𝑖)=Π_ℎ𝜃(𝑃_𝑖), 𝑃_𝑖∈Ω_ℎ, we can easily deduce 𝐵_2ℎ𝜑_ℎ;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜑_ℎ;Π_ℎ𝜃,Π^∗_ℎ𝑧_ℎ=0 (30) by the definition of 𝐵_2ℎ(⋅;⋅,Π^∗_ℎ⋅). Now we only need to bound 𝐵₂𝜑;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜑_ℎ;𝜃,Π^∗_ℎ𝑧_ℎ=_{𝑃𝑖∈Ωℎ}𝑧_ℎ𝑃_𝑖⋅_{𝑗∈Λ𝑖Γ𝑖𝑗}𝜑⋅𝜈_𝑖𝑗𝜃𝑑𝑠−_{𝑃𝑖∈Ωℎ}𝑧_ℎ𝑃_𝑖⋅_{𝑗∈Λ𝑖}𝛽_𝑖𝑗𝜑_ℎ×𝐻𝛽_𝑖𝑗𝜑_ℎ𝜃𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜑_ℎ𝜃𝑃_𝑗=_{𝑃𝑖∈Ωℎ}𝑧_ℎ𝑃_𝑖⋅_{𝑗∈Λ𝑖Γ𝑖𝑗}𝜑⋅𝜈_𝑖𝑗×𝜃−𝐻𝛽_𝑖𝑗𝜑_ℎ𝜃𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜑_ℎ𝜃𝑃_𝑗𝑑𝑠+_{𝑃𝑖∈Ωℎ}𝑧_ℎ𝑃_𝑖⋅_{𝑗∈Λ𝑖Γ𝑖𝑗}𝜑−𝜑_ℎ⋅𝜈_𝑖𝑗×𝐻𝛽𝑑𝑠_𝑖𝑗𝜑_ℎ𝜃𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜑_ℎ𝜃𝑃_𝑗. (31) We denote the last two terms on the right-hand side of (31) by 𝐼₁ and 𝐼₂, respectively. We now turn to analyze the two terms. Noting that 𝛽_𝑖𝑗=−𝛽_𝑗𝑖, we rewrite 𝐼₁ as 𝐼₁=12_{𝐾∈𝒯ℎ𝑖,𝑗∈Λ𝐾}𝑧_ℎ𝑃_𝑖−𝑧_ℎ𝑃_𝑗⋅_{Γ𝑖𝑗∩𝐾}𝜑⋅𝜈_𝑖𝑗𝐻𝛽_𝑖𝑗𝜑_ℎ𝜃−𝜃𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜑_ℎ𝜃−𝜃𝑃_𝑗𝑑𝑠, (32) here Λ_𝐾 is the set of vertex of 𝐾. From the Taylor’s Formula and the linear property of 𝑧_ℎ=(𝑧₁,𝑧₂), we obtain that 𝑧_ℎ𝑃_𝑖−𝑧_ℎ𝑃_𝑗²≤ℎ²∇𝑧₁²+∇𝑧₂²=ℎ²|𝑧|²_1,𝐾. (33)
Applying the trace inequality, we get 𝜃−𝜃𝑃_{Γ𝑖𝑗∩𝐾𝑖}𝑑𝑠≤𝐶ℎ^1/2_{Γ𝑖𝑗∩𝐾}𝜃−𝜃𝑃_𝑖²𝑑𝑠^1/2≤𝐶ℎ^1/2_{Γ𝑖𝑗∩𝐾}𝑢−𝑢𝑃_𝑖²+𝑣−𝑣𝑃_𝑖²𝑑𝑠^1/2≤𝐶ℎℎ⁻¹𝑢−𝑢𝑃_𝑖0,𝐾+𝑢−𝑢𝑃_𝑖1,𝐾²+ℎ⁻¹𝑣−𝑣𝑃_𝑖0,𝐾+𝑣−𝑣𝑃_𝑖1,𝐾^21/2≤𝐶ℎ|𝑢|²_1,𝐾+|𝑣|²_1,𝐾^1/2=𝐶ℎ|𝜃|_1,𝐾. (34) Similarly, we can deduce that 𝜃−𝜃𝑃_{Γ𝑖𝑗∩𝐾𝑗}𝑑𝑠≤𝐶ℎ|𝜃|_1,𝐾. (35) We conclude that 𝐼₁≤𝐶ℎ𝜑_∞𝑧_ℎ1|𝜃|₁. (36)
The similar argument yields the estimate 𝐼₂=12_{𝐾∈𝒯ℎ𝑖,𝑗∈Λ𝐾Γ𝑖𝑗∩𝐾}𝜑−𝜑_ℎ×𝐻𝛽⋅𝑛𝑑𝑠_𝑖𝑗𝜑_ℎ𝜃𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜑_ℎ𝜃𝑃_𝑗⋅𝑧_ℎ𝑃_𝑖−𝑧_ℎ𝑃_𝑗≤𝐶‖𝜃‖_∞𝑧_ℎ1𝜑−𝜑_ℎ+ℎ𝜑−𝜑_ℎ1. (37)
Substituting the estimates (36) and (37) into (31), we obtain 𝐵₂𝜑;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜑_ℎ;𝜃,Π^∗_ℎ𝑧_ℎ≤𝐶𝑧_ℎ1ℎ𝜑_∞|𝜃|₁+‖𝜃‖_∞×𝜑−𝜑_ℎ+ℎ𝜑−𝜑_ℎ1. (38) This yields the desired result immediately.

4. Error Bounds for Semi-Discrete Scheme

Theorem 4. Assume that 𝜃 and 𝜃_ℎ are solutions to (6) and (17), respectively. also assumes that 𝜃 is regular enough. Then there exists a positive constant 𝐶 such that 𝜃−𝜃_ℎ≤𝐶ℎ, (39) where 𝐶 depends on principally ‖𝜃₀‖₂, ‖𝜃‖_{𝐿∞((𝑊1∞)}²₎, and ‖𝜃‖_{𝐻1((𝐻2)}²₎.

Proof. We derive the following error equation from (6) and (17): 𝜕𝜃−𝜕𝑡𝜕𝜃_ℎ𝜕𝑡,Π^∗_ℎ𝑧_ℎ+𝐵₁𝜃;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵₁𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝑧_ℎ+𝐵₂𝜃;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝑧_ℎ+𝐴𝜃−𝜃_ℎ,Π^∗_ℎ𝑧_ℎ=0. (40) Let 𝜌=𝜃−Π_ℎ𝜃,𝜉=Π_ℎ𝜃−𝜃_ℎ. We rewrite the previously mentioned equation as 𝜕𝜉𝜕𝑡,Π^∗_ℎ𝑧_ℎ+𝐵₁𝜃;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵₁𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝑧_ℎ+𝐵₂𝜃;𝜃,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝑧_ℎ+𝐴𝜉,Π^∗_ℎ𝑧_ℎ=−𝜕𝜌𝜕𝑡,Π^∗_ℎ𝑧_ℎ−𝐴𝜌,Π^∗_ℎ𝑧_ℎ. (41) We choose 𝑧_ℎ=𝜉 in (41) to get 𝜕𝜉𝜕𝑡,Π^∗_ℎ𝜉+𝐴𝜉,Π^∗_ℎ𝜉=−𝜕𝜌𝜕𝑡,Π^∗_ℎ𝜉−𝐴𝜌,Π^∗_ℎ𝜉−𝐵₁𝜃;𝜃,Π^∗_ℎ𝜉−𝐵₁𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝜉−𝐵₂𝜃;𝜃,Π^∗_ℎ𝜉−𝐵_2ℎ𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝜉. (42) Using Lemmas 1, 2 and Young’s inequality, we have 12𝑑𝜉𝑑𝑡²+𝛼𝜉²₁≤𝐶𝜕𝜌𝜕𝑡²+𝜉²+ℎ²‖𝜃‖²₂+𝜀𝜉²₁+𝐵₁𝜃;𝜃,Π^∗_ℎ𝜉−𝐵₁𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝜉+𝐵₂𝜃;𝜃,Π^∗_ℎ𝜉−𝐵_2ℎ𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝜉. (43) Now we bound the last two terms on the right-hand side of (43). We need the following induction hypothesis: 1logℎ^1/2𝜉(𝑠)⟶0,ℎ⟶0,0≤𝑠<𝑡,0<𝑡≤𝑇. (44)
We know from [13] that 𝜙_∞1≤𝐶logℎ^1/2𝜙₁,∀𝜙∈𝑃_ℎ. (45) This implies that 𝜑_∞1≤𝐶logℎ^1/2𝜑₁,∀𝜑∈𝑈_ℎ. (46) Also we have the following inverse inequality: 𝜑₁≤𝐶ℎ⁻¹𝜑,∀𝜑∈𝑈_ℎ. (47)
Using (46), we get 𝐵₁𝜃;𝜃,Π^∗_ℎ𝜉−𝐵₁𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝜉≤_{𝑃𝑖∈Ωℎ}𝜉𝑃_{𝑖𝐾∗𝑖}(∇⋅𝜃)𝜃−∇⋅𝜃_ℎ𝜃_ℎ𝑑𝑥≤_{𝑃𝑖∈Ωℎ}𝜉𝑃_{𝑖𝐾∗𝑖}|(∇⋅𝜃)|𝜃−𝜃_ℎ+∇⋅𝜃−∇⋅𝜃_ℎ×𝜉+Π_ℎ𝜃‖𝑑𝑥≤𝐶‖∇⋅𝜃_∞𝜃−𝜃_ℎ𝜉+∇⋅𝜃−𝜃_ℎ×𝜉𝜉_∞+‖𝜃‖_∞≤𝐶𝜌₁+𝜉₁+1logℎ^1/2𝜉₁×𝜌₁+𝜉₁𝜉≤𝐶𝜌²₁+𝜉²+𝜉1logℎ^1/2𝜉²₁+𝜀𝜉²₁. (48) Next, we write 𝐵₂𝜃;𝜃,Π^∗_ℎ𝜉−𝐵_2ℎ𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝜉≤𝐵₂𝜃;𝜃,Π^∗_ℎ𝜉−𝐵_2ℎ𝜃_ℎ;Π_ℎ𝜃,Π^∗_ℎ𝜉+𝐵_2ℎ𝜃_ℎ;𝜉,Π^∗_ℎ𝜉=𝐷₁+𝐷₂. (49) By Choosing 𝜑=𝜃,𝜑_ℎ=𝜃_ℎ, and 𝑧_ℎ=𝜉 in Lemma 3, using (47) and the Young’s inequality, we can obtain 𝐷₁≤𝐶𝜉₁ℎ‖𝜃‖_∞|𝜃|₁+‖𝜃‖_∞×𝜃−𝜃_ℎ+ℎ𝜃−𝜃_ℎ1≤𝐶ℎ²|𝜃|²₁+𝜌²₁+𝜉²+𝜀𝜉²₁. (50) By an argument like (36) and then by (46) and (47), we have 𝐷₂=_{𝑃𝑖∈Ωℎ}𝜉𝑃_𝑖⋅_{𝑗∈Λ𝑖Γ𝑖𝑗}𝜃_ℎ⋅𝜈_𝑖𝑗×𝐻𝛽𝑑𝑠_𝑖𝑗𝜉𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜉𝑃_𝑗≤12_{𝐾∈𝑇ℎ𝑖,𝑗∈Λ𝐾}𝜉𝑃_𝑖−𝜉𝑃_𝑗×𝐻𝛽_𝑖𝑗𝜉𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜉𝑃_𝑗×_{Γ𝑖𝑗∩𝐾}𝜉⋅𝜈_𝑖𝑗+1𝑑𝑠2_{𝐾∈𝑇ℎ𝑖,𝑗∈Λ𝐾}𝜉𝑃_𝑖−𝜉𝑃_𝑗×𝐻𝛽_𝑖𝑗𝜉𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜉𝑃_𝑗×_{Γ𝑖𝑗∩𝐾}Π_ℎ𝜃⋅𝜈_𝑖𝑗𝜉𝑑𝑠≤𝐶₁𝜉+ℎ𝜉₁𝜉_∞+Π_ℎ𝜃_∞𝜉₁𝜉≤𝐶𝜉²+𝜉1logℎ^1/2𝜉²₁+𝜀𝜉²₁. (51) Substituting (50) and (51) into (49), we get 𝐵₂𝜃;𝜃,Π^∗_ℎ𝜉−𝐵_2ℎ𝜃_ℎ;𝜃_ℎ,Π^∗_ℎ𝜉≤𝐶ℎ²|𝜃|²₁+𝜌²₁+𝜉²𝜉1+𝐶logℎ^1/2𝜉²₁+2𝜀𝜉²₁. (52) Make (43), (48), and (52) together to obtain 12𝑑𝜉𝑑𝑡²+𝛼𝜉²₁≤𝐶ℎ²‖𝜃‖²₂+𝜕𝜌𝜕𝑡²+𝜌²₁+𝜉²𝜉1+𝐶logℎ^1/2𝜉²₁+4𝜀𝜉²₁. (53) Integrating the previously mentioned equation from 0 to 𝑡 and noting (44), we obtain that 12𝜉²(𝑡)−𝜉²(0+𝛼)2^𝑡₀𝜉²₁≤ℎ𝑑𝜏₂‖𝜃‖²₂+^𝑡₀𝜕𝜌𝜕𝑡²𝑑𝜏+^𝑡₀𝜌²₁+𝑑𝜏_𝑡0𝜉²𝑑𝜏 (54) for sufficiently small ℎ and 𝜀. By using Lemma 2(ii) and the Gronwall’s inequality, we have that 𝜉²(𝑡)+^𝑡₀𝜉²₁ℎ𝑑𝜏≤𝐶²‖𝜃‖²₂+^𝑡₀𝜕𝜌𝜕𝑡²𝑑𝜏+^𝑡₀𝜌²₁𝑑𝜏. (55) It follows from the interpolation estimates that 𝜉²(𝑡)+^𝑡₀𝜉²₁ℎ𝑑𝜏≤𝐶²𝜃₀²₂+^𝑇₀𝜕𝜃𝜕𝑡²₂𝑑𝜏+ℎ^4𝑇₀𝜕𝜃𝜕𝑡²₂𝑑𝜏+ℎ^2𝑇₀‖𝜃‖²₂𝑑𝜏. (56)
Now we prove the induction hypothesis (44). Noting that ‖𝜉‖(0)=0, we know that (44) holds obviously for 𝑡=0. It follows from (56) that 1logℎ^1/2𝜉𝑡)1(≤𝐶ℎlogℎ^1/2⟶0,ℎ⟶0. (57) Then (44) holds for any 𝑡∈[0,𝑇].
By (16), we have 𝜌≤𝐶ℎ²‖𝜃‖₂≤𝐶ℎ²𝜃₀₂+^𝑡₀‖𝜃‖₂𝑑𝜏. (58) From the triangle inequality, we obtain 𝜃−𝜃_ℎ≤𝐶ℎ, (59) where 𝐶 depends on ‖𝜃₀‖₂, ‖𝜃‖_{𝐻1((𝑊1∞)}²₎, and ‖𝜃‖_{𝐻1((𝐻2)}²₎. We now complete the proof of the theorem.

5. Error Bound for the Fully Discrete Scheme

Theorem 5. Assume that 𝜃 satisfies the necessary regularities and the discretization parameters obey the relation 𝜏=𝑂(ℎ). Then the error of the approximation (18) of (6) satisfies max_{0≤𝑛≤𝑇/𝜏}𝜃^𝑛−𝜃^𝑛_ℎ≤𝐶{ℎ+𝜏}, (60) where 𝐶 depends on ‖𝜃₀‖₂, ‖𝜃‖_{𝐿∞((𝑊1∞)}²₎, ‖𝜃‖_{𝐻1((𝐻2)}²₎, and ‖𝜃‖_{𝐻2((𝐿2)}²₎.

Proof. Subtract (18) from (9) to obtain that 𝜕𝜃^𝑛𝜕𝑡−𝜕_𝑡𝜃^𝑛_ℎ,Π^∗_ℎ𝑧_ℎ+𝐵₁𝜃^𝑛;𝜃^𝑛,Π^∗_ℎ𝑧_ℎ−𝐵₁𝜃_ℎ^𝑛−1;𝜃^𝑛_ℎ,Π^∗_ℎ𝑧_ℎ+𝐵₂𝜃^𝑛;𝜃^𝑛,Π^∗_ℎ𝑧_ℎ−𝐵_2ℎ𝜃_ℎ^𝑛−1;𝜃^𝑛_ℎ,Π^∗_ℎ𝑧_ℎ+𝐴𝜃^𝑛−𝜃^𝑛_ℎ,Π^∗_ℎ𝑧_ℎ=0. (61) Choose 𝑧_ℎ=𝜉^𝑛 to obtain that 𝜕_𝑡𝜉^𝑛,Π^∗_ℎ𝜉^𝑛+𝐴𝜉^𝑛,Π^∗_ℎ𝜉^𝑛=−𝜕𝜃^𝑛𝜕𝑡−𝜕_𝑡𝜃^𝑛,Π^∗_ℎ𝜉^𝑛−𝜕_𝑡𝜌^𝑛,Π^∗_ℎ𝜉^𝑛−𝐴𝜌^𝑛,Π^∗_ℎ𝜉^𝑛−𝐵₁𝜃^𝑛;𝜃^𝑛,Π^∗_ℎ𝜉^𝑛−𝐵₁𝜃_ℎ^𝑛−1;𝜃^𝑛_ℎ,Π^∗_ℎ𝜉^𝑛−𝐵₂𝜃^𝑛;𝜃^𝑛,Π^∗_ℎ𝜉^𝑛−𝐵_2ℎ𝜃_ℎ^𝑛−1;𝜃^𝑛_ℎ,Π^∗_ℎ𝜉^𝑛. (62) For the left-hand side of (62), from Lemmas 1 and 2, we have 𝜕_𝑡𝜉^𝑛,Π^∗_ℎ𝜉^𝑛=1𝜏𝜉^𝑛−𝜉^𝑛−1,Π^∗_ℎ𝜉^𝑛=1𝜉2𝜏^𝑛−𝜉^𝑛−1,Π^∗_ℎ𝜉^𝑛+𝜉^𝑛−1+1𝜉2𝜏^𝑛−𝜉^𝑛−1,Π^∗_ℎ𝜉^𝑛−𝜉^𝑛−1≥1𝜉2𝜏^𝑛2−𝜉^𝑛−12,𝐴𝜉^𝑛,Π^∗_ℎ𝜉^𝑛≥𝛼𝜉^𝑛2₁. (63) We denote terms on the right-hand side of (62) by 𝑇₁,…,𝑇₅. Then, (62) can be rewritten as 1𝜉2𝜏^𝑛2−𝜉^𝑛−12+𝛼𝜉^𝑛2₁≤𝑇₁+⋯+𝑇₅. (64)
Now we estimates the terms 𝑇₁,…,𝑇₅ one by one. From the Taylor’s formula, we have 𝜕𝜃^𝑛𝜕𝑡−𝜕_𝑡𝜃^𝑛=1𝜏^𝑡𝑛_𝑡^𝑛−1𝑡−𝑡^𝑛−1𝜕²𝜃𝜕𝑡²𝑑𝑡. (65) It follows that 𝑇₁≤𝐶𝜕𝜃^𝑛𝜕𝑡−𝜕_𝑡𝜃^𝑛𝜉^𝑛≤𝐶𝜏^𝑡𝑛_𝑡𝑛−1𝜕²𝜃𝜕𝑡²²𝑑𝑡+𝜉^𝑛2. (66) For the next two terms, we have 𝑇₂≤𝐶𝜕_𝑡𝜌^𝑛𝜉^𝑛≤𝐶𝜏^−1𝑡𝑛_𝑡𝑛−1𝜕𝜌𝜕𝑡²𝑑𝑡+𝜉^𝑛2,𝑇₃≤𝐶ℎ𝜃^𝑛₂𝜉^𝑛₁≤𝐶ℎ²𝜃^𝑛2₂+𝜀𝜉^𝑛2₁. (67)
We make the following induction hypothesis: 𝜉^𝑛−11logℎ^1/2⟶0,ℎ⟶0,1≤𝑛≤𝐿. (68)
For 𝑇₄, using the similar argument as (48) and noting (68), we deduce that 𝑇₄=𝐵₁𝜃^𝑛;𝜃^𝑛,Π^∗_ℎ𝜉^𝑛−𝐵₁𝜃_ℎ^𝑛−1;𝜃^𝑛_ℎ,Π^∗_ℎ𝜉^𝑛≤_{𝑃𝑖∈Ωℎ}𝜉^𝑛𝑃_{𝑖𝐾∗𝑖}∇⋅𝜃^𝑛𝜃^𝑛−∇⋅𝜃^𝑛_ℎ𝜃_ℎ^𝑛−1≤𝑑𝑥_{𝑃𝑖∈Ωℎ}𝜉^𝑛𝑃_{𝑖𝐾∗𝑖}∇⋅𝜃^𝑛𝜃^𝑛−𝜃_ℎ^𝑛−1+∇⋅𝜃^𝑛−∇⋅𝜃^𝑛_ℎ×𝜉^𝑛−1+Π_ℎ𝜃^𝑛−1𝑑𝑥≤𝐶∇⋅𝜃^𝑛_∞𝜃^𝑛−𝜃_ℎ^𝑛−1𝜉^𝑛+∇⋅𝜃^𝑛−𝜃^𝑛_ℎ×𝜉^𝑛_∞𝜉^𝑛−1+Π_ℎ𝜃^𝑛−1_∞𝜉^𝑛≤𝐶𝜏^𝑡𝑛_𝑡𝑛−1𝜕𝜃𝜕𝑡²𝑑𝑡^1/2+𝜌^𝑛−1+𝜉^𝑛−1×∇⋅𝜃^𝑛_∞𝜉^𝑛+𝜌^𝑛₁+𝜉^𝑛₁×1logℎ^1/2𝜉^𝑛−1𝜉^𝑛₁+𝜃^𝑛−1_∞𝜉^𝑛≤𝐶𝜏^𝑡𝑛_𝑡𝑛−1𝜕𝜃𝜕𝑡²𝑑𝑡+𝜌^𝑛−12+𝜌^𝑛2₁+𝜉^𝑛−12+𝜉^𝑛21+𝐶logℎ^1/2𝜉^𝑛−1𝜉^𝑛2₁+𝜀𝜉^𝑛2₁. (69)
Now, we write 𝑇₅=𝐵₂𝜃^𝑛;𝜃^𝑛,Π^∗_ℎ𝜉^𝑛−𝐵_2ℎ𝜃_ℎ^𝑛−1;𝜃^𝑛_ℎ,Π^∗_ℎ𝜉^𝑛≤𝐵₂𝜃^𝑛;𝜃^𝑛,Π^∗_ℎ𝜉^𝑛−𝐵_2ℎ𝜃_ℎ^𝑛−1;Π_ℎ𝜃^𝑛,Π^∗_ℎ𝜉^𝑛+𝐵_2ℎ𝜃_ℎ^𝑛−1;𝜉^𝑛,Π^∗_ℎ𝜉^𝑛=𝐸₁+𝐸₂. (70)𝐸₁ and 𝐸₂ can be handled as 𝐷₁ and 𝐷₂ in Theorem 4. Thus, we have 𝐸₁=𝐵₂𝜃^𝑛;𝜃^𝑛,Π^∗_ℎ𝜉^𝑛−𝐵_2ℎ𝜃_ℎ^𝑛−1;Π_ℎ𝜃^𝑛,Π^∗_ℎ𝜉^𝑛≤𝐶𝜃^𝑛_∞𝜉^𝑛₁ℎ𝜃^𝑛₁+𝜃^𝑛−𝜃_ℎ^𝑛−1+ℎ𝜃^𝑛−𝜃_ℎ^𝑛−1₁≤𝐶𝜃^𝑛_∞𝜉^𝑛₁×ℎ𝜃^𝑛₁+𝜌^𝑛−1+𝜉^𝑛−1+𝜏^𝑡𝑛_𝑡𝑛−1𝜕𝜃𝜕𝑡²𝑑𝑡^1/2+ℎ𝜏^𝑡𝑛_𝑡𝑛−1𝜕𝜃𝜕𝑡²₁𝑑𝑡^1/2+𝜌^𝑛−1₁+𝜉^𝑛−1₁≤𝐶𝜌^𝑛−12+𝜉^𝑛−12+𝜉^𝑛2+𝜏^𝑡𝑛_𝑡𝑛−1𝜕𝜃𝜕𝑡²₁𝜉𝑑𝑡+𝜀^𝑛−12₁+𝜀𝜉^𝑛2₁,𝐸₂=𝐵_2ℎ𝜃_ℎ^𝑛−1;𝜉^𝑛,Π^∗_ℎ𝜉^𝑛=_{𝑃𝑖∈Ωℎ}𝜉^𝑛𝑃_𝑖⋅_{𝑗∈Λ𝑖Γ𝑖𝑗}𝜃_ℎ^𝑛−1⋅𝜈_𝑖𝑗×𝐻𝛽𝑑𝑠_𝑖𝑗𝜃_ℎ^𝑛−1𝜉^𝑛𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜃_ℎ^𝑛−1𝜉^𝑛𝑃_𝑗≤12_{𝐾∈𝑇ℎ𝑖,𝑗∈Λ𝐾}𝜉^𝑛𝑃_𝑖−𝜉^𝑛𝑃_{𝑗Γ𝑖𝑗∩𝐾}𝜉^𝑛−1⋅𝜈_𝑖𝑗×𝐻𝛽𝑑𝑠_𝑖𝑗𝜃_ℎ^𝑛−1𝜉^𝑛𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜃_ℎ^𝑛−1𝜉^𝑛𝑃_𝑗+12_{𝐾∈𝑇ℎ𝑖,𝑗∈Λ𝐾}𝜉^𝑛𝑃_𝑖−𝜉^𝑛𝑃_{𝑗Γ𝑖𝑗∩𝐾}Π_ℎ𝜃^𝑛−1⋅𝜈_𝑖𝑗×𝐻𝛽𝑑𝑠_𝑖𝑗𝜃_ℎ^𝑛−1𝜉𝑃_𝑖+1−𝐻𝛽_𝑖𝑗𝜃_ℎ^𝑛−1𝜉𝑃_𝑗≤𝐶𝜉^𝑛₁𝜉^𝑛−1+ℎ𝜉^𝑛−1₁𝜉^𝑛_∞+Π_ℎ𝜃^𝑛−1_∞𝜉^𝑛₁𝜉^𝑛≤𝐶𝜉^𝑛2+𝜉^𝑛−11logℎ^1/2𝜉^𝑛2₁+𝜀𝜉^𝑛2₁. (71)
Substituting the previously mentioned estimates into (64), we get 1𝜉2𝜏^𝑛2−𝜉^𝑛−12+𝛼𝜉^𝑛2₁≤𝐶ℎ²𝜃^𝑛2₂+𝜌^𝑛2₁+𝜌^𝑛−12₁+𝜉^𝑛2+𝜉^𝑛−12+𝐶𝜏^𝑡𝑛_𝑡𝑛−1𝜕²𝜃𝜕𝑡²²+𝜕𝜃𝜕𝑡²₁𝑑𝑡+𝐶𝜏^−1𝑡𝑛_𝑡𝑛−1𝜕𝜌𝜕𝑡²1𝑑𝑡+𝐶logℎ^1/2𝜉^𝑛−1𝜉^𝑛2₁+𝜀𝜉^𝑛−12₁+4𝜀𝜉^𝑛2₁. (72)
Multiplying (72) by 2𝜏 and summing over 1≤𝑛≤𝐿, we have 𝜉^𝐿2−𝜉⁰²+2𝜏𝛼^𝐿_𝑛=1𝜉^𝑛2₁≤𝐶ℎ²𝜃^𝑛2₂+𝜏^𝐿_𝑛=0𝜌^𝑛2₁+𝜏^𝐿_𝑛=1𝜉^𝑛2+𝐶𝜏^2𝑡𝐿₀𝜕²𝜃𝜕𝑡²²+𝜕𝜃𝜕𝑡²₁𝑑𝑡+𝐶^𝑡𝐿₀𝜕𝜌𝜕𝑡²𝑑𝑡+𝐶𝜏^𝐿_𝑛=11logℎ^1/2𝜉^𝑛−1𝜉^𝑛2₁+5𝜀^𝐿_𝑛=1𝜉^𝑛2₁. (73) By choosing ℎ and 𝜀 small enough and noting Lemma 2(ii), we have 𝜉^𝐿2+𝜏^𝐿_𝑛=1𝜉^𝑛2₁≤𝐶ℎ²𝜃^𝑛2₂+𝜏^𝐿_𝑛=0𝜌^𝑛2₁+𝜏^𝐿_𝑛=1𝜉^𝑛2+𝐶𝜏^2𝑡𝐿₀𝜕²𝜃𝜕𝑡²²+𝜕𝜃𝜕𝑡²₁𝑑𝑡+𝐶^𝑡𝐿₀𝜕𝜌𝜕𝑡²𝑑𝑡. (74)
Applying the Gronwall inequality and the interpolation theory, we deduce that 𝜉^𝐿2≤𝐶𝜏^2𝑇₀𝜕²𝜃𝜕𝑡²²+𝜕𝜃𝜕𝑡²₁𝑑𝑡+𝐶ℎ²𝜃₀²₂+^𝑇₀𝜕𝜃𝜕𝑡²₂𝑑𝑡. (75)
Now we prove the induction hypothesis (68). Noting that 𝜃⁰_ℎ=Π_ℎ𝜃₀, we know that 𝜉⁰=0. From (75) and the assumption 𝜏=𝑂(ℎ), we get that 1logℎ^1/2𝜉^𝐿1≤𝐶ℎlogℎ^1/2⟶0,ℎ⟶0. (76) Thus we know that (68) holds for any 1≤𝐿≤𝑁. Using triangular inequality and the interpolation theory completes the proof.

6. Numerical Example

In this section, we will show the affectivity of our method by numerical experiments. The exact solutions to problem (1) can be obtained by employing Cole-Hopf transformation. For Ω={(𝑥₁,𝑥₂)∶0≤𝑥₁,𝑥₂≤1}, we consider the following solutions: 3𝑢=4−141+exp𝜂−4𝑥₁+4𝑥₂,3−𝑡/32𝑣=4+141+exp𝜂−4𝑥₁+4𝑥₂,−𝑡/32 (77) where 𝜂=1/𝜁. We present numerical results for the 𝐿²-norm estimates of 𝑢−𝑢_ℎ and 𝑣−𝑣_ℎ. In Tables 1 and 2, we present the numerical results for 𝜁=1 and 𝜁=0.01, respectively. In all runs, we use the uniform mesh step ℎ=Δ𝑡 and choose the time 𝑡=1. As seen in these tables, in all cases the errors decrease by a factor of about two as ℎ decreases by the factor of two. This indicates that all 𝐿²-norm error estimates are of first-order convergence, which is consistent with our theoretical analysis.

Table 1 

Numerical results for 𝜁=1.

Table 2 

Numerical results for 𝜁=0.01.

When 𝜁=1.0 and 𝜁=0.01, the figures of the exact solutions 𝑢 and the numerical solutions 𝑢_ℎ at 𝑡=1 for ℎ=1/32 are given in Figures 1, 2, 3, and 4. In order to show that our method keeps stable when 𝜁 is smaller, we also give the comparison figures of exact solution 𝑢 and numerical solution 𝑢_ℎ for 𝜁=0.001 in Figures 5 and 6. The comparison figure of numerical solution by using finite volume element method (FVEM) without upwinding is given in Figure 7, which show that the approximation produces unacceptable nonphysical oscillations.

Figure 1 

The exact solution 𝑢 when 𝜁=1 at 𝑡=1.0.

Figure 2 

The numerical solution 𝑢ℎ when 𝜁=1 at 𝑡=1.0, for ℎ=1/32.

Figure 3 

The exact solution 𝑢 when 𝜁=0.01 at 𝑡=1.0.

Figure 4 

The numerical solution of 𝑢ℎ when 𝜁=0.01 at 𝑡=1.0, for ℎ=1/32.

Figure 5 

The exact solution 𝑢 at 𝜁=0.001 at 𝑡=1.0.

Figure 6 

The numerical solution 𝑢ℎ by FVEM with upwinding when 𝜁=0.001 at 𝑡=1.

Figure 7 

The numerical solution 𝑢ℎ by FVEM without upwinding when 𝜁=0.001 at 𝑡=1.

Acknowledgments

This paper is supported by The Natural Science Foundation of China (Grant nos. 10971254 and 11171193) and the Natural Science Foundation of Shandong Province (Grant nos. ZR2009AZ003 and ZR2011AM016).