A BALANCED FINITE ELEMENT METHOD FOR SINGULARLY
c 2012 Society for Industrial and Applied Mathematics
SIAM J. NUMER. ANAL. Vol. 50, No. 5, pp. 2729–2743
A BALANCED FINITE ELEMENT METHOD FOR SINGULARLY PERTURBED REACTION-DIFFUSION PROBLEMS∗ RUNCHANG LIN† AND MARTIN STYNES‡ Abstract. Consider the singularly perturbed linear reaction-diffusion problem −ε2 Δu + bu = f in Ω ⊂ Rd , u = 0 on ∂Ω, where d ≥ 1, the domain Ω is bounded with (when d ≥ 2) Lipschitzcontinuous boundary ∂Ω, and the parameter ε satisfies 0 < ε 1. It is argued that for this type of problem, the standard energy norm v → [ε2 |v|21 + v20 ]1/2 is too weak a norm to measure adequately the errors in solutions computed by finite element methods: the multiplier ε2 gives an unbalanced norm whose different components have different orders of magnitude. A balanced and stronger norm is introduced, then for d ≥ 2 a mixed finite element method is constructed whose solution is quasioptimal in this new norm. For a problem posed on the unit square in R2 , an error bound that is uniform in ε is proved when the new method is implemented on a Shishkin mesh. Numerical results are presented to show the superiority of the new method over the standard mixed finite element method on the same mesh for this singularly perturbed problem. Key words. singularly perturbed, reaction-diffusion, mixed finite element method AMS subject classification. 65N30 DOI. 10.1137/110837784
1. Introduction. Consider the singularly perturbed linear reaction-diffusion problem (1.1a) (1.1b)
Lu := −ε2 Δu + bu = f
in Ω ⊂ Rd ,
u = 0 on ∂Ω,
where d ≥ 1, the domain Ω is bounded with (when d ≥ 2) Lipschitz-continuous boundary ∂Ω, and the parameter ε satisfies 0 < ε ≤ 1. Assume that f ∈ L2 (Ω) and ¯ with 0 < b0 ≤ b ≤ b1 for some constants b0 and b1 . Then (1.1) has a unique b ∈ C(Ω) solution u ∈ H01 (Ω). Here we have used the standard notation for the Sobolev spaces ¯ is the space of continuous functions defined on the L2 (Ω) and H01 (Ω), while C(Ω) ¯ of Ω. closure Ω Problems like (1.1) appear in certain heat transfer problems in thin domains (see [2, p. 179]), in Newton iterations of nonlinear reaction-diffusion problems, and in implicit time discretizations of parabolic reaction-diffusion problems when small step sizes are used. Furthermore, systems of equations of this type appear in several applications (see [7]), and in a future paper we shall extend our work to such systems. The singularly perturbed problem (1.1) has attracted much attention in the research literature. The properties of its solution u are discussed in [2, 5, 12, 14, 19], where it is shown that typically u exhibits sharp boundary layers near the boundary ∂Ω of Ω. Numerical methods for the approximation of u are analyzed in [2, 5, 9, 8, 11, 12, 13, 14, 18, 19, 20]. Some of these references use a finite difference approach, but here we consider only finite element methods for (1.1). ∗ Received by the editors June 17, 2011; accepted for publication (in revised form) September 10, 2012; published electronically October 23, 2012. http://www.siam.org/journals/sinum/50-5/83778.html † Department of Engineering, Mathematics, and Physics, Texas A&M International University (TAMIU), Laredo, TX 78041-1900 ([email protected]). This author was partially supported by a TAMIU University research grant. ‡ Department of Mathematics, National University of Ireland, Cork, Ireland ([email protected]).
2729
2730
RUNCHANG LIN AND MARTIN STYNES
At first sight the finite element approximation of u does not seem difficult. Define the energy norm 1/2 for each v ∈ H01 (Ω). v1,ε := ε2 |v|21 + v20 ˆ w) := ε2 (∇v, ∇w) + (bv, w) is combined with a conWhen the bilinear form B(v, forming finite element method, a standard analysis yields a quasi-optimal bound on the energy norm error u − uh 1,ε , where uh is the computed solution, and this error can be made small (independently of ε) by using, e.g., an appropriate Shishkin mesh [2, 9, 8, 12, 13], provided that u possesses sufficient regularity. Other finite element methods usually lead also to bounds on the energy norm error in the computed solution; see [10, 11, 14, 20]. But, unlike the classical case where ε = 1, for small ε the energy norm is a weak norm because (as we show in section 2) it is essentially no stronger than the L2 (Ω) norm when applied to the type of solution typically encountered in this singularly perturbed problem. A natural question then is, can one devise a finite element method whose accuracy can be bounded in a stronger norm that is better suited to (1.1)? In section 3 we show that this is indeed the case by presenting a new method that is bounded and coercive with respect to a much stronger “balanced” norm. An analysis yielding an error bound that is independent of the parameter ε is then developed in section 4 for a Shishkin mesh on a rectangular two-dimensional domain. The accuracy of the new method, measured in our balanced norm, is verified numerically in section 5. 1.1. Notation. For k = 0, 1, . . . we write H k (Ω) for the usual Sobolev space (so H 0 (Ω) = L2 (Ω)) with the associated seminorm and norm denoted by | · |k and · k . Let H01 (Ω) denote the space of functions in H 1 (Ω) whose traces vanish on ∂Ω. The L2 (Ω) inner product is denoted by (·, ·); thus (v, v) = v20 for all v ∈ L2 (Ω). This inner product notation is also used for vector-valued functions each of whose n components (in this paper, n = d or d + 1 as appropriate) n lies in L2 (Ω): that is, if v = (v1 , . . . , vn ) and w = (w1 , . . . , wn ), then (v, w) = i=1 (vi , wi ). Similarly we set n 2 1/2 |v|k = for k = 0, 1, . . . with vk defined likewise. i=1 |vi |k Set H(div; Ω) = {q ∈ [L2 (Ω)]d : ∇ · q ∈ L2 (Ω)} and H(Ω) = H(div; Ω) × H01 (Ω). In this paper, the vector-valued functions u, v, w, and f in H(Ω) have components ∇u q r 0 u= , v= , w= , and f= u v w f unless otherwise specified. Here ∇u, q, r, and 0 are vector functions that each have d components. One may think of q and r as corresponding to the gradients of v and w, respectively. Equivalently one can write, e.g., v = (q v)T , where T denotes transpose. We use C to denote a generic constant that is independent of ε and of any mesh; it can take different values at different places. 2. A priori estimates for u. In this section we motivate the new norm that will be introduced formally in section 3. Multiplying (1.1a) by u then integrating by parts over Ω yields the energy norm bound (2.1)
u1,ε ≤ C||f ||0 ,
where the constant C depends only on b0 . For small ε and smooth data, a typical solution of (1.1) has a well-behaved component whose low-order derivatives are essentially
2731
BALANCED FEM FOR REACTION-DIFFUSION
independent of ε but also contains a boundary layer term of the form e−kr(x,∂Ω)/ε , where r(x, ∂Ω) is the distance of x ∈ Rd from the boundary ∂Ω and k is some positive constant; see [5, 12, 14, 19]. The well-behaved component contributes O(ε) to ε|u|1 in u1,ε , while a short calculation shows that the layer term contributes O(ε1/2 ) to ε|u|1 . That is, u1,ε = u0 + O(ε1/2 ). Thus for small ε the energy norm of the solution u differs only slightly from the L2 norm of u; for singularly perturbed reaction-diffusion problems, the energy norm is a weak norm! If one knows only that f ∈ L2 (Ω), then the power of ε in the left-hand side of (2.1) is sharp. Consider the specific example d = 1, Ω = (0, 1), b ≡ 1, and f (x) = √ 2 sin(πx/ε). Then f 0 = 1 and an elementary calculation yields an explicit formula for u which shows that |u|1 = O(ε−1 ). Note that in this example |f |1 = O(ε−1 ). ¯ ∩ H 1 (Ω), that f L (Ω) and For the rest of this section, assume that f ∈ C(Ω) ∞ |f |1 are independent of ε, and that Ω is convex with smooth boundary ∂Ω. Then u ∈ H 2 (Ω). Let n denote the outward-pointing unit normal on ∂Ω. Multiply both sides of (1.1a) by −εΔu, then integrate over Ω. Integrating by parts and using (1.1b), we get f (∇u · n). ε3 Δu20 + ε(b∇u, ∇u) + ε(u∇b, ∇u) = ε(∇f, ∇u) − ε ∂Ω
Hence, for some C (independent of ε) one has ε
3
Δu20
+
b0 ε|u|21
b0 ε|u|21 + Cε(u20 + |f |21 ) − ε ≤ 2
f (∇u · n). ∂Ω
Invoking (2.1), we get (2.2)
ε3 Δu20 +
b0 ε|u|21 ≤ Cεf 21 − ε 2
f (∇u · n). ∂Ω
¯ define the function Let x0 ∈ ∂Ω be arbitrary but fixed. For all x ∈ Ω,
f L∞ (Ω) 1 φ(x) := b1 (x − x0 ) · n . 1 − exp b0 ε The convexity of Ω implies that φ ≥ 0 on ∂Ω. A short calculation and (1.1a) demonstrate that −ε2 Δφ+bφ ≥ f L∞(Ω) ≥ |−ε2 Δu+bu| on Ω. Thus φ is a barrier function ¯ But by construction for ±u on Ω and by a maximum principle we have |u| ≤ φ on Ω. √ φ(x0 ) = 0 = u(x0 ). Consequently |∇u(x0 ) · n| ≤ |∇φ(x0 ) · n| = b0bε1 f L∞ (Ω) . As √ x0 ∈ ∂Ω was arbitrary, it follows that ε ∂Ω f (∇u · n) ≤ b1 (meas ∂Ω)f 2L∞ (Ω) /b0 . Thus (2.2) implies that for some C one has (2.3)
ε3 Δu20 +
b0 ε|u|21 ≤ C εf 21 + f 2L∞(Ω) ≤ C f 21 + f 2L∞ (Ω) . 2
In (2.3) we have bounded ε1/2 |u|1 independently of ε. This is stronger than the energy norm inequality (2.1), which only bounded ε|u|1 . Additionally, (2.3) gives a bound on ε3/2 Δu0 . The powers of ε in the left-hand side of (2.3) are the best possible even if one assumes still more regularity of f , as can be seen by considering standard boundary layer functions such as e−x/ε when d = 1 and Ω = (0, 1).
2732
RUNCHANG LIN AND MARTIN STYNES
Finally, combine (2.1) and (2.3) to give the bound (2.4) ε3/2 Δu0 + ε1/2 |u|1 + u0 ≤ C f 1 + f L∞ (Ω) . For typical solutions u of (1.1), one finds that u0 = O(1), ε1/2 |u|1 = O(1), and ε3/2 Δu0 = O(1). Thus the weighted norm of u on the left-hand side of (2.4) is “balanced,” i.e., correctly scaled by powers of ε, in the sense that each of its components has the same order of magnitude, unlike the situation for u1,ε . The bound (2.4) was derived only to motivate the construction of our finite element method; it is not used subsequently in this paper. Remark 2.1. The unbalanced nature of the standard energy norm is implicitly exhibited in [21], where Zhang investigates the solution of one-dimensional reactiondiffusion and convection-diffusion problems by finite element methods and derives error bounds that are shown to be sharp by numerical experiments. In [21, Corollary 2.1] a discrete energy norm error bound for the reaction-diffusion problem takes the form (p is the polynomial degree of the finite elements)
C ε1/2 (N −1 ln N )p+1 + N −(p+1) , and here the factor ε1/2 can be traced to an excessive ε1/2 weighting of the H 1 component of the energy norm that is related to our discussion at the start of this section. On the other hand, the analogous error bound for the convection-diffusion problem in [21, Corollary 2.3] is C (N −1 ln N )p+1 + N −p ; the factor ε1/2 is now no longer present because in convection-diffusion problems the ε-weighted H 1 component and the L2 component of the solution are both O(1), i.e., unlike reaction-diffusion, the standard energy norm is balanced for convectiondiffusion problems. 3. A new bilinear form and finite element method. In this section we present a new bilinear form associated with (1.1) that is designed to yield coercivity and boundedness with respect to a norm inspired by (2.4). This leads directly to a new finite element method. The crucial multiplication of (1.1a) by −εΔu in the derivation of (2.3), if applied directly, would demand excessive smoothness of our finite element space, so we circumvent this difficulty by using a mixed finite element method. Thus (1.1a) is rewritten as the first-order system w = ∇u, (3.1) −ε2 ∇ · w + bu = f, which is then discretized. Our bilinear form is related to the symmetric bilinear form used in [11, section 3] but is asymmetric. For each v = (q v)T ∈ H(Ω), set √ √ ε (q − ∇v) ε (q − ∇v) ˜ 1 . (3.2) Av = and Av = (−ε∇ · q + bv) −ε2 ∇ · q + bv b Then (1.1a) is equivalent to the first-order system (3.3)
Au = f
in Ω,
2733
BALANCED FEM FOR REACTION-DIFFUSION
where f = (0 f )T as in section 1.1. We shall discretize this in the weak form ˜ = (f , Av) ˜ (Au, Av) for all v ∈ H(Ω). ˜ by 1/b is introduced to facilitate the proof The scaling of the second component of Av of coercivity in Theorem 3.1. In the special case where b = 1 and ε = 1 our method is of least-squares type (cf. [11] and the FOSLS and other approaches described in [3]), but it differs from least squares in the singularly perturbed case ε 1. Define the bilinear form B : H(Ω) × H(Ω) → R and the norm |||·|||ε by (3.4) ˜ := ε(q − ∇v, r − ∇w) + −ε2 ∇ · q + bv, 1 (−ε∇ · r + bw) B(v, w) := (Av, Aw) b and (3.5)
ε3 ε |||v|||ε = q20 + ∇v20 + b0 v20 ∇ · q20 + b1 2
1/2
for all v = (q v)T and w = (r w)T in H(Ω). If the solution u of (1.1) has typical boundary layers and is otherwise well-behaved, then on putting v = u and q = ∇u in (3.5), each term on the right-hand side is O(1) so |||·|||ε is balanced, unlike the energy norm · 1,ε . Furthermore |||·|||ε is stronger than the energy norm · 1,ε . Theorem 3.1 (boundedness and coercivity of B(·, ·) with respect to |||·|||ε ). There exists a positive constant C1 , which is independent of ε, such that (3.6a)
|B(v, w)| ≤ C1 |||v|||ε |||w|||ε for all v, w ∈ H(Ω)
and B(v, v) ≥ |||v|||2ε for all v ∈ H(Ω).
(3.6b)
Proof. To derive (3.6a), let v = (q v)T and w = (r w)T ∈ H(Ω) be arbitrary. Then ε3 |B(v, w)| ≤ ε (q0 + ∇v0 ) (r0 + ∇w0 ) + ∇ · q0 ∇ · r0 + b1 v0 w0 b0 (3.7) + ε2 |(∇ · q, w)| + ε|(v, ∇ · r)|. Integrating by parts, ε2 |(∇ · q, w)| + ε|(v, ∇ · r)| = ε2 |(q, ∇w)| + ε|(r, ∇v)| ≤ ε2 q0 ∇w0 + εr0 ∇v0 . This bound and (3.7) imply that |B(v, w)| ≤ C1 |||v|||ε |||w|||ε for some constant C1 which is independent of ε. Now we move on to the coercivity inequality. Let v = (q v)T ∈ H(Ω) be arbitrary. First, ∇·q 3 (3.8) B(v, v) = ε ∇ · q, + εq − ∇v20 − (ε + ε2 )(∇ · q, v) + (bv, v). b As v ∈ H01 (Ω), integration by parts gives
(3.9)
−(∇ · q, v) = (q · ∇v, 1) 1 = [(q + ∇v, q + ∇v) − (q − ∇v, q − ∇v)] 4 1 = q + ∇v20 − q − ∇v20 . 4
2734
RUNCHANG LIN AND MARTIN STYNES
Combining (3.8) and (3.9) yields ε + ε2 ε + ε2 ε3 q + ∇v20 + b0 v20 ∇ · q20 + ε − B(v, v) ≥ q − ∇v20 + b1 4 4 ε3 ε ≥ ∇ · q20 + q − ∇v20 + q + ∇v20 + b0 v20 , b1 4 where we used 0 < ε2 ≤ ε ≤ 1. But q − ∇v20 + q + ∇v20 = (q − ∇v, q − ∇v) + (q + ∇v, q + ∇v) = 2 q20 + ∇v20 . Thus B(v, v) ≥ |||v|||2ε , as desired. Assume that the hypotheses of Theorem 3.1 are satisfied. Let Vh ⊂ H(Ω) be any finite element space. Then our finite element approximation uh = (ph uh )T ∈ Vh of u is defined by (3.10)
˜ h ) for all vh ∈ Vh . B(uh , vh ) = (f , Av
Setting u = (∇u u)T , a standard argument yields the quasioptimality property (3.11)
|||u − uh |||ε ≤ C1 inf |||u − vh |||ε , vh ∈Vh
where C1 is the constant of Theorem 3.1. Remark 3.2. An analysis of the error |||u − uh |||ε will on a general mesh yield bounds that involve negative powers of ε. To obtain an error estimate that shows |||u − uh |||ε is small when ε < h 1, where h is the mesh diameter, one must use a mesh suited to the layers appearing in u, as in section 4 below. Furthermore, while a classical Aubin–Nitsche argument can be applied with the aim of obtaining a bound in the L2 norm that is of higher order than |||u − uh |||ε , here the small parameter ε interferes again with the analysis: one obtains only u − uh 0 ≤ Chε−1/2 |||u − uh |||ε . 4. Balanced norm error estimate on a Shishkin mesh. In this section we show how a suitably chosen mesh enables one to derive an error bound for the balanced norm error |||u − uh |||ε that is independent of ε. For this analysis one needs sharp a priori estimates on u and its derivatives together with a decomposition of u into its layer and well-behaved components. Fairly general discussions of this decomposition for domains Ω in Rd with smooth boundaries and with piecewise smooth boundaries are given in sections 2.2 and 3.2, respectively, of [19]. To keep our presentation as short as possible while taking a problem in more than one dimension, we shall consider only the special case d = 2 with Ω ≡ (0, 1)2 the unit square and use the lowest-order Raviart–Thomas space [4] to approximate functions in H(div; Ω). ¯ for some α ∈ (0, 1]. Then Assume that b and f lie in the H¨ older space C 4,α (Ω) 6,α u ∈ C (Ω). Assume also the corner compatibility conditions (4.1)
f (0, 0) = f (1, 0) = f (0, 1) = f (1, 1) = 0.
¯ see, e.g., [6]. Then u ∈ C 6,α (Ω) ∩ C 3,α (Ω); One could extend the analysis below to the case of higher-order Raviart–Thomas elements, but this would force us to assume extra corner compatibility conditions in order to bound the higher-order derivatives of u that would appear in the analysis.
BALANCED FEM FOR REACTION-DIFFUSION
2735
Notation. The edges of ∂Ω are Γ1 := {(x, 0)|0 ≤ x ≤ 1},
Γ2 := {(0, y)|0 ≤ y ≤ 1},
Γ3 := {(x, 1)|0 ≤ x ≤ 1},
Γ4 := {(1, y)|0 ≤ y ≤ 1}.
¯ as σ1 , σ2 , σ3 , σ4 , where σ1 is (0, 0) and the numbering is clockLabel the corners of Ω wise. For any measurable ω ⊂ Ω and integer m ≥ 0, we write · m,ω and | · |m,ω for the standard norm and seminorm in the Sobolev space H m (ω). A priori information about the solution u is derived in [5, 8, 13, 19]. This function can be decomposed into a sum of a smooth function, a boundary layer on each side of Ω, and a corner layer at each corner of Ω, as the following result from [13, Lemmas 1.1 and 1.2] describes. Set β = b0 /2. Lemma 4.1. The solution u of (1.1) can be decomposed as (4.2)
u=v+
4
wk +
k=1
4
zk ,
k=1
where each wk is a layer associated with the edge Γk and each zk is a layer associated ¯ one has with the corner σk . There exists a constant C such that for all (x, y) ∈ Ω m+n ∂ v(x, y) 2−m−n (4.3a) ) for 0 ≤ m + n ≤ 4, ∂xm ∂y n ≤ C(1 + ε m+n ∂ w1 (x, y) 2−m −n −βy/ε (4.3b) )ε e for 0 ≤ m + n ≤ 3, ∂xm ∂y n ≤ C(1 + ε m+n ∂ z1 (x, y) −m−n −β(x+y)/ε (4.3c) e for 0 ≤ m + n ≤ 3. ∂xm ∂y n ≤ Cε For w2 , w3 , and w4 , bounds analogous to (4.3b) can be derived. For z2 , z3 , and z4 , bounds analogous to (4.3c) can be derived. We now construct our Shishkin mesh. Let N be an even positive integer; it is the number of mesh intervals in each coordinate direction. The mesh is piecewise equidistant. The parameter λ that is used to specify where the mesh changes from coarse to fine is defined by 1 −1 , 2εβ ln N . (4.4) λ = min 4 To simplify the presentation, we make the reasonable practical assumption that ε ≤ CN −1 ,
(4.5)
as otherwise our mesh would resolve all the boundary layers and the subsequent analysis could be carried out using classical techniques. Then without loss of generality one can assume that N is so large that (4.6)
λ = 2εβ −1 ln N.
¯ = Ω11 ∪ Ω21 ∪ Ω12 ∪ Ω22 , where Partition Ω as follows (see Figure 4.1): Ω Ω11 = [λ, 1 − λ] × [λ, 1 − λ],
Ω21 = ([0, λ] ∪ [1 − λ, 1]) × [λ, 1 − λ],
Ω12 = [λ, 1 − λ] × ([0, λ] ∪ [1 − λ, 1]), Ω22 = ([0, λ] × ([0, λ] ∪ [1 − λ, 1])) ∪ ([1 − λ, 1] × ([0, λ] ∪ [1 − λ, 1])).
2736
RUNCHANG LIN AND MARTIN STYNES
Ω22 1−λ r
Ω12
Ω22
Ω21
Ω11
Ω21
λ r Ω22
Ω12
Ω r 22 1−λ
r λ
r
r r
r
Fig. 4.1. Shishkin mesh for reaction-diffusion.
Divide each of the x-intervals [0, λ] and [1 − λ, 1] into N/4 equidistant subintervals and divide [λ, 1 − λ] into N/2 equidistant subintervals. This gives a coarse mesh on [λ, 1 − λ] and a fine mesh elsewhere. Treat the y-interval [0,1] in a similar manner using λ. The final two-dimensional mesh is a tensor product of these one-dimensional piecewise equidistant Shishkin meshes; see Figure 4.1, where N = 8. This mesh divides Ω into a set Th of mesh rectangles K whose sides are parallel to the axes. The mesh is coarse on Ω11 , coarse/fine on Ω21 ∪ Ω12 , and√fine on Ω22 . The √ mesh is quasi-uniform on Ω11 and its diameter d there satisfies 2/N ≤ d ≤ 2 2/N ; on Ω12 ∪ Ω21 , the length and width of each mesh rectangle are 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 ). These quantities are used several times in our analysis. 1 ¯ to be piecewise ⊂ C(Ω) On this mesh we take the finite element space Vh,0 bilinears that vanish on the boundary ∂Ω. Let RT0,h denote the H(div; Ω)-conforming Raviart–Thomas space [4, 15, 16] of index 0 defined on Th . Set 1 ⊂ H(Ω). Vh = RT0,h × Vh,0
To carry out a satisfactory interpolation analysis on the highly anisotropic mesh Th , one uses the following sharp estimates [2]. Lemma 4.2. Let K ∈ Th be a rectangular element with sides of length hx , ky parallel to the x, y axes, respectively. Let φ ∈ H 2 (K). Then its piecewise bilinear nodal interpolant φI satisfies the bounds φ − φI 0,K ≤ C h2x φxx 0,K + ky2 φyy 0,K , (φ − φI )x 0,K ≤ C (hx φxx 0,K + ky φxy 0,K ) , (φ − φI )y 0,K ≤ C (hx φxy 0,K + ky φyy 0,K ) . By transforming to the unit square one easily obtains the standard inverse inequality ∂ψ ∂ψ (4.7) hx + ky ≤ ψ0,K for all ψ ∈ Vh1 and all K ∈ Th , ∂x 0,K ∂y 0,K where the mesh rectangle K has sides hx and ky .
2737
BALANCED FEM FOR REACTION-DIFFUSION
Let uI denote the piecewise bilinear nodal interpolant of u. Then one can prove the following L2 error estimate. Lemma 4.3 (see [13, Lemma 2.3]). There exists a constant C such that u − uI 0 ≤ CN −2 . In [13, Lemma 2.4] it is shown that ε∇(u − uI )0 ≤ CN −2 + Cε1/2 N −1 ln N . We now sharpen this bound by reducing the multiplier of ∇(u − uI )0 to ε1/2 . Lemma 4.4. There exists a constant C such that ε1/2 ∇(u − uI )0 ≤ CN −1 ln N. Proof. Lemma 4.1 gives (4.8)
ε
1/2
∇(u − u )0 = ε I
1/2
4 4 I I I (wk − wk ) + (zk − zk ) . ∇ v − v + k=1
k=1
0
Each term in this decomposition is bounded separately and (4.5) is used in several places without mentioning it. First, by Lemma 4.2 and (4.3a), 1/2 ∂ I (4.9) ε (v − v ) ≤ Cε1/2 N −1 |v|2,Ω ≤ CN −3/2 . ∂x 0 Now consider w1 . Lemma 4.2 and (4.3b) yield 2 2 ∂ w1 1/2 ∂ I 1/2 −1 ∂ w1 −1 (4.10) ε (w1 − w1 ) ≤ Cε N + ∂x∂y ≤ CN . 2 ∂x ∂x 0 0 0 Turning next to w2 , by Lemmas 4.1 and 4.2, (4.6), and (4.7) one obtains I I ∂w2 1/2 ∂(w2 − w2 ) 1/2 ∂w2 ε ≤ε + ≤ CN −3/2 , ∂x ∂x ∂x 0,Ω11 ∪Ω12 0,Ω11 ∪Ω12 0,Ω11 ∪Ω12 since |w2 | ≤ CN −2 on Ω11 ∪ Ω12 . On Ω21 ∪ Ω22 , by appealing to Lemmas 4.2 and 4.1 or inspecting [13, Lemma 2.4] one gets ε1/2 ∂(w2 − w2I )/∂x0,Ω21 ∪Ω22 ≤ CN −1 ln N . Consequently, ∂ I −1 (4.11) ε1/2 (w − w ) ln N. 2 ≤ CN ∂x 2 0 Similar results are valid for w3, w4 and the corner layer terms in (4.8). Hence, recalling (4.9)–(4.11), we have ε1/2 ∂(u − uI )/∂x0,Ω ≤ CN −1 ln N. In the same way
one can show that ε1/2 ∂(u − uI )/∂y0,Ω ≤ CN −1 ln N , and we are done. The standard local projection operator Π0,h into the Raviart–Thomas space RT0,h is defined on each mesh rectangle T = [xi−1 , xi ] × [yj−1 , yj ] as follows. On T one has RT0,h = (span{1, x}, span{1, y}). Given a function q = (q1 , q2 ) in H 1 (T ) × H 1 (T ), the first component (Π0,h q)1 of Π0,h q is defined by setting 1 (Π0,h q)1 (xk , y) = yj − yj−1
yj η=yj−1
q1 (xk , η) dη for k = i − 1, i and yj−1 ≤ y ≤ yj ,
2738
RUNCHANG LIN AND MARTIN STYNES
then interpolating linearly in x between these two values to define (Π0,h q)1 on T . The second component (Π0,h q)2 is defined similarly by interpolating linearly in y between the average values of q2 on the top and bottom of T . We shall need the following local result. Lemma 4.5. Let K ∈ Th be a rectangular element with sides of length hx , ky parallel to the x, y axes, respectively. Let Ph : L2 (K) → P0 (K) denote the local L2 projection onto a constant. Then for all q ∈ H 1 (K) and w = (w1 , w2 ) ∈ (H 1 (K))2 , there exists a constant C such that (4.12a) (4.12b)
q − Ph q0,K ≤ C(hx qx 0,K + ky qy 0,K ), (w − Π0,h w)i 0,K ≤ C[hx (wi )x 0,K + ky (wi )y 0,K ] for i = 1, 2.
Proof. See [10, Lemma 2.1] for (4.12a) and [1, Remark 4.1] for (4.12b). The operator Π0,h enjoys the “commuting diagram” property (see, e.g., [4, p. 130]) that, in the notation of Lemma 4.5, for each mesh rectangle K one has (4.13)
∇ · (Π0,h w) = Ph (∇ · w)
for all w ∈ (H 1 (K))2 ,
where Ph is defined in Lemma 4.5. Unfortunately, the standard Raviart–Thomas interpolation operator Π0,h is unsuited to singularly perturbed problems such as (1.1), for one cannot in general prove an analogue of Lemma 4.4 that bounds ε1/2 ∇u − Π0,h ∇u0 independently of ε, because the interpolation error of each boundary layer is too large on the coarse mesh rectangles that are adjacent to the fine mesh. To circumvent this difficulty we now ˜ 0,h ∇u of ∇u. define in RT0,h a modified interpolant Π 4 4 Recall from (4.2) the decomposition u = v + k=1 wk + k=1 zk . First, set ˜ 0,h ∇v = Π0,h ∇v. Π The boundary layer component w1 decays away from the side y = 0 of Ω. Let h∗ ˜ 0,h ∇w1 = Π0,h ∇w1 on each denote the fine mesh width in the Shishkin mesh. Set Π ˜ 0,h ∇w1 = (0, 0) on each mesh rectangle in mesh rectangle in [0, 1] × [0, λ − h∗ ] and Π [0, 1] × [λ, 1]. This leaves only a strip of fine mesh rectangles in [0, 1] × [λ − h∗ , λ]; on ˜ 0,h ∇w1 )1 = 0 and define (Π ˜ 0,h ∇w1 )2 to be the linear each of these rectangles set (Π ˜ interpolant in the y variable between the values of (Π0,h ∇w1 )2 on the top and bottom ˜ 0,h ∇wk of the mesh rectangle—each of these values is a constant. The interpolants Π are defined similarly for k = 2, 3, 4. ˜ 0,h ∇z1 of the corner layer z1 is defined analogously to Π ˜ 0,h ∇w1 : The interpolant Π ∗ ∗ ˜ 0,h ∇z1 = ˜ 0,h ∇z1 = Π0,h ∇z1 on each mesh rectangle in [0, λ−h ]×[0, λ−h ] and Π set Π (0, 0) on each mesh rectangle in Ω \ [0, λ) × [0, λ) , then imitate the construction of ˜ 0,h ∇w1 on each mesh rectangle in [0, λ − h∗ ] × [λ − h∗ , λ] and the construction of Π ˜ 0,h ∇w2 on each mesh rectangle in [λ − h∗ , λ] × [0, λ − h∗ ], which leaves only the Π ˜ 0,h ∇z1 is defined by linear single mesh rectangle [λ − h∗ , λ] × [λ − h∗ , λ] where Π ˜ 0,h ∇zk are defined interpolation of its values on the boundary. The interpolants Π likewise for k = 2, 3, 4. 4 ˜ ˜ ˜ 0,h ∇u = Π ˜ 0,h ∇v + 4 Π Finally, set Π k=1 0,h ∇wk + k=1 Π0,h ∇zk . Corollary 4.6. There exists a constant C such that ˜ 0,h ∇u ˜ 0,h ∇u (4.14) ε3/2 ∇ · ∇u − Π + ε1/2 ∇u − Π ≤ CN −1 ln N. 0 0 1/2 ˜ Proof. First consider ε ∇u − Π0,h ∇u 0 . For each component φ in the de˜ 0,h ∇φ = Π0,h ∇φ one can recomposition (4.2) of u, on mesh rectangles where Π peat the corresponding part of the proof of Lemma 4.4, invoking (4.12b) instead of
BALANCED FEM FOR REACTION-DIFFUSION
2739
Lemma 4.2. On the region ΩΠ (say) comprising all such rectangles, this leads to the ˜ 0,h ∇u0,ΩΠ ≤ CN −1 ln N . On the region Ωφ (say) where we set bound ε1/2 ∇u − Π ˜ 0,h ∇φ = (0, 0), one has ε1/2 ∇φ0,Ω ≤ CN −2 by the bounds of Lemma 4.1 and Π φ ˜ 0,h ∇φ0,Ω ≤ CN −2 . (4.6), so ε1/2 ∇φ − Π φ ˜ 0,h ∇φ on the region comprising the exThus it remains only to bound ∇φ − Π ˜ 0,h ∇φ was defined above by linear interpolation. ceptional mesh rectangles where Π We describe the analysis for the case φ = w1 ; the other components of u are handled similarly. The exceptional region is then R1 := [0, 1]× [λ− h∗ , λ]. By construction one ˜ 0,h ∇w1 )1 0,R1 = 0. As h∗ ≤ CεN −1 ln N , the bound (4.3b) and (4.6) yield has (Π 1/2 ˜ 0,h ∇w1 )2 (x, y)| ≤ Cε−1 N −2 for ε ∇w1 0,R1 ≤ CN −2 and the pointwise bound |(Π ∗ (x, y) ∈ R1 . But the area of R1 is h , so ˜ 0,h ∇w1 )2 0,R1 ≤ Cε−1/2 N −2 (h∗ )1/2 ≤ CN −5/2 (ln N )1/2 . ε1/2 (Π ˜ 0,h ∇u ≤ CN −1 ln N . The above arguments together show that ε1/2 ∇u − Π 0 To bound the other term in (4.14), we shall again use the decomposition (4.2) of u. First, by (4.13), Lemma 4.5, and (4.3a) we have ˜ 0,h ∇v ε3/2 ∇ · ∇v − Π = ε3/2 ∇ · (∇v − Π0,h ∇v)0 = ε3/2 Δv − Ph (Δv)0,Ω 0
≤ Cε3/2 N −1 |Δv|1 ≤ Cε1/2 N −1 . ˜ 0,h ∇w1 = Π0,h ∇w1 on the region R2 := [0, 1] × Next, consider w1 . Recall that Π ∗ [0, λ − h ]. Again appealing to (4.13), Lemma 4.5, and (4.3b), we get ˜ 0,h ∇w1 ε3/2 ∇ · ∇w1 − Π 0,R2
= ε3/2 Δw1 − Ph (Δw1 )0,R2 ∂(Δw ∂(Δw ) ) 1 1 ≤ Cε3/2 N −1 + (εN −1 ln N ) ∂x ∂y 0,R2 0,R2 ≤ CN −1 ln N. ˜ 0,h ∇w1 ≡ (0, 0), while (4.3b) and (4.6) On the region R3 := [0, 1] × [λ, 1] we have Π imply that ε3/2 Δw1 0,R3 ≤ CN −2 . Hence (4.15)
˜ 0,h ∇w1 )0,R3 ≤ CN −2 . ε3/2 ∇ · (∇w1 − Π
Now only the region R1 = [0, 1]×[λ−h∗, λ] remains to be considered. By construction ˜ 0,h ∇w1 )1 )/∂x0,R1 = 0. From (4.3b), (4.6), and h∗ ≤ CεN −1 ln N one obtains ∂((Π ˜ 0,h ∇w1 )2 (x, y)| ≤ Cε−1 N −2 for ε3/2 Δw1 0,R1 ≤ CN −2 and the pointwise bound |(Π (x, y) ∈ R1 , which implies that (4.16) ˜ 0,h ∇w1 )2 )/∂y0,R1 ≤ Cε3/2 (ε−1 N −2 /h∗ )(h∗ )1/2 ≤ CN −3/2 (ln N )−1/2 . ε3/2 ∂((Π ˜ 0,h ∇w1 )0,Ω ≤ CN −1 ln N . Combining all these bounds, we have ε3/2 ∇ · (∇w1 − Π The other layer components of u are treated similarly to w1 . 1 Define the interpolant uI of u from Vh = RT0,h × Vh,0 by uI := (Π0,h ∇u uI )T . Lemma 4.7. There exists a constant C such that |||u − uI |||ε ≤ CN −1 ln N.
2740
RUNCHANG LIN AND MARTIN STYNES
Proof. This follows from Lemma 4.3, Lemma 4.4, and Corollary 4.6. Our convergence result, which is independent of ε, is now immediate. Theorem 4.8. Consider (1.1) with Ω = (0, 1)2 . Assume that the corner com1 ¯ be the space of piecewise patibility conditions (4.1) are satisfied. Let Vh,0 ⊂ C(Ω) bilinears on the rectangular Shishkin mesh defined earlier in this section, which has N intervals in each coordinate direction. Then there exists a constant C such that the 1 , satisfies finite element solution uh defined in (3.10), with Vh = RT0,h × Vh,0 |||u − uh |||ε ≤ CN −1 ln N.
(4.17)
Proof. Combine (3.11) and Lemma 4.7. We remind the reader that each constant C is independent of ε and N . 5. Numerical results. In this section, we present numerical results for the finite element method defined in (3.10). The numerical example is taken from [13]. Other test problems gave similar results. Example 5.1. −ε2 Δu + 1 + x2 y 2 exy/2 u = f on (0, 1)2 , where f and the Dirichlet boundary conditions are chosen such that the exact solution is u(x, y) = x3 (1 + y 2 ) + sin(πx2 ) + cos(πy/2) + (x + y)[e−2x/ε + e−2(1−x)/ε + e−3y/ε + e−3(1−y)/ε ], whose properties are in agreement with the bounds of Lemma 4.1. Our earlier analysis can be extended to include inhomogeneous Dirichlet boundary conditions such as those satisfied by this solution. An approximate solution to Example 5.1 is computed using the finite element 1 on the rectangular Shishkin mesh method (3.10) with trial space Vh = RT0,h × Vh,0 of section 4. The mesh transition point λ is defined in (4.4) with β = 1/2. The stiffness matrix, the load vector, and the errors are calculated using a standard 16point Gauss–Legendre quadrature rule on each mesh element. In Tables 5.1–5.4, numerical results for a range of values of ε and N are presented. The errors appear in the lower half of the table and the rates of convergence in the Table 5.1 |||u − uh |||ε . N
ε=1
ε = 10−1
ε = 10−2
ε = 10−4
ε = 10−6
ε = 10−8
16 32 64 128 256 16 32 64 128 256 512
1.46 1.35 1.29 1.24 1.20 9.19e-1 4.62e-1 2.32e-1 1.16e-1 5.79e-2 2.90e-2
1.28 1.31 1.27 1.24 1.20 2.11e-0 1.15e-0 5.92e-1 2.98e-1 1.49e-1 7.47e-2
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.17e-0 2.30e-0 1.48e-0 8.84e-1 5.05e-1
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.16e-0 2.30e-0 1.48e-0 8.81e-1 5.03e-1
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.16e-0 2.29e-0 1.48e-0 8.81e-1 5.03e-1
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.16e-0 2.29e-0 1.48e-0 8.81e-1 5.03e-1
max
ε=10−m 2≤m≤8
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.17e-0 2.30e-0 1.48e-0 8.84e-1 5.05e-1
2741
BALANCED FEM FOR REACTION-DIFFUSION Table 5.2 ε1/2 ∇(u − uh )0 . N
ε=1
ε = 10−1
ε = 10−2
ε = 10−4
ε = 10−6
ε = 10−8
16 32 64 128 256 16 32 64 128 256 512
1.47 1.35 1.29 1.24 1.20 2.02e-1 1.01e-1 5.07e-2 2.53e-2 1.27e-2 6.33e-3
1.31 1.31 1.28 1.24 1.20 1.13e-0 6.11e-1 3.12e-1 1.57e-1 7.85e-2 3.93e-2
0.40 0.65 0.84 0.94 0.98 2.14e-0 1.77e-0 1.27e-0 8.07e-1 4.76e-1 2.71e-1
0.40 0.65 0.84 0.94 0.98 2.13e-0 1.77e-0 1.27e-0 8.05e-1 4.75e-1 2.70e-1
0.40 0.65 0.84 0.94 0.98 2.13e-0 1.77e-0 1.27e-0 8.05e-1 4.75e-1 2.70e-1
0.40 0.65 0.84 0.94 0.98 2.13e-0 1.77e-0 1.27e-0 8.05e-1 4.75e-1 2.70e-1
max
ε=10−m 2≤m≤8
0.40 0.65 0.84 0.94 0.98 2.14e-0 1.77e-0 1.27e-0 8.07e-1 4.76e-1 2.71e-1
Table 5.3 ε1/2 p − ph 0 . N
ε=1
ε = 10−1
ε = 10−2
ε = 10−4
ε = 10−6
ε = 10−8
16 32 64 128 256 16 32 64 128 256 512
1.57 1.38 1.29 1.24 1.21 8.01e-2 3.82e-2 1.89e-2 9.40e-3 4.70e-3 2.35e-3
2.54 2.45 2.11 1.67 1.35 4.07e-1 1.24e-1 3.54e-2 1.13e-2 4.46e-3 2.05e-3
0.85 1.28 1.65 1.86 1.95 1.36e-0 9.12e-1 4.74e-1 1.95e-1 6.88e-2 2.24e-2
0.84 1.28 1.65 1.86 1.95 1.34e-0 9.01e-1 4.69e-1 1.93e-1 6.82e-2 2.22e-2
0.84 1.28 1.65 1.86 1.95 1.34e-0 9.01e-1 4.69e-1 1.93e-1 6.82e-2 2.22e-2
0.84 1.28 1.65 1.86 1.95 1.34e-0 9.01e-1 4.69e-1 1.93e-1 6.82e-2 2.22e-2
max
ε=10−m 2≤m≤8
0.85 1.28 1.65 1.86 1.95 1.36e-0 9.12e-1 4.74e-1 1.95e-1 6.88e-2 2.24e-2
Table 5.4 u − uh 0 . N
ε=1
ε = 10−1
ε = 10−2
ε = 10−4
ε = 10−6
ε = 10−8
16 32 64 128 256 16 32 64 128 256 512
2.94 2.71 2.57 2.48 2.41 4.08e-3 1.02e-3 2.57e-4 6.42e-5 1.60e-5 4.01e-6
2.76 2.67 2.56 2.47 2.41 6.48e-2 1.77e-2 4.54e-3 1.14e-3 2.86e-4 7.15e-5
1.23 1.55 1.80 1.92 1.97 1.04e-1 5.82e-2 2.64e-2 1.00e-2 3.41e-3 1.09e-3
1.52 1.60 1.81 1.93 1.98 1.23e-2 6.04e-3 2.67e-3 1.01e-3 3.43e-4 1.10e-4
2.94 2.50 2.26 2.16 2.10 6.81e-3 1.71e-3 4.77e-4 1.41e-4 4.22e-5 1.26e-5
3.04 2.74 2.57 2.47 2.40 6.73e-3 1.61e-3 3.98e-4 9.95e-5 2.50e-5 6.27e-6
max
ε=10−m 2≤m≤8
1.23 1.55 1.80 1.92 1.97 1.04e-1 5.82e-2 2.64e-2 1.00e-2 3.41e-3 1.09e-3
upper half. Each convergence rate r is computed under the presumption that one has convergence of order (N −1 ln N )r , as such rates are typical for all norms on Shishkin meshes. The numerical results of Table 5.1 for the balanced norm |||·|||ε show convergence of order N −1 ln N , uniformly with respect to ε, which agrees with Theorem 4.8. Table 5.2 shows that ε1/2 ∇(u − uh )0 also exhibits convergence of order N −1 ln N . But
2742
RUNCHANG LIN AND MARTIN STYNES
in Table 5.3 we see that ε1/2 p − ph 0 attains convergence of order (N −1 ln N )2 when ε is small, which is superior to the first-order rate implied by Theorem 4.8. Finally, the L2 errors displayed in Table 5.4 show convergence of order (N −1 ln N )2 , uniformly in ε, but our current theory does not explain this. For comparison we also solved Example 5.1 on the same Shishkin meshes using other finite element methods that have been proposed for (1.1) and measured the errors in our balanced norm. This gave diverse results, as we now describe. First, consider the mixed finite element method of [10], where one rewrites (1.1) as the first-order system (3.1) then applies a standard bilinear form. The trial space is RT0,h × P0 , where P0 denotes the space of piecewise constants. Our numerical results for this method, which are not shown here, reveal that it fails to converge uniformly in ε with increasing N when errors are measured in the balanced norm ||| · ||| . On the other hand, when using the discontinuous Galerkin least-squares finite element method of [11], one obtains results similar to Table 5.1 : the error in |||u − uh |||ε is proportional to N −1 ln N and is uniform in ε. For the standard Galerkin method one observes behavior resembling Table 5.2: convergence of order N −1 ln N in the balanced norm defined by v → ε1/2 |v|1 + v0 . A recent report of Roos and Schopf [17] proves that this order of convergence is attained, but this norm is of course weaker than |||v|||ε , and furthermore, the standard Galerkin method cannot deliver the higher-order accuracy of the flux approximation that is displayed in Table 5.3. 6. Conclusion. In this paper we showed that the standard energy norm is unsuitable for singularly perturbed reaction-diffusion problems because it is essentially the same as the L2 norm. Consequently we introduced a new norm that is stronger than the energy norm and is moreover “balanced,” i.e., each component of this norm has the same order of magnitude for all values of the singular perturbation parameter ε when measuring the solution of the boundary value problem. The calculation that motivated this balanced norm also motivated a new mixed finite element method whose analysis fits naturally with the balanced norm. Theoretical convergence results (with respect to the balanced norm) were proved for this finite element method on a Shishkin mesh on a two-dimensional rectangular domain. Numerical results were presented that agree with these theoretical bounds in the balanced norm, though the rates of convergence of the L2 errors in the computed solution and its flux that were observed numerically are better than those predicted by our theory; this is a topic for future research. REFERENCES ´ n, The maximum angle condition for mixed and nonconforming [1] G. Acosta and R. G. Dura elements: Application to the Stokes equations, SIAM J. Numer. Anal., 37 (1999), pp. 18– 36. [2] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., Teubner, Stuttgart, 1999. [3] M. Berndt, T. A. Manteuffel, S. F. McCormick, and G. Starke, Analysis of first-order system least squares (FOSLS) for elliptic problems with discontinuous coefficients. I, SIAM J. Numer. Anal., 43 (2005), pp. 386–408. [4] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer Ser. Comput. Math. 15, Springer-Verlag, New York, 1991. [5] C. Clavero, J. L. Gracia, and E. O’Riordan, A parameter robust numerical method for a two dimensional reaction-diffusion problem, Math. Comp., 74 (2005), pp. 1743–1758. [6] H. Han and R. B. Kellogg, Differentiability properties of solutions of the equation −2 Δu + ru = f (x, y) in a square, SIAM J. Math. Anal., 21 (1990), pp. 394–408. ´ ska, Robust a-posteriori estimators for multilevel discretizations of [7] V. Klein and M. Peszyn reaction-diffusion systems, Int. J. Numer. Anal. Model., 8 (2011), pp. 1–27.
BALANCED FEM FOR REACTION-DIFFUSION
2743
[8] J. Li, Convergence and superconvergence analysis of finite element methods on highly nonuniform anisotropic meshes for singularly perturbed reaction-diffusion problems, Appl. Numer. Math., 36 (2001), pp. 129–154. [9] J. Li and I. M. Navon, Uniformly convergent finite element methods for singularly perturbed elliptic boundary value problems. I. Reaction-diffusion type, Comput. Math. Appl., 35 (1998), pp. 57–70. [10] J. Li and M. F. Wheeler, Uniform convergence and superconvergence of mixed finite element methods on anisotropically refined grids, SIAM J. Numer. Anal., 38 (2000), pp. 770–798. [11] R. Lin, Discontinuous discretization for least-squares formulation of singularly perturbed reaction-diffusion problems in one and two dimensions, SIAM J. Numer. Anal., 47 (2008/09), pp. 89–108. [12] T. Linß, Layer-Adapted Meshes for Reaction-Convection-Diffusion Problems, Lecture Notes in Math. 1985, Springer-Verlag, Berlin, 2010. [13] F. Liu, N. Madden, M. Stynes, and A. Zhou, A two-scale sparse grid method for a singularly perturbed reaction-diffusion problem in two dimensions, IMA J. Numer. Anal., 29 (2009), pp. 986–1007. [14] J. M. Melenk, hp-Finite Element Methods for Singular Perturbations, Lecture Notes in Math. 1796, Springer-Verlag, Berlin, 2002. [15] J.-C. N´ ed´ elec, Mixed finite elements in R3 , Numer. Math., 35 (1980), pp. 315–341. [16] P.-A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems, in Mathematical Aspects of Finite Element Methods, Lecture Notes in Math. 606, Springer, Berlin, 1977, pp. 292–315. [17] H.-G. Roos and M. Schopf, Convergence and Stability in Balanced Norms of Finite Element Methods on Shishkin Meshes for Reaction-Diffusion Problems, Technical report, Technical University of Dresden, 2012. [18] H.-G. Roos, M. Stynes, and L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, Convection-Diffusion-Reaction and Flow Problems, 2nd ed., Springer Ser. Comput. Math. 24, Springer-Verlag, Berlin, 2008. [19] G. I. Shishkin and L. P. Shishkina, Difference Methods for Singular Perturbation Problems. Chapman & Hall/CRC Monogr. Surv. Pure Appl. Math. 140, CRC Press, Boca Raton, FL, 2009. ¨ rth, Robust a posteriori error estimators for a singularly perturbed reaction-diffusion [20] R. Verfu equation, Numer. Math., 78 (1998), pp. 479–493. [21] Z. Zhang, Finite element superconvergence approximation for one-dimensional singularly perturbed problems, Numer. Methods Partial Differential Equations, 18 (2002), pp. 374–395.
SIAM J. NUMER. ANAL. Vol. 50, No. 5, pp. 2729–2743
A BALANCED FINITE ELEMENT METHOD FOR SINGULARLY PERTURBED REACTION-DIFFUSION PROBLEMS∗ RUNCHANG LIN† AND MARTIN STYNES‡ Abstract. Consider the singularly perturbed linear reaction-diffusion problem −ε2 Δu + bu = f in Ω ⊂ Rd , u = 0 on ∂Ω, where d ≥ 1, the domain Ω is bounded with (when d ≥ 2) Lipschitzcontinuous boundary ∂Ω, and the parameter ε satisfies 0 < ε 1. It is argued that for this type of problem, the standard energy norm v → [ε2 |v|21 + v20 ]1/2 is too weak a norm to measure adequately the errors in solutions computed by finite element methods: the multiplier ε2 gives an unbalanced norm whose different components have different orders of magnitude. A balanced and stronger norm is introduced, then for d ≥ 2 a mixed finite element method is constructed whose solution is quasioptimal in this new norm. For a problem posed on the unit square in R2 , an error bound that is uniform in ε is proved when the new method is implemented on a Shishkin mesh. Numerical results are presented to show the superiority of the new method over the standard mixed finite element method on the same mesh for this singularly perturbed problem. Key words. singularly perturbed, reaction-diffusion, mixed finite element method AMS subject classification. 65N30 DOI. 10.1137/110837784
1. Introduction. Consider the singularly perturbed linear reaction-diffusion problem (1.1a) (1.1b)
Lu := −ε2 Δu + bu = f
in Ω ⊂ Rd ,
u = 0 on ∂Ω,
where d ≥ 1, the domain Ω is bounded with (when d ≥ 2) Lipschitz-continuous boundary ∂Ω, and the parameter ε satisfies 0 < ε ≤ 1. Assume that f ∈ L2 (Ω) and ¯ with 0 < b0 ≤ b ≤ b1 for some constants b0 and b1 . Then (1.1) has a unique b ∈ C(Ω) solution u ∈ H01 (Ω). Here we have used the standard notation for the Sobolev spaces ¯ is the space of continuous functions defined on the L2 (Ω) and H01 (Ω), while C(Ω) ¯ of Ω. closure Ω Problems like (1.1) appear in certain heat transfer problems in thin domains (see [2, p. 179]), in Newton iterations of nonlinear reaction-diffusion problems, and in implicit time discretizations of parabolic reaction-diffusion problems when small step sizes are used. Furthermore, systems of equations of this type appear in several applications (see [7]), and in a future paper we shall extend our work to such systems. The singularly perturbed problem (1.1) has attracted much attention in the research literature. The properties of its solution u are discussed in [2, 5, 12, 14, 19], where it is shown that typically u exhibits sharp boundary layers near the boundary ∂Ω of Ω. Numerical methods for the approximation of u are analyzed in [2, 5, 9, 8, 11, 12, 13, 14, 18, 19, 20]. Some of these references use a finite difference approach, but here we consider only finite element methods for (1.1). ∗ Received by the editors June 17, 2011; accepted for publication (in revised form) September 10, 2012; published electronically October 23, 2012. http://www.siam.org/journals/sinum/50-5/83778.html † Department of Engineering, Mathematics, and Physics, Texas A&M International University (TAMIU), Laredo, TX 78041-1900 ([email protected]). This author was partially supported by a TAMIU University research grant. ‡ Department of Mathematics, National University of Ireland, Cork, Ireland ([email protected]).
2729
2730
RUNCHANG LIN AND MARTIN STYNES
At first sight the finite element approximation of u does not seem difficult. Define the energy norm 1/2 for each v ∈ H01 (Ω). v1,ε := ε2 |v|21 + v20 ˆ w) := ε2 (∇v, ∇w) + (bv, w) is combined with a conWhen the bilinear form B(v, forming finite element method, a standard analysis yields a quasi-optimal bound on the energy norm error u − uh 1,ε , where uh is the computed solution, and this error can be made small (independently of ε) by using, e.g., an appropriate Shishkin mesh [2, 9, 8, 12, 13], provided that u possesses sufficient regularity. Other finite element methods usually lead also to bounds on the energy norm error in the computed solution; see [10, 11, 14, 20]. But, unlike the classical case where ε = 1, for small ε the energy norm is a weak norm because (as we show in section 2) it is essentially no stronger than the L2 (Ω) norm when applied to the type of solution typically encountered in this singularly perturbed problem. A natural question then is, can one devise a finite element method whose accuracy can be bounded in a stronger norm that is better suited to (1.1)? In section 3 we show that this is indeed the case by presenting a new method that is bounded and coercive with respect to a much stronger “balanced” norm. An analysis yielding an error bound that is independent of the parameter ε is then developed in section 4 for a Shishkin mesh on a rectangular two-dimensional domain. The accuracy of the new method, measured in our balanced norm, is verified numerically in section 5. 1.1. Notation. For k = 0, 1, . . . we write H k (Ω) for the usual Sobolev space (so H 0 (Ω) = L2 (Ω)) with the associated seminorm and norm denoted by | · |k and · k . Let H01 (Ω) denote the space of functions in H 1 (Ω) whose traces vanish on ∂Ω. The L2 (Ω) inner product is denoted by (·, ·); thus (v, v) = v20 for all v ∈ L2 (Ω). This inner product notation is also used for vector-valued functions each of whose n components (in this paper, n = d or d + 1 as appropriate) n lies in L2 (Ω): that is, if v = (v1 , . . . , vn ) and w = (w1 , . . . , wn ), then (v, w) = i=1 (vi , wi ). Similarly we set n 2 1/2 |v|k = for k = 0, 1, . . . with vk defined likewise. i=1 |vi |k Set H(div; Ω) = {q ∈ [L2 (Ω)]d : ∇ · q ∈ L2 (Ω)} and H(Ω) = H(div; Ω) × H01 (Ω). In this paper, the vector-valued functions u, v, w, and f in H(Ω) have components ∇u q r 0 u= , v= , w= , and f= u v w f unless otherwise specified. Here ∇u, q, r, and 0 are vector functions that each have d components. One may think of q and r as corresponding to the gradients of v and w, respectively. Equivalently one can write, e.g., v = (q v)T , where T denotes transpose. We use C to denote a generic constant that is independent of ε and of any mesh; it can take different values at different places. 2. A priori estimates for u. In this section we motivate the new norm that will be introduced formally in section 3. Multiplying (1.1a) by u then integrating by parts over Ω yields the energy norm bound (2.1)
u1,ε ≤ C||f ||0 ,
where the constant C depends only on b0 . For small ε and smooth data, a typical solution of (1.1) has a well-behaved component whose low-order derivatives are essentially
2731
BALANCED FEM FOR REACTION-DIFFUSION
independent of ε but also contains a boundary layer term of the form e−kr(x,∂Ω)/ε , where r(x, ∂Ω) is the distance of x ∈ Rd from the boundary ∂Ω and k is some positive constant; see [5, 12, 14, 19]. The well-behaved component contributes O(ε) to ε|u|1 in u1,ε , while a short calculation shows that the layer term contributes O(ε1/2 ) to ε|u|1 . That is, u1,ε = u0 + O(ε1/2 ). Thus for small ε the energy norm of the solution u differs only slightly from the L2 norm of u; for singularly perturbed reaction-diffusion problems, the energy norm is a weak norm! If one knows only that f ∈ L2 (Ω), then the power of ε in the left-hand side of (2.1) is sharp. Consider the specific example d = 1, Ω = (0, 1), b ≡ 1, and f (x) = √ 2 sin(πx/ε). Then f 0 = 1 and an elementary calculation yields an explicit formula for u which shows that |u|1 = O(ε−1 ). Note that in this example |f |1 = O(ε−1 ). ¯ ∩ H 1 (Ω), that f L (Ω) and For the rest of this section, assume that f ∈ C(Ω) ∞ |f |1 are independent of ε, and that Ω is convex with smooth boundary ∂Ω. Then u ∈ H 2 (Ω). Let n denote the outward-pointing unit normal on ∂Ω. Multiply both sides of (1.1a) by −εΔu, then integrate over Ω. Integrating by parts and using (1.1b), we get f (∇u · n). ε3 Δu20 + ε(b∇u, ∇u) + ε(u∇b, ∇u) = ε(∇f, ∇u) − ε ∂Ω
Hence, for some C (independent of ε) one has ε
3
Δu20
+
b0 ε|u|21
b0 ε|u|21 + Cε(u20 + |f |21 ) − ε ≤ 2
f (∇u · n). ∂Ω
Invoking (2.1), we get (2.2)
ε3 Δu20 +
b0 ε|u|21 ≤ Cεf 21 − ε 2
f (∇u · n). ∂Ω
¯ define the function Let x0 ∈ ∂Ω be arbitrary but fixed. For all x ∈ Ω,
f L∞ (Ω) 1 φ(x) := b1 (x − x0 ) · n . 1 − exp b0 ε The convexity of Ω implies that φ ≥ 0 on ∂Ω. A short calculation and (1.1a) demonstrate that −ε2 Δφ+bφ ≥ f L∞(Ω) ≥ |−ε2 Δu+bu| on Ω. Thus φ is a barrier function ¯ But by construction for ±u on Ω and by a maximum principle we have |u| ≤ φ on Ω. √ φ(x0 ) = 0 = u(x0 ). Consequently |∇u(x0 ) · n| ≤ |∇φ(x0 ) · n| = b0bε1 f L∞ (Ω) . As √ x0 ∈ ∂Ω was arbitrary, it follows that ε ∂Ω f (∇u · n) ≤ b1 (meas ∂Ω)f 2L∞ (Ω) /b0 . Thus (2.2) implies that for some C one has (2.3)
ε3 Δu20 +
b0 ε|u|21 ≤ C εf 21 + f 2L∞(Ω) ≤ C f 21 + f 2L∞ (Ω) . 2
In (2.3) we have bounded ε1/2 |u|1 independently of ε. This is stronger than the energy norm inequality (2.1), which only bounded ε|u|1 . Additionally, (2.3) gives a bound on ε3/2 Δu0 . The powers of ε in the left-hand side of (2.3) are the best possible even if one assumes still more regularity of f , as can be seen by considering standard boundary layer functions such as e−x/ε when d = 1 and Ω = (0, 1).
2732
RUNCHANG LIN AND MARTIN STYNES
Finally, combine (2.1) and (2.3) to give the bound (2.4) ε3/2 Δu0 + ε1/2 |u|1 + u0 ≤ C f 1 + f L∞ (Ω) . For typical solutions u of (1.1), one finds that u0 = O(1), ε1/2 |u|1 = O(1), and ε3/2 Δu0 = O(1). Thus the weighted norm of u on the left-hand side of (2.4) is “balanced,” i.e., correctly scaled by powers of ε, in the sense that each of its components has the same order of magnitude, unlike the situation for u1,ε . The bound (2.4) was derived only to motivate the construction of our finite element method; it is not used subsequently in this paper. Remark 2.1. The unbalanced nature of the standard energy norm is implicitly exhibited in [21], where Zhang investigates the solution of one-dimensional reactiondiffusion and convection-diffusion problems by finite element methods and derives error bounds that are shown to be sharp by numerical experiments. In [21, Corollary 2.1] a discrete energy norm error bound for the reaction-diffusion problem takes the form (p is the polynomial degree of the finite elements)
C ε1/2 (N −1 ln N )p+1 + N −(p+1) , and here the factor ε1/2 can be traced to an excessive ε1/2 weighting of the H 1 component of the energy norm that is related to our discussion at the start of this section. On the other hand, the analogous error bound for the convection-diffusion problem in [21, Corollary 2.3] is C (N −1 ln N )p+1 + N −p ; the factor ε1/2 is now no longer present because in convection-diffusion problems the ε-weighted H 1 component and the L2 component of the solution are both O(1), i.e., unlike reaction-diffusion, the standard energy norm is balanced for convectiondiffusion problems. 3. A new bilinear form and finite element method. In this section we present a new bilinear form associated with (1.1) that is designed to yield coercivity and boundedness with respect to a norm inspired by (2.4). This leads directly to a new finite element method. The crucial multiplication of (1.1a) by −εΔu in the derivation of (2.3), if applied directly, would demand excessive smoothness of our finite element space, so we circumvent this difficulty by using a mixed finite element method. Thus (1.1a) is rewritten as the first-order system w = ∇u, (3.1) −ε2 ∇ · w + bu = f, which is then discretized. Our bilinear form is related to the symmetric bilinear form used in [11, section 3] but is asymmetric. For each v = (q v)T ∈ H(Ω), set √ √ ε (q − ∇v) ε (q − ∇v) ˜ 1 . (3.2) Av = and Av = (−ε∇ · q + bv) −ε2 ∇ · q + bv b Then (1.1a) is equivalent to the first-order system (3.3)
Au = f
in Ω,
2733
BALANCED FEM FOR REACTION-DIFFUSION
where f = (0 f )T as in section 1.1. We shall discretize this in the weak form ˜ = (f , Av) ˜ (Au, Av) for all v ∈ H(Ω). ˜ by 1/b is introduced to facilitate the proof The scaling of the second component of Av of coercivity in Theorem 3.1. In the special case where b = 1 and ε = 1 our method is of least-squares type (cf. [11] and the FOSLS and other approaches described in [3]), but it differs from least squares in the singularly perturbed case ε 1. Define the bilinear form B : H(Ω) × H(Ω) → R and the norm |||·|||ε by (3.4) ˜ := ε(q − ∇v, r − ∇w) + −ε2 ∇ · q + bv, 1 (−ε∇ · r + bw) B(v, w) := (Av, Aw) b and (3.5)
ε3 ε |||v|||ε = q20 + ∇v20 + b0 v20 ∇ · q20 + b1 2
1/2
for all v = (q v)T and w = (r w)T in H(Ω). If the solution u of (1.1) has typical boundary layers and is otherwise well-behaved, then on putting v = u and q = ∇u in (3.5), each term on the right-hand side is O(1) so |||·|||ε is balanced, unlike the energy norm · 1,ε . Furthermore |||·|||ε is stronger than the energy norm · 1,ε . Theorem 3.1 (boundedness and coercivity of B(·, ·) with respect to |||·|||ε ). There exists a positive constant C1 , which is independent of ε, such that (3.6a)
|B(v, w)| ≤ C1 |||v|||ε |||w|||ε for all v, w ∈ H(Ω)
and B(v, v) ≥ |||v|||2ε for all v ∈ H(Ω).
(3.6b)
Proof. To derive (3.6a), let v = (q v)T and w = (r w)T ∈ H(Ω) be arbitrary. Then ε3 |B(v, w)| ≤ ε (q0 + ∇v0 ) (r0 + ∇w0 ) + ∇ · q0 ∇ · r0 + b1 v0 w0 b0 (3.7) + ε2 |(∇ · q, w)| + ε|(v, ∇ · r)|. Integrating by parts, ε2 |(∇ · q, w)| + ε|(v, ∇ · r)| = ε2 |(q, ∇w)| + ε|(r, ∇v)| ≤ ε2 q0 ∇w0 + εr0 ∇v0 . This bound and (3.7) imply that |B(v, w)| ≤ C1 |||v|||ε |||w|||ε for some constant C1 which is independent of ε. Now we move on to the coercivity inequality. Let v = (q v)T ∈ H(Ω) be arbitrary. First, ∇·q 3 (3.8) B(v, v) = ε ∇ · q, + εq − ∇v20 − (ε + ε2 )(∇ · q, v) + (bv, v). b As v ∈ H01 (Ω), integration by parts gives
(3.9)
−(∇ · q, v) = (q · ∇v, 1) 1 = [(q + ∇v, q + ∇v) − (q − ∇v, q − ∇v)] 4 1 = q + ∇v20 − q − ∇v20 . 4
2734
RUNCHANG LIN AND MARTIN STYNES
Combining (3.8) and (3.9) yields ε + ε2 ε + ε2 ε3 q + ∇v20 + b0 v20 ∇ · q20 + ε − B(v, v) ≥ q − ∇v20 + b1 4 4 ε3 ε ≥ ∇ · q20 + q − ∇v20 + q + ∇v20 + b0 v20 , b1 4 where we used 0 < ε2 ≤ ε ≤ 1. But q − ∇v20 + q + ∇v20 = (q − ∇v, q − ∇v) + (q + ∇v, q + ∇v) = 2 q20 + ∇v20 . Thus B(v, v) ≥ |||v|||2ε , as desired. Assume that the hypotheses of Theorem 3.1 are satisfied. Let Vh ⊂ H(Ω) be any finite element space. Then our finite element approximation uh = (ph uh )T ∈ Vh of u is defined by (3.10)
˜ h ) for all vh ∈ Vh . B(uh , vh ) = (f , Av
Setting u = (∇u u)T , a standard argument yields the quasioptimality property (3.11)
|||u − uh |||ε ≤ C1 inf |||u − vh |||ε , vh ∈Vh
where C1 is the constant of Theorem 3.1. Remark 3.2. An analysis of the error |||u − uh |||ε will on a general mesh yield bounds that involve negative powers of ε. To obtain an error estimate that shows |||u − uh |||ε is small when ε < h 1, where h is the mesh diameter, one must use a mesh suited to the layers appearing in u, as in section 4 below. Furthermore, while a classical Aubin–Nitsche argument can be applied with the aim of obtaining a bound in the L2 norm that is of higher order than |||u − uh |||ε , here the small parameter ε interferes again with the analysis: one obtains only u − uh 0 ≤ Chε−1/2 |||u − uh |||ε . 4. Balanced norm error estimate on a Shishkin mesh. In this section we show how a suitably chosen mesh enables one to derive an error bound for the balanced norm error |||u − uh |||ε that is independent of ε. For this analysis one needs sharp a priori estimates on u and its derivatives together with a decomposition of u into its layer and well-behaved components. Fairly general discussions of this decomposition for domains Ω in Rd with smooth boundaries and with piecewise smooth boundaries are given in sections 2.2 and 3.2, respectively, of [19]. To keep our presentation as short as possible while taking a problem in more than one dimension, we shall consider only the special case d = 2 with Ω ≡ (0, 1)2 the unit square and use the lowest-order Raviart–Thomas space [4] to approximate functions in H(div; Ω). ¯ for some α ∈ (0, 1]. Then Assume that b and f lie in the H¨ older space C 4,α (Ω) 6,α u ∈ C (Ω). Assume also the corner compatibility conditions (4.1)
f (0, 0) = f (1, 0) = f (0, 1) = f (1, 1) = 0.
¯ see, e.g., [6]. Then u ∈ C 6,α (Ω) ∩ C 3,α (Ω); One could extend the analysis below to the case of higher-order Raviart–Thomas elements, but this would force us to assume extra corner compatibility conditions in order to bound the higher-order derivatives of u that would appear in the analysis.
BALANCED FEM FOR REACTION-DIFFUSION
2735
Notation. The edges of ∂Ω are Γ1 := {(x, 0)|0 ≤ x ≤ 1},
Γ2 := {(0, y)|0 ≤ y ≤ 1},
Γ3 := {(x, 1)|0 ≤ x ≤ 1},
Γ4 := {(1, y)|0 ≤ y ≤ 1}.
¯ as σ1 , σ2 , σ3 , σ4 , where σ1 is (0, 0) and the numbering is clockLabel the corners of Ω wise. For any measurable ω ⊂ Ω and integer m ≥ 0, we write · m,ω and | · |m,ω for the standard norm and seminorm in the Sobolev space H m (ω). A priori information about the solution u is derived in [5, 8, 13, 19]. This function can be decomposed into a sum of a smooth function, a boundary layer on each side of Ω, and a corner layer at each corner of Ω, as the following result from [13, Lemmas 1.1 and 1.2] describes. Set β = b0 /2. Lemma 4.1. The solution u of (1.1) can be decomposed as (4.2)
u=v+
4
wk +
k=1
4
zk ,
k=1
where each wk is a layer associated with the edge Γk and each zk is a layer associated ¯ one has with the corner σk . There exists a constant C such that for all (x, y) ∈ Ω m+n ∂ v(x, y) 2−m−n (4.3a) ) for 0 ≤ m + n ≤ 4, ∂xm ∂y n ≤ C(1 + ε m+n ∂ w1 (x, y) 2−m −n −βy/ε (4.3b) )ε e for 0 ≤ m + n ≤ 3, ∂xm ∂y n ≤ C(1 + ε m+n ∂ z1 (x, y) −m−n −β(x+y)/ε (4.3c) e for 0 ≤ m + n ≤ 3. ∂xm ∂y n ≤ Cε For w2 , w3 , and w4 , bounds analogous to (4.3b) can be derived. For z2 , z3 , and z4 , bounds analogous to (4.3c) can be derived. We now construct our Shishkin mesh. Let N be an even positive integer; it is the number of mesh intervals in each coordinate direction. The mesh is piecewise equidistant. The parameter λ that is used to specify where the mesh changes from coarse to fine is defined by 1 −1 , 2εβ ln N . (4.4) λ = min 4 To simplify the presentation, we make the reasonable practical assumption that ε ≤ CN −1 ,
(4.5)
as otherwise our mesh would resolve all the boundary layers and the subsequent analysis could be carried out using classical techniques. Then without loss of generality one can assume that N is so large that (4.6)
λ = 2εβ −1 ln N.
¯ = Ω11 ∪ Ω21 ∪ Ω12 ∪ Ω22 , where Partition Ω as follows (see Figure 4.1): Ω Ω11 = [λ, 1 − λ] × [λ, 1 − λ],
Ω21 = ([0, λ] ∪ [1 − λ, 1]) × [λ, 1 − λ],
Ω12 = [λ, 1 − λ] × ([0, λ] ∪ [1 − λ, 1]), Ω22 = ([0, λ] × ([0, λ] ∪ [1 − λ, 1])) ∪ ([1 − λ, 1] × ([0, λ] ∪ [1 − λ, 1])).
2736
RUNCHANG LIN AND MARTIN STYNES
Ω22 1−λ r
Ω12
Ω22
Ω21
Ω11
Ω21
λ r Ω22
Ω12
Ω r 22 1−λ
r λ
r
r r
r
Fig. 4.1. Shishkin mesh for reaction-diffusion.
Divide each of the x-intervals [0, λ] and [1 − λ, 1] into N/4 equidistant subintervals and divide [λ, 1 − λ] into N/2 equidistant subintervals. This gives a coarse mesh on [λ, 1 − λ] and a fine mesh elsewhere. Treat the y-interval [0,1] in a similar manner using λ. The final two-dimensional mesh is a tensor product of these one-dimensional piecewise equidistant Shishkin meshes; see Figure 4.1, where N = 8. This mesh divides Ω into a set Th of mesh rectangles K whose sides are parallel to the axes. The mesh is coarse on Ω11 , coarse/fine on Ω21 ∪ Ω12 , and√fine on Ω22 . The √ mesh is quasi-uniform on Ω11 and its diameter d there satisfies 2/N ≤ d ≤ 2 2/N ; on Ω12 ∪ Ω21 , the length and width of each mesh rectangle are 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 ). These quantities are used several times in our analysis. 1 ¯ to be piecewise ⊂ C(Ω) On this mesh we take the finite element space Vh,0 bilinears that vanish on the boundary ∂Ω. Let RT0,h denote the H(div; Ω)-conforming Raviart–Thomas space [4, 15, 16] of index 0 defined on Th . Set 1 ⊂ H(Ω). Vh = RT0,h × Vh,0
To carry out a satisfactory interpolation analysis on the highly anisotropic mesh Th , one uses the following sharp estimates [2]. Lemma 4.2. Let K ∈ Th be a rectangular element with sides of length hx , ky parallel to the x, y axes, respectively. Let φ ∈ H 2 (K). Then its piecewise bilinear nodal interpolant φI satisfies the bounds φ − φI 0,K ≤ C h2x φxx 0,K + ky2 φyy 0,K , (φ − φI )x 0,K ≤ C (hx φxx 0,K + ky φxy 0,K ) , (φ − φI )y 0,K ≤ C (hx φxy 0,K + ky φyy 0,K ) . By transforming to the unit square one easily obtains the standard inverse inequality ∂ψ ∂ψ (4.7) hx + ky ≤ ψ0,K for all ψ ∈ Vh1 and all K ∈ Th , ∂x 0,K ∂y 0,K where the mesh rectangle K has sides hx and ky .
2737
BALANCED FEM FOR REACTION-DIFFUSION
Let uI denote the piecewise bilinear nodal interpolant of u. Then one can prove the following L2 error estimate. Lemma 4.3 (see [13, Lemma 2.3]). There exists a constant C such that u − uI 0 ≤ CN −2 . In [13, Lemma 2.4] it is shown that ε∇(u − uI )0 ≤ CN −2 + Cε1/2 N −1 ln N . We now sharpen this bound by reducing the multiplier of ∇(u − uI )0 to ε1/2 . Lemma 4.4. There exists a constant C such that ε1/2 ∇(u − uI )0 ≤ CN −1 ln N. Proof. Lemma 4.1 gives (4.8)
ε
1/2
∇(u − u )0 = ε I
1/2
4 4 I I I (wk − wk ) + (zk − zk ) . ∇ v − v + k=1
k=1
0
Each term in this decomposition is bounded separately and (4.5) is used in several places without mentioning it. First, by Lemma 4.2 and (4.3a), 1/2 ∂ I (4.9) ε (v − v ) ≤ Cε1/2 N −1 |v|2,Ω ≤ CN −3/2 . ∂x 0 Now consider w1 . Lemma 4.2 and (4.3b) yield 2 2 ∂ w1 1/2 ∂ I 1/2 −1 ∂ w1 −1 (4.10) ε (w1 − w1 ) ≤ Cε N + ∂x∂y ≤ CN . 2 ∂x ∂x 0 0 0 Turning next to w2 , by Lemmas 4.1 and 4.2, (4.6), and (4.7) one obtains I I ∂w2 1/2 ∂(w2 − w2 ) 1/2 ∂w2 ε ≤ε + ≤ CN −3/2 , ∂x ∂x ∂x 0,Ω11 ∪Ω12 0,Ω11 ∪Ω12 0,Ω11 ∪Ω12 since |w2 | ≤ CN −2 on Ω11 ∪ Ω12 . On Ω21 ∪ Ω22 , by appealing to Lemmas 4.2 and 4.1 or inspecting [13, Lemma 2.4] one gets ε1/2 ∂(w2 − w2I )/∂x0,Ω21 ∪Ω22 ≤ CN −1 ln N . Consequently, ∂ I −1 (4.11) ε1/2 (w − w ) ln N. 2 ≤ CN ∂x 2 0 Similar results are valid for w3, w4 and the corner layer terms in (4.8). Hence, recalling (4.9)–(4.11), we have ε1/2 ∂(u − uI )/∂x0,Ω ≤ CN −1 ln N. In the same way
one can show that ε1/2 ∂(u − uI )/∂y0,Ω ≤ CN −1 ln N , and we are done. The standard local projection operator Π0,h into the Raviart–Thomas space RT0,h is defined on each mesh rectangle T = [xi−1 , xi ] × [yj−1 , yj ] as follows. On T one has RT0,h = (span{1, x}, span{1, y}). Given a function q = (q1 , q2 ) in H 1 (T ) × H 1 (T ), the first component (Π0,h q)1 of Π0,h q is defined by setting 1 (Π0,h q)1 (xk , y) = yj − yj−1
yj η=yj−1
q1 (xk , η) dη for k = i − 1, i and yj−1 ≤ y ≤ yj ,
2738
RUNCHANG LIN AND MARTIN STYNES
then interpolating linearly in x between these two values to define (Π0,h q)1 on T . The second component (Π0,h q)2 is defined similarly by interpolating linearly in y between the average values of q2 on the top and bottom of T . We shall need the following local result. Lemma 4.5. Let K ∈ Th be a rectangular element with sides of length hx , ky parallel to the x, y axes, respectively. Let Ph : L2 (K) → P0 (K) denote the local L2 projection onto a constant. Then for all q ∈ H 1 (K) and w = (w1 , w2 ) ∈ (H 1 (K))2 , there exists a constant C such that (4.12a) (4.12b)
q − Ph q0,K ≤ C(hx qx 0,K + ky qy 0,K ), (w − Π0,h w)i 0,K ≤ C[hx (wi )x 0,K + ky (wi )y 0,K ] for i = 1, 2.
Proof. See [10, Lemma 2.1] for (4.12a) and [1, Remark 4.1] for (4.12b). The operator Π0,h enjoys the “commuting diagram” property (see, e.g., [4, p. 130]) that, in the notation of Lemma 4.5, for each mesh rectangle K one has (4.13)
∇ · (Π0,h w) = Ph (∇ · w)
for all w ∈ (H 1 (K))2 ,
where Ph is defined in Lemma 4.5. Unfortunately, the standard Raviart–Thomas interpolation operator Π0,h is unsuited to singularly perturbed problems such as (1.1), for one cannot in general prove an analogue of Lemma 4.4 that bounds ε1/2 ∇u − Π0,h ∇u0 independently of ε, because the interpolation error of each boundary layer is too large on the coarse mesh rectangles that are adjacent to the fine mesh. To circumvent this difficulty we now ˜ 0,h ∇u of ∇u. define in RT0,h a modified interpolant Π 4 4 Recall from (4.2) the decomposition u = v + k=1 wk + k=1 zk . First, set ˜ 0,h ∇v = Π0,h ∇v. Π The boundary layer component w1 decays away from the side y = 0 of Ω. Let h∗ ˜ 0,h ∇w1 = Π0,h ∇w1 on each denote the fine mesh width in the Shishkin mesh. Set Π ˜ 0,h ∇w1 = (0, 0) on each mesh rectangle in mesh rectangle in [0, 1] × [0, λ − h∗ ] and Π [0, 1] × [λ, 1]. This leaves only a strip of fine mesh rectangles in [0, 1] × [λ − h∗ , λ]; on ˜ 0,h ∇w1 )1 = 0 and define (Π ˜ 0,h ∇w1 )2 to be the linear each of these rectangles set (Π ˜ interpolant in the y variable between the values of (Π0,h ∇w1 )2 on the top and bottom ˜ 0,h ∇wk of the mesh rectangle—each of these values is a constant. The interpolants Π are defined similarly for k = 2, 3, 4. ˜ 0,h ∇z1 of the corner layer z1 is defined analogously to Π ˜ 0,h ∇w1 : The interpolant Π ∗ ∗ ˜ 0,h ∇z1 = ˜ 0,h ∇z1 = Π0,h ∇z1 on each mesh rectangle in [0, λ−h ]×[0, λ−h ] and Π set Π (0, 0) on each mesh rectangle in Ω \ [0, λ) × [0, λ) , then imitate the construction of ˜ 0,h ∇w1 on each mesh rectangle in [0, λ − h∗ ] × [λ − h∗ , λ] and the construction of Π ˜ 0,h ∇w2 on each mesh rectangle in [λ − h∗ , λ] × [0, λ − h∗ ], which leaves only the Π ˜ 0,h ∇z1 is defined by linear single mesh rectangle [λ − h∗ , λ] × [λ − h∗ , λ] where Π ˜ 0,h ∇zk are defined interpolation of its values on the boundary. The interpolants Π likewise for k = 2, 3, 4. 4 ˜ ˜ ˜ 0,h ∇u = Π ˜ 0,h ∇v + 4 Π Finally, set Π k=1 0,h ∇wk + k=1 Π0,h ∇zk . Corollary 4.6. There exists a constant C such that ˜ 0,h ∇u ˜ 0,h ∇u (4.14) ε3/2 ∇ · ∇u − Π + ε1/2 ∇u − Π ≤ CN −1 ln N. 0 0 1/2 ˜ Proof. First consider ε ∇u − Π0,h ∇u 0 . For each component φ in the de˜ 0,h ∇φ = Π0,h ∇φ one can recomposition (4.2) of u, on mesh rectangles where Π peat the corresponding part of the proof of Lemma 4.4, invoking (4.12b) instead of
BALANCED FEM FOR REACTION-DIFFUSION
2739
Lemma 4.2. On the region ΩΠ (say) comprising all such rectangles, this leads to the ˜ 0,h ∇u0,ΩΠ ≤ CN −1 ln N . On the region Ωφ (say) where we set bound ε1/2 ∇u − Π ˜ 0,h ∇φ = (0, 0), one has ε1/2 ∇φ0,Ω ≤ CN −2 by the bounds of Lemma 4.1 and Π φ ˜ 0,h ∇φ0,Ω ≤ CN −2 . (4.6), so ε1/2 ∇φ − Π φ ˜ 0,h ∇φ on the region comprising the exThus it remains only to bound ∇φ − Π ˜ 0,h ∇φ was defined above by linear interpolation. ceptional mesh rectangles where Π We describe the analysis for the case φ = w1 ; the other components of u are handled similarly. The exceptional region is then R1 := [0, 1]× [λ− h∗ , λ]. By construction one ˜ 0,h ∇w1 )1 0,R1 = 0. As h∗ ≤ CεN −1 ln N , the bound (4.3b) and (4.6) yield has (Π 1/2 ˜ 0,h ∇w1 )2 (x, y)| ≤ Cε−1 N −2 for ε ∇w1 0,R1 ≤ CN −2 and the pointwise bound |(Π ∗ (x, y) ∈ R1 . But the area of R1 is h , so ˜ 0,h ∇w1 )2 0,R1 ≤ Cε−1/2 N −2 (h∗ )1/2 ≤ CN −5/2 (ln N )1/2 . ε1/2 (Π ˜ 0,h ∇u ≤ CN −1 ln N . The above arguments together show that ε1/2 ∇u − Π 0 To bound the other term in (4.14), we shall again use the decomposition (4.2) of u. First, by (4.13), Lemma 4.5, and (4.3a) we have ˜ 0,h ∇v ε3/2 ∇ · ∇v − Π = ε3/2 ∇ · (∇v − Π0,h ∇v)0 = ε3/2 Δv − Ph (Δv)0,Ω 0
≤ Cε3/2 N −1 |Δv|1 ≤ Cε1/2 N −1 . ˜ 0,h ∇w1 = Π0,h ∇w1 on the region R2 := [0, 1] × Next, consider w1 . Recall that Π ∗ [0, λ − h ]. Again appealing to (4.13), Lemma 4.5, and (4.3b), we get ˜ 0,h ∇w1 ε3/2 ∇ · ∇w1 − Π 0,R2
= ε3/2 Δw1 − Ph (Δw1 )0,R2 ∂(Δw ∂(Δw ) ) 1 1 ≤ Cε3/2 N −1 + (εN −1 ln N ) ∂x ∂y 0,R2 0,R2 ≤ CN −1 ln N. ˜ 0,h ∇w1 ≡ (0, 0), while (4.3b) and (4.6) On the region R3 := [0, 1] × [λ, 1] we have Π imply that ε3/2 Δw1 0,R3 ≤ CN −2 . Hence (4.15)
˜ 0,h ∇w1 )0,R3 ≤ CN −2 . ε3/2 ∇ · (∇w1 − Π
Now only the region R1 = [0, 1]×[λ−h∗, λ] remains to be considered. By construction ˜ 0,h ∇w1 )1 )/∂x0,R1 = 0. From (4.3b), (4.6), and h∗ ≤ CεN −1 ln N one obtains ∂((Π ˜ 0,h ∇w1 )2 (x, y)| ≤ Cε−1 N −2 for ε3/2 Δw1 0,R1 ≤ CN −2 and the pointwise bound |(Π (x, y) ∈ R1 , which implies that (4.16) ˜ 0,h ∇w1 )2 )/∂y0,R1 ≤ Cε3/2 (ε−1 N −2 /h∗ )(h∗ )1/2 ≤ CN −3/2 (ln N )−1/2 . ε3/2 ∂((Π ˜ 0,h ∇w1 )0,Ω ≤ CN −1 ln N . Combining all these bounds, we have ε3/2 ∇ · (∇w1 − Π The other layer components of u are treated similarly to w1 . 1 Define the interpolant uI of u from Vh = RT0,h × Vh,0 by uI := (Π0,h ∇u uI )T . Lemma 4.7. There exists a constant C such that |||u − uI |||ε ≤ CN −1 ln N.
2740
RUNCHANG LIN AND MARTIN STYNES
Proof. This follows from Lemma 4.3, Lemma 4.4, and Corollary 4.6. Our convergence result, which is independent of ε, is now immediate. Theorem 4.8. Consider (1.1) with Ω = (0, 1)2 . Assume that the corner com1 ¯ be the space of piecewise patibility conditions (4.1) are satisfied. Let Vh,0 ⊂ C(Ω) bilinears on the rectangular Shishkin mesh defined earlier in this section, which has N intervals in each coordinate direction. Then there exists a constant C such that the 1 , satisfies finite element solution uh defined in (3.10), with Vh = RT0,h × Vh,0 |||u − uh |||ε ≤ CN −1 ln N.
(4.17)
Proof. Combine (3.11) and Lemma 4.7. We remind the reader that each constant C is independent of ε and N . 5. Numerical results. In this section, we present numerical results for the finite element method defined in (3.10). The numerical example is taken from [13]. Other test problems gave similar results. Example 5.1. −ε2 Δu + 1 + x2 y 2 exy/2 u = f on (0, 1)2 , where f and the Dirichlet boundary conditions are chosen such that the exact solution is u(x, y) = x3 (1 + y 2 ) + sin(πx2 ) + cos(πy/2) + (x + y)[e−2x/ε + e−2(1−x)/ε + e−3y/ε + e−3(1−y)/ε ], whose properties are in agreement with the bounds of Lemma 4.1. Our earlier analysis can be extended to include inhomogeneous Dirichlet boundary conditions such as those satisfied by this solution. An approximate solution to Example 5.1 is computed using the finite element 1 on the rectangular Shishkin mesh method (3.10) with trial space Vh = RT0,h × Vh,0 of section 4. The mesh transition point λ is defined in (4.4) with β = 1/2. The stiffness matrix, the load vector, and the errors are calculated using a standard 16point Gauss–Legendre quadrature rule on each mesh element. In Tables 5.1–5.4, numerical results for a range of values of ε and N are presented. The errors appear in the lower half of the table and the rates of convergence in the Table 5.1 |||u − uh |||ε . N
ε=1
ε = 10−1
ε = 10−2
ε = 10−4
ε = 10−6
ε = 10−8
16 32 64 128 256 16 32 64 128 256 512
1.46 1.35 1.29 1.24 1.20 9.19e-1 4.62e-1 2.32e-1 1.16e-1 5.79e-2 2.90e-2
1.28 1.31 1.27 1.24 1.20 2.11e-0 1.15e-0 5.92e-1 2.98e-1 1.49e-1 7.47e-2
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.17e-0 2.30e-0 1.48e-0 8.84e-1 5.05e-1
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.16e-0 2.30e-0 1.48e-0 8.81e-1 5.03e-1
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.16e-0 2.29e-0 1.48e-0 8.81e-1 5.03e-1
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.16e-0 2.29e-0 1.48e-0 8.81e-1 5.03e-1
max
ε=10−m 2≤m≤8
0.40 0.62 0.81 0.93 0.97 3.80e-0 3.17e-0 2.30e-0 1.48e-0 8.84e-1 5.05e-1
2741
BALANCED FEM FOR REACTION-DIFFUSION Table 5.2 ε1/2 ∇(u − uh )0 . N
ε=1
ε = 10−1
ε = 10−2
ε = 10−4
ε = 10−6
ε = 10−8
16 32 64 128 256 16 32 64 128 256 512
1.47 1.35 1.29 1.24 1.20 2.02e-1 1.01e-1 5.07e-2 2.53e-2 1.27e-2 6.33e-3
1.31 1.31 1.28 1.24 1.20 1.13e-0 6.11e-1 3.12e-1 1.57e-1 7.85e-2 3.93e-2
0.40 0.65 0.84 0.94 0.98 2.14e-0 1.77e-0 1.27e-0 8.07e-1 4.76e-1 2.71e-1
0.40 0.65 0.84 0.94 0.98 2.13e-0 1.77e-0 1.27e-0 8.05e-1 4.75e-1 2.70e-1
0.40 0.65 0.84 0.94 0.98 2.13e-0 1.77e-0 1.27e-0 8.05e-1 4.75e-1 2.70e-1
0.40 0.65 0.84 0.94 0.98 2.13e-0 1.77e-0 1.27e-0 8.05e-1 4.75e-1 2.70e-1
max
ε=10−m 2≤m≤8
0.40 0.65 0.84 0.94 0.98 2.14e-0 1.77e-0 1.27e-0 8.07e-1 4.76e-1 2.71e-1
Table 5.3 ε1/2 p − ph 0 . N
ε=1
ε = 10−1
ε = 10−2
ε = 10−4
ε = 10−6
ε = 10−8
16 32 64 128 256 16 32 64 128 256 512
1.57 1.38 1.29 1.24 1.21 8.01e-2 3.82e-2 1.89e-2 9.40e-3 4.70e-3 2.35e-3
2.54 2.45 2.11 1.67 1.35 4.07e-1 1.24e-1 3.54e-2 1.13e-2 4.46e-3 2.05e-3
0.85 1.28 1.65 1.86 1.95 1.36e-0 9.12e-1 4.74e-1 1.95e-1 6.88e-2 2.24e-2
0.84 1.28 1.65 1.86 1.95 1.34e-0 9.01e-1 4.69e-1 1.93e-1 6.82e-2 2.22e-2
0.84 1.28 1.65 1.86 1.95 1.34e-0 9.01e-1 4.69e-1 1.93e-1 6.82e-2 2.22e-2
0.84 1.28 1.65 1.86 1.95 1.34e-0 9.01e-1 4.69e-1 1.93e-1 6.82e-2 2.22e-2
max
ε=10−m 2≤m≤8
0.85 1.28 1.65 1.86 1.95 1.36e-0 9.12e-1 4.74e-1 1.95e-1 6.88e-2 2.24e-2
Table 5.4 u − uh 0 . N
ε=1
ε = 10−1
ε = 10−2
ε = 10−4
ε = 10−6
ε = 10−8
16 32 64 128 256 16 32 64 128 256 512
2.94 2.71 2.57 2.48 2.41 4.08e-3 1.02e-3 2.57e-4 6.42e-5 1.60e-5 4.01e-6
2.76 2.67 2.56 2.47 2.41 6.48e-2 1.77e-2 4.54e-3 1.14e-3 2.86e-4 7.15e-5
1.23 1.55 1.80 1.92 1.97 1.04e-1 5.82e-2 2.64e-2 1.00e-2 3.41e-3 1.09e-3
1.52 1.60 1.81 1.93 1.98 1.23e-2 6.04e-3 2.67e-3 1.01e-3 3.43e-4 1.10e-4
2.94 2.50 2.26 2.16 2.10 6.81e-3 1.71e-3 4.77e-4 1.41e-4 4.22e-5 1.26e-5
3.04 2.74 2.57 2.47 2.40 6.73e-3 1.61e-3 3.98e-4 9.95e-5 2.50e-5 6.27e-6
max
ε=10−m 2≤m≤8
1.23 1.55 1.80 1.92 1.97 1.04e-1 5.82e-2 2.64e-2 1.00e-2 3.41e-3 1.09e-3
upper half. Each convergence rate r is computed under the presumption that one has convergence of order (N −1 ln N )r , as such rates are typical for all norms on Shishkin meshes. The numerical results of Table 5.1 for the balanced norm |||·|||ε show convergence of order N −1 ln N , uniformly with respect to ε, which agrees with Theorem 4.8. Table 5.2 shows that ε1/2 ∇(u − uh )0 also exhibits convergence of order N −1 ln N . But
2742
RUNCHANG LIN AND MARTIN STYNES
in Table 5.3 we see that ε1/2 p − ph 0 attains convergence of order (N −1 ln N )2 when ε is small, which is superior to the first-order rate implied by Theorem 4.8. Finally, the L2 errors displayed in Table 5.4 show convergence of order (N −1 ln N )2 , uniformly in ε, but our current theory does not explain this. For comparison we also solved Example 5.1 on the same Shishkin meshes using other finite element methods that have been proposed for (1.1) and measured the errors in our balanced norm. This gave diverse results, as we now describe. First, consider the mixed finite element method of [10], where one rewrites (1.1) as the first-order system (3.1) then applies a standard bilinear form. The trial space is RT0,h × P0 , where P0 denotes the space of piecewise constants. Our numerical results for this method, which are not shown here, reveal that it fails to converge uniformly in ε with increasing N when errors are measured in the balanced norm ||| · ||| . On the other hand, when using the discontinuous Galerkin least-squares finite element method of [11], one obtains results similar to Table 5.1 : the error in |||u − uh |||ε is proportional to N −1 ln N and is uniform in ε. For the standard Galerkin method one observes behavior resembling Table 5.2: convergence of order N −1 ln N in the balanced norm defined by v → ε1/2 |v|1 + v0 . A recent report of Roos and Schopf [17] proves that this order of convergence is attained, but this norm is of course weaker than |||v|||ε , and furthermore, the standard Galerkin method cannot deliver the higher-order accuracy of the flux approximation that is displayed in Table 5.3. 6. Conclusion. In this paper we showed that the standard energy norm is unsuitable for singularly perturbed reaction-diffusion problems because it is essentially the same as the L2 norm. Consequently we introduced a new norm that is stronger than the energy norm and is moreover “balanced,” i.e., each component of this norm has the same order of magnitude for all values of the singular perturbation parameter ε when measuring the solution of the boundary value problem. The calculation that motivated this balanced norm also motivated a new mixed finite element method whose analysis fits naturally with the balanced norm. Theoretical convergence results (with respect to the balanced norm) were proved for this finite element method on a Shishkin mesh on a two-dimensional rectangular domain. Numerical results were presented that agree with these theoretical bounds in the balanced norm, though the rates of convergence of the L2 errors in the computed solution and its flux that were observed numerically are better than those predicted by our theory; this is a topic for future research. REFERENCES ´ n, The maximum angle condition for mixed and nonconforming [1] G. Acosta and R. G. Dura elements: Application to the Stokes equations, SIAM J. Numer. Anal., 37 (1999), pp. 18– 36. [2] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, Adv. Numer. Math., Teubner, Stuttgart, 1999. [3] M. Berndt, T. A. Manteuffel, S. F. McCormick, and G. Starke, Analysis of first-order system least squares (FOSLS) for elliptic problems with discontinuous coefficients. I, SIAM J. Numer. Anal., 43 (2005), pp. 386–408. [4] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer Ser. Comput. Math. 15, Springer-Verlag, New York, 1991. [5] C. Clavero, J. L. Gracia, and E. O’Riordan, A parameter robust numerical method for a two dimensional reaction-diffusion problem, Math. Comp., 74 (2005), pp. 1743–1758. [6] H. Han and R. B. Kellogg, Differentiability properties of solutions of the equation −2 Δu + ru = f (x, y) in a square, SIAM J. Math. Anal., 21 (1990), pp. 394–408. ´ ska, Robust a-posteriori estimators for multilevel discretizations of [7] V. Klein and M. Peszyn reaction-diffusion systems, Int. J. Numer. Anal. Model., 8 (2011), pp. 1–27.
BALANCED FEM FOR REACTION-DIFFUSION
2743
[8] J. Li, Convergence and superconvergence analysis of finite element methods on highly nonuniform anisotropic meshes for singularly perturbed reaction-diffusion problems, Appl. Numer. Math., 36 (2001), pp. 129–154. [9] J. Li and I. M. Navon, Uniformly convergent finite element methods for singularly perturbed elliptic boundary value problems. I. Reaction-diffusion type, Comput. Math. Appl., 35 (1998), pp. 57–70. [10] J. Li and M. F. Wheeler, Uniform convergence and superconvergence of mixed finite element methods on anisotropically refined grids, SIAM J. Numer. Anal., 38 (2000), pp. 770–798. [11] R. Lin, Discontinuous discretization for least-squares formulation of singularly perturbed reaction-diffusion problems in one and two dimensions, SIAM J. Numer. Anal., 47 (2008/09), pp. 89–108. [12] T. Linß, Layer-Adapted Meshes for Reaction-Convection-Diffusion Problems, Lecture Notes in Math. 1985, Springer-Verlag, Berlin, 2010. [13] F. Liu, N. Madden, M. Stynes, and A. Zhou, A two-scale sparse grid method for a singularly perturbed reaction-diffusion problem in two dimensions, IMA J. Numer. Anal., 29 (2009), pp. 986–1007. [14] J. M. Melenk, hp-Finite Element Methods for Singular Perturbations, Lecture Notes in Math. 1796, Springer-Verlag, Berlin, 2002. [15] J.-C. N´ ed´ elec, Mixed finite elements in R3 , Numer. Math., 35 (1980), pp. 315–341. [16] P.-A. Raviart and J. M. Thomas, A mixed finite element method for 2nd order elliptic problems, in Mathematical Aspects of Finite Element Methods, Lecture Notes in Math. 606, Springer, Berlin, 1977, pp. 292–315. [17] H.-G. Roos and M. Schopf, Convergence and Stability in Balanced Norms of Finite Element Methods on Shishkin Meshes for Reaction-Diffusion Problems, Technical report, Technical University of Dresden, 2012. [18] H.-G. Roos, M. Stynes, and L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, Convection-Diffusion-Reaction and Flow Problems, 2nd ed., Springer Ser. Comput. Math. 24, Springer-Verlag, Berlin, 2008. [19] G. I. Shishkin and L. P. Shishkina, Difference Methods for Singular Perturbation Problems. Chapman & Hall/CRC Monogr. Surv. Pure Appl. Math. 140, CRC Press, Boca Raton, FL, 2009. ¨ rth, Robust a posteriori error estimators for a singularly perturbed reaction-diffusion [20] R. Verfu equation, Numer. Math., 78 (1998), pp. 479–493. [21] Z. Zhang, Finite element superconvergence approximation for one-dimensional singularly perturbed problems, Numer. Methods Partial Differential Equations, 18 (2002), pp. 374–395.
