- Email: [email protected]
A singularly perturbed problem with two parameters on a Bakhvalov-type mesh
Journal of Computational and Applied Mathematics 292 (2016) 307–319
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A singularly perturbed problem with two parameters on a Bakhvalov-type mesh Mirjana Brdar a,∗ , Helena Zarin b a
Faculty of Technology, University of Novi Sad, Bulevar cara Lazara 1, 21000 Novi Sad, Serbia
b
Department of Mathematics and Informatics, Faculty of Sciences, University of Novi Sad, Trg Dositeja Obradovića 4, 21000 Novi Sad, Serbia
article
abstract
info
Article history: Received 11 December 2014 Received in revised form 21 June 2015
A singularly perturbed problem with two small parameters is considered. On a Bakhvalovtype mesh we prove uniform convergence of a Galerkin finite element method with piecewise linear functions. Arguments in the error analysis include interpolation error bounds for a Clément quasi-interpolant as well as discretization error estimates in an energy norm. Numerical experiments support theoretical findings. © 2015 Elsevier B.V. All rights reserved.
MSC: 65L11 65L20 65L60 65L70 Keywords: Singularly perturbed problem Two small parameters Galerkin finite element method Bakhvalov-type mesh Clément quasi-interpolant
1. Introduction We consider the following two-parameter singularly perturbed boundary value problem Lu(x) := −ε1 u′′ + ε2 b(x)u′ + c (x)u = f (x), u(0) = 0,
x ∈ Ω = (0, 1),
u(1) = 0,
(1)
with b(x) ≥ b0 > 0, c ( x) −
1 2
c (x) ≥ c0 > 0,
ε2 b′ (x) ≥ γ > 0,
x ∈ Ω,
x ∈ Ω,
(2) (3)
where b, c , f are smooth functions on Ω = [0, 1], b0 , c0 , γ are constants and 0 < ε1 , ε2 ≪ 1 are small perturbation parameters. The two-parameter singularly perturbed problems appear in many areas of physics and mechanics. Among these are the mathematical models of liquid crystal materials and chemical reactions, electrical networks, control theory. Some examples can be found in [1,2].
∗
Corresponding author. E-mail addresses: [email protected] (M. Brdar), [email protected] (H. Zarin).
http://dx.doi.org/10.1016/j.cam.2015.07.011 0377-0427/© 2015 Elsevier B.V. All rights reserved.
308
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
In case ε2 = 0, the problem (1)–(3) belongs to the class of reaction–diffusion differential equations whose solution has 1/2 1/2 two boundary layers of the widths O (ε1 | ln ε1 |). For ε2 = 1, problem (1)–(3) reduces to the convection–diffusion case with a solution exhibiting a layer at x = 1 of the width O (ε1 | ln ε1 |). There is a vast literature on numerical solving of these two problems classes, see for instance [3,4] and references therein. Here we are interested in the purely two-parameter case 0 < ε1 , ε2 ≪ 1 from [5], where a relationship between the perturbation parameters defines the layer structure at x = 0 and x = 1. The main goal in the construction of numerical methods for singularly perturbed problems is acquiring their uniform convergence (robustness) with respect to the perturbation parameters. For the one-dimensional problem (1)–(3), the authors in [6–10] consider various robust finite difference schemes on piecewise uniform meshes of Shishkin type. In [11], an analysis of a simple upwind scheme on arbitrary meshes has been presented. Bakhvalov-type meshes appear in [12–14] p+1/2 where for a special ε2 = ε1 , p > 0, uniform convergence has been proved for higher-order finite difference methods. The streamline-diffusion finite element method on the Shishkin mesh for (1)–(3) has been presented in [5] where authors prove (almost) second order of convergence in the maximum norm. For this method, a robust a posteriori error estimate in the maximum norm can be found in [15]. In this paper we analyze a linear finite element method on a mesh of Bakhvalov-type for the given problem. The first optimal result of convergence in the energy norm but for a one-dimensional convection–diffusion problem has been presented in [16], where an error analysis includes estimates for a quasi-interpolant. Here we adopt the approach from [16] and for the two-parameter singularly perturbed problem prove the first order of uniform convergence in the energy norm for the Galerkin method on the Bakhvalov-type mesh. The paper is organized as follows. In the next section we present properties of the exact solution as well as its derivatives, with special emphasis on the layer structure. In Section 3 we give a weak formulation and the Galerkin approach to the resolution of our problem. The layer-adapted mesh has been presented in Section 4, while in Section 5 we prove interpolation error bounds for the Clément approximation operator. Section 6 contains the main result on the uniform error bound in the energy norm. The last section is devoted to numerical experiments that illustrate our theoretical results. Notation 1. Throughout the paper, C will denote a generic positive constant independent of the perturbation parameters ε1 , ε2 and the number of degrees of freedom N. For a set D ⊂ R, we use the standard notation for Banach spaces Lp (D), Sobolev spaces W k,p (D), H k (D) = W k,2 (D), norms ∥ · ∥Lp (D) and seminorms | · |H k (D) . 2. Properties of the continuous problem In order to describe the layers in the solution of (1)–(3), the characteristic equation is introduced
−ε1 r 2 (x) + ε2 b(x)r (x) + c (x) = 0. This equation defines two continuous functions r0 , r1 : [0, 1] → R with r0 (x) < 0 and r1 (x) > 0. Let
µ0 = − max r0 (x), x∈[0,1]
µ1 = min r1 (x), x∈[0,1]
that is
µ0,1 = min
∓ε2 b(x) +
ε22 b2 (x) + 4ε1 c (x)
.
2ε1 Different properties of µ0 and µ1 can be derived, while here we emphasize the following, [17], x∈[0,1]
µ0 ≤ µ1 ,
1/2
1 max{µ− 0 , ε1 µ1 } ≤ C (ε2 + ε1 ),
−1/2
(4)
1/2 2
ε2 µ0 ≤ b0 ∥c ∥L∞ (Ω ) , ε2 (ε1 µ1 ) ≤ Cε . (5) The values µ0 and µ1 determine the decay of the boundary layers. Since |r0 | < r1 , the layer at x = 1 is stronger than the layer at x = 0. The following assertion from [3] provides more details on the behavior of the solution of (1)–(3) and its −1
derivatives. Theorem 1 ([3]). Let b, c , f ∈ C q (Ω ) for some q ≥ 1, and let p, κ ∈ (0, 1) be arbitrary. Assume that q∥b′ ∥L∞ (Ω ) ε2 ≤ κ(1 − p). Then
|u(k) (x)| ≤ C (1 + µk0 e−pµ0 x + µk1 e−pµ1 (1−x) ),
x ∈ Ω,
for 0 ≤ k ≤ q. Now, from [18] the solution u to the problem (1)–(3) has the representation u = S + E0 + E1 , with
|S (k) (x)| ≤ C , |E0(k) (x)| ≤ C µk0 e−pµ0 x , for x ∈ Ω and 0 ≤ k ≤ q.
|E1(k) (x)| ≤ C µk1 e−pµ1 (1−x) ,
(6)
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
309
Fig. 1. The exact solution (37) for ε1 = 10−3 , ε2 = 10−1 (left) and ε1 = 10−3 , ε2 = 10−6 (right).
Fig. 2. The exact solution (37) for ε1 = 10−12 , ε2 = 10−1 (left) and ε1 = 10−12 , ε2 = 10−12 (right).
Depending on the relation between the perturbation parameters, there are two extremal situations, cf. [3,5], −1/2
(i) when ε22 ≪ ε1 ≪ 1, then µ0 and µ1 behave like ε1 and the layers are similar to the reaction–diffusion case (see right graphs on Figs. 1 and 2), (ii) when ε1 ≪ ε22 ≪ 1, then µ0 = O (ε2−1 ) and µ1 = O (ε2 ε1−1 ); now µ1 is much larger than µ0 and the boundary layer in the vicinity of x = 1 is stronger than the boundary layer at x = 0 (see left graphs on Figs. 1 and 2). In other cases of the relations between the perturbation parameters, µ0 and µ1 vary between the previously given values. Notice that 1 ≪ µ0 ≤ µ1 . 3. Weak formulation and the Galerkin approach For the problem (1)–(3), the standard weak formulation is: find u ∈ H01 (Ω ) such that a(u, v) = (f , v),
∀v ∈ H01 (Ω ),
(7)
where (·, ·) denotes the standard scalar product in L (Ω ), and a(·, ·) is the bilinear form 2
a(u, v) := ε1 (u′ , v ′ ) + ε2 (bu′ , v) + (cu, v),
u, v ∈ H01 (Ω ).
Let N + 1, N ∈ N, be the number of mesh points 0 = x0 < x1 < · · · < xN −1 < xN = 1
(8)
that are arbitrarily located at the interval Ω . With V we denote the space of continuous, piecewise linear functions with the standard basis {φi : i = 0, 1, . . . , N } of hat functions, i.e. φi (xj ) = δij , i, j = 0, 1, . . . , N. Let V0N = V N ∩ H01 (Ω ). The Galerkin method is characterized by: find uN ∈ V0N such that N
a(uN , v N ) = (f , v N ),
∀v N ∈ V0N .
(9)
Under the assumption (3), the bilinear form a is coercive with respect to the energy norm
∥v∥2E := ε1 |v|2H 1 (Ω ) + ∥v∥2L2 (Ω ) . Hence, the problems (7) and (9) have unique solutions.
(10)
310
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
The Galerkin method (9) generates a finite difference scheme that can be written as 1
LN ui =
xi+1
h¯ i
f φi dx,
i = 1, 2, . . . , N − 1,
xi−1
with u0 = uN = 0, where ui ≈ u(xi ), hi = xi − xi−1 , i = 1, 2, . . . , N, are the local mesh step sizes and h¯ i = (hi + hi+1 )/2, i = 1, 2, . . . , N − 1. The discrete operator takes form LN ui = −ε1
D + ui − D − ui h¯ i
ui + 1 − ui ui − ui − 1 1 + ε2 βi+1 − βi−1 + (γi−1 ui−1 + γi ui + γi+1 ui+1 ) , h¯ i h¯ i h¯ i
where D+ and D− are the usual forward and backward difference operators and
βi−1 =
xi
bφi′−1 φi dx,
βi+1 =
xi
bφi′+1 φi dx,
xi
xi−1
γi−1 =
xi+1
c φi−1 φi dx,
γi =
xi+1
c φi2 dx,
γi+1 =
xi+1
c φi+1 φi dx.
xi
xi−1
xi−1
In the sequel we shall prove uniform convergence of the Galerkin method (9) on a specially chosen layer-adapted mesh. The error bound will be derived in the energy norm (10) starting from the triangle inequality applied to the error decomposition u − uN = η + χ , where uN = S N + E0N + E1N , η = ηS + η0 + η1 , χ = χS + χ0 + χ1 and
ηS = S − S I , π
χS = S I − S N , π
ηj = Ej − Ej + κ ,
χj = Ej −
N j
(11) EjN
−κ , N j
j = 0, 1.
(12)
Here S I ∈ V N is the standard interpolant of the regular solution component S. The functions Ejπ ∈ V N , j = 0, 1, are special
interpolants of the layer components Ej that will be introduced in Section 5, and κjN ∈ V N are certain correction terms that allow ηj (0) = ηj (1) = 0, j = 0, 1 (the necessity for the new interpolant in the analysis that includes layer components is already described in [16]). Hence, the key ingredient in our error analysis is the solution decomposition from Theorem 1 that also plays an important role in the construction of the mesh used for the discretization. 4. Bakhvalov-type layer-adapted mesh Let N ∈ N, N ≥ 8, be divisible by 4 and
σj =
τ pµj
ln µj <
1 4
,
j = 0, 1,
(13)
where τ ≥ 1 is a user-chosen parameter and p ∈ (0, 1) is the parameter from Theorem 1. Since 1 ≪ µ0 ≤ µ1 , throughout the paper we assume 1 −1 −θ µ− , 1 ≤ µ0 ≤ CN
for some 0 < θ ≤ 1,
(14)
which in practical implementations does not represent a restriction. The solution u of the problem (1)–(3) exhibits two boundary layers in the vicinity of x = 0 and x = 1 of the widths 1 −1 O (µ− 0 ln µ0 ) and O (µ1 ln µ1 ), respectively. Therefore, we define a mesh on Ω such that it is equidistant on Ω c with N /2 mesh subintervals, and gradually divided on Ω 0 and Ω 1 with N /4 mesh subintervals, where
Ω0 = (0, σ0 ),
Ωc = (σ0 , 1 − σ1 ),
Ω1 = (1 − σ1 , 1).
Due to the existence of two boundary layers, we use two mesh generating functions ϕ0 and ϕ1 , with
1 ϕ0 (t ) = − ln 1 − 4(1 − µ− 0 )t , 1 ϕ1 (t ) = − ln 1 − 4(1 − µ− 1 )(1 − t ) . The functions ϕ0 and ϕ1 are piecewise continuously differentiable and satisfy the following conditions
ϕ0 (0) = 0, ϕ0 (1/4) = ln µ0 , ϕ1 (3/4) = ln µ1 , ϕ1 (1) = 0.
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
311
The mesh points xi , i = 0, 1, . . . , N, are defined with
τ ϕ0 (ti ), p µ0 1 (1 − σ0 − σ1 ), xi = σ0 + 2 ti − 4 τ 1 − ϕ1 (ti ), p µ1
i = 0, 1, . . . , i= i=
N N
,
,
4
,
+ 1, . . . ,
4 4 3N 3N 4
N
4
3N 4
,
(15)
+ 1, . . . , N ,
where ti = i/N , i = 0, 1, . . . , N. The transition points, i.e. the points where the mesh switches from fine to coarse and vice versa, are chosen to be x N = σ0 , x 3N = 1 − σ1 , in order to have 4
4
exp(−pµ0 x N ) = µ0
−τ
and
4
exp(−pµ1 (1 − x 3N )) = µ−τ 1 . 4
With τ ≥ 1 we control the smallness of layer functions in the vicinity of the transition points. Remark 1. Following [13], instead of (14) we can introduce a stronger assumption 1
µj
ln µj ≤ CN −1 ,
j = 0, 1,
emanating from the practice that in the case of an equidistant mesh there are no mesh points placed inside the layer region. Then for small values of perturbation parameters we have (14) with θ = 1. Remark 2. There exists a threshold value µmin = µmin (τ , p) such that (13) is fulfilled when µ1 ≥ µ0 ≥ µmin > 1. For example, when p = 0.5 and τ = 1, 2, 3, 4, 5, we get µmin ≈ 27, 68, 114, 164, 215, respectively. Since µj ≫ 1, j = 0, 1, the inequality (13) is not a limiting condition. Nevertheless, as our numerical experiments from the last section show, the condition (13) influences the admissible range for the perturbation parameters for which the mesh is layer-adapted. For example, for the first test problem when b(x) = c (x) = 1, we get ε1 ≤ 10−5 , ε2 ≤ 10−3 , while in the second test problem with b(x) = 3 − 2x2 , c (x) = 1, we obtain ε1 ≤ 10−6 , ε2 ≤ 10−3 . As already remarked, the mesh (15) is dense in the layer regions. Moreover, it can be easily shown that hi > hi−1 ,
for i = 2, 3, . . . , N /4,
(16)
hi−1 > hi ,
for i = 3N /4 + 2, . . . , N .
(17)
In the subsequent analysis we shall also need the following properties of the mesh step sizes. Lemma 1. The mesh sizes hi = xi − xi−1 , i = 1, 2, . . . , N, satisfy
−1 C µ0 , hi ≤ CN −1 , −1 C µ1 ,
i = 1, 2, . . . , N /4 − 1, i = N /4, . . . , 3N /4 + 1, i = 3N /4 + 2, . . . , N .
Proof. The functions ϕ0 and ϕ0′ are both strictly increasing on [0, 1/4]. Now, if i = 1, 2, . . . , N /4 − 1, by the mean value theorem there exists ξi ∈ (ti−1 , ti ) such that hi =
τ −1 ′ τ −1 ′ τ −1 4 1 N ϕ0 (ξi ) < N ϕ0 (t N −1 ) < N = C µ− 0 . 4 p µ0 p µ0 p µ0 1 − 4t N −1 4
For i = N /4 and ξi ∈ (t N −1 , t N ), 4
hN = 4
4
1 τ −1 4(1 − µ− τ −1 ′ 0 ) N ϕ0 (ξi ) ≤ N ≤ CN −1 . 1 1 pµ0 p µ0 1 − 4(1 − µ− 0 )4
Trivially, hi = 2N −1 (1 − σ0 − σ1 ) ≤ CN −1 ,
i=
N 4
+ 1, . . . ,
3N 4
.
The mesh step sizes on the interval Ω 1 can be analogously estimated.
Notice that hi ≤ CN −1 for all indices i. Now using xi−1 = xi − hi and the previous lemma, we easily get the following
312
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
Corollary 1. For the mesh points (15) we have e−µ0 xi−1 ≤ C e−µ0 xi ,
i = 1, 2, . . . ,
e−µ1 xi−1 ≤ C e−µ1 xi ,
i=
3N 4
N 4
− 1,
+ 2, . . . , N .
It is important to notice that from (16) it follows h N > h N −1 4
4
1 τ 4N −1 (1 − µ− 1 −1/2 0 ) = ln 1 + ≥ C µ− . 0 N 1 p µ0 1 − 4(1 − µ− 0 )t N −1 4
1 −1/2 Similarly, hi ≥ C µ− for i = 3N /4 + 1, 3N /4 + 2. Also, hi ≥ CN −1 for i = N /4 + 1, . . . , 3N /4. 1 N
Lemma 2. The mesh step sizes hi = xi − xi−1 , i = 1, 2, . . . , N, satisfy hi +1 hi hi −1 hi
≤ C,
i = 1, 2, . . . ,
≤ C,
i=
3N 4
N 4
− 2,
+ 3, . . . , N .
Proof. For each i = 1, 2, . . . , N /4 − 2 there exists ξi ∈ (ti−1 , ti ) such that hi =
τ −1 ′ N ϕ0 (ξi ). p µ0
Now hi +1 hi
=
1 1 − 4(1 − µ− 0 )ξi 1 1 − 4(1 − µ− 0 )ξi+1
8N
≤ 1+
−1
≤1+
(1 − µ0 )
1 8N −1 (1 − µ− 0 ) 1 1 − 4(1 − µ− 0 )ξi+1
−1
1 µ0 + 4N −1 (1 − µ− 0 )
−1
≤ C.
The other inequality can be derived analogously.
5. Interpolation error estimates for Clément approximation operator In the error analysis that includes the terms from (11)–(12), beside using the usual nodal interpolation operator from V N , we introduce a quasi-interpolant of Clément-type from [19] enriched with certain stability properties (more on this operator can be found in [20]). We first cite [16] and sketch the construction of the quasi-interpolant. For every mesh point xi , i = 0, 1, . . . , N, of a general regular mesh (8), by ∆i we denote a macroelement (xi−1 , xi+1 ), i = 1, . . . , N − 1, with ∆0 = (x0 , x1 ), ∆N = (xN −1 , xN ). Let P1 (∆i ) be the space of linear polynomials on ∆i . To a given function v ∈ L1 (Ω ) we define the local L2 -projection pi ∈ P1 (∆i ) by
∆i
pi s dx =
∆i
v s dx,
∀s ∈ P1 (∆i ).
Now, the quasi-interpolant v π ∈ V N of the function v ∈ L1 (Ω ) is
v π (x) :=
N
pi (xi )φi (x),
i=0
where φi represent hat functions. Notice that in general v π (xi ) ̸= v(xi ), i = 0, 1, . . . , N. For such a constructed operator, for a function v ∈ W k,p (∆) we have the following estimates
∥v − v π ∥Lp (e) ≤ CH k ∥Dk v∥Lp (∆) ,
p ∈ [1, ∞], k = 1, 2,
(18)
∥v π ∥Lp (e) ≤ C ∥v∥Lp (∆) , p ∈ [1, ∞]. (19) Here e = (xi−1 , xi ) denotes an element and ∆ is the union of the two corresponding macroelements ∆i−1 and ∆i , with the diameter H. We shall also use the inequality
|v − v π |H 1 (e) ≤ CH k−1 |v|H k (∆) , k = 1, 2, (20) which is valid only for elements e = (xi−1 , xi ), i = 1, 2, . . . , N /4 − 2 and i = 3N /4 + 3, . . . , N, according to [20] and Lemma 2.
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
313
Lemma 3. On the Bakhvalov-type mesh (15), for τ ≥ α/2, α > 0, we have 2
min{hi+1 µ0 , 1}e− α pµ0 xi ≤ CN −1 ,
i = 1, 2, . . . ,
2
min{hi µ1 , 1}e− α pµ1 (1−xi ) ≤ CN −1 ,
i=
3N
N 4
− 1,
+ 1, . . . , N − 1.
4
Proof. The idea for the proof is similar to [21,16]. First notice that for i = 1, 2, . . . , N /4 − 1 we have
min{hi+1 µ0 , 1}e
− α2 pµ0 xi
= min
1 4τ N −1 (1 − µ− 0 )
p(1 − 4(1 − µ0 )ξi+1 ) −1
,1
1 1 − 4(1 − µ− 0 )ti
2ατ
for some ξi+1 ∈ (ti , ti+1 ). From the assumption τ ≥ α/2, we get
min{hi+1 µ0 , 1}e
− α2 pµ0 xi
≤ min
4τ N −1 (1 − µ0−1 ) p(1 − 4(1 − µ0 )ti+1 ) −1
= min
1 4 τ N − 1 ( 1 − µ− 0 ) 1 p(1 − 4(1 − µ− 0 )ti+1 )
,1
1 1 − 4(1 − µ− 0 )ti
,1
1 1 −1 1 − 4(1 − µ− (1 − µ− 0 )ti+1 + 4N 0 ) .
The last expression can be written as
τ ω1 min , 1 (ω1 + ω2 ), p ω2
(21)
where 1 −1 ω1 = 4N −1 (1 − µ− , 0 ) ≤ CN
1 ω2 = 1 − 4(1 − µ− 0 )ti+1 .
If in (21) the minimum is 1, then
min
τ ω1 τ ω1 , 1 (ω1 + ω2 ) = ω1 + ω2 ≤ ω1 + ≤ CN −1 , p ω2 p
otherwise
τ ω1 ω1 τ ω1 , 1 (ω1 + ω2 ) = + 1 ≤ CN −1 . min p ω2 p ω2
The analysis for other values of the index i can be derived similarly.
Since the quasi-interpolant in general does not satisfy boundary conditions, we introduce correction functions κjN ∈ V N , j = 0, 1, defined as
κjN (xi ) =
(Ejπ − Ej )(xi ), 0,
i = 0, N , i = 1, 2, . . . , N − 1.
(22)
Now, recalling the earlier notation for ηj in (12), we can start with the interpolation error analysis. Lemma 4. Let the assumptions of Theorem 1 hold. If κjN are correction functions from (22) and Ejπ are Clément quasi-interpolants
of the solution components Ej , j = 0, 1, respectively, then the functions ηj = Ej − Ejπ + κjN on the Bakhvalov-type mesh (15) with
τ ≥ max{5/2, 2θ −1 } satisfy −1/2
∥ηj ∥L2 (Ω ) ≤ C µj
N −2 ,
j = 0, 1.
Proof. First we estimate the layer component η0 . Let e = (xi−1 , xi ), i = 2, 3, . . . , N /4 − 1, be some subinterval on the fine mesh and ∆ = (xi−2 , xi+1 ) ⊂ [0, σ0 ] an appropriate macroelement. The triangle inequality, stability estimate (19) and the properties of E0 in (6) together yield 1 −2pµ0 xi−2 ∥E0 − E0π ∥2L2 (e) ≤ C ∥E0 ∥2L2 (∆) ≤ C µ− . 0 e
(23)
On the other hand, if we use (18), (16) and estimate the integral by the maximum of the positive and monotonically decreasing integrand times the measure of the domain, we get
∥E0 − E0π ∥2L2 (e) ≤ CH 4 ∥E0′′ ∥2L2 (∆) ≤ Ch4i+1 µ40 ≤ Ch5i+1 µ40 e−2pµ0 xi−2 .
xi+1
e−2pµ0 x dx
xi−2
(24)
314
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
Combining both estimates (23) and (24) with Corollary 1, we have
5 2 1 1 −5 min{hi+1 µ0 , 1}e− 5 pµ0 xi ≤ C µ− ∥E0 − E0π ∥2L2 (e) ≤ C µ− . 0 0 N
(25)
The last inequality follows from Lemma 3 with α = 5. When e = (x0 , x1 ), analogous arguments yield (25) with i = 1. Also 1/2
1/2
∥κ0N ∥L2 (x0 ,x1 ) ≤ Ch1 |E0 (x0 ) − E0π (x0 )| ≤ Ch1 ∥E0 − E0π ∥L∞ (x0 ,x1 ) 5/2
1/2
≤ Ch1 (h1 + h2 )2 ∥E0′′ ∥L∞ (x0 ,x2 ) ≤ Ch1 µ20 −1/2
≤ C µ0
N −5/2 ,
(26)
1 −1 since h1 ≤ C µ− . Now 0 N N
N
−1 4
2 0 L2 (x i−1 ,xi )
∥η ∥
≤C
i =1
−1 4
∥E0 − E0π ∥2L2 (x
i−1 ,xi )
+ C ∥κ0N ∥2L2 (x
0 ,x 1 )
i=1 1 −4 ≤ C µ− . 0 N
(27)
For i ≥ N /4, we again use the triangle inequality, (19), (6) and (14) in order to get N
∥E0 − E0π ∥2L2 (x
i−1 ,xi )
≤C
i= N 4
N −1
∥E0 ∥2L2 (x
i−2 ,xi+1 )
+ C ∥E0 ∥2L2 (x
N −2 ,xN )
i= N 4
xN
≤C xN 4
1 e−2pµ0 x dx ≤ C µ− 0 e
−2pµ0 x N 4
−2
−2
2τ −1 2τ 1 1 1 ≤ C µ− µ0 + N − 1 = C µ− 1 − 4(1 − µ− 0 0 0 )t N −2 4
1 −2θτ ≤ C µ− . 0 N
Since κ0N is exponentially small on (xN −1 , xN ), the last inequality together with (27) and −θ τ ≤ −2 finally imply −1/2
∥η0 ∥L2 (Ω ) ≤ C µ0
N −2 .
Similar arguments lead to the interpolation error bound for the solution component E1 .
Lemma 5. Let the assumptions of Theorem 1 hold. If κjN are correction functions from (22) and Ejπ are Clément quasi-interpolants of the solution components Ej , j = 0, 1, respectively, then the functions ηj = Ej − Ejπ + κjN on the Bakhvalov-type mesh (15) with
τ ≥ 2θ −1 satisfy
1/2
|ηj |H 1 (Ω ) ≤ C µj N −1 ,
j = 0, 1.
Proof. Similarly as in the previous lemma, we start with the solution component E0 . First let e = (xi−1 , xi ), when i = 2, 3, . . . , N /4 − 2. Analogously to (23) we derive
|E0 − E0π |2H 1 (e) ≤ C µ0 e−2pµ0 xi−2 ,
(28)
where we have used (20) with k = 1. If we apply again (20) but with k = 2, then we have
|E0 − E0π |2H 1 (e) ≤ Ch3i+1 µ40 e−2pµ0 xi−2 and
3 2 |E0 − E0π |2H 1 (e) ≤ C µ0 min{hi+1 µ0 , 1}e− 3 pµ0 xi ≤ C µ0 N −3 , which follows from Lemma 3 with α = 3. The same estimate holds for the element e = (x0 , x1 ). Summation for i = 1, 2, . . . , N /4 − 2 yields N
−2 4 i =1
|E0 − E0π |2H 1 (x
i−1 ,xi )
≤ C µ0 N − 2 .
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
315
The correction function can be bounded with the inverse inequality and (26). Hence
|κ0N |H 1 (x0 ,x1 ) ≤
C h1
3/2
1/2
∥κ0N ∥L2 (x0 ,x1 ) ≤ Ch1 µ20 ≤ C µ0 N −3/2
and consequently N
−2 4
|η0 |2H 1 (x
i−1 ,xi )
≤ C µ0 N − 2 .
(29)
i=1
For i = N /4 − 1, . . . , 3N /4, the triangle inequality, the inverse inequality as well as the stability estimate (19) result in
|E0 − E0π |H 1 (e) ≤ |E0 |H 1 (e) +
1 hi
∥E0π ∥L2 (e) ≤ |E0 |H 1 (e) +
C hi
∥E0 ∥L2 (∆) .
Based on (6), Corollary 1, lower bounds for hi , (4) and (14), we have 3N 4
|E0 − E0π |2H 1 (x
i−1 ,xi )
1 2 −1 −1 2 τ ≤ C (µ0 + µ− ) ≤ C µ0 N 2(1−θ τ ) . 0 N )(µ0 + N
i= N −1 4
For i > 3N /4, it is sufficient to use the smallness of the component E0 and the correction function κ0N on Ω1 . Finally, we conclude 1/2
1/2
1/2
|η0 |H 1 (Ω ) ≤ C µ0 N −1 + C µ0 N 1−θτ ≤ C µ0 N −1 , since 1 − θ τ ≤ −1. The estimate for |η1 |H 1 (Ω ) can be similarly derived.
The following assertion represents the main result on the interpolation error. Theorem 2. Let the assumptions of Theorem 1 hold. On the Bakhvalov-type mesh (15) with τ ≥ max{5/2, 2θ −1 }, the functions ηS = S − S I and ηj = Ej − Ejπ + κjN can be estimated with 1/2
∥ηS ∥E ≤ C ε1 N −1 + CN −2 , ∥ηj ∥E ≤ C (ε2 + ε
1/2 1/2 −1 N 1
)
(30)
,
j = 0, 1,
(31)
where S I ∈ V N is the standard interpolant of S , Ejπ are quasi-interpolants of Ej , j = 0, 1, and κjN , j = 0, 1, are functions from (22). Proof. The bounds for ∥ηj ∥E are a direct consequence of the previous two lemmata, (4) and −1/2
∥ηj ∥E ≤ C (ε1 µ1 )1/2 N −1 + C µ0
N −2 ,
j = 0, 1.
The estimate for the energy norm of ηS = S − S I is already known (see for instance [3]), and follows from
∥ηS ∥L2 (Ω ) ≤ CN −2 and |ηS |H 1 (Ω ) ≤ CN −1 , which completes the proof.
6. Convergence in the energy norm In the sequel we shall first derive the estimate for ∥χ∥E , which together with Theorem 2 will provide us with the final error bound. Theorem 3. Let the assumptions of Theorems 1 and 2 hold. The discretization error for the problem (1)–(3) on the Bakhvalov-type mesh (15) can be estimated with 1/2
∥χ∥E ≤ C ε2 N −1 + CN −2 . Proof. The bilinear form a(·, ·) is coercive with respect to the energy norm. Let α = min{1, γ }. Coercivity and Galerkin orthogonality imply
α∥χ∥2E ≤ a(χ , χ ) = −a(η, χ ) = −ε2 (bη′ , χ ) − (c η, χ ).
(32)
316
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
The last term in (32) can be directly estimated applying the Cauchy–Schwarz inequality, Lemma 4 and the L2 -norm bound on ηS in the proof of Theorem 2 and (4). Therefore
|(c η, χ )| ≤ C ∥η∥L2 (Ω ) ∥χ∥L2 (Ω ) −1/2
≤ C (N −2 + µ0
N −2 )∥χ∥E ≤ CN −2 ∥χ∥E .
(33)
The remaining term in (32) we analyze according to the splitting η = (ηS + η0 ) + η1 . First we get
|ε2 (b(ηS′ + η0′ ), χ )| ≤ C ε2 (|ηS |H 1 (Ω ) + |η0 |H 1 (Ω ) )∥χ∥L2 (Ω ) 1/2
1/2
≤ C ε2 (N −1 + µ0 N −1 )∥χ∥E ≤ C ε2 N −1 ∥χ ∥E .
(34)
Here we have used the Cauchy–Schwarz inequality, Lemma 5 and the H 1 -seminorm bound on ηS in the proof of Theorem 2 and (5). These arguments are not sufficient for estimating the expression with η1 . Instead we apply integration by parts and since η1 satisfies homogeneous Dirichlet boundary conditions, we can derive
|ε2 (bη1′ , χ )| = ε2 |(η1 , b′ χ + bχ ′ )| ≤ C ε2 ∥η1 ∥L2 (Ω ) (∥χ∥L2 (Ω ) + |χ|H 1 (Ω ) ) −1/2
≤ C ε2 µ1 −1/2
Due to ε2 (ε1 µ1 )
≤ Cε
1/2 2
−1/2
N − 2 ( 1 + ε1
)∥χ ∥E .
from (5), we have
1/2
|ε2 (bη1′ , χ )| ≤ C ε2 N −2 ∥χ∥E . Collecting (33)–(35) into (32) and dividing the inequality with ∥χ∥E , we complete the proof.
(35)
The final error estimate in the last theorem represents a direct consequence of Theorems 2, 3 and (4). Theorem 4. Let uN be the Galerkin approximation of the solution u of (1)–(3) on the Bakhvalov-type mesh (15). If the assumptions of Theorem 2 are satisfied, then 1/2
∥u − uN ∥E ≤ C (ε2 + ε1 )1/2 N −1 + CN −2 .
(36)
Remark 3. Under the assumption ε1 ≤ CN −1 , the author in [16] proved ∥u − uN ∥E ≤ CN −1 for the convection–diffusion problem. We can consider (36) as the generalization of the aforementioned result, since for ε2 = 1 we place a uniform mesh 1 −1 on [0, 1 − σ1 ] and use µ− . 1 ≤ C ε1 ≤ CN Remark 4. In the reaction–diffusion case, i.e. ε2 = 0, the linear Galerkin finite element approximation uN on an appropriate Shishkin mesh satisfies, [22], 1/4
∥u − uN ∥E ≤ C ε1 N −1 ln N + CN −2 , where N + 1 is the total number of nodes. The error estimate (36) remains valid in this setting even though Theorem 4 cannot be applied. Also, we would like to draw attention to a result in [23] for the two-dimensional two-parameter singularly perturbed problem and the Galerkin FEM on the Shishkin mesh. The energy norm error estimate is similar to (36), i.e. 1/2
∥u − uN ∥E ≤ C (ε2 + ε1 )1/2 N −1 ln N + CN −2 , where again a logarithmic factor stems from the choice of transition points. Here, N + 1 is the number of nodes located in each direction. Remark 5. As numerical tests from the next section suggest, a generalization with higher order polynomials should be possible. Let V N ,r = {v ∈ C (Ω ) : v|(xi−1 ,xi ) ∈ Pr (xi−1 , xi ), i = 1, 2, . . . , N } be the finite-dimensional space of element-wise polynomials of the highest degree r ≥ 1. For a finite element approximation uN ∈ V N ,r ∩ H01 (Ω ) we conjecture an error bound 1/2
∥u − uN ∥E ≤ C (ε2 + ε1 )1/2 N −r + CN −(r +1) . The proof on the layer-adapted mesh (15) with τ ≥ max{r + 3/2, 2r θ −1 } requires an approximation operator that exhibits properties similar to (18)–(20). In [24] such an operator has been constructed as a mapping from Lp (Ω ) onto the space of spline functions of the order r + 1. If for an element e = (xi−1 , xi ) we define ∆ = (xi−r −1 , xi+r ), then estimates
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
317
of the type (18) and (19) can be used with k ≤ r + 1, cf. [24, Chapter 5, Theorem 4.4 and Theorem 4.5]. Moreover, if ∆ ⊂ [0, xN /4−1 ] ∪ [x3N /4+1 , 1], then from [24, Chapter 7, Theorem 7.4] we have H 1 -stability as well as the estimate that corresponds to (20), again with k ≤ r + 1. These arguments should be sufficient to prove an analogue of Theorem 2 on the interpolation error. At the end, in the analysis of the discretization error we have the additional term ε1 (η′ , χ ′ ) in (32) that can be estimated using the Cauchy–Schwarz inequality. 7. Numerical experiments In this section we present numerical results conducted in Mathematica 8 in order to verify theoretical estimates. For each test problem, we have a range R(ε1 , ε2 ) of perturbation parameters for which the condition (13) is fulfilled and the mesh is purely of the Bakhvalov-type. In general, this set of parameters constrains the evolution from a smooth (ε1 = 1) to a layered solution (ε1 ≪ 1). Hence, we have also carried out tests on a wider set of parameters, including even ε1 = ε2 = 1. When σj ≥ 1/4, we simply take a uniform mesh on the corresponding part of the domain. In all our experiments we have covered the cases when ε22 ≤ C ε1 , which resembles to the reaction–diffusion problem, as well as the case when ε22 ≥ C ε1 , that can be considered similar to the convection–diffusion problem. Here the condition for ε2 in Theorem 1 is always fulfilled. All integrals are approximated using the 5-point Gauss– Legendre quadrature rule. For the mesh parameters in the transition point (13) we choose p = 0.5 and τ = 4, since θ = 1/2. Varying these values influences on R(ε1 , ε2 ) such that increasing quotient τ /p reduces this set, while the rate of convergence changes insignificantly. The first test problem is
−ε1 u′′ + ε2 u′ + u = cos π x,
x ∈ Ω , u(0) = u(1) = 0,
with the exact solution u(x) = a cos π x + b sin π x + Ae−µ0 x + Be−µ1 (1−x) ,
(37)
where a=
ε1 π 2 + 1 , ε π + (ε1 π 2 + 1)2 2 2
A = −a
2
1 + e−µ1 1 − e−µ0 −µ1
,
b=
B=a
ε2 π , ε22 π 2 + (ε1 π 2 + 1)2
1 + e−µ0 1 − e−µ0 −µ1
,
µ0,1 =
∓ε2 +
ε22 + 4ε1
2ε1
.
For different values of perturbation parameters, the graph of the exact solution is depicted in Figs. 1 and 2. As already noted in Remark 2, for this problem the admissible set of perturbation parameters is
R(ε1 , ε2 ) = {(ε1 , ε2 ) : 0 < ε1 ≤ 10−5 , 0 < ε2 ≤ 10−3 }. For a fixed ε1 , ε2 and N, we have evaluated the energy norm of the error eNε1 ,ε2 = ∥u − uN ∥E , where u is the exact solution
given by (37) and uN represents its numerical approximation. In Table 1 we present the results for eNR =
max ε1 =10−5 ,...,10−12 ε2 =10−3 ,...,10−12
eNε1 ,ε2
where perturbation parameters take values from R(ε1 , ε2 ). The corresponding rate of convergence pNR =
ln eNR − ln e2N R
ln 2 in the third column of Table 1 confirms the first order of convergence of the Galerkin FEM with linear elements when applied to the Bakhvalov-type mesh (15). The last two columns of Table 1 show the values for eN =
max ε1 =1,10−1 ,...,10−12 ε2 =1,10−1 ,...,10−12
eNε1 ,ε2
and the corresponding rate of convergence pN =
ln eN − ln e2N
,
ln 2 which is in compliance with the main bound (36). We have also tested Galerkin FEM on the mesh (15), but with quadratic (r = 2) and cubic (r = 3) elements. The results in Table 2 show r-order uniform convergence of the method that has been indicated in Remark 5. To illustrate theoretical findings, in particular Theorem 4, we also consider the second test problem
−ε1 u′′ (x) + ε2 (3 − 2x2 )u′ (x) + u(x) = (1 + x)2 ,
x ∈ Ω , u(0) = u(1) = 0,
whose exact solution is not available. To estimate the energy error we use the following variant of the double mesh principle from [25].
318
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319 Table 1 Galerkin FEM with linear elements for the first test problem. N
eNR
pNR
eN
pN
24 25 26 27 28 29 210 211 212
3.591(−2) 1.763(−2) 8.790(−3) 4.396(−3) 2.197(−3) 1.096(−3) 5.469(−4) 2.730(−4) 1.364(−4)
1.027 1.004 1.000 1.001 1.003 1.003 1.002 1.001 –
2.229(−1) 1.612(−1) 9.618(−2) 5.096(−2) 2.588(−2) 1.299(−2) 6.503(−3) 3.252(−3) 1.626(−3)
0.468 0.745 0.916 0.977 0.994 0.999 1.000 1.000 –
Table 2 Galerkin FEM with quadratic and cubic elements for the first test problem. r =2
N
24 25 26 27 28 29 210 211 212
r =3
eNR
pNR
eN
pN
eNR
pNR
eN
pN
7.310(−3) 1.880(−3) 4.734(−4) 1.186(−4) 2.967(−5) 7.423(−6) 1.857(−6) 4.643(−7) 1.161(−7)
1.959 1.990 1.997 1.999 1.999 1.999 2.000 2.000 –
1.242(−1) 5.874(−2) 1.987(−2) 5.504(−3) 1.416(−3) 3.565(−4) 8.929(−5) 2.233(−5) 5.584(−6)
1.080 1.564 1.852 1.959 1.989 1.997 1.999 2.000 –
2.227(−3) 3.001(−4) 3.824(−5) 4.808(−6) 6.047(−7) 7.665(−8) 9.778(−9) 1.242(−9) 1.565(−10)
2.891 2.972 2.992 2.991 2.980 2.971 2.977 2.989 –
5.490(−2) 1.548(−2) 3.285(−3) 5.192(−4) 6.966(−5) 8.872(−6) 1.114(−6) 1.395(−7) 1.743(−8)
1.827 2.236 2.662 2.898 2.973 2.993 2.998 3.000 –
Table 3 Galerkin FEM with linear elements for the second test problem. N
dNR
N rR
dN
rN
24 25 26 27 28 29 210 211 212
6.045(−2) 3.049(−2) 1.527(−2) 7.640(−3) 3.820(−3) 1.910(−3) 9.550(−4) 4.775(−4) 2.388(−4)
0.987 0.997 0.999 1.000 1.000 1.000 1.000 1.000 –
8.845(−1) 6.089(−1) 3.424(−1) 1.768(−1) 8.915(−2) 4.467(−2) 2.235(−2) 1.117(−2) 5.587(−3)
0.539 0.830 0.953 0.988 0.997 0.999 1.000 1.000 –
For a fixed ε1 , ε2 and N, let u2N be the discrete solution obtained on a new mesh consisting from the mesh points (15) and its midpoints xi+1/2 = (xi+1 + xi )/2, i = 0, . . . , N − 1. Let dNε1 ,ε2 denote the energy norm ∥uN − u2N ∥E . In Table 3 we present the results for dNR =
max ε1 =10−6 ,...,10−12 ε2 =10−3 ,...,10−12
dNε1 ,ε2
and dN =
max ε1 =1,10−1 ,...,10−12 ε2 =1,10−1 ,...,10−12
dNε1 ,ε2 ,
N together with the corresponding rates of convergence rR and r N . Notice that now for the admissible set of perturbation parameters we have
R(ε1 , ε2 ) = {(ε1 , ε2 ) : 0 < ε1 ≤ 10−6 , 0 < ε2 ≤ 10−3 }. The results from Table 3 again clearly show parameter uniform convergence of the first order in the energy norm for both dNR and dN . Further research in this direction should include the error analysis in the maximum norm, application on other types of layer-adapted meshes, cf. [26,27], or different finite element methods. Moreover, recent literature on numerical solving of singularly perturbed problems addresses the weakness of the standard ε1 -weighted energy norm ∥ · ∥E and introduces balanced norms, cf. [28–31]. Since the boundary layer functions have the following properties
max{∥e
−pµ0 x
∥E , ∥e
−pµ1 (1−x)
∥E } ≤
1/4
ε22 ≤ C ε1 ,
1/4
ε1 ≤ C ε22 ,
C ε1 ,
C ε2 ,
it is necessary to further investigate the sharpness of the error bound from Theorem 4.
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
319
Acknowledgments The work has been supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia, Projects 174030 and III44006. The authors would like to thank the anonymous reviewers for their valuable comments and suggestions which improved the quality of the paper. References [1] A.H. Nayfeh, Introduction to Perturbation Techniques, John Wiley and Sons, New York, NY, USA, 1993. [2] R.E. O’Malley Jr., Singular Perturbation Methods for Ordinary Differential Equations, in: Applied Mathematical Sciences, Springer, New York, NY, USA, 1991, p. viii+225. [3] T. Linss, Layer-Adapted Meshes for Reaction–Convection–Diffusion Problems, Springer-Verlag, Berlin Heidelberg, 2010. [4] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, in: Springer Series in Computational Mathematics, vol. 24, Springer–Verlag, Berlin Heidelberg, 2008. [5] H.-G. Roos, Z. Uzelac, The SDFEM for a convection–diffusion problem with two small parameters, Comput. Methods Appl. Math. 3 (3) (2003) 443–458. [6] J.L. Gracia, E. O’Riordan, M.L. Pickett, A parameter robust higher order numerical method for a singularly perturbed two-parameter problem, Appl. Numer. Math. 56 (7) (2006) 962–980. [7] T. Linss, H.-G. Roos, Analysis of a finite-difference scheme for a singularly perturbed problem with two small parameters, J. Math. Anal. Appl. 289 (2004) 355–366. [8] S. Natesan, J.L. Gracia, C. Clavero, Singularly perturbed boundary-value problems with two small parameters — a defect correction approach, in: Proceedings of the International Conference on Boundary and Interior Layers — Computational and Asymptotic Methods, BAIL 2004, ONERA-Centre de Toulouse, France, pp. 1–6. [9] E. O’Riordan, M.L. Pickett, G.I. Shishkin, Singularly perturbed problems modelling reaction–convection–diffusion processes, Comput. Methods Appl. Math. 3 (3) (2003) 424–442. [10] K. Surla, Z. Uzelac, Lj. Teofanov, The discrete minimum principle for quadratic spline discretization of a singularly perturbed problem, Math. Comput. Simul. 79 (8) (2009) 2490–2505. [11] T. Linss, Layer-adapted meshes for one-dimensional reaction–convection–diffusion problems, J. Numer. Math. 12 (3) (2004) 193–205. [12] Dj. Herceg, Fourth-order finite-difference method for boundary value problems with two small parameters, Appl. Math. Comput. 218 (2) (2011) 616–627. [13] R. Vulanović, A higher-order scheme for quasilinear boundary value problems with two small parameters, Computing 67 (4) (2001) 287–303. [14] R. Vulanović, Special meshes and higher-order schemes for singularly perturbed boundary values problems, Novi Sad J. Math. 31 (1) (2001) 1–7. [15] T. Linss, A posteriori error estimation for a singularly perturbed problem with two small parameters, Int. J. Numer. Anal. Model. 7 (3) (2010) 491–506. [16] H.-G. Roos, Error estimates for linear finite elements on Bakhvalov-type meshes, Appl. Math. 51 (1) (2006) 63–72. [17] Lj. Teofanov, H.-G. Roos, An elliptic singularly perturbed problem with two parameters I: Solution decomposition, J. Comput. Appl. Math. 206 (2007) 1082–1097. [18] T. Linss, The necessity of Shishkin decompositions, Appl. Math. Lett. 14 (2001) 891–896. [19] P. Clement, Approximation by finite element functions using local regularization, RAIRO Anal. Numer. R-2 (1975) 77–84. [20] M. Dobrowolski, Finite Elemente, in: University Textbook of Würzburg, 1998. [21] N. Kopteva, M. Stynes, Approximation of derivatives in a convection–diffusion two-point boundary value problem, Appl. Numer. Math. 39 (1) (2001) 47–60. [22] H.-G. Roos, Robust numerical methods for singularly perturbed differential equations: A survey covering 2008–2012, ISRN Appl. Math. 2012 (2012) 30. [23] Lj. Teofanov, H.-G. Roos, An elliptic singularly perturbed problem with two parameters II: Robust finite element solution, J. Comput. Appl. Math. 212 (2) (2008) 374–389. [24] R.A. DeVore, G.G. Lorentz, Constructive Approximation, Springer-Verlag, Berlin Heidelberg, 1993. [25] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman and Hall/CRC Press, Boca Raton, FL, 2000. [26] R.G. Duran, A.L. Lombardi, Finite element approximation of convection diffusion problems using graded meshes, Appl. Numer. Math. 56 (2006) 1314–1325. [27] H.-G. Roos, T. Skalický, A comparison of the finite element method on Shishkin and Gartland-type meshes for convection–diffusion problems, CWI Q. 10 (1997) 277–300. [28] S. Franz, H.-G. Roos, Error estimation in a balanced norm for a convection–diffusion problem with two different boundary layers, Calcolo 51 (2014) 423–440. [29] R. Lin, M. Stynes, A balanced finite element method for singularly perturbed reaction–diffusion problems, SIAM J. Numer. Anal. 50 (5) (2012) 2729–2743. [30] H.-G. Roos, M. Schopf, Error estimation in energy norms: Is it necessary to fit the mesh to boundary layers, in: Numerical Analysis and Its Applications, in: Lecture Notes in Computer Science Volume, vol. 8236, 2013, pp. 95–109. [31] H.-G. Roos, M. Schopf, Convergence and stability in balanced norms of finite element methods on Shishkin meshes for reaction–diffusion problems, ZAMM Z. Angew. Math. Mech. (2014) 1–15.
Contents lists available at ScienceDirect
Journal of Computational and Applied Mathematics journal homepage: www.elsevier.com/locate/cam
A singularly perturbed problem with two parameters on a Bakhvalov-type mesh Mirjana Brdar a,∗ , Helena Zarin b a
Faculty of Technology, University of Novi Sad, Bulevar cara Lazara 1, 21000 Novi Sad, Serbia
b
Department of Mathematics and Informatics, Faculty of Sciences, University of Novi Sad, Trg Dositeja Obradovića 4, 21000 Novi Sad, Serbia
article
abstract
info
Article history: Received 11 December 2014 Received in revised form 21 June 2015
A singularly perturbed problem with two small parameters is considered. On a Bakhvalovtype mesh we prove uniform convergence of a Galerkin finite element method with piecewise linear functions. Arguments in the error analysis include interpolation error bounds for a Clément quasi-interpolant as well as discretization error estimates in an energy norm. Numerical experiments support theoretical findings. © 2015 Elsevier B.V. All rights reserved.
MSC: 65L11 65L20 65L60 65L70 Keywords: Singularly perturbed problem Two small parameters Galerkin finite element method Bakhvalov-type mesh Clément quasi-interpolant
1. Introduction We consider the following two-parameter singularly perturbed boundary value problem Lu(x) := −ε1 u′′ + ε2 b(x)u′ + c (x)u = f (x), u(0) = 0,
x ∈ Ω = (0, 1),
u(1) = 0,
(1)
with b(x) ≥ b0 > 0, c ( x) −
1 2
c (x) ≥ c0 > 0,
ε2 b′ (x) ≥ γ > 0,
x ∈ Ω,
x ∈ Ω,
(2) (3)
where b, c , f are smooth functions on Ω = [0, 1], b0 , c0 , γ are constants and 0 < ε1 , ε2 ≪ 1 are small perturbation parameters. The two-parameter singularly perturbed problems appear in many areas of physics and mechanics. Among these are the mathematical models of liquid crystal materials and chemical reactions, electrical networks, control theory. Some examples can be found in [1,2].
∗
Corresponding author. E-mail addresses: [email protected] (M. Brdar), [email protected] (H. Zarin).
http://dx.doi.org/10.1016/j.cam.2015.07.011 0377-0427/© 2015 Elsevier B.V. All rights reserved.
308
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
In case ε2 = 0, the problem (1)–(3) belongs to the class of reaction–diffusion differential equations whose solution has 1/2 1/2 two boundary layers of the widths O (ε1 | ln ε1 |). For ε2 = 1, problem (1)–(3) reduces to the convection–diffusion case with a solution exhibiting a layer at x = 1 of the width O (ε1 | ln ε1 |). There is a vast literature on numerical solving of these two problems classes, see for instance [3,4] and references therein. Here we are interested in the purely two-parameter case 0 < ε1 , ε2 ≪ 1 from [5], where a relationship between the perturbation parameters defines the layer structure at x = 0 and x = 1. The main goal in the construction of numerical methods for singularly perturbed problems is acquiring their uniform convergence (robustness) with respect to the perturbation parameters. For the one-dimensional problem (1)–(3), the authors in [6–10] consider various robust finite difference schemes on piecewise uniform meshes of Shishkin type. In [11], an analysis of a simple upwind scheme on arbitrary meshes has been presented. Bakhvalov-type meshes appear in [12–14] p+1/2 where for a special ε2 = ε1 , p > 0, uniform convergence has been proved for higher-order finite difference methods. The streamline-diffusion finite element method on the Shishkin mesh for (1)–(3) has been presented in [5] where authors prove (almost) second order of convergence in the maximum norm. For this method, a robust a posteriori error estimate in the maximum norm can be found in [15]. In this paper we analyze a linear finite element method on a mesh of Bakhvalov-type for the given problem. The first optimal result of convergence in the energy norm but for a one-dimensional convection–diffusion problem has been presented in [16], where an error analysis includes estimates for a quasi-interpolant. Here we adopt the approach from [16] and for the two-parameter singularly perturbed problem prove the first order of uniform convergence in the energy norm for the Galerkin method on the Bakhvalov-type mesh. The paper is organized as follows. In the next section we present properties of the exact solution as well as its derivatives, with special emphasis on the layer structure. In Section 3 we give a weak formulation and the Galerkin approach to the resolution of our problem. The layer-adapted mesh has been presented in Section 4, while in Section 5 we prove interpolation error bounds for the Clément approximation operator. Section 6 contains the main result on the uniform error bound in the energy norm. The last section is devoted to numerical experiments that illustrate our theoretical results. Notation 1. Throughout the paper, C will denote a generic positive constant independent of the perturbation parameters ε1 , ε2 and the number of degrees of freedom N. For a set D ⊂ R, we use the standard notation for Banach spaces Lp (D), Sobolev spaces W k,p (D), H k (D) = W k,2 (D), norms ∥ · ∥Lp (D) and seminorms | · |H k (D) . 2. Properties of the continuous problem In order to describe the layers in the solution of (1)–(3), the characteristic equation is introduced
−ε1 r 2 (x) + ε2 b(x)r (x) + c (x) = 0. This equation defines two continuous functions r0 , r1 : [0, 1] → R with r0 (x) < 0 and r1 (x) > 0. Let
µ0 = − max r0 (x), x∈[0,1]
µ1 = min r1 (x), x∈[0,1]
that is
µ0,1 = min
∓ε2 b(x) +
ε22 b2 (x) + 4ε1 c (x)
.
2ε1 Different properties of µ0 and µ1 can be derived, while here we emphasize the following, [17], x∈[0,1]
µ0 ≤ µ1 ,
1/2
1 max{µ− 0 , ε1 µ1 } ≤ C (ε2 + ε1 ),
−1/2
(4)
1/2 2
ε2 µ0 ≤ b0 ∥c ∥L∞ (Ω ) , ε2 (ε1 µ1 ) ≤ Cε . (5) The values µ0 and µ1 determine the decay of the boundary layers. Since |r0 | < r1 , the layer at x = 1 is stronger than the layer at x = 0. The following assertion from [3] provides more details on the behavior of the solution of (1)–(3) and its −1
derivatives. Theorem 1 ([3]). Let b, c , f ∈ C q (Ω ) for some q ≥ 1, and let p, κ ∈ (0, 1) be arbitrary. Assume that q∥b′ ∥L∞ (Ω ) ε2 ≤ κ(1 − p). Then
|u(k) (x)| ≤ C (1 + µk0 e−pµ0 x + µk1 e−pµ1 (1−x) ),
x ∈ Ω,
for 0 ≤ k ≤ q. Now, from [18] the solution u to the problem (1)–(3) has the representation u = S + E0 + E1 , with
|S (k) (x)| ≤ C , |E0(k) (x)| ≤ C µk0 e−pµ0 x , for x ∈ Ω and 0 ≤ k ≤ q.
|E1(k) (x)| ≤ C µk1 e−pµ1 (1−x) ,
(6)
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
309
Fig. 1. The exact solution (37) for ε1 = 10−3 , ε2 = 10−1 (left) and ε1 = 10−3 , ε2 = 10−6 (right).
Fig. 2. The exact solution (37) for ε1 = 10−12 , ε2 = 10−1 (left) and ε1 = 10−12 , ε2 = 10−12 (right).
Depending on the relation between the perturbation parameters, there are two extremal situations, cf. [3,5], −1/2
(i) when ε22 ≪ ε1 ≪ 1, then µ0 and µ1 behave like ε1 and the layers are similar to the reaction–diffusion case (see right graphs on Figs. 1 and 2), (ii) when ε1 ≪ ε22 ≪ 1, then µ0 = O (ε2−1 ) and µ1 = O (ε2 ε1−1 ); now µ1 is much larger than µ0 and the boundary layer in the vicinity of x = 1 is stronger than the boundary layer at x = 0 (see left graphs on Figs. 1 and 2). In other cases of the relations between the perturbation parameters, µ0 and µ1 vary between the previously given values. Notice that 1 ≪ µ0 ≤ µ1 . 3. Weak formulation and the Galerkin approach For the problem (1)–(3), the standard weak formulation is: find u ∈ H01 (Ω ) such that a(u, v) = (f , v),
∀v ∈ H01 (Ω ),
(7)
where (·, ·) denotes the standard scalar product in L (Ω ), and a(·, ·) is the bilinear form 2
a(u, v) := ε1 (u′ , v ′ ) + ε2 (bu′ , v) + (cu, v),
u, v ∈ H01 (Ω ).
Let N + 1, N ∈ N, be the number of mesh points 0 = x0 < x1 < · · · < xN −1 < xN = 1
(8)
that are arbitrarily located at the interval Ω . With V we denote the space of continuous, piecewise linear functions with the standard basis {φi : i = 0, 1, . . . , N } of hat functions, i.e. φi (xj ) = δij , i, j = 0, 1, . . . , N. Let V0N = V N ∩ H01 (Ω ). The Galerkin method is characterized by: find uN ∈ V0N such that N
a(uN , v N ) = (f , v N ),
∀v N ∈ V0N .
(9)
Under the assumption (3), the bilinear form a is coercive with respect to the energy norm
∥v∥2E := ε1 |v|2H 1 (Ω ) + ∥v∥2L2 (Ω ) . Hence, the problems (7) and (9) have unique solutions.
(10)
310
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
The Galerkin method (9) generates a finite difference scheme that can be written as 1
LN ui =
xi+1
h¯ i
f φi dx,
i = 1, 2, . . . , N − 1,
xi−1
with u0 = uN = 0, where ui ≈ u(xi ), hi = xi − xi−1 , i = 1, 2, . . . , N, are the local mesh step sizes and h¯ i = (hi + hi+1 )/2, i = 1, 2, . . . , N − 1. The discrete operator takes form LN ui = −ε1
D + ui − D − ui h¯ i
ui + 1 − ui ui − ui − 1 1 + ε2 βi+1 − βi−1 + (γi−1 ui−1 + γi ui + γi+1 ui+1 ) , h¯ i h¯ i h¯ i
where D+ and D− are the usual forward and backward difference operators and
βi−1 =
xi
bφi′−1 φi dx,
βi+1 =
xi
bφi′+1 φi dx,
xi
xi−1
γi−1 =
xi+1
c φi−1 φi dx,
γi =
xi+1
c φi2 dx,
γi+1 =
xi+1
c φi+1 φi dx.
xi
xi−1
xi−1
In the sequel we shall prove uniform convergence of the Galerkin method (9) on a specially chosen layer-adapted mesh. The error bound will be derived in the energy norm (10) starting from the triangle inequality applied to the error decomposition u − uN = η + χ , where uN = S N + E0N + E1N , η = ηS + η0 + η1 , χ = χS + χ0 + χ1 and
ηS = S − S I , π
χS = S I − S N , π
ηj = Ej − Ej + κ ,
χj = Ej −
N j
(11) EjN
−κ , N j
j = 0, 1.
(12)
Here S I ∈ V N is the standard interpolant of the regular solution component S. The functions Ejπ ∈ V N , j = 0, 1, are special
interpolants of the layer components Ej that will be introduced in Section 5, and κjN ∈ V N are certain correction terms that allow ηj (0) = ηj (1) = 0, j = 0, 1 (the necessity for the new interpolant in the analysis that includes layer components is already described in [16]). Hence, the key ingredient in our error analysis is the solution decomposition from Theorem 1 that also plays an important role in the construction of the mesh used for the discretization. 4. Bakhvalov-type layer-adapted mesh Let N ∈ N, N ≥ 8, be divisible by 4 and
σj =
τ pµj
ln µj <
1 4
,
j = 0, 1,
(13)
where τ ≥ 1 is a user-chosen parameter and p ∈ (0, 1) is the parameter from Theorem 1. Since 1 ≪ µ0 ≤ µ1 , throughout the paper we assume 1 −1 −θ µ− , 1 ≤ µ0 ≤ CN
for some 0 < θ ≤ 1,
(14)
which in practical implementations does not represent a restriction. The solution u of the problem (1)–(3) exhibits two boundary layers in the vicinity of x = 0 and x = 1 of the widths 1 −1 O (µ− 0 ln µ0 ) and O (µ1 ln µ1 ), respectively. Therefore, we define a mesh on Ω such that it is equidistant on Ω c with N /2 mesh subintervals, and gradually divided on Ω 0 and Ω 1 with N /4 mesh subintervals, where
Ω0 = (0, σ0 ),
Ωc = (σ0 , 1 − σ1 ),
Ω1 = (1 − σ1 , 1).
Due to the existence of two boundary layers, we use two mesh generating functions ϕ0 and ϕ1 , with
1 ϕ0 (t ) = − ln 1 − 4(1 − µ− 0 )t , 1 ϕ1 (t ) = − ln 1 − 4(1 − µ− 1 )(1 − t ) . The functions ϕ0 and ϕ1 are piecewise continuously differentiable and satisfy the following conditions
ϕ0 (0) = 0, ϕ0 (1/4) = ln µ0 , ϕ1 (3/4) = ln µ1 , ϕ1 (1) = 0.
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
311
The mesh points xi , i = 0, 1, . . . , N, are defined with
τ ϕ0 (ti ), p µ0 1 (1 − σ0 − σ1 ), xi = σ0 + 2 ti − 4 τ 1 − ϕ1 (ti ), p µ1
i = 0, 1, . . . , i= i=
N N
,
,
4
,
+ 1, . . . ,
4 4 3N 3N 4
N
4
3N 4
,
(15)
+ 1, . . . , N ,
where ti = i/N , i = 0, 1, . . . , N. The transition points, i.e. the points where the mesh switches from fine to coarse and vice versa, are chosen to be x N = σ0 , x 3N = 1 − σ1 , in order to have 4
4
exp(−pµ0 x N ) = µ0
−τ
and
4
exp(−pµ1 (1 − x 3N )) = µ−τ 1 . 4
With τ ≥ 1 we control the smallness of layer functions in the vicinity of the transition points. Remark 1. Following [13], instead of (14) we can introduce a stronger assumption 1
µj
ln µj ≤ CN −1 ,
j = 0, 1,
emanating from the practice that in the case of an equidistant mesh there are no mesh points placed inside the layer region. Then for small values of perturbation parameters we have (14) with θ = 1. Remark 2. There exists a threshold value µmin = µmin (τ , p) such that (13) is fulfilled when µ1 ≥ µ0 ≥ µmin > 1. For example, when p = 0.5 and τ = 1, 2, 3, 4, 5, we get µmin ≈ 27, 68, 114, 164, 215, respectively. Since µj ≫ 1, j = 0, 1, the inequality (13) is not a limiting condition. Nevertheless, as our numerical experiments from the last section show, the condition (13) influences the admissible range for the perturbation parameters for which the mesh is layer-adapted. For example, for the first test problem when b(x) = c (x) = 1, we get ε1 ≤ 10−5 , ε2 ≤ 10−3 , while in the second test problem with b(x) = 3 − 2x2 , c (x) = 1, we obtain ε1 ≤ 10−6 , ε2 ≤ 10−3 . As already remarked, the mesh (15) is dense in the layer regions. Moreover, it can be easily shown that hi > hi−1 ,
for i = 2, 3, . . . , N /4,
(16)
hi−1 > hi ,
for i = 3N /4 + 2, . . . , N .
(17)
In the subsequent analysis we shall also need the following properties of the mesh step sizes. Lemma 1. The mesh sizes hi = xi − xi−1 , i = 1, 2, . . . , N, satisfy
−1 C µ0 , hi ≤ CN −1 , −1 C µ1 ,
i = 1, 2, . . . , N /4 − 1, i = N /4, . . . , 3N /4 + 1, i = 3N /4 + 2, . . . , N .
Proof. The functions ϕ0 and ϕ0′ are both strictly increasing on [0, 1/4]. Now, if i = 1, 2, . . . , N /4 − 1, by the mean value theorem there exists ξi ∈ (ti−1 , ti ) such that hi =
τ −1 ′ τ −1 ′ τ −1 4 1 N ϕ0 (ξi ) < N ϕ0 (t N −1 ) < N = C µ− 0 . 4 p µ0 p µ0 p µ0 1 − 4t N −1 4
For i = N /4 and ξi ∈ (t N −1 , t N ), 4
hN = 4
4
1 τ −1 4(1 − µ− τ −1 ′ 0 ) N ϕ0 (ξi ) ≤ N ≤ CN −1 . 1 1 pµ0 p µ0 1 − 4(1 − µ− 0 )4
Trivially, hi = 2N −1 (1 − σ0 − σ1 ) ≤ CN −1 ,
i=
N 4
+ 1, . . . ,
3N 4
.
The mesh step sizes on the interval Ω 1 can be analogously estimated.
Notice that hi ≤ CN −1 for all indices i. Now using xi−1 = xi − hi and the previous lemma, we easily get the following
312
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
Corollary 1. For the mesh points (15) we have e−µ0 xi−1 ≤ C e−µ0 xi ,
i = 1, 2, . . . ,
e−µ1 xi−1 ≤ C e−µ1 xi ,
i=
3N 4
N 4
− 1,
+ 2, . . . , N .
It is important to notice that from (16) it follows h N > h N −1 4
4
1 τ 4N −1 (1 − µ− 1 −1/2 0 ) = ln 1 + ≥ C µ− . 0 N 1 p µ0 1 − 4(1 − µ− 0 )t N −1 4
1 −1/2 Similarly, hi ≥ C µ− for i = 3N /4 + 1, 3N /4 + 2. Also, hi ≥ CN −1 for i = N /4 + 1, . . . , 3N /4. 1 N
Lemma 2. The mesh step sizes hi = xi − xi−1 , i = 1, 2, . . . , N, satisfy hi +1 hi hi −1 hi
≤ C,
i = 1, 2, . . . ,
≤ C,
i=
3N 4
N 4
− 2,
+ 3, . . . , N .
Proof. For each i = 1, 2, . . . , N /4 − 2 there exists ξi ∈ (ti−1 , ti ) such that hi =
τ −1 ′ N ϕ0 (ξi ). p µ0
Now hi +1 hi
=
1 1 − 4(1 − µ− 0 )ξi 1 1 − 4(1 − µ− 0 )ξi+1
8N
≤ 1+
−1
≤1+
(1 − µ0 )
1 8N −1 (1 − µ− 0 ) 1 1 − 4(1 − µ− 0 )ξi+1
−1
1 µ0 + 4N −1 (1 − µ− 0 )
−1
≤ C.
The other inequality can be derived analogously.
5. Interpolation error estimates for Clément approximation operator In the error analysis that includes the terms from (11)–(12), beside using the usual nodal interpolation operator from V N , we introduce a quasi-interpolant of Clément-type from [19] enriched with certain stability properties (more on this operator can be found in [20]). We first cite [16] and sketch the construction of the quasi-interpolant. For every mesh point xi , i = 0, 1, . . . , N, of a general regular mesh (8), by ∆i we denote a macroelement (xi−1 , xi+1 ), i = 1, . . . , N − 1, with ∆0 = (x0 , x1 ), ∆N = (xN −1 , xN ). Let P1 (∆i ) be the space of linear polynomials on ∆i . To a given function v ∈ L1 (Ω ) we define the local L2 -projection pi ∈ P1 (∆i ) by
∆i
pi s dx =
∆i
v s dx,
∀s ∈ P1 (∆i ).
Now, the quasi-interpolant v π ∈ V N of the function v ∈ L1 (Ω ) is
v π (x) :=
N
pi (xi )φi (x),
i=0
where φi represent hat functions. Notice that in general v π (xi ) ̸= v(xi ), i = 0, 1, . . . , N. For such a constructed operator, for a function v ∈ W k,p (∆) we have the following estimates
∥v − v π ∥Lp (e) ≤ CH k ∥Dk v∥Lp (∆) ,
p ∈ [1, ∞], k = 1, 2,
(18)
∥v π ∥Lp (e) ≤ C ∥v∥Lp (∆) , p ∈ [1, ∞]. (19) Here e = (xi−1 , xi ) denotes an element and ∆ is the union of the two corresponding macroelements ∆i−1 and ∆i , with the diameter H. We shall also use the inequality
|v − v π |H 1 (e) ≤ CH k−1 |v|H k (∆) , k = 1, 2, (20) which is valid only for elements e = (xi−1 , xi ), i = 1, 2, . . . , N /4 − 2 and i = 3N /4 + 3, . . . , N, according to [20] and Lemma 2.
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
313
Lemma 3. On the Bakhvalov-type mesh (15), for τ ≥ α/2, α > 0, we have 2
min{hi+1 µ0 , 1}e− α pµ0 xi ≤ CN −1 ,
i = 1, 2, . . . ,
2
min{hi µ1 , 1}e− α pµ1 (1−xi ) ≤ CN −1 ,
i=
3N
N 4
− 1,
+ 1, . . . , N − 1.
4
Proof. The idea for the proof is similar to [21,16]. First notice that for i = 1, 2, . . . , N /4 − 1 we have
min{hi+1 µ0 , 1}e
− α2 pµ0 xi
= min
1 4τ N −1 (1 − µ− 0 )
p(1 − 4(1 − µ0 )ξi+1 ) −1
,1
1 1 − 4(1 − µ− 0 )ti
2ατ
for some ξi+1 ∈ (ti , ti+1 ). From the assumption τ ≥ α/2, we get
min{hi+1 µ0 , 1}e
− α2 pµ0 xi
≤ min
4τ N −1 (1 − µ0−1 ) p(1 − 4(1 − µ0 )ti+1 ) −1
= min
1 4 τ N − 1 ( 1 − µ− 0 ) 1 p(1 − 4(1 − µ− 0 )ti+1 )
,1
1 1 − 4(1 − µ− 0 )ti
,1
1 1 −1 1 − 4(1 − µ− (1 − µ− 0 )ti+1 + 4N 0 ) .
The last expression can be written as
τ ω1 min , 1 (ω1 + ω2 ), p ω2
(21)
where 1 −1 ω1 = 4N −1 (1 − µ− , 0 ) ≤ CN
1 ω2 = 1 − 4(1 − µ− 0 )ti+1 .
If in (21) the minimum is 1, then
min
τ ω1 τ ω1 , 1 (ω1 + ω2 ) = ω1 + ω2 ≤ ω1 + ≤ CN −1 , p ω2 p
otherwise
τ ω1 ω1 τ ω1 , 1 (ω1 + ω2 ) = + 1 ≤ CN −1 . min p ω2 p ω2
The analysis for other values of the index i can be derived similarly.
Since the quasi-interpolant in general does not satisfy boundary conditions, we introduce correction functions κjN ∈ V N , j = 0, 1, defined as
κjN (xi ) =
(Ejπ − Ej )(xi ), 0,
i = 0, N , i = 1, 2, . . . , N − 1.
(22)
Now, recalling the earlier notation for ηj in (12), we can start with the interpolation error analysis. Lemma 4. Let the assumptions of Theorem 1 hold. If κjN are correction functions from (22) and Ejπ are Clément quasi-interpolants
of the solution components Ej , j = 0, 1, respectively, then the functions ηj = Ej − Ejπ + κjN on the Bakhvalov-type mesh (15) with
τ ≥ max{5/2, 2θ −1 } satisfy −1/2
∥ηj ∥L2 (Ω ) ≤ C µj
N −2 ,
j = 0, 1.
Proof. First we estimate the layer component η0 . Let e = (xi−1 , xi ), i = 2, 3, . . . , N /4 − 1, be some subinterval on the fine mesh and ∆ = (xi−2 , xi+1 ) ⊂ [0, σ0 ] an appropriate macroelement. The triangle inequality, stability estimate (19) and the properties of E0 in (6) together yield 1 −2pµ0 xi−2 ∥E0 − E0π ∥2L2 (e) ≤ C ∥E0 ∥2L2 (∆) ≤ C µ− . 0 e
(23)
On the other hand, if we use (18), (16) and estimate the integral by the maximum of the positive and monotonically decreasing integrand times the measure of the domain, we get
∥E0 − E0π ∥2L2 (e) ≤ CH 4 ∥E0′′ ∥2L2 (∆) ≤ Ch4i+1 µ40 ≤ Ch5i+1 µ40 e−2pµ0 xi−2 .
xi+1
e−2pµ0 x dx
xi−2
(24)
314
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
Combining both estimates (23) and (24) with Corollary 1, we have
5 2 1 1 −5 min{hi+1 µ0 , 1}e− 5 pµ0 xi ≤ C µ− ∥E0 − E0π ∥2L2 (e) ≤ C µ− . 0 0 N
(25)
The last inequality follows from Lemma 3 with α = 5. When e = (x0 , x1 ), analogous arguments yield (25) with i = 1. Also 1/2
1/2
∥κ0N ∥L2 (x0 ,x1 ) ≤ Ch1 |E0 (x0 ) − E0π (x0 )| ≤ Ch1 ∥E0 − E0π ∥L∞ (x0 ,x1 ) 5/2
1/2
≤ Ch1 (h1 + h2 )2 ∥E0′′ ∥L∞ (x0 ,x2 ) ≤ Ch1 µ20 −1/2
≤ C µ0
N −5/2 ,
(26)
1 −1 since h1 ≤ C µ− . Now 0 N N
N
−1 4
2 0 L2 (x i−1 ,xi )
∥η ∥
≤C
i =1
−1 4
∥E0 − E0π ∥2L2 (x
i−1 ,xi )
+ C ∥κ0N ∥2L2 (x
0 ,x 1 )
i=1 1 −4 ≤ C µ− . 0 N
(27)
For i ≥ N /4, we again use the triangle inequality, (19), (6) and (14) in order to get N
∥E0 − E0π ∥2L2 (x
i−1 ,xi )
≤C
i= N 4
N −1
∥E0 ∥2L2 (x
i−2 ,xi+1 )
+ C ∥E0 ∥2L2 (x
N −2 ,xN )
i= N 4
xN
≤C xN 4
1 e−2pµ0 x dx ≤ C µ− 0 e
−2pµ0 x N 4
−2
−2
2τ −1 2τ 1 1 1 ≤ C µ− µ0 + N − 1 = C µ− 1 − 4(1 − µ− 0 0 0 )t N −2 4
1 −2θτ ≤ C µ− . 0 N
Since κ0N is exponentially small on (xN −1 , xN ), the last inequality together with (27) and −θ τ ≤ −2 finally imply −1/2
∥η0 ∥L2 (Ω ) ≤ C µ0
N −2 .
Similar arguments lead to the interpolation error bound for the solution component E1 .
Lemma 5. Let the assumptions of Theorem 1 hold. If κjN are correction functions from (22) and Ejπ are Clément quasi-interpolants of the solution components Ej , j = 0, 1, respectively, then the functions ηj = Ej − Ejπ + κjN on the Bakhvalov-type mesh (15) with
τ ≥ 2θ −1 satisfy
1/2
|ηj |H 1 (Ω ) ≤ C µj N −1 ,
j = 0, 1.
Proof. Similarly as in the previous lemma, we start with the solution component E0 . First let e = (xi−1 , xi ), when i = 2, 3, . . . , N /4 − 2. Analogously to (23) we derive
|E0 − E0π |2H 1 (e) ≤ C µ0 e−2pµ0 xi−2 ,
(28)
where we have used (20) with k = 1. If we apply again (20) but with k = 2, then we have
|E0 − E0π |2H 1 (e) ≤ Ch3i+1 µ40 e−2pµ0 xi−2 and
3 2 |E0 − E0π |2H 1 (e) ≤ C µ0 min{hi+1 µ0 , 1}e− 3 pµ0 xi ≤ C µ0 N −3 , which follows from Lemma 3 with α = 3. The same estimate holds for the element e = (x0 , x1 ). Summation for i = 1, 2, . . . , N /4 − 2 yields N
−2 4 i =1
|E0 − E0π |2H 1 (x
i−1 ,xi )
≤ C µ0 N − 2 .
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
315
The correction function can be bounded with the inverse inequality and (26). Hence
|κ0N |H 1 (x0 ,x1 ) ≤
C h1
3/2
1/2
∥κ0N ∥L2 (x0 ,x1 ) ≤ Ch1 µ20 ≤ C µ0 N −3/2
and consequently N
−2 4
|η0 |2H 1 (x
i−1 ,xi )
≤ C µ0 N − 2 .
(29)
i=1
For i = N /4 − 1, . . . , 3N /4, the triangle inequality, the inverse inequality as well as the stability estimate (19) result in
|E0 − E0π |H 1 (e) ≤ |E0 |H 1 (e) +
1 hi
∥E0π ∥L2 (e) ≤ |E0 |H 1 (e) +
C hi
∥E0 ∥L2 (∆) .
Based on (6), Corollary 1, lower bounds for hi , (4) and (14), we have 3N 4
|E0 − E0π |2H 1 (x
i−1 ,xi )
1 2 −1 −1 2 τ ≤ C (µ0 + µ− ) ≤ C µ0 N 2(1−θ τ ) . 0 N )(µ0 + N
i= N −1 4
For i > 3N /4, it is sufficient to use the smallness of the component E0 and the correction function κ0N on Ω1 . Finally, we conclude 1/2
1/2
1/2
|η0 |H 1 (Ω ) ≤ C µ0 N −1 + C µ0 N 1−θτ ≤ C µ0 N −1 , since 1 − θ τ ≤ −1. The estimate for |η1 |H 1 (Ω ) can be similarly derived.
The following assertion represents the main result on the interpolation error. Theorem 2. Let the assumptions of Theorem 1 hold. On the Bakhvalov-type mesh (15) with τ ≥ max{5/2, 2θ −1 }, the functions ηS = S − S I and ηj = Ej − Ejπ + κjN can be estimated with 1/2
∥ηS ∥E ≤ C ε1 N −1 + CN −2 , ∥ηj ∥E ≤ C (ε2 + ε
1/2 1/2 −1 N 1
)
(30)
,
j = 0, 1,
(31)
where S I ∈ V N is the standard interpolant of S , Ejπ are quasi-interpolants of Ej , j = 0, 1, and κjN , j = 0, 1, are functions from (22). Proof. The bounds for ∥ηj ∥E are a direct consequence of the previous two lemmata, (4) and −1/2
∥ηj ∥E ≤ C (ε1 µ1 )1/2 N −1 + C µ0
N −2 ,
j = 0, 1.
The estimate for the energy norm of ηS = S − S I is already known (see for instance [3]), and follows from
∥ηS ∥L2 (Ω ) ≤ CN −2 and |ηS |H 1 (Ω ) ≤ CN −1 , which completes the proof.
6. Convergence in the energy norm In the sequel we shall first derive the estimate for ∥χ∥E , which together with Theorem 2 will provide us with the final error bound. Theorem 3. Let the assumptions of Theorems 1 and 2 hold. The discretization error for the problem (1)–(3) on the Bakhvalov-type mesh (15) can be estimated with 1/2
∥χ∥E ≤ C ε2 N −1 + CN −2 . Proof. The bilinear form a(·, ·) is coercive with respect to the energy norm. Let α = min{1, γ }. Coercivity and Galerkin orthogonality imply
α∥χ∥2E ≤ a(χ , χ ) = −a(η, χ ) = −ε2 (bη′ , χ ) − (c η, χ ).
(32)
316
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
The last term in (32) can be directly estimated applying the Cauchy–Schwarz inequality, Lemma 4 and the L2 -norm bound on ηS in the proof of Theorem 2 and (4). Therefore
|(c η, χ )| ≤ C ∥η∥L2 (Ω ) ∥χ∥L2 (Ω ) −1/2
≤ C (N −2 + µ0
N −2 )∥χ∥E ≤ CN −2 ∥χ∥E .
(33)
The remaining term in (32) we analyze according to the splitting η = (ηS + η0 ) + η1 . First we get
|ε2 (b(ηS′ + η0′ ), χ )| ≤ C ε2 (|ηS |H 1 (Ω ) + |η0 |H 1 (Ω ) )∥χ∥L2 (Ω ) 1/2
1/2
≤ C ε2 (N −1 + µ0 N −1 )∥χ∥E ≤ C ε2 N −1 ∥χ ∥E .
(34)
Here we have used the Cauchy–Schwarz inequality, Lemma 5 and the H 1 -seminorm bound on ηS in the proof of Theorem 2 and (5). These arguments are not sufficient for estimating the expression with η1 . Instead we apply integration by parts and since η1 satisfies homogeneous Dirichlet boundary conditions, we can derive
|ε2 (bη1′ , χ )| = ε2 |(η1 , b′ χ + bχ ′ )| ≤ C ε2 ∥η1 ∥L2 (Ω ) (∥χ∥L2 (Ω ) + |χ|H 1 (Ω ) ) −1/2
≤ C ε2 µ1 −1/2
Due to ε2 (ε1 µ1 )
≤ Cε
1/2 2
−1/2
N − 2 ( 1 + ε1
)∥χ ∥E .
from (5), we have
1/2
|ε2 (bη1′ , χ )| ≤ C ε2 N −2 ∥χ∥E . Collecting (33)–(35) into (32) and dividing the inequality with ∥χ∥E , we complete the proof.
(35)
The final error estimate in the last theorem represents a direct consequence of Theorems 2, 3 and (4). Theorem 4. Let uN be the Galerkin approximation of the solution u of (1)–(3) on the Bakhvalov-type mesh (15). If the assumptions of Theorem 2 are satisfied, then 1/2
∥u − uN ∥E ≤ C (ε2 + ε1 )1/2 N −1 + CN −2 .
(36)
Remark 3. Under the assumption ε1 ≤ CN −1 , the author in [16] proved ∥u − uN ∥E ≤ CN −1 for the convection–diffusion problem. We can consider (36) as the generalization of the aforementioned result, since for ε2 = 1 we place a uniform mesh 1 −1 on [0, 1 − σ1 ] and use µ− . 1 ≤ C ε1 ≤ CN Remark 4. In the reaction–diffusion case, i.e. ε2 = 0, the linear Galerkin finite element approximation uN on an appropriate Shishkin mesh satisfies, [22], 1/4
∥u − uN ∥E ≤ C ε1 N −1 ln N + CN −2 , where N + 1 is the total number of nodes. The error estimate (36) remains valid in this setting even though Theorem 4 cannot be applied. Also, we would like to draw attention to a result in [23] for the two-dimensional two-parameter singularly perturbed problem and the Galerkin FEM on the Shishkin mesh. The energy norm error estimate is similar to (36), i.e. 1/2
∥u − uN ∥E ≤ C (ε2 + ε1 )1/2 N −1 ln N + CN −2 , where again a logarithmic factor stems from the choice of transition points. Here, N + 1 is the number of nodes located in each direction. Remark 5. As numerical tests from the next section suggest, a generalization with higher order polynomials should be possible. Let V N ,r = {v ∈ C (Ω ) : v|(xi−1 ,xi ) ∈ Pr (xi−1 , xi ), i = 1, 2, . . . , N } be the finite-dimensional space of element-wise polynomials of the highest degree r ≥ 1. For a finite element approximation uN ∈ V N ,r ∩ H01 (Ω ) we conjecture an error bound 1/2
∥u − uN ∥E ≤ C (ε2 + ε1 )1/2 N −r + CN −(r +1) . The proof on the layer-adapted mesh (15) with τ ≥ max{r + 3/2, 2r θ −1 } requires an approximation operator that exhibits properties similar to (18)–(20). In [24] such an operator has been constructed as a mapping from Lp (Ω ) onto the space of spline functions of the order r + 1. If for an element e = (xi−1 , xi ) we define ∆ = (xi−r −1 , xi+r ), then estimates
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
317
of the type (18) and (19) can be used with k ≤ r + 1, cf. [24, Chapter 5, Theorem 4.4 and Theorem 4.5]. Moreover, if ∆ ⊂ [0, xN /4−1 ] ∪ [x3N /4+1 , 1], then from [24, Chapter 7, Theorem 7.4] we have H 1 -stability as well as the estimate that corresponds to (20), again with k ≤ r + 1. These arguments should be sufficient to prove an analogue of Theorem 2 on the interpolation error. At the end, in the analysis of the discretization error we have the additional term ε1 (η′ , χ ′ ) in (32) that can be estimated using the Cauchy–Schwarz inequality. 7. Numerical experiments In this section we present numerical results conducted in Mathematica 8 in order to verify theoretical estimates. For each test problem, we have a range R(ε1 , ε2 ) of perturbation parameters for which the condition (13) is fulfilled and the mesh is purely of the Bakhvalov-type. In general, this set of parameters constrains the evolution from a smooth (ε1 = 1) to a layered solution (ε1 ≪ 1). Hence, we have also carried out tests on a wider set of parameters, including even ε1 = ε2 = 1. When σj ≥ 1/4, we simply take a uniform mesh on the corresponding part of the domain. In all our experiments we have covered the cases when ε22 ≤ C ε1 , which resembles to the reaction–diffusion problem, as well as the case when ε22 ≥ C ε1 , that can be considered similar to the convection–diffusion problem. Here the condition for ε2 in Theorem 1 is always fulfilled. All integrals are approximated using the 5-point Gauss– Legendre quadrature rule. For the mesh parameters in the transition point (13) we choose p = 0.5 and τ = 4, since θ = 1/2. Varying these values influences on R(ε1 , ε2 ) such that increasing quotient τ /p reduces this set, while the rate of convergence changes insignificantly. The first test problem is
−ε1 u′′ + ε2 u′ + u = cos π x,
x ∈ Ω , u(0) = u(1) = 0,
with the exact solution u(x) = a cos π x + b sin π x + Ae−µ0 x + Be−µ1 (1−x) ,
(37)
where a=
ε1 π 2 + 1 , ε π + (ε1 π 2 + 1)2 2 2
A = −a
2
1 + e−µ1 1 − e−µ0 −µ1
,
b=
B=a
ε2 π , ε22 π 2 + (ε1 π 2 + 1)2
1 + e−µ0 1 − e−µ0 −µ1
,
µ0,1 =
∓ε2 +
ε22 + 4ε1
2ε1
.
For different values of perturbation parameters, the graph of the exact solution is depicted in Figs. 1 and 2. As already noted in Remark 2, for this problem the admissible set of perturbation parameters is
R(ε1 , ε2 ) = {(ε1 , ε2 ) : 0 < ε1 ≤ 10−5 , 0 < ε2 ≤ 10−3 }. For a fixed ε1 , ε2 and N, we have evaluated the energy norm of the error eNε1 ,ε2 = ∥u − uN ∥E , where u is the exact solution
given by (37) and uN represents its numerical approximation. In Table 1 we present the results for eNR =
max ε1 =10−5 ,...,10−12 ε2 =10−3 ,...,10−12
eNε1 ,ε2
where perturbation parameters take values from R(ε1 , ε2 ). The corresponding rate of convergence pNR =
ln eNR − ln e2N R
ln 2 in the third column of Table 1 confirms the first order of convergence of the Galerkin FEM with linear elements when applied to the Bakhvalov-type mesh (15). The last two columns of Table 1 show the values for eN =
max ε1 =1,10−1 ,...,10−12 ε2 =1,10−1 ,...,10−12
eNε1 ,ε2
and the corresponding rate of convergence pN =
ln eN − ln e2N
,
ln 2 which is in compliance with the main bound (36). We have also tested Galerkin FEM on the mesh (15), but with quadratic (r = 2) and cubic (r = 3) elements. The results in Table 2 show r-order uniform convergence of the method that has been indicated in Remark 5. To illustrate theoretical findings, in particular Theorem 4, we also consider the second test problem
−ε1 u′′ (x) + ε2 (3 − 2x2 )u′ (x) + u(x) = (1 + x)2 ,
x ∈ Ω , u(0) = u(1) = 0,
whose exact solution is not available. To estimate the energy error we use the following variant of the double mesh principle from [25].
318
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319 Table 1 Galerkin FEM with linear elements for the first test problem. N
eNR
pNR
eN
pN
24 25 26 27 28 29 210 211 212
3.591(−2) 1.763(−2) 8.790(−3) 4.396(−3) 2.197(−3) 1.096(−3) 5.469(−4) 2.730(−4) 1.364(−4)
1.027 1.004 1.000 1.001 1.003 1.003 1.002 1.001 –
2.229(−1) 1.612(−1) 9.618(−2) 5.096(−2) 2.588(−2) 1.299(−2) 6.503(−3) 3.252(−3) 1.626(−3)
0.468 0.745 0.916 0.977 0.994 0.999 1.000 1.000 –
Table 2 Galerkin FEM with quadratic and cubic elements for the first test problem. r =2
N
24 25 26 27 28 29 210 211 212
r =3
eNR
pNR
eN
pN
eNR
pNR
eN
pN
7.310(−3) 1.880(−3) 4.734(−4) 1.186(−4) 2.967(−5) 7.423(−6) 1.857(−6) 4.643(−7) 1.161(−7)
1.959 1.990 1.997 1.999 1.999 1.999 2.000 2.000 –
1.242(−1) 5.874(−2) 1.987(−2) 5.504(−3) 1.416(−3) 3.565(−4) 8.929(−5) 2.233(−5) 5.584(−6)
1.080 1.564 1.852 1.959 1.989 1.997 1.999 2.000 –
2.227(−3) 3.001(−4) 3.824(−5) 4.808(−6) 6.047(−7) 7.665(−8) 9.778(−9) 1.242(−9) 1.565(−10)
2.891 2.972 2.992 2.991 2.980 2.971 2.977 2.989 –
5.490(−2) 1.548(−2) 3.285(−3) 5.192(−4) 6.966(−5) 8.872(−6) 1.114(−6) 1.395(−7) 1.743(−8)
1.827 2.236 2.662 2.898 2.973 2.993 2.998 3.000 –
Table 3 Galerkin FEM with linear elements for the second test problem. N
dNR
N rR
dN
rN
24 25 26 27 28 29 210 211 212
6.045(−2) 3.049(−2) 1.527(−2) 7.640(−3) 3.820(−3) 1.910(−3) 9.550(−4) 4.775(−4) 2.388(−4)
0.987 0.997 0.999 1.000 1.000 1.000 1.000 1.000 –
8.845(−1) 6.089(−1) 3.424(−1) 1.768(−1) 8.915(−2) 4.467(−2) 2.235(−2) 1.117(−2) 5.587(−3)
0.539 0.830 0.953 0.988 0.997 0.999 1.000 1.000 –
For a fixed ε1 , ε2 and N, let u2N be the discrete solution obtained on a new mesh consisting from the mesh points (15) and its midpoints xi+1/2 = (xi+1 + xi )/2, i = 0, . . . , N − 1. Let dNε1 ,ε2 denote the energy norm ∥uN − u2N ∥E . In Table 3 we present the results for dNR =
max ε1 =10−6 ,...,10−12 ε2 =10−3 ,...,10−12
dNε1 ,ε2
and dN =
max ε1 =1,10−1 ,...,10−12 ε2 =1,10−1 ,...,10−12
dNε1 ,ε2 ,
N together with the corresponding rates of convergence rR and r N . Notice that now for the admissible set of perturbation parameters we have
R(ε1 , ε2 ) = {(ε1 , ε2 ) : 0 < ε1 ≤ 10−6 , 0 < ε2 ≤ 10−3 }. The results from Table 3 again clearly show parameter uniform convergence of the first order in the energy norm for both dNR and dN . Further research in this direction should include the error analysis in the maximum norm, application on other types of layer-adapted meshes, cf. [26,27], or different finite element methods. Moreover, recent literature on numerical solving of singularly perturbed problems addresses the weakness of the standard ε1 -weighted energy norm ∥ · ∥E and introduces balanced norms, cf. [28–31]. Since the boundary layer functions have the following properties
max{∥e
−pµ0 x
∥E , ∥e
−pµ1 (1−x)
∥E } ≤
1/4
ε22 ≤ C ε1 ,
1/4
ε1 ≤ C ε22 ,
C ε1 ,
C ε2 ,
it is necessary to further investigate the sharpness of the error bound from Theorem 4.
M. Brdar, H. Zarin / Journal of Computational and Applied Mathematics 292 (2016) 307–319
319
Acknowledgments The work has been supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia, Projects 174030 and III44006. The authors would like to thank the anonymous reviewers for their valuable comments and suggestions which improved the quality of the paper. References [1] A.H. Nayfeh, Introduction to Perturbation Techniques, John Wiley and Sons, New York, NY, USA, 1993. [2] R.E. O’Malley Jr., Singular Perturbation Methods for Ordinary Differential Equations, in: Applied Mathematical Sciences, Springer, New York, NY, USA, 1991, p. viii+225. [3] T. Linss, Layer-Adapted Meshes for Reaction–Convection–Diffusion Problems, Springer-Verlag, Berlin Heidelberg, 2010. [4] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, in: Springer Series in Computational Mathematics, vol. 24, Springer–Verlag, Berlin Heidelberg, 2008. [5] H.-G. Roos, Z. Uzelac, The SDFEM for a convection–diffusion problem with two small parameters, Comput. Methods Appl. Math. 3 (3) (2003) 443–458. [6] J.L. Gracia, E. O’Riordan, M.L. Pickett, A parameter robust higher order numerical method for a singularly perturbed two-parameter problem, Appl. Numer. Math. 56 (7) (2006) 962–980. [7] T. Linss, H.-G. Roos, Analysis of a finite-difference scheme for a singularly perturbed problem with two small parameters, J. Math. Anal. Appl. 289 (2004) 355–366. [8] S. Natesan, J.L. Gracia, C. Clavero, Singularly perturbed boundary-value problems with two small parameters — a defect correction approach, in: Proceedings of the International Conference on Boundary and Interior Layers — Computational and Asymptotic Methods, BAIL 2004, ONERA-Centre de Toulouse, France, pp. 1–6. [9] E. O’Riordan, M.L. Pickett, G.I. Shishkin, Singularly perturbed problems modelling reaction–convection–diffusion processes, Comput. Methods Appl. Math. 3 (3) (2003) 424–442. [10] K. Surla, Z. Uzelac, Lj. Teofanov, The discrete minimum principle for quadratic spline discretization of a singularly perturbed problem, Math. Comput. Simul. 79 (8) (2009) 2490–2505. [11] T. Linss, Layer-adapted meshes for one-dimensional reaction–convection–diffusion problems, J. Numer. Math. 12 (3) (2004) 193–205. [12] Dj. Herceg, Fourth-order finite-difference method for boundary value problems with two small parameters, Appl. Math. Comput. 218 (2) (2011) 616–627. [13] R. Vulanović, A higher-order scheme for quasilinear boundary value problems with two small parameters, Computing 67 (4) (2001) 287–303. [14] R. Vulanović, Special meshes and higher-order schemes for singularly perturbed boundary values problems, Novi Sad J. Math. 31 (1) (2001) 1–7. [15] T. Linss, A posteriori error estimation for a singularly perturbed problem with two small parameters, Int. J. Numer. Anal. Model. 7 (3) (2010) 491–506. [16] H.-G. Roos, Error estimates for linear finite elements on Bakhvalov-type meshes, Appl. Math. 51 (1) (2006) 63–72. [17] Lj. Teofanov, H.-G. Roos, An elliptic singularly perturbed problem with two parameters I: Solution decomposition, J. Comput. Appl. Math. 206 (2007) 1082–1097. [18] T. Linss, The necessity of Shishkin decompositions, Appl. Math. Lett. 14 (2001) 891–896. [19] P. Clement, Approximation by finite element functions using local regularization, RAIRO Anal. Numer. R-2 (1975) 77–84. [20] M. Dobrowolski, Finite Elemente, in: University Textbook of Würzburg, 1998. [21] N. Kopteva, M. Stynes, Approximation of derivatives in a convection–diffusion two-point boundary value problem, Appl. Numer. Math. 39 (1) (2001) 47–60. [22] H.-G. Roos, Robust numerical methods for singularly perturbed differential equations: A survey covering 2008–2012, ISRN Appl. Math. 2012 (2012) 30. [23] Lj. Teofanov, H.-G. Roos, An elliptic singularly perturbed problem with two parameters II: Robust finite element solution, J. Comput. Appl. Math. 212 (2) (2008) 374–389. [24] R.A. DeVore, G.G. Lorentz, Constructive Approximation, Springer-Verlag, Berlin Heidelberg, 1993. [25] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman and Hall/CRC Press, Boca Raton, FL, 2000. [26] R.G. Duran, A.L. Lombardi, Finite element approximation of convection diffusion problems using graded meshes, Appl. Numer. Math. 56 (2006) 1314–1325. [27] H.-G. Roos, T. Skalický, A comparison of the finite element method on Shishkin and Gartland-type meshes for convection–diffusion problems, CWI Q. 10 (1997) 277–300. [28] S. Franz, H.-G. Roos, Error estimation in a balanced norm for a convection–diffusion problem with two different boundary layers, Calcolo 51 (2014) 423–440. [29] R. Lin, M. Stynes, A balanced finite element method for singularly perturbed reaction–diffusion problems, SIAM J. Numer. Anal. 50 (5) (2012) 2729–2743. [30] H.-G. Roos, M. Schopf, Error estimation in energy norms: Is it necessary to fit the mesh to boundary layers, in: Numerical Analysis and Its Applications, in: Lecture Notes in Computer Science Volume, vol. 8236, 2013, pp. 95–109. [31] H.-G. Roos, M. Schopf, Convergence and stability in balanced norms of finite element methods on Shishkin meshes for reaction–diffusion problems, ZAMM Z. Angew. Math. Mech. (2014) 1–15.