- Email: [email protected]
Using rectangular Qp elements in the SDFEM for a convection–diffusion problem with a boundary layer
Applied Numerical Mathematics 58 (2008) 1789–1802 www.elsevier.com/locate/apnum
Using rectangular Qp elements in the SDFEM for a convection–diffusion problem with a boundary layer ✩ Martin Stynes a,∗ , Lutz Tobiska b a National University of Ireland, Cork, Ireland b Institut für Analysis und Numerik, Otto-von-Guericke Universität Magdeburg, Postfach 4120, D-39016 Magdeburg, Germany
Available online 23 November 2007
Abstract The streamline diffusion finite element method (SDFEM; the method is also known as SUPG) is applied to a convection– diffusion problem posed on the unit square whose solution has exponential boundary layers. A rectangular Shishkin mesh is used. The trial functions in the SDFEM are piecewise polynomials that lie in the space Qp , i.e., are tensor products of polynomials of degree p in one variable, where p > 1. The error bound IN u − uN SD CN −(p+1/2) is proved; here uN is the computed SDFEM solution, IN u is chosen in the finite element space to be a special approximant of the true solution u, and ·SD is the streamline-diffusion norm. This result is compared with previously known results for the case p = 1. The error bound is a superclose result; uN can be enhanced using local postprocessing to yield a modified solution u˜ N for which u − u˜ N SD CN −(p+1/2) . © 2007 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 65N15; 65N30 Keywords: Convection–diffusion; Streamline-diffusion method; Shishkin mesh
1. Introduction We consider the singularly perturbed boundary value problem Lu := −εu + b · ∇u + cu = f
on Ω = (0, 1)2 ,
u = 0 on ∂Ω,
(1.1)
where the parameter ε lies in the interval (0, 1]. Assume that the function b(x, y) = (b1 (x, y), b2 (x, y)) with ¯ and b1 (x, y) > β1 > 0 and b2 (x, y) > β2 > 0, c(x, y) 0 on Ω, b(x, y) ¯ (1.2) c(x, y) − div c0 > 0 on Ω, 2 where β1 , β2 and c0 are some constants. We assume that the functions b, c and f are sufficiently smooth. These hypotheses ensure that (1.1) has a unique solution in H01 (Ω) ∩ H 2 (Ω) for all f ∈ L2 (Ω). Note that for sufficiently small ε, the other hypotheses imply that (1.2) can always be ensured by the simple change of variable v(x, y) = e−σ x u(x, y) when σ is chosen suitably. ✩
We dedicate this paper to our colleague and friend Pieter Hemker on the occasion of his 65th birthday.
* Corresponding author.
E-mail addresses: [email protected] (M. Stynes), [email protected] (L. Tobiska). 0168-9274/$30.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2007.11.004
1790
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
The convection–diffusion problem (1.1) is a well-known 2-dimensional testbed for numerical methods for advection-dominated transport phenomena; it exhibits exponential boundary layer behaviour at the sides x = 1 and y = 1 of Ω that is typical of such problems. Thus layer-adapted meshes are usually used to solve (1.1). A piecewiseuniform Shishkin mesh can be chosen a priori when one has some understanding of the structure of these layers [3,8,10]. This mesh has O(N 2 ) points in a rectangular grid that is refined near the sides x = 1 and y = 1. In [12] we showed that for piecewise bilinears (the space Q1 ), on writing uI for the piecewise bilinear nodal interpolant to u, then the error uI − uN L2 (Ω) is second-order convergent (up to a logarithmic factor), uniformly in ε; this is the optimal rate of convergence. What about higher-order elements (the space Qp , where p > 1)? The present paper is the first to analyse higher-order finite element methods on Shishkin meshes for convection–diffusion problems. Errors will be estimated in both the streamline diffusion and L2 norms. Does one again attain the optimal order of convergence in L2 ? We answer this question and discover that p = 1 is atypical. The extension of the analysis of [12] to Qp is not routine: for example, instead of a standard interpolant from Qp , one should use the vertices– edges–element–approximant of Lin [5]. The non-standard anisotropic interpolation properties of this approximant are studied in Lemma 5 below. The organization of the paper is as follows. In Section 2 we describe the Shishkin mesh and the SDFEM. A decomposition of the solution u and some approximation-theoretical facts that we shall need later are presented in Section 3. The convergence of the SDFEM is analysed in Section 4. Notation. Throughout the paper C will denote a generic positive constant that is independent of ε and the mesh. Subscripted constants such as C1 are also independent of ε and the mesh, but have a fixed value. We use the standard Sobolev spaces W k,p (D), H k (D) = W k,2 (D), H0k (D), Lp (D) = W 0,p (D) for non-negative integers k and 1 p ∞, and write (·,·)D for the L2 (D) inner product. Here D is any measurable subset of Ω. Then |·|k,D and ·k,D are the usual Sobolev seminorm and norm on H k (D). When D = Ω, we drop D from the notation for simplicity. 2. The Shishkin mesh and the SDFEM Shishkin meshes [8–10] are piecewise uniform meshes, constructed a priori, that are refined inside layers. Let N be an even positive integer. We let λx and λy denote two mesh transition parameters that will be used to specify where the mesh changes from coarse to fine; these are defined by 1 (p + 1)ε 1 (p + 1)ε λx = min , ln N and λy = min , ln N , 2 β1 2 β2 where our trial space, which is defined below, comprises functions that are piecewise in Qp for some integer p > 1. In fact we assume that λx = ((p + 1)ε/(β1 )) ln N and λy = ((p + 1)ε/(β2 )) ln N , as otherwise N −1 is much smaller than ε (which is very unlikely in practice) and one can then simplify our analysis of the method. In these formulae, if the multiplier of ε/β1 and ε/β2 is chosen too small then the final order of convergence obtained will decrease. For the case p = 1 this dependence is investigated both theoretically and numerically in [11]. We divide Ω as in Fig. 1: Ω¯ = Ω11 ∪ Ω21 ∪ Ω12 ∪ Ω22 , where Ω11 = [0, 1 − λx ] × [0, 1 − λy ], Ω21 = [1 − λx , 1] × [0, 1 − λy ], Ω12 = [0, 1 − λx ] × [1 − λy , 1], Ω22 = [1 − λx , 1] × [1 − λy , 1]. ¯ i, j = 0, . . . , N} are the rectangular lattice defined by The mesh points Ω N = {(xi , yj ) ∈ Ω: for i = 0, . . . , N/2, 2i(1 − λx )/N xi = 1 − 2(N − i)λx /N for i = N/2 + 1, . . . , N, and
yj =
2j (1 − λy )/N 1 − 2(N − j )λy /N
for j = 0, . . . , N/2, for j = N/2 + 1, . . . , N.
Our mesh is constructed by drawing lines parallel to the coordinate axes through these mesh points, so it is a tensor product of two one-dimensional piecewise uniform meshes. This divides Ω into a set TN of mesh rectangles K whose sides are parallel to the axes—see Fig. 1. The mesh is coarse on Ω11 , coarse/fine on Ω21 ∪ Ω12 , and fine on Ω22 .
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1791
Fig. 1. Division of Ω and Shishkin mesh for T8 .
√ √ The mesh is quasiuniform on Ω11 and its diameter d there satisfies 2/N d 2 2/N ; on Ω12 ∪ Ω21 , each mesh rectangle has dimensions O(N −1 ) by O(εN −1 ln N ); and on Ω22 each rectangle is O(εN −1 ln N ) by O(εN −1 ln N ). We shall use these properties several times in our analysis. Given a mesh rectangle K, its dimensions are written as hx,K by hy,K and its barycentre is denoted by (xK , yK ). We now describe the SDFEM on this rectangular mesh. Let p be a fixed integer with p > 1. (In [12] we analysed the case p = 1.) Our trial space is ¯ v|∂Ω = 0 and v|K ∈ Qp (K) ∀K ∈ TN , V N = v ∈ C(Ω): where Qp (K) = span{x i y j : 0 i, j p}. Define the bilinear form BSD (·,·) by BSD (w, v) = BGal (w, v) + BStab (w, v), where BGal (w, v) = ε(∇w, ∇v) + (b · ∇w + cw, v), δK (−εw + b · ∇w + cw, b · ∇v)K , BStab (w, v) =
(2.1) (2.2)
K∈TN
where H˜ 1 (Ω) denotes the set of functions in H 1 (Ω) that lie in H 2 (K) for each K, for all and δK 0 is a user-chosen piecewise constant parameter. The SDFEM is defined as follows: Find uN ∈ V N such that BSD (uN , v N ) = (f, v N ) + δK (f, b · ∇v N )K ∀v N ∈ V N . (2.3) (w, v) ∈ H˜ 1 (Ω) × H 1 (Ω),
K∈TN
We clearly have the Galerkin orthogonality property BSD (u − uN , v N ) = 0 ∀v N ∈ V N .
(2.4)
The choice of δK is discussed in, e.g., [10, §III.3.2.1]. For each K ∈ TN , set hK = min{hx,K , hy,K }. Let μ be a constant such that the inverse inequality N v N 0,K μh−1 K ∇v 0,K
∀v N ∈ V N , ∀K ⊂ Ω11
is valid. Similarly to [10, p. 233], we set ⎧ if K ⊂ Ω11 and ε N −1 , ⎨ C1 N −1 −1 −2 δK = C1 ε N if K ⊂ Ω11 and ε > N −1 , ⎩ 0 otherwise, where the positive constant C1 is chosen (independently of ε) such that h2K c0 1 ∀K ⊂ Ω11 . , 0 δK min 2 (maxK |c(x, y)|)2 μ2 ε
(2.5)
(2.6)
1792
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
Then the argument of [10, p. 231] shows that the inequality 1 BSD (v N , v N ) v N 2SD ∀v N ∈ V N 2 holds, where ·SD is given by 1/2 2 2 2 vSD = ε|v|1 + δK b · ∇v0,K + c0 v0
(2.7)
∀v ∈ H 1 (Ω).
(2.8)
K∈TN
It follows that (2.3) has a unique solution uN . We shall also use the ε-weighted energy norm v1,ε = (ε|v|21 + c0 v20 )1/2 for all v ∈ H 1 (Ω). This norm is clearly weaker than ·SD . 3. Solution decomposition, Lin’s approximant We now assume that u can be decomposed in a precise way that is typical of the behaviour observed in solutions of (1.1) when interior layers are absent. Assumption 1. Assume that u = S + E21 + E12 + E22 , where there exists a constant C such that for all (x, y) ∈ Ω we have i+j ∂ S ∂x i ∂y j (x, y) C for 0 i + j p + 1, i+j ∂ E21 Cε −i e−β1 (1−x)/ε for 0 i + j p + 2, (x, y) ∂x i ∂y j i+j ∂ E12 Cε −j e−β2 (1−y)/ε for 0 i + j p + 2, (x, y) ∂x i ∂y j i+j ∂ E22 −(i+j ) −(β1 (1−x)+β2 (1−y))/ε e for 0 i + j p + 2. ∂x i ∂y j (x, y) Cε
(3.1)
(3.2) (3.3) (3.4) (3.5)
Furthermore, assume that S ∈ H p+2 (Ω) with Sp+2 C.
(3.6)
Here S is the smooth part of u, E21 is an exponential boundary layer along the side x = 1 of Ω, E12 is an exponential boundary layer along the side y = 1, while E22 is an exponential corner layer at (1, 1). Remark 2. In [7] a proof is given that under certain compatibility conditions on the data f , the bounds (3.2)–(3.6) of Assumption 1 are valid if p = 1. The extension of these results to our case p > 1 is non-trivial but an inspection of the argument in [7] convinces one that it can be done by placing a sufficient number of compatibility conditions on f . The number of these conditions increases rapidly with p, which is restrictive. Nevertheless in the case of a problem that did not satisfy all these assumptions, one could combine the argument of the current paper with cut-off functions to obtain bounds on the error in the computed solution in regions that were sufficiently far from all layers induced by inadequate compatibility of the data, much like the analysis in [4]. ¯ we shall follow the approach of Lin [5] in defining a non-standard “vertices–edges– Given any function v ∈ C(Ω), element approximant” IK v ∈ Qp on each mesh rectangle K ∈ TN . We use this approximant because one then has, compared with standard interpolants, sharper bounds (see Lemma 4 and Remark 11 below) for the errors that appear in various terms in the finite element analysis. Lin’s approximant is based on nodal values of v at the corners of K together with moments of v on the edges of K and on K itself. It is initially defined on the reference domain Kˆ = (−1, +1)2 , whose vertices and edges we
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1793
ˆ The approximation operator Iˆ : C(K) ˆ → Qp (K) ˆ denote by aˆ i and eˆi respectively for i = 1, . . . , 4. Let v(·,·) ˆ ∈ C(K). is defined by requiring ˆ aˆ i ) for i = 1, . . . , 4, (Iˆv)( ˆ aˆ i ) = v( (Iˆv)q ˆ dγ = vq ˆ dγ ∀q ∈ Pp−2 (eˆi ) for i = 1, . . . , 4, eˆi
eˆi
(Iˆv)q ˆ dξ dη =
Kˆ
vq ˆ dξ dη
(3.7a) (3.7b)
ˆ ∀q ∈ Qp−2 (K);
(3.7c)
Kˆ
here Pp−2 (eˆi ) is the space of polynomials of degree at most p − 2 in the single variable whose axis is parallel to the edge eˆi . Lemma 3. The approximation operator Iˆ is uniquely defined. Proof. There is a total of (p + 1)2 conditions on Iˆvˆ and we have (p + 1)2 degrees of freedom for Iˆvˆ ∈ Qp . Thus it is enough to show that taking all degrees of freedom equal to zero implies Iˆvˆ = 0. ˆ the function rˆ = Iˆv| Assume that all degrees of freedom are equal to zero. Then for any edge eˆ of K, ˆ eˆ is a 2 polynomial of degree p in one variable on (−1, 1) that is L -orthogonal to the space Pp−2 (e) ˆ of polynomials of p degree p − 2. The L2 -orthogonality of the Legendre polynomials {Lm }m=0 implies that rˆ = ap−1 Lp−1 + ap Lp with unknown constants ap−1 and ap . In addition, rˆ must vanish at the endpoints of e, ˆ i.e., rˆ (±1) = 0, which implies that rˆ ≡ 0. Applying this argument for each edge of Kˆ we see that Iˆvˆ vanishes on the boundary of Kˆ and therefore (Iˆv)(ξ, ˆ η) = (1 − ξ 2 )(1 − η2 )q(ξ, ˆ η)
ˆ with qˆ ∈ Qp−2 (K).
But our assumption that all degrees of freedom equal zero requires in particular that in Iˆvˆ ≡ 0. 2
ˆ
ˆ qˆ dξ Kˆ (I v)
dη = 0, which results
Using the transformation x = xK + hx,K ξ/2, y = yK + hy,K η/2 to map from Kˆ to an arbitrary K ∈ TN , one obtains the corresponding approximation operator IK : C(K) → Qp (K). Then a continuous global approximation ¯ → V N is defined by setting operator IN : C(Ω) (IN v)|K := IK (v|K )
∀K ∈ TN .
Lemma 4. Let K ∈ TN . Let w ∈ H p+2 (K) and let IK w ∈ Qp (K) be its vertices–edges–element approximant. Then for each v N ∈ Qp (K) we have p+2 w v N 0,K (IK w − w)x v N dx dy Chp+1 ∂ x x y,K p+1 ∂x∂y 0,K K
and
p+2 w (IK w − w)y v N dx dy Chp+1 ∂ y x,K p+1 ∂x ∂y K
Proof. For each mesh rectangle K, we set hy,K 2 1 2 . FK (y) = (y − yK ) − 2 2 Now, from the proof of [5, Lemma 4] it follows that
0,K
vyN 0,K .
1794
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
(IK w − w)x vxN dx dy = K
(−2)p (2p + 1)! +
(−2)p
FK (y)FK (y) p
∂ p+2 w(x, y) ∂ p+1 v N (x, yK ) dx dy ∂x∂y p ∂x∂y p+1
K
p
FK (y)
(2p)!
∂ p+2 w(x, y) ∂ p v N (x, yK ) dx dy. ∂x∂y p+1 ∂x∂y p−1
K
Clearly p F (y)
K
h2y,K
p
8
hy,K , and FK (y) 2
and for arbitrary functions v N ∈ Qp (K) one has the identities ∂ p+1 v N (x, yK ) ∂ p+1 v N (x, y) = ∂x∂y p ∂x∂y p and ∂ p v N (x, yK ) ∂ p v N (x, y) ∂ p+1 v N (x, y) = − (y − yK ) p−1 p−1 ∂x∂y p ∂x∂y ∂x∂y as well as the inverse inequalities n+1 N N ∂v v n ∂ C hy,K ∂x n ∂x∂y 0,K 0,K
for n = p − 1, p.
Invoking these above, we get the first statement of lemma. The second follows analogously.
2
To obtain satisfactory approximation error estimates on Ω12 ∪ Ω21 , where the Shishkin mesh is highly anisotropic, we shall use some technical results of Apel and Dobrowolski [1]. Let r be an integer with 1 r p. Define the anisotropic space i+j +1 ∂ vˆ r ˆ 1 ˆ < ∞, 0 i + j r Hξ (K) = vˆ ∈ L (K): i+1 j ∂ξ ∂η 0,Kˆ ˆ similarly. and define Hηr (K) Lemma 5. Let r be an integer with 1 r p. There exists a constant C independent of vˆ such that for the vertex– edges–element approximant Iˆvˆ Dξ (vˆ − Iˆv) ˆ ∩ Hξr (K). ˆ ˆ 0,Kˆ C|Dξ v| ˆ r,Kˆ ∀vˆ ∈ C(K) (3.8) ˆ are replaced by Dη and Hηr (K), ˆ respectively. An analogous inequality holds if Dξ and Hξr (K) ˆ → R by Proof. Define r + r 2 continuous linear functionals Fi : H r (K) vˆ → ξ vˆ dξ dη for = 0, . . . , r − 1, Kˆ
vˆ →
ξ j ηk Kˆ
∂ vˆ dξ dη ∂η
for j = 0, . . . , r − 1 and k = 0, . . . , r − 1.
Then, recalling the definition of Iˆ, it is straightforward to verify that ˆ ∩ Hξr (K) ˆ for i = 1, . . . , r + r 2 ˆ = 0 ∀vˆ ∈ C(K) Fi Dξ (Iˆvˆ − v) and also that ˆ = 0 for i = 1, . . . , r + r 2 Fi (Dξ v)
ˆ implies Dξ vˆ = 0 ∀vˆ ∈ Qr (K).
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1795
Now inequality (3.8) follows immediately from [1, Lemma 3] by setting (in the notation of that paper) Q = {(i, j ): 0 i, j r − 1}, P = {(i, j ): 0 i, j r}, and γ = (1, 0). For the analogous inequality involving Dη , modify the above linear functionals by interchanging ξ with η and replacing the η-derivative by the ξ -derivative, and take γ = (0, 1) in [1, Lemma 3]. 2 The anisotropic estimate of Lemma 5 implies the following sharp bounds. Lemma 6. Let r be an integer with 1 r p. Let K ∈ TN . Then there exists a constant C such that the vertices– edges–element approximant IN v satisfies r+1 j i ∂ v , (v − IN v)x C h h x,K y,K 0,K i+1 j ∂x ∂y 0,K i+j =r r+1 ∂ v j i (v − IN v)y C hx,K hy,K ∂x i ∂y j +1 0,K 0,K i+j =r
for all v ∈ H r+1 (K). Proof. Using the transformation x = xK +
hx,K ξ, 2
y = yK +
hy,K η, 2
ˆ apply the anisotropic result of Lemma 5, then map back to the we map the element K onto the reference element K, original element K. Thus 2|K|1/2 |K|1/2 (v − IN v)x (vˆ − Iˆv) ˆ ξ 0,Kˆ C |Dξ v| ˆ r,Kˆ 0,K hx,K hx,K i+j +1 i+j +1 ∂ ∂ v v |K|1/2 j i C C hx,K hy,K i+1 j i+1 h ∂ξ ∂η ∂x ∂y j ˆ x,K i+j =r
0,K
The proof of the second inequality is along the same lines.
i+j =r
.
0,K
2
As well as Lemma 6, we have the following more standard anisotropic approximation error estimate. Lemma 7. Let K ∈ TN . Let q ∈ [1, ∞]. Then there exists a constant C such that the vertices–edges–element approximant IN v satisfies p+1 ∂ v j i hx,K hy,K v − IN vLq (K) C ∂x i ∂y j q L (K) i+j =p+1
for all v ∈ W p+1,q (K). Proof. Lemma 4 of [1] implies that we can invoke [1, Lemma 3] with Q = P = {(i, j ): 0 i, j p} and γ = ˆ (0, 0). This bounds v − IN vLq (K) ˆ , and the desired result follows by mapping between K and K as in the proof of Lemma 6. 2 For the Lin approximant one has the stability bound IN vL∞ (K) CvL∞ (K)
(3.9)
¯ where the constant C is independent of v and K ∈ Tn . This can be proved by mapping (3.7) to the for all v ∈ C(K), ˆ reference element K. From Lemma 7 one can deduce a pointwise bound on the error in the approximant.
1796
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
Lemma 8. Let Assumption 1 hold. Then there exists a constant C such that −(p+1) if (x, y) ∈ Ω11 , (u − IN u)(x, y) CN C(N −1 ln N )p+1 otherwise. Proof. Using the decomposition (3.1), write u − IN u = (S − IN S) + (E12 − IN E12 ) + (E21 − IN E21 ) + (E22 − IN E22 ). By Lemma 7 (with q = ∞) and Assumption 1, (S − IN S)(x, y) CN −(p+1) S p+1,∞ CN −(p+1) W Similarly, on Ω12 ∪ Ω22 one gets (E12 − IN E12 )(x, y) C
(3.10)
∀(x, y) ∈ Ω.
N −i (εN −1 ln N )j ε −j C(N −1 ln N )p+1 ;
i+j =p+1
on Ω11 ∪ Ω21 , the choice of transition point 1 − λy in the Shishkin mesh and (3.9) yield IN E12 (x, y) C max E12 (x, y) CN −(p+1) . Ω11 ∪Ω21
The functions E21 and E22 are handled similarly. Putting the various bounds into (3.10) gives the desired result.
2
The next lemma collects several results that involve the approximants of the various terms in the decomposition (3.1). Lemma 9. Let Assumption 1 hold. Let E denote any one of the functions E12 , E21 or E22 . Then there exists a constant C such that the following approximation error estimates hold: S − IN S0,Ω CN −(p+1) , E0,Ω11 Cε
1/2
N
−(p+1)
(3.11) (3.12)
,
εEL1 (Ω11 ) + ∇EL1 (Ω11 ) CN
−(p+1)
(3.13)
,
E − IN E0,Ω C(N −1 ln N )p+1 ,
(3.14) −1 1/2
IN E12 0,Ω11 ∪Ω21 + IN E21 0,Ω11 ∪Ω12 C(ε + N ) 1/2 −(p+1) N . IN E22 0,Ω12 ∪Ω21 C ε(ε + N −1 ) ln N
N
−(p+1)
(3.15)
,
(3.16)
Proof. Imitate [12, Lemma 3.2] to prove (3.11)–(3.14). For (3.15) we shall give only the argument to bound IN E12 0,Ω11 ∪Ω21 as IN E21 0,Ω11 ∪Ω12 is handled analogously. Using (3.9) and the decay property (3.4) of E12 , IN E12 Ω11 ∪Ω21 =
N/2
yj
1
2 IN E12 (x, y) dx dy
1/2
j =1 y=yj −1 x=0
C
N/2
yj
1/2
1 e
−2β2 (1−yj )/ε
dx dy
j =1 y=yj −1 x=0
C
N/2−1 y j +1 1
e
−2β2 (1−y)/ε
y N/2
dx dy +
j =1 y=yj x=0
1/2 C (ε + N −1 )N −2(p+1) , by the definition of λy . Similarly, using (3.5) instead of (3.4),
1/2
1 e
−2β2 λy /ε
dx dy
y=yN/2−1 x=0
(3.17)
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
IN E22 Ω21 C
N/2−1 y j +1
1 e
−2β2 (1−y)/ε
y N/2
e
−2β2 λy /ε
dx dy
y=yN/2−1 x=1−λx
1/2 C ε(ε + N −1 )N −2(p+1) ln N . The bound on IN E22 Ω12 is analogous.
1/2
1
dx dy +
j =1 y=yj x=1−λ x
1797
(3.18)
2
4. Error bound for IN u − uN In the analysis of this section, we shall try to give the sharpest possible bound for each term in order to reveal the source of the essential difference between the order of accuracy attained by Q1 elements and that attained by Qp elements when p > 1. The extreme nonuniformity of the mesh forces us to employ different analytical techniques on the different regions of Fig. 1. The proof of our main result (Theorem 15) requires a bound on BSD (uN − IN u, uN − IN u). To achieve this, write BSD = BGal + BStab , then bounds on the individual terms in each of these bilinear forms are proved in Lemmas 10–13. Lemma 10. Let Assumption 1 hold. Then there exists a constant C such that for all w N ∈ V N we have ε∇(u − IN u), ∇w N C N −(p+1) + ε 1/2 (N −1 ln N )p+1 w N 1,ε .
(4.1)
Proof. Let w N ∈ V N be arbitrary but fixed. As in (3.10), we write u − IN u = (S − IN S) + (E12 − IN E12 ) + (E21 − IN E21 ) + (E22 − IN E22 ). Then consider each corresponding term in (ε∇(u − IN u), ∇w N ). First, Lemma 4 implies that ε(S − IN S)x , w N CεN −(p+1) Sp+2 w N 0 Cε 1/2 N −(p+1) w N 1,ε , x x since Sp+2 C by Assumption 1. Next, consider ε(E12 − IN E12 )x , w N ε(E12 − IN E12 )x , w N x x Ω
12 ∪Ω22
(4.2)
+ ε(E12 − IN E12 )x , w N
x Ω11 ∪Ω21
.
The calculation on Ω12 ∪ Ω22 resembles that for S: ε(E12 − IN E12 )x , w N
x Ω12 ∪Ω22
Cε(εN −1 ln N )p+1
ε
−2(p+1) −2β2 (1−y)/ε
e
1/2 dx dy
wxN 0
Ω12 ∪Ω22
Cε(N
−1
ln N )
On Ω11 ∪ Ω21 we use Lemma 6 with r = 1: ε(E12 − IN E12 )x , w N x Ω11 ∪Ω21 CεN −1 (E12 )xx 0,Ω ∪Ω + (E12 )xy 0,Ω 11
CεN
−1
21
p+1
w N 1,ε .
11 ∪Ω21
(4.3)
N wx 0
1/2 1−λ y −2 −2β2 (1−y)/ε ε e dy wxN 0 CN −(p+2) w N 1,ε ,
(4.4)
y=0
by the choice of transition point 1 − λy . The E21 integral is also broken into two regions as ε(E21 − IN E21 )x , w N ε(E21 − IN E21 )x , w N x x Ω
21 ∪Ω22
+ ε(E21 − IN E21 )x , w N
For the first term on the right-hand side here, we again use Lemma 4 to get
x Ω11 ∪Ω12
.
1798
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
ε(E21 − IN E21 )x , w N x Ω21 ∪Ω22 p+2 ∂ E21 CεN −(p+1) wxN 0,Ω21 ∪Ω22 ∂x∂y p+1 0,Ω21 ∪Ω22 1/2 −(p+1) −2 −2β1 (1−x)/ε CεN ε e dx dy wxN 0 CN −(p+1) w N 1,ε .
(4.5)
Ω21 ∪Ω22
The other term will be handled differently for ε N −1 and ε < N −1 , respectively. If ε N −1 we use Lemma 6 and the choice of transition point 1 − λx to get N ε(E21 − IN E21 )x , w N (E21 )xy CεN −1 (E21 )xx wx 0 + x Ω ∪Ω 0,Ω ∪Ω 0,Ω ∪Ω 11
12
11
12
11
12
Cε 1/2 N −1 (ε −3/2 + ε −1/2 )N −(p+1) ε 1/2 wxN 0 CN −(p+1) w N 1,ε . For ε < N −1 we have ε(E21 − IN E21 )x , w N
x Ω11 ∪Ω12
ε (E21 )x
0,Ω11 ∪Ω12
+ (IN E21 )x 0,Ω
11 ∪Ω12
N wx Ω11 ∪Ω12 .
Assumption 1 and the choice of transition point 1 − λx imply that 1/2 1−λ x 1/2 1/2 −2 −2β1 (1−x)/ε (E21 )x ε Cε ε e dx CN −(p+1) . 0,Ω11 ∪Ω12
x=0
Invoking an inverse inequality and (3.15), we have ε 1/2 (IN E21 )x Cε 1/2 NIN E21 Ω Ω11 ∪Ω12
11 ∪Ω12
CN −(p+1)
(4.6)
since ε < N −1 . Putting these inequalities together, CN −(p+1) w N 1,ε . ε(E21 − IN E21 )x , w N
(4.7)
x Ω11 ∪Ω12
The E22 term is handled like the other two exponential layer terms. First, + ε(E22 − IN E22 )x , w N ε(E22 − IN E22 )x , w N ε(E22 − IN E22 )x , w N x x Ω x Ω\Ω . 22
22
Then, similarly to (4.3), Lemma 4 yields
p+2 E22 Cε(εN −1 ln N )p+1 ∂ wxN 0 x Ω12 ∪Ω22 ∂x∂y p+1 0,Ω12 ∪Ω22 1/2 Cε(εN −1 ln N )p+1 ε −2(p+2) e−2β1 (1−x)/ε e−2β2 (1−y)/ε dx dy wxN 0
ε(E22 − IN E22 )x , w N
Ω12 ∪Ω22
Cε
1/2
(N
−1
p+1
ln N )
w N 1,ε .
(4.8)
N −1
< N −1
N −1 ,
and ε separately. For ε apply Lemma 6 with r = 1 to On Ω11 ∪ Ω21 we consider the cases ε get ε(E22 − IN E22 )x , w N ε 1/2 N −1 (E22 )xx + (E22 )xy 0,Ω ∪Ω w N 1,ε x Ω ∪Ω 0,Ω ∪Ω 11
21
11
ε For ε
< N −1
we use the splitting ε(E22 − IN E22 )x , w N x Ω
11 ∪Ω21
1/2
−1
(εN )
N
−(p+1)
ε (E22 )x
21
w 1,ε ε
0,Ω11 ∪Ω21
N
11
1/2
N
−(p+1)
+ (IN E22 )x 0,Ω
w 1,ε .
11 ∪Ω21
N wx 0 ;
Assumption 1 and the choice of λy imply that (E22 )x 0,Ω
1−λ y 1 11 ∪Ω21
C
1/2 ε
y=0 x=0
−2 −2β1 (1−x)/ε −2β2 (1−y)/ε
e
e
dy
21
N
CN −(p+1) ,
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1799
while a calculation like (4.6) yields (IN E22 )x 0,Ω11 CN −(p+1) , and (IN E22 )x C(εN −1 ln N )−1 IN E22 0,Ω21 C(ε ln N )−1/2 N −(p+1/2) , 0,Ω 21
where we used (3.16) and an inverse estimate that can be proved in the standard way. Hence CN −(p+1/2) (ln−1/2 N )w N 1,ε . ε(E22 − IN E22 )x , w N x Ω\Ω 22
(4.9)
Adding (4.2)–(4.9) and similar estimates for ε|((E − IN E)y , wyN )| completes the proof. Here E stands for E12 , E21 , and E22 . 2 Remark 11. If instead of the improved approximation estimates of Lemma 4 we simply apply Cauchy–Schwarz inequalities and the anisotropic bounds of Lemmas 6 (with r = p) and 7, then the conclusion of Lemma 10 becomes ε∇(u − IN u), ∇w N C(N −1 ln N )p w N 1,ε . (4.10) Now we can bound BGal (IN u − u, w N ). Lemma 12. Let Assumption 1 hold. Then there exists a constant C such that for all w N ∈ V N we have BGal (IN u − u, w N ) CN −(p+1/2) w N SD . Proof. Now
BGal (IN u − u, w N ) = ε(∇IN u − u, ∇w N ) + b · ∇(IN u − u) + c(IN u − u), w N .
Here b · ∇(IN u − u) + c(IN u − u), w N = IN u − u, (c − div b)w N − b · ∇w N CIN u − u0 w N 0 + IN u − u0,Ω11 b · ∇w N 0,Ω11 + IN u − uL∞ (Ω\Ω11 ) b · ∇w N L1 (Ω\Ω11 ) C(N −1 ln N )p+1 w N 0 + CN −(p+1) b · ∇w N 0,Ω11 + C(N −1 ln N )p+1 (ε 1/2 ln1/2 N )b · ∇w N 0,Ω\Ω11 ,
(4.11)
where we invoked Lemma 8 and applied the Cauchy–Schwarz inequality on the domain Ω \ Ω11 , whose measure is bounded by Cε ln N . Now consider separately the cases ε N −1 and ε < N −1 when bounding the second term in (4.11). For ε N −1 , one obtains CN −(p+1) b · ∇w N 0,Ω11 CN −(p+1/2) ε 1/2 |w N |1,Ω11 CN −(p+1/2) w N SD , while in the case ε < N −1 (for which δK = C1 N −1 ) we estimate 1/2 CN −(p+1) b · ∇w N 0,Ω11 CN −(p+1/2) δK b · ∇w N 20,K CN −(p+1/2) w N SD . K⊂Ω11
Hence from (4.11) we obtain b · ∇(IN u − u) + c(IN u − u), w N C (N −1 ln N )p+1 + N −(p+1/2) + (N −1 ln N )p+1 ln1/2 N w N SD CN −(p+1/2) w N SD . Combine this inequality with the result of Lemma 10 to complete the proof.
(4.12) 2
Lemma 13. Let Assumption 1 hold. Then there exists a constant C such that for all w N ∈ V N we have BStab (IN u − u, w N ) CN −(p+1/2) w N SD .
1800
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
Proof. Let w N ∈ V N be arbitrary but fixed. The contributions from the smooth and layer components of u to |BStab (IN u − u, w N )| will be analysed separately. Writing E for E12 , E21 or E22 , we have BStab (IN E − E, w N ) = δK −ε(IN E − E) + b · ∇(IN E − E) + c(IN E − E), b · ∇w N K . K⊂Ω11
Let ε
N −1
and δK C1 N −1 . Then by standard approximation estimates we have 1/2 1/2 BStab (IN E − E, w N ) C δK (εN −(p−1) + N −p + N −(p+1) )Ep+1,K δK b · ∇w N 0,K so εδK = C1
N −2
K⊂Ω11
C(ε 1/2 N −p + N −(p+1/2) )Ep+1,Ω11 w N SD . Using Assumption 1 and the choices of the transition points we have 1/2 −2(p+1) −2β1 (1−x)/ε E21 p+1,Ω11 C ε e dx dy Cε −1/2 (εN )−p N −1 Cε −1/2 N −1 . Ω11
Analogously, we obtain E12 p+1,Ω11 Cε −1/2 N −1 and E22 p+1,Ω11 Cε −1 N −(p+3) Cε −1/2 N −(p+5/2) . Thus for ε N −1 we get BStab (IN E − E, w N ) C N −(p+1) + (εN )−1/2 N −(p+1) w N SD CN −(p+1) w N SD . < N −1 ,
(4.13)
N −1 ,
so δK = C1 and use the splitting Now we consider the case ε δK −ε(IN E − E) + b · ∇(IN E − E) + c(IN E − E), b · ∇w N K K⊂Ω11
C1 N −1 εEL1 (Ω11 ) + ∇EL1 (Ω11 ) b · ∇w N L∞ (Ω11 ) + C1 N −1 ε (IN E)0,Ω + ∇(IN E)0,Ω + IN E − E0,Ω11 b · ∇w N 0,Ω11 . 11
11
Here εEL1 (Ω11 ) + ∇EL1 (Ω11 ) CN −(p+1) by (3.13). Inverse estimates and (3.15) yield ε (IN E)0,Ω + ∇(IN E)0,Ω C(εN + 1)∇(IN E)0,Ω CN −(p+1/2) , 11
11
11
and from (3.14) one has E − IN E0,Ω Taking b˜ to be the piecewise constant approximation of b that equals b at the barycentre of each K ∈ Tn ∩ Ω11 , ˜ · ∇w N ∞ b · ∇w N L∞ (K) (b − b) + b˜ · ∇w N L∞ (K) L (K) C(N −1 ln N )p+1 .
CN −1 ∇w N L∞ (K) + CN b˜ · ∇w N 0,K C∇w N 0,K + CN (b˜ − b) · ∇w N 0,K + b · ∇w N 0,K CN b · ∇w N 0,Ω11 + w N 0,Ω11 , ˜ CN −1 , inverse inequalities, and the equivalence of norms on finite-dimensional spaces after where we used |b − b| ˆ Consequently mapping b˜ · ∇w N to the reference element K. BStab (IN E − E, w N ) CN −1 N −(p+1) N b · ∇w N 0,Ω + w N 0,Ω + N −(p+1/2) b · ∇w N 0,Ω 11 11 11 CN −(p+1/2) w N SD .
(4.14)
Turning now to the smooth part S of u, εδK (IN S − S), b · ∇w N K + δK b · ∇(S − S I ), b · ∇w N K BStab (IN S − S, w N ) = K⊂Ω11
+
K⊂Ω11
δK c(S − S I ), b · ∇w N K .
K⊂Ω11
(4.15)
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1801
Using the Cauchy–Schwarz inequality and a standard interpolation bound [2, Theorem 3.1.5] we get N εδK (IN S − S), b · ∇w K ε 1/2 δK (IN S − S)0,Ω ε 1/2 w N 1,Ω11 11
K⊂Ω11
CN −(p+1/2) w N SD , )1/2
CN −1
(4.16)
CN −1
since (εδK and δK from the definition of δK . For the second term in (4.15), one has I N Cδ 1/2 N −p wN SD CN −(p+1/2) wN SD . b · ∇(S − S δ ), b · ∇w K K K
(4.17)
K⊂Ω11
Finally, the third term in (4.15) is bounded by I N Cδ 1/2 N −(p+1) w N SD CN −(p+3/2) wN SD . c(S − S δ ), b · ∇w K K K
(4.18)
K⊂Ω11
Combining (4.13)–(4.18), the conclusion of the lemma follows.
2
Remark 14. The argument used in [12] to squeeze a higher order of convergence from K⊂Ω11 εδK ((IN S − S), b · ∇w N )K cannot be used here, since it leads to terms containing jumps in (IN S) across element boundaries. In the case p = 1 of [12] these terms were all zero.
Furthermore, when p > 1 we do not have any analogue of Lemma 4 for mixed-derivative terms such as K (IK w − w)x wyN dx dy. This is another difference from the case p = 1 of [12]. It is now straightforward to prove our main result. Theorem 15. Let Assumption 1 hold. Then the SDFEM solution uN satisfies uN − IN uSD CN −(p+1/2) .
(4.19)
Proof. By (2.7) and (2.4), we have 1 N u − IN u2SD BSD (uN − IN u, uN − IN u) = BSD (u − IN u, uN − IN u) 2 = BGal (u − IN u, uN − IN u) + BStab (u − IN u, uN − IN u). Now invoke Lemmas 12 and 13 to complete the proof.
2
5. Conclusions In this paper higher-order finite element methods (Qp with p > 1) on Shishkin meshes for convection–diffusion problems have been analysed for the first time. The main result, Theorem 15, uses a special approximant of Lin to yield an order of convergence 1/2 higher than can be obtained by standard approaches. Nevertheless this theorem is inferior to the case p = 1 considered in [12], where a gain of an entire order of convergence over the standard analysis is achieved. What is the reason for this inconsistency? It turns out that when p = 1, for each K ∈ TN an identity of Lin [6,14] expresses (b · ∇(IN u − u), w N )K in terms of values of u, w N and their derivatives on the boundary of K, with a higher-order remainder. Then the lower-order terms for neighbouring K can be combined to deduce a higher-order estimate. The discussion in [13, Remark 20] reveals that in the one-dimensional case when Ω is an interval, analogous identities exist also for p > 1, but now the lower-order terms for neighbouring K do not combine in a productive way. This is the fundamental reason for the unexpected result that the case p = 1 is exceptional. The gain of a half order of convergence in Theorem 15 can be exploited in practice. It shows that IN u is “superclose” to u, i.e. that uN − IN uSD exhibits a higher order of convergence than uN − uSD does in general. Consequently a straightforward generalization of certain techniques from [12,13] works here: an inexpensive local postprocessing of uN will yield a piecewise Qp+1 solution u˜ N such that u˜ N − uSD CN −(p+1/2) .
1802
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
References [1] T. Apel, M. Dobrowolski, Anisotropic interpolation with applications to the finite element method, Computing 47 (1992) 277–293. [2] P. Ciarlet, The Finite Element Method for Elliptic Problems, SIAM, Philadelphia, PA, 2002. [3] P. Farrell, A. Hegarty, J. Miller, E. O’Riordan, G. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman & Hall/CRC, Boca Raton, FL, 2000. [4] C. Johnson, A. Schatz, L. Wahlbin, Crosswind smear and pointwise errors in the streamline diffusion finite element methods, Math. Comp. 49 (1987) 25–38. [5] Q. Lin, N. Yan, A. Zhou, A rectangle test for interpolated finite elements, in: Proc. Syst. Sci. Eng., Great Wall (H.K.) Culture Publish Co., 1991, pp. 217–229. [6] T. Linß, Uniform superconvergence of a Galerkin finite element method on Shishkin-type meshes, Numer. Methods Partial Differential Equations 16 (2000) 426–440. [7] T. Linß, M. Stynes, Asymptotic analysis and Shishkin-type decomposition for an elliptic convection–diffusion problem, J. Math. Anal. Appl. 261 (2001) 604–632. [8] J. Miller, E. O’Riordan, G. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [9] H. Roos, Layer-adapted grids for singular perturbation problems, Z. Angew. Math. Mech. 78 (1998) 291–309. [10] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, Berlin, 1996. [11] M. Stynes, L. Tobiska, A finite difference analysis of a streamline diffusion method on a Shishkin mesh, Numer. Algorithms 18 (1998) 337–360. [12] M. Stynes, L. Tobiska, The SDFEM for a convection–diffusion problem with a boundary layer: optimal error analysis and enhancement of accuracy, SIAM J. Numer. Anal. 41 (2003) 1620–1642. [13] L. Tobiska, Analysis of a new stabilized higher order finite element method for advection–diffusion equations, Comput. Methods Appl. Mech. Engrg. 196 (2006) 538–550. [14] Z. Zhang, Finite element superconvergence on Shishkin mesh for 2-d convection–diffusion problems, Math. Comp. 72 (2003) 1147–1177.
Using rectangular Qp elements in the SDFEM for a convection–diffusion problem with a boundary layer ✩ Martin Stynes a,∗ , Lutz Tobiska b a National University of Ireland, Cork, Ireland b Institut für Analysis und Numerik, Otto-von-Guericke Universität Magdeburg, Postfach 4120, D-39016 Magdeburg, Germany
Available online 23 November 2007
Abstract The streamline diffusion finite element method (SDFEM; the method is also known as SUPG) is applied to a convection– diffusion problem posed on the unit square whose solution has exponential boundary layers. A rectangular Shishkin mesh is used. The trial functions in the SDFEM are piecewise polynomials that lie in the space Qp , i.e., are tensor products of polynomials of degree p in one variable, where p > 1. The error bound IN u − uN SD CN −(p+1/2) is proved; here uN is the computed SDFEM solution, IN u is chosen in the finite element space to be a special approximant of the true solution u, and ·SD is the streamline-diffusion norm. This result is compared with previously known results for the case p = 1. The error bound is a superclose result; uN can be enhanced using local postprocessing to yield a modified solution u˜ N for which u − u˜ N SD CN −(p+1/2) . © 2007 IMACS. Published by Elsevier B.V. All rights reserved. MSC: 65N15; 65N30 Keywords: Convection–diffusion; Streamline-diffusion method; Shishkin mesh
1. Introduction We consider the singularly perturbed boundary value problem Lu := −εu + b · ∇u + cu = f
on Ω = (0, 1)2 ,
u = 0 on ∂Ω,
(1.1)
where the parameter ε lies in the interval (0, 1]. Assume that the function b(x, y) = (b1 (x, y), b2 (x, y)) with ¯ and b1 (x, y) > β1 > 0 and b2 (x, y) > β2 > 0, c(x, y) 0 on Ω, b(x, y) ¯ (1.2) c(x, y) − div c0 > 0 on Ω, 2 where β1 , β2 and c0 are some constants. We assume that the functions b, c and f are sufficiently smooth. These hypotheses ensure that (1.1) has a unique solution in H01 (Ω) ∩ H 2 (Ω) for all f ∈ L2 (Ω). Note that for sufficiently small ε, the other hypotheses imply that (1.2) can always be ensured by the simple change of variable v(x, y) = e−σ x u(x, y) when σ is chosen suitably. ✩
We dedicate this paper to our colleague and friend Pieter Hemker on the occasion of his 65th birthday.
* Corresponding author.
E-mail addresses: [email protected] (M. Stynes), [email protected] (L. Tobiska). 0168-9274/$30.00 © 2007 IMACS. Published by Elsevier B.V. All rights reserved. doi:10.1016/j.apnum.2007.11.004
1790
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
The convection–diffusion problem (1.1) is a well-known 2-dimensional testbed for numerical methods for advection-dominated transport phenomena; it exhibits exponential boundary layer behaviour at the sides x = 1 and y = 1 of Ω that is typical of such problems. Thus layer-adapted meshes are usually used to solve (1.1). A piecewiseuniform Shishkin mesh can be chosen a priori when one has some understanding of the structure of these layers [3,8,10]. This mesh has O(N 2 ) points in a rectangular grid that is refined near the sides x = 1 and y = 1. In [12] we showed that for piecewise bilinears (the space Q1 ), on writing uI for the piecewise bilinear nodal interpolant to u, then the error uI − uN L2 (Ω) is second-order convergent (up to a logarithmic factor), uniformly in ε; this is the optimal rate of convergence. What about higher-order elements (the space Qp , where p > 1)? The present paper is the first to analyse higher-order finite element methods on Shishkin meshes for convection–diffusion problems. Errors will be estimated in both the streamline diffusion and L2 norms. Does one again attain the optimal order of convergence in L2 ? We answer this question and discover that p = 1 is atypical. The extension of the analysis of [12] to Qp is not routine: for example, instead of a standard interpolant from Qp , one should use the vertices– edges–element–approximant of Lin [5]. The non-standard anisotropic interpolation properties of this approximant are studied in Lemma 5 below. The organization of the paper is as follows. In Section 2 we describe the Shishkin mesh and the SDFEM. A decomposition of the solution u and some approximation-theoretical facts that we shall need later are presented in Section 3. The convergence of the SDFEM is analysed in Section 4. Notation. Throughout the paper C will denote a generic positive constant that is independent of ε and the mesh. Subscripted constants such as C1 are also independent of ε and the mesh, but have a fixed value. We use the standard Sobolev spaces W k,p (D), H k (D) = W k,2 (D), H0k (D), Lp (D) = W 0,p (D) for non-negative integers k and 1 p ∞, and write (·,·)D for the L2 (D) inner product. Here D is any measurable subset of Ω. Then |·|k,D and ·k,D are the usual Sobolev seminorm and norm on H k (D). When D = Ω, we drop D from the notation for simplicity. 2. The Shishkin mesh and the SDFEM Shishkin meshes [8–10] are piecewise uniform meshes, constructed a priori, that are refined inside layers. Let N be an even positive integer. We let λx and λy denote two mesh transition parameters that will be used to specify where the mesh changes from coarse to fine; these are defined by 1 (p + 1)ε 1 (p + 1)ε λx = min , ln N and λy = min , ln N , 2 β1 2 β2 where our trial space, which is defined below, comprises functions that are piecewise in Qp for some integer p > 1. In fact we assume that λx = ((p + 1)ε/(β1 )) ln N and λy = ((p + 1)ε/(β2 )) ln N , as otherwise N −1 is much smaller than ε (which is very unlikely in practice) and one can then simplify our analysis of the method. In these formulae, if the multiplier of ε/β1 and ε/β2 is chosen too small then the final order of convergence obtained will decrease. For the case p = 1 this dependence is investigated both theoretically and numerically in [11]. We divide Ω as in Fig. 1: Ω¯ = Ω11 ∪ Ω21 ∪ Ω12 ∪ Ω22 , where Ω11 = [0, 1 − λx ] × [0, 1 − λy ], Ω21 = [1 − λx , 1] × [0, 1 − λy ], Ω12 = [0, 1 − λx ] × [1 − λy , 1], Ω22 = [1 − λx , 1] × [1 − λy , 1]. ¯ i, j = 0, . . . , N} are the rectangular lattice defined by The mesh points Ω N = {(xi , yj ) ∈ Ω: for i = 0, . . . , N/2, 2i(1 − λx )/N xi = 1 − 2(N − i)λx /N for i = N/2 + 1, . . . , N, and
yj =
2j (1 − λy )/N 1 − 2(N − j )λy /N
for j = 0, . . . , N/2, for j = N/2 + 1, . . . , N.
Our mesh is constructed by drawing lines parallel to the coordinate axes through these mesh points, so it is a tensor product of two one-dimensional piecewise uniform meshes. This divides Ω into a set TN of mesh rectangles K whose sides are parallel to the axes—see Fig. 1. The mesh is coarse on Ω11 , coarse/fine on Ω21 ∪ Ω12 , and fine on Ω22 .
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1791
Fig. 1. Division of Ω and Shishkin mesh for T8 .
√ √ The mesh is quasiuniform on Ω11 and its diameter d there satisfies 2/N d 2 2/N ; on Ω12 ∪ Ω21 , each mesh rectangle has dimensions O(N −1 ) by O(εN −1 ln N ); and on Ω22 each rectangle is O(εN −1 ln N ) by O(εN −1 ln N ). We shall use these properties several times in our analysis. Given a mesh rectangle K, its dimensions are written as hx,K by hy,K and its barycentre is denoted by (xK , yK ). We now describe the SDFEM on this rectangular mesh. Let p be a fixed integer with p > 1. (In [12] we analysed the case p = 1.) Our trial space is ¯ v|∂Ω = 0 and v|K ∈ Qp (K) ∀K ∈ TN , V N = v ∈ C(Ω): where Qp (K) = span{x i y j : 0 i, j p}. Define the bilinear form BSD (·,·) by BSD (w, v) = BGal (w, v) + BStab (w, v), where BGal (w, v) = ε(∇w, ∇v) + (b · ∇w + cw, v), δK (−εw + b · ∇w + cw, b · ∇v)K , BStab (w, v) =
(2.1) (2.2)
K∈TN
where H˜ 1 (Ω) denotes the set of functions in H 1 (Ω) that lie in H 2 (K) for each K, for all and δK 0 is a user-chosen piecewise constant parameter. The SDFEM is defined as follows: Find uN ∈ V N such that BSD (uN , v N ) = (f, v N ) + δK (f, b · ∇v N )K ∀v N ∈ V N . (2.3) (w, v) ∈ H˜ 1 (Ω) × H 1 (Ω),
K∈TN
We clearly have the Galerkin orthogonality property BSD (u − uN , v N ) = 0 ∀v N ∈ V N .
(2.4)
The choice of δK is discussed in, e.g., [10, §III.3.2.1]. For each K ∈ TN , set hK = min{hx,K , hy,K }. Let μ be a constant such that the inverse inequality N v N 0,K μh−1 K ∇v 0,K
∀v N ∈ V N , ∀K ⊂ Ω11
is valid. Similarly to [10, p. 233], we set ⎧ if K ⊂ Ω11 and ε N −1 , ⎨ C1 N −1 −1 −2 δK = C1 ε N if K ⊂ Ω11 and ε > N −1 , ⎩ 0 otherwise, where the positive constant C1 is chosen (independently of ε) such that h2K c0 1 ∀K ⊂ Ω11 . , 0 δK min 2 (maxK |c(x, y)|)2 μ2 ε
(2.5)
(2.6)
1792
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
Then the argument of [10, p. 231] shows that the inequality 1 BSD (v N , v N ) v N 2SD ∀v N ∈ V N 2 holds, where ·SD is given by 1/2 2 2 2 vSD = ε|v|1 + δK b · ∇v0,K + c0 v0
(2.7)
∀v ∈ H 1 (Ω).
(2.8)
K∈TN
It follows that (2.3) has a unique solution uN . We shall also use the ε-weighted energy norm v1,ε = (ε|v|21 + c0 v20 )1/2 for all v ∈ H 1 (Ω). This norm is clearly weaker than ·SD . 3. Solution decomposition, Lin’s approximant We now assume that u can be decomposed in a precise way that is typical of the behaviour observed in solutions of (1.1) when interior layers are absent. Assumption 1. Assume that u = S + E21 + E12 + E22 , where there exists a constant C such that for all (x, y) ∈ Ω we have i+j ∂ S ∂x i ∂y j (x, y) C for 0 i + j p + 1, i+j ∂ E21 Cε −i e−β1 (1−x)/ε for 0 i + j p + 2, (x, y) ∂x i ∂y j i+j ∂ E12 Cε −j e−β2 (1−y)/ε for 0 i + j p + 2, (x, y) ∂x i ∂y j i+j ∂ E22 −(i+j ) −(β1 (1−x)+β2 (1−y))/ε e for 0 i + j p + 2. ∂x i ∂y j (x, y) Cε
(3.1)
(3.2) (3.3) (3.4) (3.5)
Furthermore, assume that S ∈ H p+2 (Ω) with Sp+2 C.
(3.6)
Here S is the smooth part of u, E21 is an exponential boundary layer along the side x = 1 of Ω, E12 is an exponential boundary layer along the side y = 1, while E22 is an exponential corner layer at (1, 1). Remark 2. In [7] a proof is given that under certain compatibility conditions on the data f , the bounds (3.2)–(3.6) of Assumption 1 are valid if p = 1. The extension of these results to our case p > 1 is non-trivial but an inspection of the argument in [7] convinces one that it can be done by placing a sufficient number of compatibility conditions on f . The number of these conditions increases rapidly with p, which is restrictive. Nevertheless in the case of a problem that did not satisfy all these assumptions, one could combine the argument of the current paper with cut-off functions to obtain bounds on the error in the computed solution in regions that were sufficiently far from all layers induced by inadequate compatibility of the data, much like the analysis in [4]. ¯ we shall follow the approach of Lin [5] in defining a non-standard “vertices–edges– Given any function v ∈ C(Ω), element approximant” IK v ∈ Qp on each mesh rectangle K ∈ TN . We use this approximant because one then has, compared with standard interpolants, sharper bounds (see Lemma 4 and Remark 11 below) for the errors that appear in various terms in the finite element analysis. Lin’s approximant is based on nodal values of v at the corners of K together with moments of v on the edges of K and on K itself. It is initially defined on the reference domain Kˆ = (−1, +1)2 , whose vertices and edges we
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1793
ˆ The approximation operator Iˆ : C(K) ˆ → Qp (K) ˆ denote by aˆ i and eˆi respectively for i = 1, . . . , 4. Let v(·,·) ˆ ∈ C(K). is defined by requiring ˆ aˆ i ) for i = 1, . . . , 4, (Iˆv)( ˆ aˆ i ) = v( (Iˆv)q ˆ dγ = vq ˆ dγ ∀q ∈ Pp−2 (eˆi ) for i = 1, . . . , 4, eˆi
eˆi
(Iˆv)q ˆ dξ dη =
Kˆ
vq ˆ dξ dη
(3.7a) (3.7b)
ˆ ∀q ∈ Qp−2 (K);
(3.7c)
Kˆ
here Pp−2 (eˆi ) is the space of polynomials of degree at most p − 2 in the single variable whose axis is parallel to the edge eˆi . Lemma 3. The approximation operator Iˆ is uniquely defined. Proof. There is a total of (p + 1)2 conditions on Iˆvˆ and we have (p + 1)2 degrees of freedom for Iˆvˆ ∈ Qp . Thus it is enough to show that taking all degrees of freedom equal to zero implies Iˆvˆ = 0. ˆ the function rˆ = Iˆv| Assume that all degrees of freedom are equal to zero. Then for any edge eˆ of K, ˆ eˆ is a 2 polynomial of degree p in one variable on (−1, 1) that is L -orthogonal to the space Pp−2 (e) ˆ of polynomials of p degree p − 2. The L2 -orthogonality of the Legendre polynomials {Lm }m=0 implies that rˆ = ap−1 Lp−1 + ap Lp with unknown constants ap−1 and ap . In addition, rˆ must vanish at the endpoints of e, ˆ i.e., rˆ (±1) = 0, which implies that rˆ ≡ 0. Applying this argument for each edge of Kˆ we see that Iˆvˆ vanishes on the boundary of Kˆ and therefore (Iˆv)(ξ, ˆ η) = (1 − ξ 2 )(1 − η2 )q(ξ, ˆ η)
ˆ with qˆ ∈ Qp−2 (K).
But our assumption that all degrees of freedom equal zero requires in particular that in Iˆvˆ ≡ 0. 2
ˆ
ˆ qˆ dξ Kˆ (I v)
dη = 0, which results
Using the transformation x = xK + hx,K ξ/2, y = yK + hy,K η/2 to map from Kˆ to an arbitrary K ∈ TN , one obtains the corresponding approximation operator IK : C(K) → Qp (K). Then a continuous global approximation ¯ → V N is defined by setting operator IN : C(Ω) (IN v)|K := IK (v|K )
∀K ∈ TN .
Lemma 4. Let K ∈ TN . Let w ∈ H p+2 (K) and let IK w ∈ Qp (K) be its vertices–edges–element approximant. Then for each v N ∈ Qp (K) we have p+2 w v N 0,K (IK w − w)x v N dx dy Chp+1 ∂ x x y,K p+1 ∂x∂y 0,K K
and
p+2 w (IK w − w)y v N dx dy Chp+1 ∂ y x,K p+1 ∂x ∂y K
Proof. For each mesh rectangle K, we set hy,K 2 1 2 . FK (y) = (y − yK ) − 2 2 Now, from the proof of [5, Lemma 4] it follows that
0,K
vyN 0,K .
1794
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
(IK w − w)x vxN dx dy = K
(−2)p (2p + 1)! +
(−2)p
FK (y)FK (y) p
∂ p+2 w(x, y) ∂ p+1 v N (x, yK ) dx dy ∂x∂y p ∂x∂y p+1
K
p
FK (y)
(2p)!
∂ p+2 w(x, y) ∂ p v N (x, yK ) dx dy. ∂x∂y p+1 ∂x∂y p−1
K
Clearly p F (y)
K
h2y,K
p
8
hy,K , and FK (y) 2
and for arbitrary functions v N ∈ Qp (K) one has the identities ∂ p+1 v N (x, yK ) ∂ p+1 v N (x, y) = ∂x∂y p ∂x∂y p and ∂ p v N (x, yK ) ∂ p v N (x, y) ∂ p+1 v N (x, y) = − (y − yK ) p−1 p−1 ∂x∂y p ∂x∂y ∂x∂y as well as the inverse inequalities n+1 N N ∂v v n ∂ C hy,K ∂x n ∂x∂y 0,K 0,K
for n = p − 1, p.
Invoking these above, we get the first statement of lemma. The second follows analogously.
2
To obtain satisfactory approximation error estimates on Ω12 ∪ Ω21 , where the Shishkin mesh is highly anisotropic, we shall use some technical results of Apel and Dobrowolski [1]. Let r be an integer with 1 r p. Define the anisotropic space i+j +1 ∂ vˆ r ˆ 1 ˆ < ∞, 0 i + j r Hξ (K) = vˆ ∈ L (K): i+1 j ∂ξ ∂η 0,Kˆ ˆ similarly. and define Hηr (K) Lemma 5. Let r be an integer with 1 r p. There exists a constant C independent of vˆ such that for the vertex– edges–element approximant Iˆvˆ Dξ (vˆ − Iˆv) ˆ ∩ Hξr (K). ˆ ˆ 0,Kˆ C|Dξ v| ˆ r,Kˆ ∀vˆ ∈ C(K) (3.8) ˆ are replaced by Dη and Hηr (K), ˆ respectively. An analogous inequality holds if Dξ and Hξr (K) ˆ → R by Proof. Define r + r 2 continuous linear functionals Fi : H r (K) vˆ → ξ vˆ dξ dη for = 0, . . . , r − 1, Kˆ
vˆ →
ξ j ηk Kˆ
∂ vˆ dξ dη ∂η
for j = 0, . . . , r − 1 and k = 0, . . . , r − 1.
Then, recalling the definition of Iˆ, it is straightforward to verify that ˆ ∩ Hξr (K) ˆ for i = 1, . . . , r + r 2 ˆ = 0 ∀vˆ ∈ C(K) Fi Dξ (Iˆvˆ − v) and also that ˆ = 0 for i = 1, . . . , r + r 2 Fi (Dξ v)
ˆ implies Dξ vˆ = 0 ∀vˆ ∈ Qr (K).
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1795
Now inequality (3.8) follows immediately from [1, Lemma 3] by setting (in the notation of that paper) Q = {(i, j ): 0 i, j r − 1}, P = {(i, j ): 0 i, j r}, and γ = (1, 0). For the analogous inequality involving Dη , modify the above linear functionals by interchanging ξ with η and replacing the η-derivative by the ξ -derivative, and take γ = (0, 1) in [1, Lemma 3]. 2 The anisotropic estimate of Lemma 5 implies the following sharp bounds. Lemma 6. Let r be an integer with 1 r p. Let K ∈ TN . Then there exists a constant C such that the vertices– edges–element approximant IN v satisfies r+1 j i ∂ v , (v − IN v)x C h h x,K y,K 0,K i+1 j ∂x ∂y 0,K i+j =r r+1 ∂ v j i (v − IN v)y C hx,K hy,K ∂x i ∂y j +1 0,K 0,K i+j =r
for all v ∈ H r+1 (K). Proof. Using the transformation x = xK +
hx,K ξ, 2
y = yK +
hy,K η, 2
ˆ apply the anisotropic result of Lemma 5, then map back to the we map the element K onto the reference element K, original element K. Thus 2|K|1/2 |K|1/2 (v − IN v)x (vˆ − Iˆv) ˆ ξ 0,Kˆ C |Dξ v| ˆ r,Kˆ 0,K hx,K hx,K i+j +1 i+j +1 ∂ ∂ v v |K|1/2 j i C C hx,K hy,K i+1 j i+1 h ∂ξ ∂η ∂x ∂y j ˆ x,K i+j =r
0,K
The proof of the second inequality is along the same lines.
i+j =r
.
0,K
2
As well as Lemma 6, we have the following more standard anisotropic approximation error estimate. Lemma 7. Let K ∈ TN . Let q ∈ [1, ∞]. Then there exists a constant C such that the vertices–edges–element approximant IN v satisfies p+1 ∂ v j i hx,K hy,K v − IN vLq (K) C ∂x i ∂y j q L (K) i+j =p+1
for all v ∈ W p+1,q (K). Proof. Lemma 4 of [1] implies that we can invoke [1, Lemma 3] with Q = P = {(i, j ): 0 i, j p} and γ = ˆ (0, 0). This bounds v − IN vLq (K) ˆ , and the desired result follows by mapping between K and K as in the proof of Lemma 6. 2 For the Lin approximant one has the stability bound IN vL∞ (K) CvL∞ (K)
(3.9)
¯ where the constant C is independent of v and K ∈ Tn . This can be proved by mapping (3.7) to the for all v ∈ C(K), ˆ reference element K. From Lemma 7 one can deduce a pointwise bound on the error in the approximant.
1796
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
Lemma 8. Let Assumption 1 hold. Then there exists a constant C such that −(p+1) if (x, y) ∈ Ω11 , (u − IN u)(x, y) CN C(N −1 ln N )p+1 otherwise. Proof. Using the decomposition (3.1), write u − IN u = (S − IN S) + (E12 − IN E12 ) + (E21 − IN E21 ) + (E22 − IN E22 ). By Lemma 7 (with q = ∞) and Assumption 1, (S − IN S)(x, y) CN −(p+1) S p+1,∞ CN −(p+1) W Similarly, on Ω12 ∪ Ω22 one gets (E12 − IN E12 )(x, y) C
(3.10)
∀(x, y) ∈ Ω.
N −i (εN −1 ln N )j ε −j C(N −1 ln N )p+1 ;
i+j =p+1
on Ω11 ∪ Ω21 , the choice of transition point 1 − λy in the Shishkin mesh and (3.9) yield IN E12 (x, y) C max E12 (x, y) CN −(p+1) . Ω11 ∪Ω21
The functions E21 and E22 are handled similarly. Putting the various bounds into (3.10) gives the desired result.
2
The next lemma collects several results that involve the approximants of the various terms in the decomposition (3.1). Lemma 9. Let Assumption 1 hold. Let E denote any one of the functions E12 , E21 or E22 . Then there exists a constant C such that the following approximation error estimates hold: S − IN S0,Ω CN −(p+1) , E0,Ω11 Cε
1/2
N
−(p+1)
(3.11) (3.12)
,
εEL1 (Ω11 ) + ∇EL1 (Ω11 ) CN
−(p+1)
(3.13)
,
E − IN E0,Ω C(N −1 ln N )p+1 ,
(3.14) −1 1/2
IN E12 0,Ω11 ∪Ω21 + IN E21 0,Ω11 ∪Ω12 C(ε + N ) 1/2 −(p+1) N . IN E22 0,Ω12 ∪Ω21 C ε(ε + N −1 ) ln N
N
−(p+1)
(3.15)
,
(3.16)
Proof. Imitate [12, Lemma 3.2] to prove (3.11)–(3.14). For (3.15) we shall give only the argument to bound IN E12 0,Ω11 ∪Ω21 as IN E21 0,Ω11 ∪Ω12 is handled analogously. Using (3.9) and the decay property (3.4) of E12 , IN E12 Ω11 ∪Ω21 =
N/2
yj
1
2 IN E12 (x, y) dx dy
1/2
j =1 y=yj −1 x=0
C
N/2
yj
1/2
1 e
−2β2 (1−yj )/ε
dx dy
j =1 y=yj −1 x=0
C
N/2−1 y j +1 1
e
−2β2 (1−y)/ε
y N/2
dx dy +
j =1 y=yj x=0
1/2 C (ε + N −1 )N −2(p+1) , by the definition of λy . Similarly, using (3.5) instead of (3.4),
1/2
1 e
−2β2 λy /ε
dx dy
y=yN/2−1 x=0
(3.17)
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
IN E22 Ω21 C
N/2−1 y j +1
1 e
−2β2 (1−y)/ε
y N/2
e
−2β2 λy /ε
dx dy
y=yN/2−1 x=1−λx
1/2 C ε(ε + N −1 )N −2(p+1) ln N . The bound on IN E22 Ω12 is analogous.
1/2
1
dx dy +
j =1 y=yj x=1−λ x
1797
(3.18)
2
4. Error bound for IN u − uN In the analysis of this section, we shall try to give the sharpest possible bound for each term in order to reveal the source of the essential difference between the order of accuracy attained by Q1 elements and that attained by Qp elements when p > 1. The extreme nonuniformity of the mesh forces us to employ different analytical techniques on the different regions of Fig. 1. The proof of our main result (Theorem 15) requires a bound on BSD (uN − IN u, uN − IN u). To achieve this, write BSD = BGal + BStab , then bounds on the individual terms in each of these bilinear forms are proved in Lemmas 10–13. Lemma 10. Let Assumption 1 hold. Then there exists a constant C such that for all w N ∈ V N we have ε∇(u − IN u), ∇w N C N −(p+1) + ε 1/2 (N −1 ln N )p+1 w N 1,ε .
(4.1)
Proof. Let w N ∈ V N be arbitrary but fixed. As in (3.10), we write u − IN u = (S − IN S) + (E12 − IN E12 ) + (E21 − IN E21 ) + (E22 − IN E22 ). Then consider each corresponding term in (ε∇(u − IN u), ∇w N ). First, Lemma 4 implies that ε(S − IN S)x , w N CεN −(p+1) Sp+2 w N 0 Cε 1/2 N −(p+1) w N 1,ε , x x since Sp+2 C by Assumption 1. Next, consider ε(E12 − IN E12 )x , w N ε(E12 − IN E12 )x , w N x x Ω
12 ∪Ω22
(4.2)
+ ε(E12 − IN E12 )x , w N
x Ω11 ∪Ω21
.
The calculation on Ω12 ∪ Ω22 resembles that for S: ε(E12 − IN E12 )x , w N
x Ω12 ∪Ω22
Cε(εN −1 ln N )p+1
ε
−2(p+1) −2β2 (1−y)/ε
e
1/2 dx dy
wxN 0
Ω12 ∪Ω22
Cε(N
−1
ln N )
On Ω11 ∪ Ω21 we use Lemma 6 with r = 1: ε(E12 − IN E12 )x , w N x Ω11 ∪Ω21 CεN −1 (E12 )xx 0,Ω ∪Ω + (E12 )xy 0,Ω 11
CεN
−1
21
p+1
w N 1,ε .
11 ∪Ω21
(4.3)
N wx 0
1/2 1−λ y −2 −2β2 (1−y)/ε ε e dy wxN 0 CN −(p+2) w N 1,ε ,
(4.4)
y=0
by the choice of transition point 1 − λy . The E21 integral is also broken into two regions as ε(E21 − IN E21 )x , w N ε(E21 − IN E21 )x , w N x x Ω
21 ∪Ω22
+ ε(E21 − IN E21 )x , w N
For the first term on the right-hand side here, we again use Lemma 4 to get
x Ω11 ∪Ω12
.
1798
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
ε(E21 − IN E21 )x , w N x Ω21 ∪Ω22 p+2 ∂ E21 CεN −(p+1) wxN 0,Ω21 ∪Ω22 ∂x∂y p+1 0,Ω21 ∪Ω22 1/2 −(p+1) −2 −2β1 (1−x)/ε CεN ε e dx dy wxN 0 CN −(p+1) w N 1,ε .
(4.5)
Ω21 ∪Ω22
The other term will be handled differently for ε N −1 and ε < N −1 , respectively. If ε N −1 we use Lemma 6 and the choice of transition point 1 − λx to get N ε(E21 − IN E21 )x , w N (E21 )xy CεN −1 (E21 )xx wx 0 + x Ω ∪Ω 0,Ω ∪Ω 0,Ω ∪Ω 11
12
11
12
11
12
Cε 1/2 N −1 (ε −3/2 + ε −1/2 )N −(p+1) ε 1/2 wxN 0 CN −(p+1) w N 1,ε . For ε < N −1 we have ε(E21 − IN E21 )x , w N
x Ω11 ∪Ω12
ε (E21 )x
0,Ω11 ∪Ω12
+ (IN E21 )x 0,Ω
11 ∪Ω12
N wx Ω11 ∪Ω12 .
Assumption 1 and the choice of transition point 1 − λx imply that 1/2 1−λ x 1/2 1/2 −2 −2β1 (1−x)/ε (E21 )x ε Cε ε e dx CN −(p+1) . 0,Ω11 ∪Ω12
x=0
Invoking an inverse inequality and (3.15), we have ε 1/2 (IN E21 )x Cε 1/2 NIN E21 Ω Ω11 ∪Ω12
11 ∪Ω12
CN −(p+1)
(4.6)
since ε < N −1 . Putting these inequalities together, CN −(p+1) w N 1,ε . ε(E21 − IN E21 )x , w N
(4.7)
x Ω11 ∪Ω12
The E22 term is handled like the other two exponential layer terms. First, + ε(E22 − IN E22 )x , w N ε(E22 − IN E22 )x , w N ε(E22 − IN E22 )x , w N x x Ω x Ω\Ω . 22
22
Then, similarly to (4.3), Lemma 4 yields
p+2 E22 Cε(εN −1 ln N )p+1 ∂ wxN 0 x Ω12 ∪Ω22 ∂x∂y p+1 0,Ω12 ∪Ω22 1/2 Cε(εN −1 ln N )p+1 ε −2(p+2) e−2β1 (1−x)/ε e−2β2 (1−y)/ε dx dy wxN 0
ε(E22 − IN E22 )x , w N
Ω12 ∪Ω22
Cε
1/2
(N
−1
p+1
ln N )
w N 1,ε .
(4.8)
N −1
< N −1
N −1 ,
and ε separately. For ε apply Lemma 6 with r = 1 to On Ω11 ∪ Ω21 we consider the cases ε get ε(E22 − IN E22 )x , w N ε 1/2 N −1 (E22 )xx + (E22 )xy 0,Ω ∪Ω w N 1,ε x Ω ∪Ω 0,Ω ∪Ω 11
21
11
ε For ε
< N −1
we use the splitting ε(E22 − IN E22 )x , w N x Ω
11 ∪Ω21
1/2
−1
(εN )
N
−(p+1)
ε (E22 )x
21
w 1,ε ε
0,Ω11 ∪Ω21
N
11
1/2
N
−(p+1)
+ (IN E22 )x 0,Ω
w 1,ε .
11 ∪Ω21
N wx 0 ;
Assumption 1 and the choice of λy imply that (E22 )x 0,Ω
1−λ y 1 11 ∪Ω21
C
1/2 ε
y=0 x=0
−2 −2β1 (1−x)/ε −2β2 (1−y)/ε
e
e
dy
21
N
CN −(p+1) ,
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1799
while a calculation like (4.6) yields (IN E22 )x 0,Ω11 CN −(p+1) , and (IN E22 )x C(εN −1 ln N )−1 IN E22 0,Ω21 C(ε ln N )−1/2 N −(p+1/2) , 0,Ω 21
where we used (3.16) and an inverse estimate that can be proved in the standard way. Hence CN −(p+1/2) (ln−1/2 N )w N 1,ε . ε(E22 − IN E22 )x , w N x Ω\Ω 22
(4.9)
Adding (4.2)–(4.9) and similar estimates for ε|((E − IN E)y , wyN )| completes the proof. Here E stands for E12 , E21 , and E22 . 2 Remark 11. If instead of the improved approximation estimates of Lemma 4 we simply apply Cauchy–Schwarz inequalities and the anisotropic bounds of Lemmas 6 (with r = p) and 7, then the conclusion of Lemma 10 becomes ε∇(u − IN u), ∇w N C(N −1 ln N )p w N 1,ε . (4.10) Now we can bound BGal (IN u − u, w N ). Lemma 12. Let Assumption 1 hold. Then there exists a constant C such that for all w N ∈ V N we have BGal (IN u − u, w N ) CN −(p+1/2) w N SD . Proof. Now
BGal (IN u − u, w N ) = ε(∇IN u − u, ∇w N ) + b · ∇(IN u − u) + c(IN u − u), w N .
Here b · ∇(IN u − u) + c(IN u − u), w N = IN u − u, (c − div b)w N − b · ∇w N CIN u − u0 w N 0 + IN u − u0,Ω11 b · ∇w N 0,Ω11 + IN u − uL∞ (Ω\Ω11 ) b · ∇w N L1 (Ω\Ω11 ) C(N −1 ln N )p+1 w N 0 + CN −(p+1) b · ∇w N 0,Ω11 + C(N −1 ln N )p+1 (ε 1/2 ln1/2 N )b · ∇w N 0,Ω\Ω11 ,
(4.11)
where we invoked Lemma 8 and applied the Cauchy–Schwarz inequality on the domain Ω \ Ω11 , whose measure is bounded by Cε ln N . Now consider separately the cases ε N −1 and ε < N −1 when bounding the second term in (4.11). For ε N −1 , one obtains CN −(p+1) b · ∇w N 0,Ω11 CN −(p+1/2) ε 1/2 |w N |1,Ω11 CN −(p+1/2) w N SD , while in the case ε < N −1 (for which δK = C1 N −1 ) we estimate 1/2 CN −(p+1) b · ∇w N 0,Ω11 CN −(p+1/2) δK b · ∇w N 20,K CN −(p+1/2) w N SD . K⊂Ω11
Hence from (4.11) we obtain b · ∇(IN u − u) + c(IN u − u), w N C (N −1 ln N )p+1 + N −(p+1/2) + (N −1 ln N )p+1 ln1/2 N w N SD CN −(p+1/2) w N SD . Combine this inequality with the result of Lemma 10 to complete the proof.
(4.12) 2
Lemma 13. Let Assumption 1 hold. Then there exists a constant C such that for all w N ∈ V N we have BStab (IN u − u, w N ) CN −(p+1/2) w N SD .
1800
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
Proof. Let w N ∈ V N be arbitrary but fixed. The contributions from the smooth and layer components of u to |BStab (IN u − u, w N )| will be analysed separately. Writing E for E12 , E21 or E22 , we have BStab (IN E − E, w N ) = δK −ε(IN E − E) + b · ∇(IN E − E) + c(IN E − E), b · ∇w N K . K⊂Ω11
Let ε
N −1
and δK C1 N −1 . Then by standard approximation estimates we have 1/2 1/2 BStab (IN E − E, w N ) C δK (εN −(p−1) + N −p + N −(p+1) )Ep+1,K δK b · ∇w N 0,K so εδK = C1
N −2
K⊂Ω11
C(ε 1/2 N −p + N −(p+1/2) )Ep+1,Ω11 w N SD . Using Assumption 1 and the choices of the transition points we have 1/2 −2(p+1) −2β1 (1−x)/ε E21 p+1,Ω11 C ε e dx dy Cε −1/2 (εN )−p N −1 Cε −1/2 N −1 . Ω11
Analogously, we obtain E12 p+1,Ω11 Cε −1/2 N −1 and E22 p+1,Ω11 Cε −1 N −(p+3) Cε −1/2 N −(p+5/2) . Thus for ε N −1 we get BStab (IN E − E, w N ) C N −(p+1) + (εN )−1/2 N −(p+1) w N SD CN −(p+1) w N SD . < N −1 ,
(4.13)
N −1 ,
so δK = C1 and use the splitting Now we consider the case ε δK −ε(IN E − E) + b · ∇(IN E − E) + c(IN E − E), b · ∇w N K K⊂Ω11
C1 N −1 εEL1 (Ω11 ) + ∇EL1 (Ω11 ) b · ∇w N L∞ (Ω11 ) + C1 N −1 ε (IN E)0,Ω + ∇(IN E)0,Ω + IN E − E0,Ω11 b · ∇w N 0,Ω11 . 11
11
Here εEL1 (Ω11 ) + ∇EL1 (Ω11 ) CN −(p+1) by (3.13). Inverse estimates and (3.15) yield ε (IN E)0,Ω + ∇(IN E)0,Ω C(εN + 1)∇(IN E)0,Ω CN −(p+1/2) , 11
11
11
and from (3.14) one has E − IN E0,Ω Taking b˜ to be the piecewise constant approximation of b that equals b at the barycentre of each K ∈ Tn ∩ Ω11 , ˜ · ∇w N ∞ b · ∇w N L∞ (K) (b − b) + b˜ · ∇w N L∞ (K) L (K) C(N −1 ln N )p+1 .
CN −1 ∇w N L∞ (K) + CN b˜ · ∇w N 0,K C∇w N 0,K + CN (b˜ − b) · ∇w N 0,K + b · ∇w N 0,K CN b · ∇w N 0,Ω11 + w N 0,Ω11 , ˜ CN −1 , inverse inequalities, and the equivalence of norms on finite-dimensional spaces after where we used |b − b| ˆ Consequently mapping b˜ · ∇w N to the reference element K. BStab (IN E − E, w N ) CN −1 N −(p+1) N b · ∇w N 0,Ω + w N 0,Ω + N −(p+1/2) b · ∇w N 0,Ω 11 11 11 CN −(p+1/2) w N SD .
(4.14)
Turning now to the smooth part S of u, εδK (IN S − S), b · ∇w N K + δK b · ∇(S − S I ), b · ∇w N K BStab (IN S − S, w N ) = K⊂Ω11
+
K⊂Ω11
δK c(S − S I ), b · ∇w N K .
K⊂Ω11
(4.15)
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
1801
Using the Cauchy–Schwarz inequality and a standard interpolation bound [2, Theorem 3.1.5] we get N εδK (IN S − S), b · ∇w K ε 1/2 δK (IN S − S)0,Ω ε 1/2 w N 1,Ω11 11
K⊂Ω11
CN −(p+1/2) w N SD , )1/2
CN −1
(4.16)
CN −1
since (εδK and δK from the definition of δK . For the second term in (4.15), one has I N Cδ 1/2 N −p wN SD CN −(p+1/2) wN SD . b · ∇(S − S δ ), b · ∇w K K K
(4.17)
K⊂Ω11
Finally, the third term in (4.15) is bounded by I N Cδ 1/2 N −(p+1) w N SD CN −(p+3/2) wN SD . c(S − S δ ), b · ∇w K K K
(4.18)
K⊂Ω11
Combining (4.13)–(4.18), the conclusion of the lemma follows.
2
Remark 14. The argument used in [12] to squeeze a higher order of convergence from K⊂Ω11 εδK ((IN S − S), b · ∇w N )K cannot be used here, since it leads to terms containing jumps in (IN S) across element boundaries. In the case p = 1 of [12] these terms were all zero.
Furthermore, when p > 1 we do not have any analogue of Lemma 4 for mixed-derivative terms such as K (IK w − w)x wyN dx dy. This is another difference from the case p = 1 of [12]. It is now straightforward to prove our main result. Theorem 15. Let Assumption 1 hold. Then the SDFEM solution uN satisfies uN − IN uSD CN −(p+1/2) .
(4.19)
Proof. By (2.7) and (2.4), we have 1 N u − IN u2SD BSD (uN − IN u, uN − IN u) = BSD (u − IN u, uN − IN u) 2 = BGal (u − IN u, uN − IN u) + BStab (u − IN u, uN − IN u). Now invoke Lemmas 12 and 13 to complete the proof.
2
5. Conclusions In this paper higher-order finite element methods (Qp with p > 1) on Shishkin meshes for convection–diffusion problems have been analysed for the first time. The main result, Theorem 15, uses a special approximant of Lin to yield an order of convergence 1/2 higher than can be obtained by standard approaches. Nevertheless this theorem is inferior to the case p = 1 considered in [12], where a gain of an entire order of convergence over the standard analysis is achieved. What is the reason for this inconsistency? It turns out that when p = 1, for each K ∈ TN an identity of Lin [6,14] expresses (b · ∇(IN u − u), w N )K in terms of values of u, w N and their derivatives on the boundary of K, with a higher-order remainder. Then the lower-order terms for neighbouring K can be combined to deduce a higher-order estimate. The discussion in [13, Remark 20] reveals that in the one-dimensional case when Ω is an interval, analogous identities exist also for p > 1, but now the lower-order terms for neighbouring K do not combine in a productive way. This is the fundamental reason for the unexpected result that the case p = 1 is exceptional. The gain of a half order of convergence in Theorem 15 can be exploited in practice. It shows that IN u is “superclose” to u, i.e. that uN − IN uSD exhibits a higher order of convergence than uN − uSD does in general. Consequently a straightforward generalization of certain techniques from [12,13] works here: an inexpensive local postprocessing of uN will yield a piecewise Qp+1 solution u˜ N such that u˜ N − uSD CN −(p+1/2) .
1802
M. Stynes, L. Tobiska / Applied Numerical Mathematics 58 (2008) 1789–1802
References [1] T. Apel, M. Dobrowolski, Anisotropic interpolation with applications to the finite element method, Computing 47 (1992) 277–293. [2] P. Ciarlet, The Finite Element Method for Elliptic Problems, SIAM, Philadelphia, PA, 2002. [3] P. Farrell, A. Hegarty, J. Miller, E. O’Riordan, G. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman & Hall/CRC, Boca Raton, FL, 2000. [4] C. Johnson, A. Schatz, L. Wahlbin, Crosswind smear and pointwise errors in the streamline diffusion finite element methods, Math. Comp. 49 (1987) 25–38. [5] Q. Lin, N. Yan, A. Zhou, A rectangle test for interpolated finite elements, in: Proc. Syst. Sci. Eng., Great Wall (H.K.) Culture Publish Co., 1991, pp. 217–229. [6] T. Linß, Uniform superconvergence of a Galerkin finite element method on Shishkin-type meshes, Numer. Methods Partial Differential Equations 16 (2000) 426–440. [7] T. Linß, M. Stynes, Asymptotic analysis and Shishkin-type decomposition for an elliptic convection–diffusion problem, J. Math. Anal. Appl. 261 (2001) 604–632. [8] J. Miller, E. O’Riordan, G. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [9] H. Roos, Layer-adapted grids for singular perturbation problems, Z. Angew. Math. Mech. 78 (1998) 291–309. [10] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, Berlin, 1996. [11] M. Stynes, L. Tobiska, A finite difference analysis of a streamline diffusion method on a Shishkin mesh, Numer. Algorithms 18 (1998) 337–360. [12] M. Stynes, L. Tobiska, The SDFEM for a convection–diffusion problem with a boundary layer: optimal error analysis and enhancement of accuracy, SIAM J. Numer. Anal. 41 (2003) 1620–1642. [13] L. Tobiska, Analysis of a new stabilized higher order finite element method for advection–diffusion equations, Comput. Methods Appl. Mech. Engrg. 196 (2006) 538–550. [14] Z. Zhang, Finite element superconvergence on Shishkin mesh for 2-d convection–diffusion problems, Math. Comp. 72 (2003) 1147–1177.