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

Authors: Qing Yang
Journal: Abstract and Applied Analysis
 (2013)

---

## 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 inL2-norm. Numerical experiment is presented finally to validate the theoretical analysis.


---

## Body

## 1. Introduction

We consider the following two-dimensional Burgers equation [1–3]:
(1)(a)∂u∂t+u∂u∂x1+v∂u∂x2=ζΔu,x=(x1,x2)∈Ω,t∈J=(0,T],(b)∂v∂t+u∂v∂x1+v∂v∂x2=ζΔv,(x,t)∈Ω×J,(c)u(x,0)=ϕ(x),v(x,0)=ψ(x),x∈Ω,(d)u=g1,v=g2,(x,t)∈∂Ω×J,
for the unknown functions u and v in a bounded spatial domain Ω⊂ℝ2, over a time interval [0,T]. 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:(2)L2(Ω)={f:∫Ω‍|f|2dx<∞},∥f∥=[∫Ω‍|f|2dx]1/2,L∞(Ω)={f:esssupΩ|f|<∞},∥f∥∞=esssupΩ|f|,Hm(Ω)={f:∂|α|f∂xα∈L2(Ω),|α|≤m},m≥0,∥f∥m=[∑|α|≤m‍∥∂|α|f∂xα∥2]1/2,W∞m(Ω)={f:∂|α|f∂xα∈L∞(Ω),|α|≤m},m≥0.∥f∥m,∞=max|α|≤m∥∂|α|f∂xα∥L∞,
In particular, H0(Ω)=L2(Ω),  W∞0(Ω)=L∞(Ω). Let [a,b]⊂[0,T] and let X be any of the spaces just defined. If f(x,t) represents functions on Ω×[a,b], we set
(3)Hm(a,b;X)={f:∫ab‍∥∂αf∂tα(·,t)∥X2dt<∞,α≤m},∥f∥Hm(a,b;X)=[∑α=0m‍∫ab‍∥∂αf∂tα(·,t)∥X2dt]1/2,m≥0,W∞m(a,b;X)={f:esssup[a,b]∥∂αf∂tα(·,t)∥X<∞,α≤m},∥f∥W∞m(a,b;X)=max0≤α≤messsup[a,b]∥∂αf∂tα(·,t)∥X,m≥0,L2(a,b;X)=H0(a,b;X),L∞(a,b;X)=W∞0(a,b;X).
If [a,b]=[0,T], we drop it from the notation. We also drop Ω; thus, we write L∞(W∞1) for L∞(0,T;W∞1(Ω)).Ifw=(w1,w2) is a vector function, we say that w∈X2 if w1∈X and w2∈X.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 L2-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 byC and Ci(i=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 w=(w1,w2):ℝ2→ℝ2 is a matrix, and the divergence of a matrix function A=(aij)1≤i,j≤2:ℝ2→ℝ2×2 is a vector
(4)∇w=(∂wi∂xj)1≤i,j≤2,∇·A=(∑j=12‍∂a1j∂xj,∑j=12‍∂a2j∂xj).
Consequently, we have for a vector function w=(w1,w2)(5)Δw=∇·∇w=(Δw1,Δw2).
Let θ=(u,v), θ0(x)=(u0(x),v0(x)), and let g=(g1,g2); then the system (1) can be written as the following vector form:
(6)(a)∂θ∂t+θ·∇θ-ζΔθ=0,(x,t)∈Ω×J,(b)θ(x,0)=θ0(x),x∈Ω,(c)θ(x,t)=g,(x,t)∈∂Ω×J,
where
(7)θ·∇θ=u∂θ∂x1+v∂θ∂x2.Let𝒯h={K} be a triangulation of the domain Ω, and as usual, we assume the triangles K to be shape regular. Denote by Ω-h={Pi} the set of the vertices of all the triangles K, and let Ωh=Ω-h∖∂Ω. For a given triangulation 𝒯h, we construct a dual mesh 𝒯h* whose elements are called control volumes. Each triangle K∈𝒯h can be divided into three subdomains by connecting an inner point of the triangle to the midpoints of the three edges. Around each Pi∈Ωh, we associate a control volume Ki*=KPi*, which consists of the union of subregions having Pi as a vertex. For a vertex Pi∈∂Ω, we can define its control volume in a similar way. Then we define the dual partition 𝒯h*={KPi*,Pi∈Ω-h} 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 g1=g2=0. We define the following finite dimensional spaces:
(8)Ph={ωh∈H01(Ω),ωh|K∈𝒫1(K),K∈Th},Yh={φh∈L2(Ω),φh|Ki*∈𝒫0(K*),Pi∈Ωh;φh|Ki*=0,Pi∈∂Ωφh∈L2(Ω),φh|Ki*∈𝒫0(K*)},Uh=Ph2,Vh=Yh2,
where 𝒫l(V)(l=0,1) denotes the set of polynomials on V with a degree of not more than l.Multiplying (6a) by test functionz∈Vh and integrating by parts yield
(9)(∂θ∂t,z)+B1(θ;θ,z)+B2(θ;θ,z)+A(θ,z)=0,∀z∈Vh,
where
(10)B1(φ;w,z)=-∑Pi∈Ωh‍z(Pi)·∫Ki*‍(∇·w)φdx,B2(φ;w,z)=∑Pi∈Ωh‍z(Pi)·∫∂Ki*‍(φ·ν)wds,A(w,z)=∑Pi∈Ωh‍z(Pi)·∫∂Ki*‍ζ(ν·∇w)ds,
here ν is the unit outward normal vector of ∂Ki*.Now we approximateB2(φ;w,z) by using the upwind technique.LetΛi={j:Pj is adjoint with Pi}. Assuming that j∈Λi, let Γij=∂Ki*∩∂Kj* and γij is the length of Γij. Denote by νij the unit outward normal vector of Γij when Γij is regarded as the boundary of Ki*. Define
(11)βij(φ)=∫Γij‍φ·νijds.
Let
(12)βij+(φ)=max(βij(φ),0),βij-(φ)=max(-βij(φ),0),∫∂Ki*‍(φ·ν)wds≈∑j∈Λi‍{βij+(φ)w(Pi)-βij-(φ)w(Pj)}.
The upwind discretization of the nonlinear term B2(φ;w,z) is defined by the form
(13)B2h(φ;w,z)=∑Pi∈Ωh‍∑j∈Λi‍{βij+(φ)w(Pi)-βij-(φ)w(Pj)}·z(Pi).
Using the heaviside function
(14)H(r)={1,r≥0,0,r<0,
we can write B2h(φ;w,z) as
(15)B2h(φ;w,z)=∑Pi∈Ωh‍∑j∈Λi‍βij(φ)×{H(βij(φ))w(Pi)+(1-H(βij(φ)))w(Pj)}·z(Pi).Introduce the interpolation operatorsΠh:H01(Ω)→Ph and Πh*:Ph→Yh, respectively. For w=(w1,w2), define Πhw=(Πhw1,Πhw2) and Πh*w=(Πh*w1,Πh*w2). Assuming that w∈H2(Ω)2, we can easily get the following interpolation estimates:
(16)∥w-Πhw∥s≤h2-s∥w∥2,s=0,1.The semi-discrete upwind finite volume scheme of (6) is as follows: find θh:[0,T]→Uh such that
(17)(∂θh∂t,Πh*zh)+B1(θh;θh,Πh*zh)+B2h(θh;θh,Πh*zh)+A(θh,Πh*zh)=0,∀zh∈Uh,θh(x,0)=θ0h(x),
where θ0h(x) is the interpolation projection of θ0, that is, θ0h(x)=Πhθ0.Partition[0,T] into 0=t0<t1<⋯<tN=T, with τn=tn-tn-1. Our analysis is valid for variable time steps, but we drop the superscript from τ for convenience. For functions f on Ω×J, we write fn(x) for f(x,tn). By approximating ∂θh/∂t at the time t=tn with the backward difference ∂tθhn=(θhn-θhn-1)/τ, we define the fully discrete upwind finite volume scheme for (6) as follows: find θhn∈Uh, such that
(18)(∂tθhn,Πh*zh)+B1h(θhn-1;θhn,Πh*zh)+B2h(θhn-1;θhn,Πh*zh)+A(θhn,Πh*zh)=0,n≥1,∀zh∈Uh,θh0=θ0h.


## 3. Some Lemmas

Now we present several lemmas. Letwh=(w1,w2)∈Uh,w~h=(w~1,w~2)∈Uh.Lemma 1.
(
i
)   Πh* is a self-adjoint operator, that is,
(19)(wh,Πh*w~h)=(w~h,Πh*wh),∀wh,w~h∈Uh.
(
ii
) Let |∥wh∥|=(wh,Πh*wh)1/2. Then, for some positive constants C1 and C2 that are independent of h,
(20)C1∥wh∥≤|∥wh∥|≤C2∥wh∥,∀wh∈Uh.Proof.
It is easy to know that(21)(wh,Πh*w~h)=(w1,Πh*w~1)+(w2,Πh*w~2),|∥wh∥|2=(wh,Πh*wh)=(w1,Πh*w1)+(w2,Πh*w2).
From [4] we know that for w1, w2, w~1, w~2∈Ph,
(22)(w1,Πh*w~1)=(w~1,Πh*w1),(w2,Πh*w~2)=(w~2,Πh*w2),C12∥w1∥2≤(w1,Πh*w1)≤C22∥w1∥2,C12∥w2∥2≤(w2,Πh*w2)≤C22∥w2∥2,
where C1 and C2 are some positive constants that are independent of h. Thus we obtain (19) and (20) immediately.Lemma 2.
For the bilinear formA(·,Πh*·), one has the following conclusions:(i)
Forwh,w~h∈Uh, one has
(23)A(wh,Πh*w~h)=A(w~h,Πh*wh).(ii)
There exists a positive constantC such that
(24)|A(w-Πhw,Πh*w~h)|≤Ch∥w∥2∥w~h∥1,∀w∈H2(Ω)2,w~h∈Uh.(iii)
There exists a positive constantα such that
(25)A(wh,Πh*wh)≥α∥wh∥12,∀wh∈Uh.Proof.
Forϕ,ψ∈Ph, define the bilinear form
(26)a(ϕ,Πh*ψ)=-∑Pi∈Ωh‍ψ(Pi)∫∂Ki*‍ζ∇ϕ·νds.
Then, we get
(27)A(wh,Πh*w~h)=-∑Pi∈Ωh‍w~h(Pi)·∫∂Ki*‍ζ(ν·∇wh)ds=-∑Pi∈Ωh‍w~1(Pi)∫∂K*‍ζ∇w1·νds-∑Pi∈Ωh‍w~2(Pi)∫∂Ki*‍ζ∇w2·νds=a(w1,Πh*w~1)+a(w2,Πh*w~2),A(w-Πhw,Πh*w~h)=a(w1-Πhw1,Πh*w~1)+a(w2-Πhw2,Πh*w~2).
By combining the above results and the corresponding conclusions for a(·,Πh*·) in [4], we can obtain (23)–(25).Lemma 3.
Forφ∈(W∞0(Ω))2,θ∈(H01(Ω))2,φh∈Uh, and zh∈Uh, one has
(28)|B2(φ;θ,Πh*zh)-B2h(φh;Πhθ,Πh*zh)|≤|zh|1{h∥φ∥∞|θ|1+∥θ∥∞(∥φ-φh∥+h|φ-φh|1)}.Proof.
First we have(29)B2(φ;θ,Πh*zh)-B2h(φh;Πhθ,Πh*zh)=B2(φ;θ,Πh*zh)-B2h(φh;θ,Πh*zh)+B2h(φh;θ,Πh*zh)-B2h(φh;Πhθ,Πh*zh).
Noting that θ(Pi)=Πhθ(Pi), Pi∈Ωh, we can easily deduce
(30)B2h(φh;θ,Πh*zh)-B2h(φh;Πhθ,Πh*zh)=0
by the definition of B2h(·;·,Πh*·). Now we only need to bound
(31)B2(φ;θ,Πh*zh)-B2h(φh;θ,Πh*zh)=∑Pi∈Ωh‍zh(Pi)·∑j∈Λi‍∫Γij‍(φ·νij)θds-∑Pi∈Ωh‍zh(Pi)·∑j∈Λi‍βij(φh)×[H(βij(φh))θ(Pi)+(1-H(βij(φh)))θ(Pj)]=∑Pi∈Ωh‍zh(Pi)·∑j∈Λi‍∫Γij‍(φ·νij)×{θ-[H(βij(φh))θ(Pi)+(1-H(βij(φh)))θ(Pj)]}ds+∑Pi∈Ωh‍zh(Pi)·∑j∈Λi‍∫Γij‍(φ-φh)·νijds×[H(βij(φh))θ(Pi)+(1-H(βij(φh)))θ(Pj)].
We denote the last two terms on the right-hand side of (31) by I1 and I2, respectively. We now turn to analyze the two terms. Noting that βij=-βji, we rewrite I1 as
(32)I1=12∑K∈𝒯h‍∑i,j∈ΛK‍[zh(Pi)-zh(Pj)]·∫Γij∩K‍(φ·νij)[H(βij(φh))(θ-θ(Pi))∫Γij∩K‍(φ·νij)+(1-H(βij(φh)))(θ-θ(Pj))]ds,
here ΛK is the set of vertex of K. From the Taylor’s Formula and the linear property of zh=(z1,z2), we obtain that
(33)|zh(Pi)-zh(Pj)|2≤h2(|∇z1|2+|∇z2|2)=h2|z|1,K2.
Applying the trace inequality, we get(34)∫Γij∩K‍|θ-θ(Pi)|ds≤Ch1/2{∫Γij∩K‍|θ-θ(Pi)|2ds}1/2≤Ch1/2{∫Γij∩K‍[|u-u(Pi)|2+|v-v(Pi)|2]ds}1/2≤Ch{(h-1|u-u(Pi)|0,K+|u-u(Pi)|1,K)2+(h-1|v-v(Pi)|0,K+|v-v(Pi)|1,K)2}1/2≤Ch(|u|1,K2+|v|1,K2)1/2=Ch|θ|1,K.
Similarly, we can deduce that
(35)∫Γij∩K‍|θ-θ(Pj)|ds≤Ch|θ|1,K.
We conclude that
(36)|I1|≤Ch∥φ∥∞|zh|1|θ|1.
The similar argument yields the estimate(37)|I2|=|12∑K∈𝒯h‍∑i,j∈ΛK‍∫Γij∩K‍[(φ-φh)·n]ds×[H(βij(φh))θ(Pi)+(1-H(βij(φh)))θ(Pj)]·[zh(Pi)-zh(Pj)]12∑K∈𝒯h‍∑i,j∈ΛK‍∫Γij∩K‍|≤C∥θ∥∞|zh|1(∥φ-φh∥+h|φ-φh|1).
Substituting the estimates (36) and (37) into (31), we obtain
(38)|B2(φ;θ,Πh*zh)-B2h(φh;θ,Πh*zh)|≤C|zh|1{h∥φ∥∞|θ|1+∥θ∥∞×(∥φ-φh∥+h|φ-φh|1)}.
This yields the desired result immediately.


## 4. Error Bounds for Semi-Discrete Scheme

Theorem 4.
Assume thatθ and θh are solutions to (6) and (17), respectively. also assumes that θ is regular enough. Then there exists a positive constant C such that
(39)∥θ-θh∥≤Ch,
where C depends on principally ∥θ0∥2, ∥θ∥L∞((W∞1)2), and ∥θ∥H1((H2)2).Proof.
We derive the following error equation from (6) and (17):
(40)(∂θ∂t-∂θh∂t,Πh*zh)+B1(θ;θ,Πh*zh)-B1(θh;θh,Πh*zh)+B2(θ;θ,Πh*zh)-B2h(θh;θh,Πh*zh)+A(θ-θh,Πh*zh)=0.
Let ρ=θ-Πhθ,ξ=Πhθ-θh. We rewrite the previously mentioned equation as
(41)(∂ξ∂t,Πh*zh)+B1(θ;θ,Πh*zh)-B1(θh;θh,Πh*zh)+B2(θ;θ,Πh*zh)-B2h(θh;θh,Πh*zh)+A(ξ,Πh*zh)=-(∂ρ∂t,Πh*zh)-A(ρ,Πh*zh).
We choose zh=ξ in (41) to get
(42)(∂ξ∂t,Πh*ξ)+A(ξ,Πh*ξ)=-(∂ρ∂t,Πh*ξ)-A(ρ,Πh*ξ)-[B1(θ;θ,Πh*ξ)-B1(θh;θh,Πh*ξ)]-[B2(θ;θ,Πh*ξ)-B2h(θh;θh,Πh*ξ)].
Using Lemmas 1, 2 and Young’s inequality, we have
(43)12ddt|∥ξ∥|2+α∥ξ∥12≤C(∥∂ρ∂t∥2+∥ξ∥2+h2∥θ∥22)+ε∥ξ∥12+|B1(θ;θ,Πh*ξ)-B1(θh;θh,Πh*ξ)|+|B2(θ;θ,Πh*ξ)-B2h(θh;θh,Πh*ξ)|.
Now we bound the last two terms on the right-hand side of (43). We need the following induction hypothesis:
(44)(log1h)1/2∥ξ∥(s)⟶0,h⟶0,0≤s<t,0<t≤T.
We know from [13] that
(45)∥ϕ∥∞≤C(log1h)1/2∥ϕ∥1,∀ϕ∈Ph.
This implies that
(46)∥φ∥∞≤C(log1h)1/2∥φ∥1,∀φ∈Uh.
Also we have the following inverse inequality:
(47)∥φ∥1≤Ch-1∥φ∥,∀φ∈Uh.
Using (46), we get
(48)|B1(θ;θ,Πh*ξ)-B1(θh;θh,Πh*ξ)|≤∑Pi∈Ωh‍|ξ(Pi)|∫Ki*‍|(∇·θ)θ-(∇·θh)θh|dx≤∑Pi∈Ωh‍|ξ(Pi)|∫Ki*‍{|(∇·θ)||θ-θh|+|∇·θ-∇·θh|≤∑Pi∈Ωh‍|ξ(Pi)|∫Ki*‍×(|ξ|+|Πhθ|)}dx≤C{∥∇·θ∥∞∥θ-θh∥∥ξ∥+∥∇·(θ-θh)∥≤C×∥ξ∥(∥ξ∥∞+∥θ∥∞)}≤C{∥ρ∥1+∥ξ∥1+(log1h)1/2∥ξ∥1≤C×(∥ρ∥1+∥ξ∥1)(log1h)1/2}∥ξ∥≤C{∥ρ∥12+∥ξ∥2+∥ξ∥(log1h)1/2∥ξ∥12}+ε∥ξ∥12.
Next, we write
(49)|B2(θ;θ,Πh*ξ)-B2h(θh;θh,Πh*ξ)|≤|B2(θ;θ,Πh*ξ)-B2h(θh;Πhθ,Πh*ξ)|+|B2h(θh;ξ,Πh*ξ)|=D1+D2.
By Choosing φ=θ,φh=θh, and zh=ξ in Lemma 3, using (47) and the Young’s inequality, we can obtain
(50)D1≤C|ξ|1{h∥θ∥∞|θ|1+∥θ∥∞D1≤C|ξ|1×(∥θ-θh∥+h|θ-θh|1)}≤C{h2|θ|12+∥ρ∥12+∥ξ∥2}+ε∥ξ∥12.
By an argument like (36) and then by (46) and (47), we have
(51)D2=|∑Pi∈Ωh‍ξ(Pi)·∑j∈Λi‍∫Γij‍(θh·νij)ds×[H(βij)ξ(Pi)+(1-H(βij))ξ(Pj)]∑Pi∈Ωh‍|≤12∑K∈Th‍∑i,j∈ΛK‍|ξ(Pi)-ξ(Pj)|×|H(βij)ξ(Pi)+(1-H(βij))ξ(Pj)|×∫Γij∩K‍|ξ·νij|ds+12∑K∈Th‍∑i,j∈ΛK‍|ξ(Pi)-ξ(Pj)|×|H(βij)ξ(Pi)+(1-H(βij))ξ(Pj)|×∫Γij∩K‍|Πhθ·νij|ds≤C{∥ξ∥1(∥ξ∥+h|ξ|1)∥ξ∥∞+∥Πhθ∥∞∥ξ∥1∥ξ∥}≤C{∥ξ∥2+∥ξ∥(log1h)1/2∥ξ∥12}+ε∥ξ∥12.
Substituting (50) and (51) into (49), we get
(52)|B2(θ;θ,Πh*ξ)-B2h(θh;θh,Πh*ξ)|≤C{h2|θ|12+∥ρ∥12+∥ξ∥2}+C∥ξ∥(log1h)1/2∥ξ∥12+2ε∥ξ∥12.
Make (43), (48), and (52) together to obtain
(53)12ddt|∥ξ∥|2+α∥ξ∥12≤C{h2∥θ∥22+∥∂ρ∂t∥2+∥ρ∥12+∥ξ∥2}+C∥ξ∥(log1h)1/2∥ξ∥12+4ε∥ξ∥12.
Integrating the previously mentioned equation from 0 to t and noting (44), we obtain that
(54)12{|∥ξ∥|2(t)-|∥ξ∥|2(0)}+α2∫0t‍∥ξ∥12dτ≤{h2∥θ∥22+∫0t‍∥∂ρ∂t∥2dτ+∫0t‍∥ρ∥12dτ+∫0t‍∥ξ∥2dτ{h2∥θ∥22+∫0t‍∥∂ρ∂t∥2dτ+∫0t‍∥ρ∥12dτ}
for sufficiently small h and ε. By using Lemma 2(ii) and the Gronwall’s inequality, we have that
(55)∥ξ∥2(t)+∫0t‍∥ξ∥12dτ≤C{h2∥θ∥22+∫0t‍∥∂ρ∂t∥2dτ+∫0t‍∥ρ∥12dτ}.
It follows from the interpolation estimates that
(56)∥ξ∥2(t)+∫0t‍∥ξ∥12dτ≤C{h2(∥θ0∥22+∫0T‍∥∂θ∂t∥22dτ)+h4∫0T‍∥∂θ∂t∥22dτ+h2∫0T‍∥θ∥22dτ}.
Now we prove the induction hypothesis (44). Noting that ∥ξ∥(0)=0, we know that (44) holds obviously for t=0. It follows from (56) that
(57)(log1h)1/2∥ξ∥(t)≤Ch(log1h)1/2⟶0,h⟶0.
Then (44) holds for any t∈[0,T].
By (16), we have
(58)∥ρ∥≤Ch2∥θ∥2≤Ch2{∥θ0∥2+∫0t‍∥θ∥2dτ}.
From the triangle inequality, we obtain
(59)∥θ-θh∥≤Ch,
where C depends on ∥θ0∥2, ∥θ∥H1((W∞1)2), and ∥θ∥H1((H2)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 τ=O(h). Then the error of the approximation (18) of (6) satisfies
(60)max0≤n≤T/τ∥θn-θhn∥≤C{h+τ},
where C depends on ∥θ0∥2, ∥θ∥L∞((W∞1)2), ∥θ∥H1((H2)2), and ∥θ∥H2((L2)2).Proof.
Subtract (18) from (9) to obtain that
(61)(∂θn∂t-∂tθhn,Πh*zh)+B1(θn;θn,Πh*zh)-B1(θhn-1;θhn,Πh*zh)+B2(θn;θn,Πh*zh)-B2h(θhn-1;θhn,Πh*zh)+A(θn-θhn,Πh*zh)=0.
Choose zh=ξn to obtain that
(62)(∂tξn,Πh*ξn)+A(ξn,Πh*ξn)=-(∂θn∂t-∂tθn,Πh*ξn)-(∂tρn,Πh*ξn)-A(ρn,Πh*ξn)-[B1(θn;θn,Πh*ξn)-B1(θhn-1;θhn,Πh*ξn)]-[B2(θn;θn,Πh*ξn)-B2h(θhn-1;θhn,Πh*ξn)].
For the left-hand side of (62), from Lemmas 1 and 2, we have
(63)(∂tξn,Πh*ξn)=1τ(ξn-ξn-1,Πh*ξn)=12τ(ξn-ξn-1,Πh*(ξn+ξn-1))+12τ(ξn-ξn-1,Πh*(ξn-ξn-1))≥12τ(|∥ξn∥|2-|∥ξn-1∥|2),A(ξn,Πh*ξn)≥α∥ξn∥12.
We denote terms on the right-hand side of (62) by T1,…,T5. Then, (62) can be rewritten as
(64)12τ(|∥ξn∥|2-|∥ξn-1∥|2)+α∥ξn∥12≤T1+⋯+T5.
Now we estimates the termsT1,…,T5 one by one. From the Taylor’s formula, we have
(65)∂θn∂t-∂tθn=1τ∫tn-1tn‍(t-tn-1)∂2θ∂t2dt.
It follows that
(66)|T1|≤C∥∂θn∂t-∂tθn∥∥ξn∥≤C{τ∫tn-1tn‍∥∂2θ∂t2∥2dt+∥ξn∥2}.
For the next two terms, we have
(67)|T2|≤C∥∂tρn∥∥ξn∥≤C{τ-1∫tn-1tn‍∥∂ρ∂t∥2dt+∥ξn∥2},|T3|≤Ch∥θn∥2∥ξn∥1≤Ch2∥θn∥22+ε∥ξn∥12.
We make the following induction hypothesis:(68)∥ξn-1∥(log1h)1/2⟶0,h⟶0,1≤n≤L.
ForT4, using the similar argument as (48) and noting (68), we deduce that
(69)|T4|=|B1(θn;θn,Πh*ξn)-B1(θhn-1;θhn,Πh*ξn)|≤∑Pi∈Ωh‍|ξn(Pi)|∫Ki*‍|(∇·θn)θn-(∇·θhn)θhn-1|dx≤∑Pi∈Ωh‍|ξn(Pi)|∫Ki*{|∇·θn||θn-θhn-1|+|∇·θn-∇·θhn|‍≤∑Pi∈Ωh‍|ξn(Pi)|×(|ξn-1|+|Πhθn-1|)|θn-θhn-1|}dx≤C{∥∇·θn∥∞∥θn-θhn-1∥∥ξn∥+∥∇·(θn-θhn)∥×(∥ξn∥∞∥ξn-1∥+∥Πhθn-1∥∞∥ξn∥)∥θn-θhn-1∥∥∇·θn∥∞}≤C{[(τ∫tn-1tn‍∥∂θ∂t∥2dt)1/2+∥ρn-1∥+∥ξn-1∥]×∥∇·θn∥∞∥ξn∥+(∥ρn∥1+∥ξn∥1)×[(log1h)1/2∥ξn-1∥∥ξn∥1+∥θn-1∥∞∥ξn∥][(τ∫tn-1tn‍∥∂θ∂t∥2dt)1/2+∥ρn-1∥+∥ξn-1∥]}≤C{τ∫tn-1tn‍∥∂θ∂t∥2dt+∥ρn-1∥2+∥ρn∥12+∥ξn-1∥2+∥ξn∥2∫tn-1tn‍∥∂θ∂t∥2}+C(log1h)1/2∥ξn-1∥∥ξn∥12+ε∥ξn∥12.
Now, we write(70)|T5|=|B2(θn;θn,Πh*ξn)-B2h(θhn-1;θhn,Πh*ξn)|≤|B2(θn;θn,Πh*ξn)-B2h(θhn-1;Πhθn,Πh*ξn)|+|B2h(θhn-1;ξn,Πh*ξn)|=E1+E2.E1 and E2 can be handled as D1 and D2 in Theorem 4. Thus, we have
(71)E1=|B2(θn;θn,Πh*ξn)-B2h(θhn-1;Πhθn,Πh*ξn)|≤C∥θn∥∞|ξn|1[h|θn|1+(∥θn-θhn-1∥+h|θn-θhn-1|1)]≤C∥θn∥∞|ξn|1×{(τ∫tn-1tn‍∥∂θ∂t∥2dt)1/2h|θn|1+∥ρn-1∥+∥ξn-1∥+(τ∫tn-1tn‍∥∂θ∂t∥2dt)1/2+h[(τ∫tn-1tn∥∂θ∂t∥12‍dt)1/2+∥ρn-1∥1+∥ξn-1∥1]}≤C{∥ρn-1∥2+∥ξn-1∥2+∥ξn∥2+τ∫tn-1tn‍∥∂θ∂t∥12dt}+ε∥ξn-1∥12+ε∥ξn∥12,E2=|B2h(θhn-1;ξn,Πh*ξn)|=|∑Pi∈Ωh‍ξn(Pi)·∑j∈Λi‍∫Γij‍(θhn-1·νij)ds×[H(βij(θhn-1))ξn(Pi)+(1-H(βij(θhn-1)))ξn(Pj)(1-H(βij(θhn-1)))]∑Pi∈Ωh‍|≤12∑K∈Th‍∑i,j∈ΛK‍|ξn(Pi)-ξn(Pj)|∫Γij∩K‍|ξn-1·νij|ds×|H(βij(θhn-1))ξn(Pi)+(1-H(βij(θhn-1)))ξn(Pj)|+12∑K∈Th‍∑i,j∈ΛK‍|ξn(Pi)-ξn(Pj)|∫Γij∩K‍|Πhθn-1·νij|ds×|H(βij(θhn-1))ξ(Pi)+(1-H(βij(θhn-1)))ξ(Pj)|≤C{∥ξn∥1(∥ξn-1∥+h|ξn-1|1)∥ξn∥∞+∥Πhθn-1∥∞∥ξn∥1∥ξn∥}≤C{∥ξn∥2+∥ξn-1∥(log1h)1/2∥ξn∥12}+ε∥ξn∥12.
Substituting the previously mentioned estimates into (64), we get
(72)12τ(|∥ξn∥|2-|∥ξn-1∥|2)+α∥ξn∥12≤C{h2∥θn∥22+∥ρn∥12+∥ρn-1∥12+∥ξn∥2+∥ξn-1∥2}+Cτ∫tn-1tn‍(∥∂2θ∂t2∥2+∥∂θ∂t∥12)dt+Cτ-1∫tn-1tn‍∥∂ρ∂t∥2dt+C(log1h)1/2∥ξn-1∥∥ξn∥12+ε∥ξn-1∥12+4ε∥ξn∥12.
Multiplying (72) by 2τ and summing over 1≤n≤L, we have
(73)|∥ξL∥|2-|∥ξ0∥|2+2τα∑n=1L‍∥ξn∥12≤C{h2∥θn∥22+τ∑n=0L‍∥ρn∥12+τ∑n=1L‍∥ξn∥2}+Cτ2∫0tL‍(∥∂2θ∂t2∥2+∥∂θ∂t∥12)dt+C∫0tL‍∥∂ρ∂t∥2dt+Cτ∑n=1L‍(log1h)1/2∥ξn-1∥∥ξn∥12+5ε∑n=1L‍∥ξn∥12.
By choosing h and ε small enough and noting Lemma 2(ii), we have
(74)∥ξL∥2+τ∑n=1L‍∥ξn∥12≤C{h2∥θn∥22+τ∑n=0L‍∥ρn∥12}+τ∑n=1L‍∥ξn∥2+Cτ2∫0tL‍(∥∂2θ∂t2∥2+∥∂θ∂t∥12)dt+C∫0tL‍∥∂ρ∂t∥2dt.
Applying the Gronwall inequality and the interpolation theory, we deduce that(75)∥ξL∥2≤Cτ2∫0T‍(∥∂2θ∂t2∥2+∥∂θ∂t∥12)dt+Ch2{∥θ0∥22+∫0T‍∥∂θ∂t∥22dt}.
Now we prove the induction hypothesis (68). Noting that θh0=Πhθ0, we know that ξ0=0. From (75) and the assumption τ=O(h), we get that
(76)(log1h)1/2∥ξL∥≤Ch(log1h)1/2⟶0,h⟶0.
Thus we know that (68) holds for any 1≤L≤N. 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 Ω={(x1,x2):0≤x1,x2≤1}, we consider the following solutions:
(77)u=34-14(1+exp(η(-4x1+4x2-t)/32)),v=34+14(1+exp(η(-4x1+4x2-t)/32)),
where η=1/ζ. We present numerical results for the L2-norm estimates of u-uh and v-vh. 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 h=Δt and choose the time t=1. As seen in these tables, in all cases the errors decrease by a factor of about two as h decreases by the factor of two. This indicates that all L2-norm error estimates are of first-order convergence, which is consistent with our theoretical analysis.Table 1
Numerical results forζ=1.
h
1
/
8
1
/
16
1
/
32
1
/
64
∥
u
-
u
h
∥
h
1.18416
e
-
007
5.33942
e
-
008
2.52582
e
-
008
1.22755
e
-
008
Rate
1.15
1.08
1.05
∥
v
-
v
h
∥
h
6.73307
e
-
008
2.36565
e
-
008
9.49449
e
-
009
4.21298
e
-
009
Rate
1.51
1.32
1.17Table 2
Numerical results forζ=0.01.
h
1
/
8
1
/
16
1
/
32
1
/
64
∥
u
-
u
h
∥
h
7.57108
e
-
003
3.81036
e
-
003
1.92640
e
-
003
9.71683
e
-
004
Rate
0.99
0.98
0.99
∥
v
-
v
h
∥
h
7.57108
e
-
003
3.81036
e
-
003
1.92640
e
-
003
9.71683
e
-
004
Rate
0.99
0.98
0.99Whenζ=1.0 and ζ=0.01, the figures of the exact solutions u and the numerical solutions uh at t=1 for h=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 u and numerical solution uh 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 solutionu when ζ=1 at t=1.0.Figure 2
The numerical solutionuh when ζ=1 at t=1.0, for h=1/32.Figure 3
The exact solutionu when ζ=0.01 at t=1.0.Figure 4
The numerical solution ofuh when ζ=0.01 at t=1.0, for h=1/32.Figure 5
The exact solutionu at ζ=0.001 at t=1.0.Figure 6
The numerical solutionuh by FVEM with upwinding when ζ=0.001 at t=1.Figure 7
The numerical solutionu~h by FVEM without upwinding when ζ=0.001 at t=1.



---

*Source: 351619-2013-03-13.xml*