FINITE ELEMENT SUPERCONVERGENCE ON SHISHKIN MESH ...

MATHEMATICS OF COMPUTATION Volume 72, Number 243, Pages 1147–1177 S 0025-5718(03)01486-8 Article electronically published on February 3, 2003

FINITE ELEMENT SUPERCONVERGENCE ON SHISHKIN MESH FOR 2-D CONVECTION-DIFFUSION PROBLEMS ZHIMIN ZHANG

Abstract. In this work, the bilinear finite element method on a Shishkin mesh for convection-diffusion problems is analyzed in the two-dimensional setting. A superconvergence rate O(N −2 ln2 N + N −1.5 ln N ) in a discrete -weighted energy norm is established under certain regularity assumptions. This convergence rate is uniformly valid with respect to the singular perturbation parameter . Numerical tests indicate that the rate O(N −2 ln2 N ) is sharp for the boundary layer terms. As a by-product, an -uniform convergence of the same order is obtained for the L2 -norm. Furthermore, under the same regularity assumption, an -uniform convergence of order N −3/2 ln5/2 N + N −1 ln1/2 N in the L∞ norm is proved for some mesh points in the boundary layer region.

1. Introduction There has been extensive research in numerical solutions of singular perturbation problems because of the practical importance of these problems (for example, the Navier-Stokes equations at high Reynolds number). One of the typical behaviors of singularly perturbed problems is the boundary layer phenomenon: the solution varies rapidly within very thin layer regions near the boundary. Most of the traditional numerical methods fail to catch the rapid change of the solution in boundary layers, and this failure in turn pollutes the numerical approximation on the whole domain. See [18] and [22]. Many methods have been developed to overcome the numerical difficulty caused by boundary layers. The reader is referred to three 1996 books [13, 14, 16] for the significant progress that has been made in this field, and articles [2, 4, 7, 8, 11, 12, 15, 18, 19, 20, 21, 24, 25] for more information. A realistic approach in practice may be starting with a certain up-winding scheme, such as the streamline-diffusion method, followed by an adaptive procedure to refine the mesh, eventually resolving the boundary layer, and maybe locating some possible internal layers. Then a question arises naturally: Is there any superconvergence phenomenon when the boundary layer is successfully resolved? The current work intends to answer this question for a specific situation. Received by the editor July 19, 2000 and, in revised form, December 10, 2001. 2000 Mathematics Subject Classification. Primary 65N30, 65N15. Key words and phrases. Convection, diffusion, singularly perturbed, boundary layer, Shishkin mesh, finite element method. This research was partially supported by the National Science Foundation grants DMS0074301, DMS-0079743, and INT-0196139. c

2003 American Mathematical Society

1147

License or copyright restrictions may apply to redistribution; see http://www.ams.org/journal-terms-of-use

1148

ZHIMIN ZHANG

We shall analyze the standard finite element method combined with one kind of local refinement strategy, namely, the Shishkin mesh. Roughly speaking, the Shishkin mesh is a piecewise uniform mesh with an anisotropic mesh of high ratio in the boundary layer region. The analysis in this paper shows that superconvergence is uniformly valid with respect to the singular perturbation parameter for the bilinear finite element method with the Shishkin mesh for our model problem. This finding is consistent with the symmetry theory [17] in the finite element superconvergence, since for a piecewise uniform mesh there are indeed many symmetries. However, we are not able to apply the symmetry theory directly to convectiondiffusion equations because of the use of highly anisotropic meshes. For general theory and new developments of finite element superconvergence, the reader is referred to the recent books [1], [9], [23], and the conference proceedings [6]. Recently, Li and Wheeler have obtained a superconvergence result for the lowest Raviart-Thomas rectangular element in approximating singularly perturbed reaction-diffusion equations in a mixed formulation [8]. By a local postprocessing, the authors are able to prove an O(N −2 ) convergence rate for the gradient. However, we have not seen any superconvergence result for convection-diffusion equations (which is more difficult) in the displacement formulation. In the current work, we consider the standard finite element method for a convection-diffusion model problem, (1.1)

− ∆u + β~ · ∇u + cu =

f

in Ω = (0, 1) × (0, 1),

(1.2)

u =

0

on ∂Ω,

~ y) = (β1 (x, y), β2 (x, y)) ≥ (α, α) > (0, 0), where is a small positive number, β(x, ¯ c(x, y) ≥ 0 for all (x, y) ∈ Ω, and 1 ~ y) ≥ c0 > 0 c(x, y) − divβ(x, 2 ~ c, and f are sufficiently smooth. These with constants α and c0 . We assume that β, assumptions guarantee that (1.1)–(1.2) has a unique solution in H 2 (Ω) ∩ H01 (Ω) for all f ∈ L2 (Ω). Note that when is sufficiently small, condition (1.3) can be ensured by the other hypotheses and a transformation u = et(x+y) v with a positive constant t that satisfies 1 ~ y) ≥ c0 . c(x, y) + (β1 (x, y) + β2 (x, y))t − 2t2 − divβ(x, 2 Indeed, it is the case in which is very small that we are interested in. With the above assumption, the solution of (1.1)–(1.2) typically has boundary layers of width O( ln 1 ) at the outflow boundary x = 1 and y = 1. With some further assumptions, it is possible to characterize the boundary layers more precisely (see the regularity result in the next section). Our main concern here is superconvergence in a discrete -weighted energy norm k·k,N (see (2.6)) in the presence of exponential boundary layers. We shall establish an error bound of order N −2 ln2 N + N −3/2 in the discrete -weighted energy norm under certain regularity assumptions. For the one-dimensional case, see a recent work of the author [24]. As a consequence of the superconvergence result, we obtain convergence of the same order in the L2 -norm and pointwise convergence of order N −3/2 ln5/2 N + N −1 ln1/2 N at some mesh points inside the boundary layer under the same regularity assumption. These results are all uniformly valid with respect (1.3)


FINITE ELEMENT SUPERCONVERGENCE

1149

to . Furthermore, numerical tests indicate that the estimate N −2 ln2 N is sharp. It is worth pointing out that the error bounds obtained here are different from the error bounds obtained by Zhou [25] in that the Sobolev norms (kuk2 or kuk3 ) of the solution do not appear in the bounding constants. Recently, Melenk and Schwab have done some work on the p and the hp finite element methods for singularly perturbed problems in the two-dimensional setting. Their mesh design follows earlier work of Schwab and Suri in the one-dimensional reaction-diffusion problem [19], namely, the mesh size κp in the exponential boundary layer region is adopted. Here p is the polynomial degree in the finite element space and κ is a user-supplied constant. In [11], a robust exponential convergence rate is established for the reaction-diffusion equation under the analytic assumption on the input data. In [2], similar results are obtained for the dominant components (the smooth part and the layer part) of convection-diffusion problems. So far, a complete regularity analysis on the convection-diffusion equation seems lacking, although the counterpart results for reaction-diffusion problems are relatively rich [3, 5, 12]. Here is the outline of the article. After this brief introduction we introduce the method in Section 2. Section 3 serves as a preliminary to the analysis. In Section 4, we establish all ingredients for the proof of our main theorems, and in Section 5, we present and prove the main theorems. Finally, some numerical results are presented in Section 6. Throughout the article, the standard notation for the Sobolev spaces and norms will be used; and generic constants C, Ci are independent of and N . An index will be attached to indicate an inner product or a norm on a subdomain, for example, (·, ·)Ωx and k · kΩy . 2. The finite element method on a Shishkin mesh The regularity result. Regularity is a very complicated issue, and most of the known results are for reaction-diffusion equations. See [3], [5], [12], and [16]. Regarding convection-diffusion equations, the reader is referred to [20] and [10]. Here we adopt the result from the latter. Define the operator Li , i = 0, 1, by c ∂v ∂ i β2 ∂i Li v := +v i . ∂y ∂xi β1 ∂x β1 ~ and c be smooth, and let f ∈ C 4,1 (Ω) ¯ satisfy the compatibility Lemma 2.1. Let β conditions f (0, 0) = f (0, 1) = f (1, 0) = f (1, 1) = 0, and f f (0, 0) = (0, 0), β1 y β2 x f f f − L0 (0, 0) = (0, 0), β1 x β1 β 2 xx y f f f f f − L0 − L0 (0, 0) = (0, 0), − 2L1 β1 xx β1 x β1 β1 β 2 xxx y ! f f (0, 0). (0, 0) = β1 β2 β2 xx β1 yy


1150

ZHIMIN ZHANG

Then the boundary problem (1.1)–(1.2) has a classical solution u ∈ C 3,1 (Ω) which can be decomposed into u=u ¯ + w0 + w1 + w2 , where for all (x, y) ∈ Ω we have i+j ∂ u ¯ (2.1) ∂xi ∂y j (x, y) ≤ C for 0 ≤ i + j ≤ 2 and i+j ∂ w1 (2.2) (x, y) ∂xi ∂y j i+j ∂ w2 (2.3) (x, y) ∂xi ∂y j i+j ∂ w0 (2.4) ∂xi ∂y j (x, y)

≤ C−i e−α(1−x)/ , ≤ C−j e−α(1−y)/ , ≤ C−(i+j) e−α(1−x)/ e−α(1−y)/

~ c and f . for 0 ≤ i + j ≤ 3. Here the constant C depends on various norms of β, See [10, Theorem 5.1] for details. Note that when ∂ i+j f (0, 0) = 0, ∂xi ∂y j

1 ≤ i + j ≤ 3,

then the last four compatibility conditions of Lemma 2.1 are satisfied. The Shishkin mesh. Define the transition parameter 1 κ τ = min( , ln N ) 2 α with κ = 2.5, and divide Ω into four subdomains Ω0 = (0, 1 − τ )2 ,

Ωx = (1 − τ, 1) × (0, 1 − τ ),

Ωy = (0, 1 − τ ) × (1 − τ, 1),

Ωxy = (1 − τ, 1)2 .

Each subdomain is then decomposed into N × N (N ≥ 2) uniform rectangles (see Figure 1). Therefore, there are (2N + 1)2 nodes (xi , yj ), i, j = 0, 1, 2, . . . , 2N , and 4N 2 elements Ωij = (xi−1 , xi ) × (yj−1 , yj ), We denote H=

1−τ , N

h=

j = 1, 2, . . . , 2N. τ . N

2.5 ln N , since otherwise N −1 is much In the later analysis, we assume that τ = α less than and the traditional finite element analysis can be applied. For small , the Shishkin mesh is highly graded with ratio of H/h = O(−1 ). It is neither regular nor quasi-uniform. The parameter τ is selected so as to deal with the singular behavior of the boundary layer functions w1 , w2 , and w0 . In the boundary layer region, the small mesh size compensates for the sharp change of the solution. We see that 2.5 ln N h = . α N



1151

L ?

Ωy

Ωxy

S

S

Ωx

Ω0 S

Figure 1. The Shishkin mesh Outside the boundary layer, the exponential decay of w1 , w2 , and w0 dominate: Z 1−τ −2.5 e−α(1−x)/ dx ≤ e−ατ / = eln N = . α α αN 2.5 0 In the analysis, these facts are used repeatedly. Remark 2.1. In the literature, κ = 2 is widely used in determining the transition point for the Shishkin mesh. Our numerical results reveal the same convergent rates for κ = 1.5, κ = 2, and κ = 2.5. However, κ = 2 has a better error distribution than the other nearby numbers (see Section 6). For technical reasons, we use κ = 2.5 in our analysis. Variational formulation. The weak formulation of the model problem (1.1)– (1.2) reads: Find u ∈ H01 (Ω) such that B (u, v) = f (v),

∀v ∈ H01 (Ω),


1152

ZHIMIN ZHANG

where Z ~ · ∇u, v) + (cu, v), B (u, v) = (∇u, ∇v) + (β

f (v) = (f, v) =

f vdxdy. Ω

We define an energy norm k · k by kvk2 = k∇vk2 + kvk2 = |v|2 + kvk2 , where k · k is the L2 -norm. We have, from integration by parts and applying (1.3), 1 ~ v) ≥ |v|2 + c0 kvk2 ≥ min(1, c0 )kvk2 . (2.5) B (v, v) = (∇v, ∇v) + ((c − divβ)v, 2 Let VN ⊂ H01 (Ω) be the C 0 bilinear finite element space on the Shishkin mesh; we look for uN ∈ VN such that B (uN , v) = f (v),

∀v ∈ VN .

We define a discrete energy norm k · k,N by (2.6)

kvk2,N = |v|2,N + kvk2

with |v|2,N =

X

4hK ~K |∇v(xK , yK )|2 .

K

Here K = (xK − hK , xK + hK ) × (yK − ~K , yK + ~K ) is an element (see Figure 2). For the Shishkin mesh, 2hK , 2~K are either h or H. Main task and difficulties. The main task is to establish the approximability of the bilinear finite element space to functions with exponential terms of arbitrarily large parameters in the energy norm as well as in the discrete energy norm (2.6). There are two difficulties: (i) The bilinear form B does not satisfy the uniform stability |B (u, v)| ≤ Ckuk kvk for a constant C independent of , although it does satisfy the coercivity condition (2.5). (ii) The bilinear interpolant uI of the solution u cannot be uniformly bounded by u in either the L2 -norm or the H 1 -norm as kuI k ≤ Ckuk,

k∇uI k ≤ Ck∇uk.

for a constant C independent of . However, all the error bounds must be -uniform. The standard finite element analysis cannot produce the expected result, and the situation is further complicated by the superconvergent consideration. In this work we shall use a different framework to overcome these difficulties. Furthermore, integral identities developed in the 90’s (see the Appendix) are used to prove superconvergence. The analysis is very delicate.



1153

l3 ~K hK

l4

l2

(xK , yK )

l1 Figure 2. Geometry of the element K

3. Preliminaries On an individual rectangular element K (see Figure 2), v ∈ VN is defined as v(xK + hK ξ, yK + ~K η) = +

vK v1K (1 − ξ)(1 − η) + 2 (1 + ξ)(1 − η) 4 4

vK v3K (1 + ξ)(1 + η) + 4 (1 − ξ)(1 + η) = vˆ(ξ, η), 4 4

b = [−1, 1]2 , (ξ, η) ∈ K

where v1K = v(xK − hK , yK − ~K ),

v2K = v(xK + hK , yK − ~K ),

v3K = v(xK + hK , yK + ~K ),

v4K = v(xK − hK , yK + ~K ).

As a preliminary, we first introduce some inequalities for v ∈ VN that will be used in the analysis. Their proofs are straightforward calculations, and hence are omitted. There are general results for most of these inequalities; however, the results here provide specific information about the bounding constants which may not appear elsewhere. Imbedding inequalities: Z ! Z Z ! Z Z Z 9 9 2 2 2 + v dxdy, + v 2 dxdy. v dy ≤ v dx ≤ (3.1) K K K h ~ K K lK l K l l K 2 4 1 3 Inverse inequalities: Z 2 Z Z Z 2 ∂v ∂v 9 9 2 (3.2) dxdy ≤ 2 v dxdy, dxdy ≤ 2 v 2 dxdy; ∂x hK K ∂y ~K K K K (3.3) Z 2 Z 2 Z 2 2 ∂ v ∂v ∂v 3 3 dxdy ≤ 2 dxdy, or ≤ 2 dxdy. ∂x∂y 2~ ∂x 2h ∂y K K K K K Stability inequality: (3.4)

kvk,N ≤ kvk .


1154

ZHIMIN ZHANG

Discrete inequalities: Z (3.5) v 2 dxdy ≤ hK ~K [(v1K )2 + (v2K )2 + (v3K )2 + (v4K )2 ], K Z h H 2 + [(v1K )2 + (v2K )2 + (v3K )2 + (v4K )2 ]. (3.6) |∇v| dxdy ≤ h H K In this article, we shall frequently use the bilinear interpolation wI of a given function w. We start from two identities which again can be derived through simple 1 (Ω), calculation. When w ∈ W∞ (3.7) (3.8)

∂w ∂w (x, yK − ~K ) + (x, yK + ~K ) dx, ∂x ∂x xK −hK Z yK +~K I ∂w 1 ∂w ∂w (xK , yK ) = (xK + hK , y) + (xK − hK , y) dy; ∂y 4~K yK −~K ∂y ∂y 1 ∂wI (xK , yK ) = ∂x 4hK

Z

xK +hK

3 (Ω), we have and if w ∈ W∞

∂ (wI − w)(xK , yK ) ∂x Z hK 2 3 ~2K ∂ 3 w t ∂ w ∂3w 1 + − t~K 2 = (xK + s1 t, yK − s1 ~K ) 4hK −hK 2 ∂x3 ∂x ∂y 2 ∂x∂y 2 2 3 ~2K ∂ 3 w t ∂ w ∂3w + (3.9) + + t~K 2 (xK + s2 t, yK + s2 ~K ) dt, 2 ∂x3 ∂x ∂y 2 ∂x∂y 2 ∂ (wI − w)(xK , yK ) ∂y Z ~K 2 hK ∂ 3 w ∂3w t2 ∂ 3 w 1 − h t + (xK − s3 hK , yK + s3 t) = K 4~K −~K 2 ∂x2 ∂y ∂x∂y 2 2 ∂y 3 2 hK ∂ 3 w t2 ∂ 3 w ∂3w + hK (3.10) + + (xK + s4 hK , yK + s4 t) dt, 2 ∂x2 ∂y ∂x∂y 2 2 ∂y 3 where 0 < si < 1 for i = 1, 2, 3, 4. Also, (3.11)

kw − wI k2K ≤ C(kh2K

2 ∂2w 2 ∂ 2w 2 2 ∂ w 2 k k + kh ~ + k~ k ), K K K K K ∂x2 ∂x∂y ∂y 2 K

where C is a constant independent of hK , ~K , and w. Finally, we list some inequalities regarding the exponential boundary layer functions which will be frequently used in the next section. (3.12) (3.13)

2hK e

−2α(1−xK )/

Z

xK +hK