- Email: [email protected]
Finite element approximation of convection–diffusion problems using an exponentially graded mesh
Computers and Mathematics with Applications 72 (2016) 1532–1540
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Finite element approximation of convection–diffusion problems using an exponentially graded mesh Christos Xenophontos a,∗ , Sebastian Franz b , Lars Ludwig c a
Department of Mathematics and Statistics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus
b
Technische Universität Dresden, Institut für Numerische Mathematik, 01062 Dresden, Germany
c
Institut für Wissenschaftliches Rechnen, Technische Universität Dresden, 01062 Dresden, Germany
article
info
Article history: Received 3 May 2016 Received in revised form 5 July 2016 Accepted 11 July 2016 Available online 9 August 2016 Keywords: Boundary layers Convection–diffusion Exponential mesh Uniform optimal convergence
abstract We present the analysis of an h version Finite Element Method for the approximation of the solution to convection–diffusion problems. The method uses piece-wise polynomials of degree p ≥ 1, defined on an exponentially graded mesh, optimally constructed for the approximation of exponential layers. We consider a model convection–diffusion problem, posed on the unit square and establish robust, optimal convergence rates in the energy and in the maximum norm. We also present the results of some numerical computations that illustrate our theoretical findings and compare the proposed method with others found in the literature. © 2016 Elsevier Ltd. All rights reserved.
1. Introduction Singularly perturbed problems have been studied for a number of decades (see [1,2] and the references therein) as they arise in a variety of applications. One family of such problems are convection-dominated equations, in which the following two difficulties arise: stability of discretizations and the presence of boundary layers in the solution. In the context of the Finite Element Method (FEM), the robust approximation of boundary layers requires either the use of the h version on non-uniform, layer adapted meshes (such as the Shishkin [3] or Bakhvalov [4] mesh), or the use of the high order p and hp versions on appropriately designed (variable) meshes [5]. In both cases, the a priori knowledge of the position of the layers is taken into account and mesh-degree combinations can be chosen for which uniform error estimates can be established (see, e.g., [6]). Among the various layer adapted meshes that have been proposed in the literature over the years, one finds the non-uniform, layer adapted, exponentially graded mesh, which was constructed by optimizing certain upper bounds on the error (see [7] for details). The finite element analysis on this mesh appears in [8] where robust, optimal rates of convergence were proved for one-dimensional reaction–diffusion and convection–diffusion problems. The purpose of this article is to extend the results of [8] to a two-dimensional convection–diffusion problem posed on a square, under sufficient regularity assumptions (which exclude the presence of corner singularities in the solution). The proposed method not only yields the optimal rate of convergence, but it does so without the assumption ε ≤ CN −1 , as is done in most other methods; here ε is the singular perturbation parameter, N is the number of meshpoints and C is a positive constant independent of ε and N. In addition, in numerical tests it outperforms the equally optimal Bakhvalov–Shishkin mesh (see [8] for one-dimensional results and Section 5 ahead).
∗
Corresponding author. E-mail address: [email protected] (C. Xenophontos).
http://dx.doi.org/10.1016/j.camwa.2016.07.008 0898-1221/© 2016 Elsevier Ltd. All rights reserved.
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
1533
The rest of the paper is organized as follows: in Section 2 we present the model problem and its regularity. The discretization using the exponentially graded mesh is presented in Section 3 and in Section 4 we present our main result of uniform, optimal convergence. Section 5 shows the results of some numerical computations that illustrate the theoretical findings and Section 6 gives our conclusions. With Ω ⊂ R2 a bounded open set with Lipschitz boundary ∂ Ω and measure |Ω |, we will denote by C k (Ω ) the space of continuous functions on Ω with continuous derivatives up to order k and by C k,α (Ω ) we will denote the space of functions whose derivatives up to order k are Hölder continuous with exponent α . We will use the usual Sobolev spaces W k,m (Ω ) of functions on Ω with 0, 1, 2, . . . , k generalized derivatives in Lm (Ω ), equipped with the norm and seminorm ∥·∥k,m,Ω and |·|k,m,Ω , respectively. When m = 2, we will write H k (Ω ) instead of W k,2 (Ω ), and for the norm and seminorm, we will write ∥·∥k,Ω and |·|k,Ω , respectively. We will also use the space H01 (Ω ) = u ∈ H 1 (Ω ) : u|∂ Ω = 0 .
The norm of the space L∞ (Ω ) of essentially bounded functions is denoted by ∥ · ∥∞,Ω . Finally, the notation ‘‘a . b’’ means ‘‘a ≤ Cb’’ with C being a generic positive constant, independent of any discretization or singular perturbation parameters. 2. The model problem and its regularity We consider the following convection–diffusion boundary value problem (BVP): Find u ∈ C 2 (Ω ) ∩ C Ω such that
−ε 1u + bT ∇ u + cu = f in Ω = (0, 1)2 , u = 0 on ∂ Ω ,
(1a) (1b)
where bT = [b1 (x, y), b2 (x, y)] ≥ [β1 , β2 ] > 0
on Ω ,
ε ∈ (0, 1] is a small parameter, b1 , b2 , c , f are sufficiently smooth functions and β1 , β2 are constants. We will assume, as usual, that c−
1 2
div b > 0,
c ≥ 0 on Ω
and our problem will have a unique (weak) solution u ∈ H01 (Ω ) for all f ∈ L2 (Ω ). The analysis will be performed under the assumption that the solution does not contain any corner singularities. This will happen if certain compatibility conditions are imposed on the data (see, e.g., [9]). Under such an assumption, the solution u will admit a decomposition (see, e.g., [6]) into a smooth part and a boundary layer part. Bounds on the derivatives of each part are available only in certain cases (see, e.g., [10]). To avoid corner singularities and focus on the boundary layers, we will make the following assumption: Assumption 1. The BVP (1) has a classical solution u ∈ C q+1,α Ω , q ≥ 2 for some α ∈ (0, 1), which can be decomposed as
u = S + E1 + E2 + E12 ,
(A1)
where for all (x, y) ∈ Ω there holds
i+j ∂ S . 1, ( x , y ) ∂ xi ∂ yj i+j ∂ E1 −i −β1 (1−x)/ε , ∂ xi ∂ yj (x, y) . ε e i+j ∂ E2 . ε −j e−β2 (1−y)/ε ( x , y ) ∂ xi ∂ yj
(A2)
(A3)
(A4)
and
i+j ∂ E12 −(i+j) −(β1 (1−x)+β2 (1−y))/ε e , ∂ xi ∂ yj (x, y) . ε
(A5)
for 0 ≤ i + j ≤ q. The above was proven for the case q = 2, Ω = (0, 1)2 in [11] and the precise compatibility of the data was given. See also [10] for a similar decomposition that includes corner singularities. Here, we assume that our problem has a unique solution that satisfies Assumption 1 for q ≥ 2.
1534
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
The variational formulation of (1) reads: Find u ∈ H01 (Ω ) such that B(u, v) = F (v)
∀ v ∈ H01 (Ω ) ,
(2)
where, with ⟨·, ·⟩Ω the usual L2 (Ω ) inner-product, B(u, v) = ε ⟨∇ u, ∇v⟩Ω + bT ∇ u + cu, v
F (v) = ⟨f , v⟩Ω .
Ω
,
Associated with the above problem is the energy norm
1/2 ∥v∥E := ε ∥∇v∥20,Ω + ∥v∥20,Ω , v ∈ H 1 (Ω ) . The bilinear form B(·, ·) is coercive with respect to the energy norm, i.e. B(v, v) & ∥v∥2E
(3)
∀ v ∈ H01 (Ω ) .
(4)
However, it is not uniformly bounded, with respect to ε , and as a result the finite element analysis requires non-standard techniques (see the proof of Theorem 9 ahead). The discrete version of (2) reads: find uh ∈ Vh ⊂ H01 (Ω ) such that B(uh , v) = F (v)
∀ v ∈ Vh ⊂ H01 (Ω ) ,
(5)
where Vh is a finite dimensional subspace of H01 (Ω ) that will be defined shortly. There holds B(u − uh , v) = 0 ∀ v ∈ Vh ⊂ H01 (Ω ) .
(6)
3. Discretization using an exponential mesh We now describe the design of the mesh, beginning with the construction in one-dimension; the mesh on Ω will then be obtained in a tensor product way. We point out that in view of (A3)–(A5), the ‘challenge’ lies in approximating the onedimensional boundary layer function e−β(1−x)/ε ,
x ∈ [0, 1], ε ∈ (0, 1], β ∈ R+ .
3.1. Description of the mesh in one-dimension
N
Let N > 2 be even and divide the interval I := (0, 1) into N subintervals Ij , using the nodal points xj j=0 and set hj = Ij = xj − xj−1 , j = 1, . . . , N. Let Pp (I ) denote the space of polynomials of degree ≤ p on I ⊂ R. Since N is a multiple of 2, we split the interval I into
[0, σ ],
[σ , 1]
(σ will be defined shortly) and on [0, σ ] we choose an equidistant mesh with N /2 + 1 elements. For the other subinterval the mesh will be exponentially graded with N /2 − 1 elements, using a mesh generating function φ that satisfies φ(0) = 0. First, we define the points x∗k =
ε (p + 1)φ(k/N ), β
k = 0, . . . , N /2 − 1,
where
φ(t ) = − ln 1 − 2Cp,ε t ,
t ∈ [0, 1/2 − 1/N ],
with
Cp,ε = 1 − exp −
β (p + 1)ε
∈ R+ .
(7)
N
We then set σ = 1 − x∗N /2−1 and the mesh points xj j=0 on I are given by (cf. [8])
∗ xN /2−1 2j , xj = 1+N ∗ 1 − x N −j ,
j = 0, . . . , N /2 + 1
(8)
j = N /2 + 1, . . . , N .
An example is shown in Fig. 1. It was shown in [8] that, under the assumption βε (p + 1) ln(N − 2) < 1, which merely requires that ε is not large, there holds hj .
ε (p + 1)N −1 , β
j = N /2 + 1, . . . , N
(9)
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
1535
Fig. 1. Example of the exponential mesh in one-dimension.
Fig. 2. Left: Subdivision of Ω . Right: Example of the exponential mesh on Ω .
and e
−β x∗N /2−1 /ε
. N −(p+1) .
(10)
3.2. Two-dimensional mesh and interpolation results
N
We follow the one-dimensional construction described above and define the mesh points {xi }Ni=0 , yj j=0 as in (8), using β = β1 for the points in the x-direction and β = β2 for the points in the y-direction. We set h(i x) = xi − xi−1 , h(j y) = yj − yj−1 , i, j = 1, . . . , N and we define the two-dimensional exponential mesh ∆ =
xi , y j
N
i,j=0
on Ω , using these nodal
points. We denote each element in the mesh by Ωk , k = 1, . . . , (N + 1) . An example of the mesh is shown in Fig. 2. We then define the finite dimensional subspace Vh that appears in (5) as 2
Vh = u ∈ H01 (Ω ) : u|Ωk ∈ Qp (Ωk ) , k = 1, . . . , (N + 1)2 ,
where
Qp (Ω ) = span xi yj : x, y ∈ Ω , i, j = 0, . . . , p .
Eqs. (9), (10) hold in ‘both directions’, namely: (x)
(y)
max{hi , hj } . (p + 1)ε N −1 ,
max e
−β1 1−x∗N /2−1 /ε
,e
−β2
j = N /2 + 1, . . . , N ,
1−y∗ /ε N /2−1
. N −(p+1) .
(11) (12)
Next, we note that the mesh ‘divides’ the domain Ω into four regions as follows:
Ω11 = 0, 1 − x∗N /2−1 × 0, 1 − y∗N /2−1 , Ω12 = 0, 1 − x∗N /2−1 × 1 − y∗N /2−1 , 1 ,
Ω21 = 1 − x∗N /2−1 , 1 × 0, 1 − y∗N /2−1 Ω22 = 1 − x∗N /2−1 , 1 × 1 − y∗N /2−1 , 1 .
(See Fig. 2.) By (12) and (A3)–(A5), we have that the boundary layer components of the solution satisfy max{|E1 (x, y)|, |E2 (x, y)|, |E12 (x, y)|} . N −(p+1)
∀ (x, y) ∈ Ω11 .
(13)
We next define the two-dimensional interpolant we will use, usually referred to as a ‘nodal-edge-element’ interpolant (see, e.g., [12]). To this end, let K = (−1, 1)2 be the reference element with vertices ai and edges ei , i = 1, . . . , 4, and let Pp (a, b) denote the space of polynomials on (a, b) of degree ≤ q. Definition 2. Let v(·, ·) ∈ C ( K ). The approximation operator I : C ( K ) → Qp ( K ) is defined as follows:
I v ( ai ) = v ( ai ) , i = 1, . . . , 4, I v q= v q ∀ q ∈ Pp−2 ( ei ), i = 1, . . . , 4, ei ei I v q= v q ∀ q ∈ Qp−2 ( K ), i = 1, . . . , 4. K
K
1536
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
It was shown in [13] that the above approximation operator I is uniquely defined. Now, using the transformation (y) (x) x = xK + hK ξ /2, y = yK + hK η/2 we map the reference element K to any element K ∈ ∆ and correspondingly we obtain the approximation operator IK : C (K ) → Qp (K ). The global, continuous approximation operator IN : C (Ω ) → Vh is defined element-wise by
(IN v)|K = IK ( v|K ) ∀ K ∈ ∆. Moreover, the following stability bound holds [13]:
∥IN v∥∞,K . ∥v∥∞,K .
(14)
We close this section with some results from [13] pertaining to the interpolation operator of Definition 2. Lemma 3. Let p ≥ 1 be an integer, let K ∈ ∆ be an element in the mesh, with |K | = hx hy , and let v ∈ H p+1 (K ). Then, with IN being the interpolation operator of Definition 2, there holds
p+1 ∂ u i j ∂ (v − IN v) . , hx hy i+1 j ∂x ∂ x ∂ y 0,K 0 ,K i+j=p p+1 ∂ u i j ∂ (v − IN v) . hx hy i j +1 . ∂y ∂ x ∂ y 0,K 0 ,K i+j=p Proof. This was shown in [13, Lemma 6].
Lemma 4. Let q ∈ [1, ∞], let K ∈ ∆ be an element in the mesh, with |K | = hx hy , and let v ∈ W p+1,q (K ). Then, with IN being the interpolation operator of Definition 2, there holds
∥v − IN v∥Lq (K ) .
i+j=p+1
p+1 ∂ u ∂ xi ∂ yj q
hix hjy
L (K )
.
Proof. This was established in [13, Lemma 7], but the proof is based on [14].
4. Error estimates I Recall the decomposition u = S + E1 + E2 + E12 and let uI = S I + E1I + E2I + E12 be the interpolant of u based on Definition 2, with the obvious notation. We have the following lemma.
Lemma 5. Under Assumption 1, we have
I max E1I ∞,Ω ∪Ω , E2I ∞,Ω ∪Ω , E12 ∞,Ω \Ω22 11 12 11 21
. N −(p+1) .
Proof. For Ω11 , the desired bounds follow immediately from (13). For Ω12 , we use (14) and (A3)–(A5) to get for, e.g., E1I
I E
1 ∞,Ω12
. ∥E1 ∥∞,Ω12 . e−β1 (1−x)/ε ∞,Ω
12
. e−β1 (1−xN /2−1 )/ε . N −(p+1) ,
by (10). The other two estimates are shown analogously.
The next result gives bounds on the interpolation error. Theorem 6. Let u satisfy Assumption 1 and let uI be its interpolant as in Definition 2. Then
u − uI
∞,Ω
. N −(p+1) .
I . By (A2) and Lemma 4, we have Proof. We again use the decompositions u = S + E1 + E2 + E12 , uI = S I + E1I + E2I + E12
S − S I
∞,Ω
. N −(p+1) .
The estimation of the remaining terms is quite adapted meshes and it follows the arguments found in, standard on layer e.g. [13]; we will only highlight the steps for E1 − E1I (x, y). For (x, y) ∈ Ω11 ∪ Ω12 , we use (A3), (13) and (14), to get
E1 − E I
1 ∞,Ω11 ∪Ω12
. ∥E1 ∥∞,Ω11 ∪Ω12 . N −(p+1) . (x)
(y)
For an element Ωk ∈ Ω21 there holds hk . βε (p + 1)N −1 , hk . N −1 , so that Lemma 4, (A3) and (11) give E1 − E1I ∞,Ω . 1 21 N −(p+1) . For Ωk ∈ Ω22 we use Lemma 4, (A3) and the fact that x ∈ [1 − x∗N /2−1 , 1], y ∈ [1 − y∗N /2−1 , 1] along with (11), to
get E1 − E1I ∞,Ω . N −(p+1) . Combining these estimates we have E1 − E1I ∞,Ω . N −(p+1) .
k
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
1537
The following lemma gives bounds on the interpolation error in the derivatives. Lemma 7. Let u satisfy Assumption 1 and let uI be its interpolant as in Definition 2. Then
u − uI
1,Ω
. ε −1/2 N −p .
Proof. Again the arguments are standard and follow [13], so we will give the main steps for only ∂∂x u − uI 0,Ω . With I we have from Lemma 3 and (A2), u = S + E1 + E2 + E12 , uI = S I + E1I + E2I + E12
∂ I ∂x S − S
. N −p .
0,Ω
The remaining terms will be dealt with one at a time,starting with E1 . For Ωk ⊂ Ω21 ∪ Ω22 , there holds |Ωk | = hx hy , hx . ε (p + 1)N −1 , hy . N −1 and Lemma 3 and (A3) give ∂∂x E1I − E1 0,Ω . ε 1/2 N −(p+1) , so that summing over all elements β 1 Ωk ⊂ Ω21 ∪ Ω22 we get ∂∂x E1I − E1 0,Ω
k
21 ∪Ω22
. ε 1/2 N −p . For Ωk ⊂ Ω12 ∪ Ω11 , with |Ωk | = hx hy . N −2 , we have by the
triangle inequality, (14), an inverse inequality, (A3) and (13)
∂ I E − E 1 ∂x 1
0,Ωk
≤ ε −1/2 N −(p+1) ,
so that summing over all elements Ωk ⊂ Ω12 ∪ Ω11 gives the desired result. We now look at E2 : following the same steps as above, we obtain ∂∂x E2I − E2 0,Ω ∪Ω . ε 1/2 N −p . For Ωk ⊂ Ω21 ∪ Ω11 , we have by Lemma 3 (in Ω21 ) and standard 12
interpolation estimates (in Ω11 )
∂ I E − E 2 ∂x 2
.N
−1
0,Ωk
22
2 ∂ 2 E2 ∂ E2 + . ε −1/2 N −(p+1) , ∂ x2 ∂ x∂ y 0,Ωk 0,Ωk
where also used. Summing over all Ωk ⊂ Ω12 ∪ Ω11 and combining the above estimates we get (A4) and (11) were ∂ E I − E2 . ε −1/2 N −p . We finally consider E12 . In Ω12 ∪ Ω11 we use the same technique as for E2 and in Ω22 2 ∂x 0,Ω we use theanisotropic estimate of Lemma 3. In order to estimate the E12 term in Ω I 21 , weuse the triangle inequality and (14) to get ∂∂x E12 − E12 0,Ω . ε −1/2 N −p . Combining all of the above gives ∂∂x u − uI 0,Ω . ε −1/2 N −p . 21
Theorem 6 and Lemma 7 immediately yield the following. Corollary 8. Let u be the solution of (2) and let uI be its interpolant as in Definition 2. Then
u − uI . N − p . E We are now in a position to present our main result. Theorem 9. Let the solution u of (2) satisfy Assumption 1 and let uh , the solution of (5), be its finite element approximation based on the exponential mesh. Then
∥ u − uh ∥ E . N − p . Proof. As already mentioned in Section 2, the analysis does not follow from standard techniques. Following [15], we have from (4) and (6),
I u − uh 2 . B uI − uh , uI − uh = B uI − u, uI − uh + B u − uh , uI − uh E = B uI − u, uI − uh , where uI is the interpolant of u. Now, B uI − u, uI − uh . uI − uE uI − uh E + uI − u0,Ω
T I b ∇ u − uh 0,Ω
11 + uI − u∞,Ω \Ω bT ∇ uI − uh L1 (Ω \Ω 11
11
An inverse estimate gives
T I b ∇ u − uh . N uI − uh 0,Ω , 0 ,Ω 11
11
while the Cauchy–Schwarz inequality yields
T I b ∇ u − uh 1 L ( Ω \Ω
11 )
. ε 1/2 ∇ uI − uh 0,Ω \Ω . 11
11
. )
1538
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
Fig. 3. Energy norm convergence of the method using the eXp mesh and a given solution.
Combining these estimates, we have
I u − uh . uI − u + N uI − u + uI − u∞,Ω \Ω . N −p , 0 ,Ω E E 11
11
by Theorem 6 and Lemma 7. The proof is completed by noticing that
∥ u − uh ∥ E ≤ u − uI E + uI − uh E and using Corollary 8.
5. Numerical results We illustrate the theoretical findings by considering the BVP (1) with bT = [1, 1] , c = 1 and two different choices for the right hand side f . For the first one we calculate f such that the exact solution is given by
u(x, y) =
sin(π x/2) −
e−(1−x)/ε − e−1/ε 1 − e−1/ε
y−
e−2(1−y)/ε − e−2/ε
1 − e−2/ε
.
We use piecewise polynomials of degree p = 1, 2, 3, defined on the exponentially graded mesh (cf. Fig. 2) with N points in each direction and we compute the error in the energy norm between the exact and numerical solution. The performance of the method, for various values of ε , is shown in Fig. 3, in which the percentage relative error is plotted versus N, in a log–log scale. The straight lines have slope −p, just as Theorem 9 predicts and moreover they coincide for all ε , verifying the robustness of the method. We also show, in Fig. 4, a comparison between the proposed exponential (eXp) mesh and two other meshes from the literature: the Bakhvalov–Shishkin (B–S) mesh (see, e.g., [6] and the references therein) and another optimally constructed, layer adapted mesh with robust convergence rate from [16], which we label RTU. We compare the three methods for ε = 10−6 and p = 1, 2, 3—other values of ε gave the same results. As can be seen, all methods yield the optimal rate of convergence, but the eXp mesh performs slightly better, perhaps due to the constant being smaller in the error estimate. For our second example we set the right hand side f = 1. In this case the exact solution is unknown and we use a numerically computed reference solution uref as a substitute for the exact solution in computing the errors. For computing uref , a polynomial degree of 4 and an exponentially refined mesh with N = 512 is used. Furthermore, the exact solution of this problem does not fulfil Assumption 1, as the basic compatibility conditions f (0, 0) = f (1, 0) = f (0, 1) = f (1, 1) = 0, given in [10], are not fulfilled. Therefore the exact solution has corner singularities and reduced regularity. Fig. 5 shows the convergence behaviour for several values of ε and polynomial degrees p = 1, 2, 3. Clearly the uniformity in ε is observable for p = 1 and p = 2. The optimal convergence rate of order one is still attained for p = 1, while for p = 2 we already see a reduction from two to 1.8. For p = 3 the situation is more drastic as now the convergence rate decreases for ε becoming smaller. For ε = 10−8 we observe only a rate of 1.33. 6. Conclusions We presented the finite element analysis for a two-dimensional exponentially graded mesh, appropriate for singularly perturbed problems with exponential layers. By considering a model convection–diffusion problem, we proved robust,
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
1539
Fig. 4. Energy norm convergence: comparison of the eXp, B–S and RTU meshes.
Fig. 5. Energy norm convergence of the method using the eXp mesh and a f ≡ 1.
optimal convergence in the energy and maximum norms. Based on numerical experiments, it appears that the proposed mesh outperforms the commonly used layer adapted Bakhvalov–Shishkin mesh, hence it presents itself as an attractive alternative. References [1] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [2] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, second ed., in: Springer Series in Computational Mathematics, vol. 24, Springer-Verlag, Berlin, 2008, Convection–diffusion-reaction and flow problems. [3] G.I. Shishkin, Grid approximation of singularly perturbed boundary value problems with a regular boundary layer, Sov. J. Numer. Anal. Math. Model. 4 (1989) 397–417. [4] N.S. Bakhvalov, Towards optimization of methods for solving boundary value problems in the presence of boundary layers, Zh. Vychisl. Mat. Mat. Fiz. 9 (1969) 841–859 (in Russian). [5] C. Schwab, M. Suri, The p and hp versions of the finite element method for problems with boundary layers, Math. Comp. 65 (1996) 1403–1429. [6] T. Linß, Layer-Adapted Meshes for Reaction-Convection-Diffusion Problems, in: Lecture Notes in Mathematics, vol. 1985, Springer-Verlag, 2010. [7] C. Xenophontos, Optimal mesh design for the finite element approximation of reaction–diffusion problems, Internat. J. Numer. Methods Engrg. 53 (2002) 929–943. [8] P. Cnstantinou, C. Xenophontos, Finite element analysis of an exponentially graded mesh for singularly perturbed problems, Comput. Methods Appl. Math. 15 (2) (2015) 135–143. [9] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman Publishing, 1985. [10] R.B. Kellogg, M. Stynes, Corner singularities and boundary layers in a simple convection–diffusion problem, J. Differential Equations 213 (2005) 81–120. [11] T. Linß, Analysis of a Galerkin finite element method on a Bakhvalov–Shishkin mesh for a linear convection–diffusion problem, IMA J. Numer. Anal. 20 (2000) 621–632.
1540
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
[12] Q. Lin, N. Yan, A. Zhou, A rectangle test for interpolated finite elements, in: Proc. Syst. Sci. Eng., Culture Publishing Co., Great Wall (HK), 1991, pp. 217–229. [13] M. Stynes, L. Tobiska, Using rectangular Qp elements in the SDFEM for a convection–diffusion problem with a boundary layer, Appl. Numer. Math. 58 (2008) 1789–1802. [14] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, in: Advances in Numerical Mathematics, B. G. Teubner, Stuttgart, 1999. [15] E. O’Riordan, M. Stynes, A uniformly convergent Galerkin method on a Shishkin mesh for a convection–diffusion problem, J. Math. Anal. Appl. 214 (1997) 36–54. [16] H.-G. Roos, Lj. Teofanov, Z. Uzelac, Graded meshes for higher order FEM, J. Comput. Math. 33 (2015) 1–16.
Contents lists available at ScienceDirect
Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa
Finite element approximation of convection–diffusion problems using an exponentially graded mesh Christos Xenophontos a,∗ , Sebastian Franz b , Lars Ludwig c a
Department of Mathematics and Statistics, University of Cyprus, P.O. Box 20537, 1678 Nicosia, Cyprus
b
Technische Universität Dresden, Institut für Numerische Mathematik, 01062 Dresden, Germany
c
Institut für Wissenschaftliches Rechnen, Technische Universität Dresden, 01062 Dresden, Germany
article
info
Article history: Received 3 May 2016 Received in revised form 5 July 2016 Accepted 11 July 2016 Available online 9 August 2016 Keywords: Boundary layers Convection–diffusion Exponential mesh Uniform optimal convergence
abstract We present the analysis of an h version Finite Element Method for the approximation of the solution to convection–diffusion problems. The method uses piece-wise polynomials of degree p ≥ 1, defined on an exponentially graded mesh, optimally constructed for the approximation of exponential layers. We consider a model convection–diffusion problem, posed on the unit square and establish robust, optimal convergence rates in the energy and in the maximum norm. We also present the results of some numerical computations that illustrate our theoretical findings and compare the proposed method with others found in the literature. © 2016 Elsevier Ltd. All rights reserved.
1. Introduction Singularly perturbed problems have been studied for a number of decades (see [1,2] and the references therein) as they arise in a variety of applications. One family of such problems are convection-dominated equations, in which the following two difficulties arise: stability of discretizations and the presence of boundary layers in the solution. In the context of the Finite Element Method (FEM), the robust approximation of boundary layers requires either the use of the h version on non-uniform, layer adapted meshes (such as the Shishkin [3] or Bakhvalov [4] mesh), or the use of the high order p and hp versions on appropriately designed (variable) meshes [5]. In both cases, the a priori knowledge of the position of the layers is taken into account and mesh-degree combinations can be chosen for which uniform error estimates can be established (see, e.g., [6]). Among the various layer adapted meshes that have been proposed in the literature over the years, one finds the non-uniform, layer adapted, exponentially graded mesh, which was constructed by optimizing certain upper bounds on the error (see [7] for details). The finite element analysis on this mesh appears in [8] where robust, optimal rates of convergence were proved for one-dimensional reaction–diffusion and convection–diffusion problems. The purpose of this article is to extend the results of [8] to a two-dimensional convection–diffusion problem posed on a square, under sufficient regularity assumptions (which exclude the presence of corner singularities in the solution). The proposed method not only yields the optimal rate of convergence, but it does so without the assumption ε ≤ CN −1 , as is done in most other methods; here ε is the singular perturbation parameter, N is the number of meshpoints and C is a positive constant independent of ε and N. In addition, in numerical tests it outperforms the equally optimal Bakhvalov–Shishkin mesh (see [8] for one-dimensional results and Section 5 ahead).
∗
Corresponding author. E-mail address: [email protected] (C. Xenophontos).
http://dx.doi.org/10.1016/j.camwa.2016.07.008 0898-1221/© 2016 Elsevier Ltd. All rights reserved.
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
1533
The rest of the paper is organized as follows: in Section 2 we present the model problem and its regularity. The discretization using the exponentially graded mesh is presented in Section 3 and in Section 4 we present our main result of uniform, optimal convergence. Section 5 shows the results of some numerical computations that illustrate the theoretical findings and Section 6 gives our conclusions. With Ω ⊂ R2 a bounded open set with Lipschitz boundary ∂ Ω and measure |Ω |, we will denote by C k (Ω ) the space of continuous functions on Ω with continuous derivatives up to order k and by C k,α (Ω ) we will denote the space of functions whose derivatives up to order k are Hölder continuous with exponent α . We will use the usual Sobolev spaces W k,m (Ω ) of functions on Ω with 0, 1, 2, . . . , k generalized derivatives in Lm (Ω ), equipped with the norm and seminorm ∥·∥k,m,Ω and |·|k,m,Ω , respectively. When m = 2, we will write H k (Ω ) instead of W k,2 (Ω ), and for the norm and seminorm, we will write ∥·∥k,Ω and |·|k,Ω , respectively. We will also use the space H01 (Ω ) = u ∈ H 1 (Ω ) : u|∂ Ω = 0 .
The norm of the space L∞ (Ω ) of essentially bounded functions is denoted by ∥ · ∥∞,Ω . Finally, the notation ‘‘a . b’’ means ‘‘a ≤ Cb’’ with C being a generic positive constant, independent of any discretization or singular perturbation parameters. 2. The model problem and its regularity We consider the following convection–diffusion boundary value problem (BVP): Find u ∈ C 2 (Ω ) ∩ C Ω such that
−ε 1u + bT ∇ u + cu = f in Ω = (0, 1)2 , u = 0 on ∂ Ω ,
(1a) (1b)
where bT = [b1 (x, y), b2 (x, y)] ≥ [β1 , β2 ] > 0
on Ω ,
ε ∈ (0, 1] is a small parameter, b1 , b2 , c , f are sufficiently smooth functions and β1 , β2 are constants. We will assume, as usual, that c−
1 2
div b > 0,
c ≥ 0 on Ω
and our problem will have a unique (weak) solution u ∈ H01 (Ω ) for all f ∈ L2 (Ω ). The analysis will be performed under the assumption that the solution does not contain any corner singularities. This will happen if certain compatibility conditions are imposed on the data (see, e.g., [9]). Under such an assumption, the solution u will admit a decomposition (see, e.g., [6]) into a smooth part and a boundary layer part. Bounds on the derivatives of each part are available only in certain cases (see, e.g., [10]). To avoid corner singularities and focus on the boundary layers, we will make the following assumption: Assumption 1. The BVP (1) has a classical solution u ∈ C q+1,α Ω , q ≥ 2 for some α ∈ (0, 1), which can be decomposed as
u = S + E1 + E2 + E12 ,
(A1)
where for all (x, y) ∈ Ω there holds
i+j ∂ S . 1, ( x , y ) ∂ xi ∂ yj i+j ∂ E1 −i −β1 (1−x)/ε , ∂ xi ∂ yj (x, y) . ε e i+j ∂ E2 . ε −j e−β2 (1−y)/ε ( x , y ) ∂ xi ∂ yj
(A2)
(A3)
(A4)
and
i+j ∂ E12 −(i+j) −(β1 (1−x)+β2 (1−y))/ε e , ∂ xi ∂ yj (x, y) . ε
(A5)
for 0 ≤ i + j ≤ q. The above was proven for the case q = 2, Ω = (0, 1)2 in [11] and the precise compatibility of the data was given. See also [10] for a similar decomposition that includes corner singularities. Here, we assume that our problem has a unique solution that satisfies Assumption 1 for q ≥ 2.
1534
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
The variational formulation of (1) reads: Find u ∈ H01 (Ω ) such that B(u, v) = F (v)
∀ v ∈ H01 (Ω ) ,
(2)
where, with ⟨·, ·⟩Ω the usual L2 (Ω ) inner-product, B(u, v) = ε ⟨∇ u, ∇v⟩Ω + bT ∇ u + cu, v
F (v) = ⟨f , v⟩Ω .
Ω
,
Associated with the above problem is the energy norm
1/2 ∥v∥E := ε ∥∇v∥20,Ω + ∥v∥20,Ω , v ∈ H 1 (Ω ) . The bilinear form B(·, ·) is coercive with respect to the energy norm, i.e. B(v, v) & ∥v∥2E
(3)
∀ v ∈ H01 (Ω ) .
(4)
However, it is not uniformly bounded, with respect to ε , and as a result the finite element analysis requires non-standard techniques (see the proof of Theorem 9 ahead). The discrete version of (2) reads: find uh ∈ Vh ⊂ H01 (Ω ) such that B(uh , v) = F (v)
∀ v ∈ Vh ⊂ H01 (Ω ) ,
(5)
where Vh is a finite dimensional subspace of H01 (Ω ) that will be defined shortly. There holds B(u − uh , v) = 0 ∀ v ∈ Vh ⊂ H01 (Ω ) .
(6)
3. Discretization using an exponential mesh We now describe the design of the mesh, beginning with the construction in one-dimension; the mesh on Ω will then be obtained in a tensor product way. We point out that in view of (A3)–(A5), the ‘challenge’ lies in approximating the onedimensional boundary layer function e−β(1−x)/ε ,
x ∈ [0, 1], ε ∈ (0, 1], β ∈ R+ .
3.1. Description of the mesh in one-dimension
N
Let N > 2 be even and divide the interval I := (0, 1) into N subintervals Ij , using the nodal points xj j=0 and set hj = Ij = xj − xj−1 , j = 1, . . . , N. Let Pp (I ) denote the space of polynomials of degree ≤ p on I ⊂ R. Since N is a multiple of 2, we split the interval I into
[0, σ ],
[σ , 1]
(σ will be defined shortly) and on [0, σ ] we choose an equidistant mesh with N /2 + 1 elements. For the other subinterval the mesh will be exponentially graded with N /2 − 1 elements, using a mesh generating function φ that satisfies φ(0) = 0. First, we define the points x∗k =
ε (p + 1)φ(k/N ), β
k = 0, . . . , N /2 − 1,
where
φ(t ) = − ln 1 − 2Cp,ε t ,
t ∈ [0, 1/2 − 1/N ],
with
Cp,ε = 1 − exp −
β (p + 1)ε
∈ R+ .
(7)
N
We then set σ = 1 − x∗N /2−1 and the mesh points xj j=0 on I are given by (cf. [8])
∗ xN /2−1 2j , xj = 1+N ∗ 1 − x N −j ,
j = 0, . . . , N /2 + 1
(8)
j = N /2 + 1, . . . , N .
An example is shown in Fig. 1. It was shown in [8] that, under the assumption βε (p + 1) ln(N − 2) < 1, which merely requires that ε is not large, there holds hj .
ε (p + 1)N −1 , β
j = N /2 + 1, . . . , N
(9)
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
1535
Fig. 1. Example of the exponential mesh in one-dimension.
Fig. 2. Left: Subdivision of Ω . Right: Example of the exponential mesh on Ω .
and e
−β x∗N /2−1 /ε
. N −(p+1) .
(10)
3.2. Two-dimensional mesh and interpolation results
N
We follow the one-dimensional construction described above and define the mesh points {xi }Ni=0 , yj j=0 as in (8), using β = β1 for the points in the x-direction and β = β2 for the points in the y-direction. We set h(i x) = xi − xi−1 , h(j y) = yj − yj−1 , i, j = 1, . . . , N and we define the two-dimensional exponential mesh ∆ =
xi , y j
N
i,j=0
on Ω , using these nodal
points. We denote each element in the mesh by Ωk , k = 1, . . . , (N + 1) . An example of the mesh is shown in Fig. 2. We then define the finite dimensional subspace Vh that appears in (5) as 2
Vh = u ∈ H01 (Ω ) : u|Ωk ∈ Qp (Ωk ) , k = 1, . . . , (N + 1)2 ,
where
Qp (Ω ) = span xi yj : x, y ∈ Ω , i, j = 0, . . . , p .
Eqs. (9), (10) hold in ‘both directions’, namely: (x)
(y)
max{hi , hj } . (p + 1)ε N −1 ,
max e
−β1 1−x∗N /2−1 /ε
,e
−β2
j = N /2 + 1, . . . , N ,
1−y∗ /ε N /2−1
. N −(p+1) .
(11) (12)
Next, we note that the mesh ‘divides’ the domain Ω into four regions as follows:
Ω11 = 0, 1 − x∗N /2−1 × 0, 1 − y∗N /2−1 , Ω12 = 0, 1 − x∗N /2−1 × 1 − y∗N /2−1 , 1 ,
Ω21 = 1 − x∗N /2−1 , 1 × 0, 1 − y∗N /2−1 Ω22 = 1 − x∗N /2−1 , 1 × 1 − y∗N /2−1 , 1 .
(See Fig. 2.) By (12) and (A3)–(A5), we have that the boundary layer components of the solution satisfy max{|E1 (x, y)|, |E2 (x, y)|, |E12 (x, y)|} . N −(p+1)
∀ (x, y) ∈ Ω11 .
(13)
We next define the two-dimensional interpolant we will use, usually referred to as a ‘nodal-edge-element’ interpolant (see, e.g., [12]). To this end, let K = (−1, 1)2 be the reference element with vertices ai and edges ei , i = 1, . . . , 4, and let Pp (a, b) denote the space of polynomials on (a, b) of degree ≤ q. Definition 2. Let v(·, ·) ∈ C ( K ). The approximation operator I : C ( K ) → Qp ( K ) is defined as follows:
I v ( ai ) = v ( ai ) , i = 1, . . . , 4, I v q= v q ∀ q ∈ Pp−2 ( ei ), i = 1, . . . , 4, ei ei I v q= v q ∀ q ∈ Qp−2 ( K ), i = 1, . . . , 4. K
K
1536
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
It was shown in [13] that the above approximation operator I is uniquely defined. Now, using the transformation (y) (x) x = xK + hK ξ /2, y = yK + hK η/2 we map the reference element K to any element K ∈ ∆ and correspondingly we obtain the approximation operator IK : C (K ) → Qp (K ). The global, continuous approximation operator IN : C (Ω ) → Vh is defined element-wise by
(IN v)|K = IK ( v|K ) ∀ K ∈ ∆. Moreover, the following stability bound holds [13]:
∥IN v∥∞,K . ∥v∥∞,K .
(14)
We close this section with some results from [13] pertaining to the interpolation operator of Definition 2. Lemma 3. Let p ≥ 1 be an integer, let K ∈ ∆ be an element in the mesh, with |K | = hx hy , and let v ∈ H p+1 (K ). Then, with IN being the interpolation operator of Definition 2, there holds
p+1 ∂ u i j ∂ (v − IN v) . , hx hy i+1 j ∂x ∂ x ∂ y 0,K 0 ,K i+j=p p+1 ∂ u i j ∂ (v − IN v) . hx hy i j +1 . ∂y ∂ x ∂ y 0,K 0 ,K i+j=p Proof. This was shown in [13, Lemma 6].
Lemma 4. Let q ∈ [1, ∞], let K ∈ ∆ be an element in the mesh, with |K | = hx hy , and let v ∈ W p+1,q (K ). Then, with IN being the interpolation operator of Definition 2, there holds
∥v − IN v∥Lq (K ) .
i+j=p+1
p+1 ∂ u ∂ xi ∂ yj q
hix hjy
L (K )
.
Proof. This was established in [13, Lemma 7], but the proof is based on [14].
4. Error estimates I Recall the decomposition u = S + E1 + E2 + E12 and let uI = S I + E1I + E2I + E12 be the interpolant of u based on Definition 2, with the obvious notation. We have the following lemma.
Lemma 5. Under Assumption 1, we have
I max E1I ∞,Ω ∪Ω , E2I ∞,Ω ∪Ω , E12 ∞,Ω \Ω22 11 12 11 21
. N −(p+1) .
Proof. For Ω11 , the desired bounds follow immediately from (13). For Ω12 , we use (14) and (A3)–(A5) to get for, e.g., E1I
I E
1 ∞,Ω12
. ∥E1 ∥∞,Ω12 . e−β1 (1−x)/ε ∞,Ω
12
. e−β1 (1−xN /2−1 )/ε . N −(p+1) ,
by (10). The other two estimates are shown analogously.
The next result gives bounds on the interpolation error. Theorem 6. Let u satisfy Assumption 1 and let uI be its interpolant as in Definition 2. Then
u − uI
∞,Ω
. N −(p+1) .
I . By (A2) and Lemma 4, we have Proof. We again use the decompositions u = S + E1 + E2 + E12 , uI = S I + E1I + E2I + E12
S − S I
∞,Ω
. N −(p+1) .
The estimation of the remaining terms is quite adapted meshes and it follows the arguments found in, standard on layer e.g. [13]; we will only highlight the steps for E1 − E1I (x, y). For (x, y) ∈ Ω11 ∪ Ω12 , we use (A3), (13) and (14), to get
E1 − E I
1 ∞,Ω11 ∪Ω12
. ∥E1 ∥∞,Ω11 ∪Ω12 . N −(p+1) . (x)
(y)
For an element Ωk ∈ Ω21 there holds hk . βε (p + 1)N −1 , hk . N −1 , so that Lemma 4, (A3) and (11) give E1 − E1I ∞,Ω . 1 21 N −(p+1) . For Ωk ∈ Ω22 we use Lemma 4, (A3) and the fact that x ∈ [1 − x∗N /2−1 , 1], y ∈ [1 − y∗N /2−1 , 1] along with (11), to
get E1 − E1I ∞,Ω . N −(p+1) . Combining these estimates we have E1 − E1I ∞,Ω . N −(p+1) .
k
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
1537
The following lemma gives bounds on the interpolation error in the derivatives. Lemma 7. Let u satisfy Assumption 1 and let uI be its interpolant as in Definition 2. Then
u − uI
1,Ω
. ε −1/2 N −p .
Proof. Again the arguments are standard and follow [13], so we will give the main steps for only ∂∂x u − uI 0,Ω . With I we have from Lemma 3 and (A2), u = S + E1 + E2 + E12 , uI = S I + E1I + E2I + E12
∂ I ∂x S − S
. N −p .
0,Ω
The remaining terms will be dealt with one at a time,starting with E1 . For Ωk ⊂ Ω21 ∪ Ω22 , there holds |Ωk | = hx hy , hx . ε (p + 1)N −1 , hy . N −1 and Lemma 3 and (A3) give ∂∂x E1I − E1 0,Ω . ε 1/2 N −(p+1) , so that summing over all elements β 1 Ωk ⊂ Ω21 ∪ Ω22 we get ∂∂x E1I − E1 0,Ω
k
21 ∪Ω22
. ε 1/2 N −p . For Ωk ⊂ Ω12 ∪ Ω11 , with |Ωk | = hx hy . N −2 , we have by the
triangle inequality, (14), an inverse inequality, (A3) and (13)
∂ I E − E 1 ∂x 1
0,Ωk
≤ ε −1/2 N −(p+1) ,
so that summing over all elements Ωk ⊂ Ω12 ∪ Ω11 gives the desired result. We now look at E2 : following the same steps as above, we obtain ∂∂x E2I − E2 0,Ω ∪Ω . ε 1/2 N −p . For Ωk ⊂ Ω21 ∪ Ω11 , we have by Lemma 3 (in Ω21 ) and standard 12
interpolation estimates (in Ω11 )
∂ I E − E 2 ∂x 2
.N
−1
0,Ωk
22
2 ∂ 2 E2 ∂ E2 + . ε −1/2 N −(p+1) , ∂ x2 ∂ x∂ y 0,Ωk 0,Ωk
where also used. Summing over all Ωk ⊂ Ω12 ∪ Ω11 and combining the above estimates we get (A4) and (11) were ∂ E I − E2 . ε −1/2 N −p . We finally consider E12 . In Ω12 ∪ Ω11 we use the same technique as for E2 and in Ω22 2 ∂x 0,Ω we use theanisotropic estimate of Lemma 3. In order to estimate the E12 term in Ω I 21 , weuse the triangle inequality and (14) to get ∂∂x E12 − E12 0,Ω . ε −1/2 N −p . Combining all of the above gives ∂∂x u − uI 0,Ω . ε −1/2 N −p . 21
Theorem 6 and Lemma 7 immediately yield the following. Corollary 8. Let u be the solution of (2) and let uI be its interpolant as in Definition 2. Then
u − uI . N − p . E We are now in a position to present our main result. Theorem 9. Let the solution u of (2) satisfy Assumption 1 and let uh , the solution of (5), be its finite element approximation based on the exponential mesh. Then
∥ u − uh ∥ E . N − p . Proof. As already mentioned in Section 2, the analysis does not follow from standard techniques. Following [15], we have from (4) and (6),
I u − uh 2 . B uI − uh , uI − uh = B uI − u, uI − uh + B u − uh , uI − uh E = B uI − u, uI − uh , where uI is the interpolant of u. Now, B uI − u, uI − uh . uI − uE uI − uh E + uI − u0,Ω
T I b ∇ u − uh 0,Ω
11 + uI − u∞,Ω \Ω bT ∇ uI − uh L1 (Ω \Ω 11
11
An inverse estimate gives
T I b ∇ u − uh . N uI − uh 0,Ω , 0 ,Ω 11
11
while the Cauchy–Schwarz inequality yields
T I b ∇ u − uh 1 L ( Ω \Ω
11 )
. ε 1/2 ∇ uI − uh 0,Ω \Ω . 11
11
. )
1538
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
Fig. 3. Energy norm convergence of the method using the eXp mesh and a given solution.
Combining these estimates, we have
I u − uh . uI − u + N uI − u + uI − u∞,Ω \Ω . N −p , 0 ,Ω E E 11
11
by Theorem 6 and Lemma 7. The proof is completed by noticing that
∥ u − uh ∥ E ≤ u − uI E + uI − uh E and using Corollary 8.
5. Numerical results We illustrate the theoretical findings by considering the BVP (1) with bT = [1, 1] , c = 1 and two different choices for the right hand side f . For the first one we calculate f such that the exact solution is given by
u(x, y) =
sin(π x/2) −
e−(1−x)/ε − e−1/ε 1 − e−1/ε
y−
e−2(1−y)/ε − e−2/ε
1 − e−2/ε
.
We use piecewise polynomials of degree p = 1, 2, 3, defined on the exponentially graded mesh (cf. Fig. 2) with N points in each direction and we compute the error in the energy norm between the exact and numerical solution. The performance of the method, for various values of ε , is shown in Fig. 3, in which the percentage relative error is plotted versus N, in a log–log scale. The straight lines have slope −p, just as Theorem 9 predicts and moreover they coincide for all ε , verifying the robustness of the method. We also show, in Fig. 4, a comparison between the proposed exponential (eXp) mesh and two other meshes from the literature: the Bakhvalov–Shishkin (B–S) mesh (see, e.g., [6] and the references therein) and another optimally constructed, layer adapted mesh with robust convergence rate from [16], which we label RTU. We compare the three methods for ε = 10−6 and p = 1, 2, 3—other values of ε gave the same results. As can be seen, all methods yield the optimal rate of convergence, but the eXp mesh performs slightly better, perhaps due to the constant being smaller in the error estimate. For our second example we set the right hand side f = 1. In this case the exact solution is unknown and we use a numerically computed reference solution uref as a substitute for the exact solution in computing the errors. For computing uref , a polynomial degree of 4 and an exponentially refined mesh with N = 512 is used. Furthermore, the exact solution of this problem does not fulfil Assumption 1, as the basic compatibility conditions f (0, 0) = f (1, 0) = f (0, 1) = f (1, 1) = 0, given in [10], are not fulfilled. Therefore the exact solution has corner singularities and reduced regularity. Fig. 5 shows the convergence behaviour for several values of ε and polynomial degrees p = 1, 2, 3. Clearly the uniformity in ε is observable for p = 1 and p = 2. The optimal convergence rate of order one is still attained for p = 1, while for p = 2 we already see a reduction from two to 1.8. For p = 3 the situation is more drastic as now the convergence rate decreases for ε becoming smaller. For ε = 10−8 we observe only a rate of 1.33. 6. Conclusions We presented the finite element analysis for a two-dimensional exponentially graded mesh, appropriate for singularly perturbed problems with exponential layers. By considering a model convection–diffusion problem, we proved robust,
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
1539
Fig. 4. Energy norm convergence: comparison of the eXp, B–S and RTU meshes.
Fig. 5. Energy norm convergence of the method using the eXp mesh and a f ≡ 1.
optimal convergence in the energy and maximum norms. Based on numerical experiments, it appears that the proposed mesh outperforms the commonly used layer adapted Bakhvalov–Shishkin mesh, hence it presents itself as an attractive alternative. References [1] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [2] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, second ed., in: Springer Series in Computational Mathematics, vol. 24, Springer-Verlag, Berlin, 2008, Convection–diffusion-reaction and flow problems. [3] G.I. Shishkin, Grid approximation of singularly perturbed boundary value problems with a regular boundary layer, Sov. J. Numer. Anal. Math. Model. 4 (1989) 397–417. [4] N.S. Bakhvalov, Towards optimization of methods for solving boundary value problems in the presence of boundary layers, Zh. Vychisl. Mat. Mat. Fiz. 9 (1969) 841–859 (in Russian). [5] C. Schwab, M. Suri, The p and hp versions of the finite element method for problems with boundary layers, Math. Comp. 65 (1996) 1403–1429. [6] T. Linß, Layer-Adapted Meshes for Reaction-Convection-Diffusion Problems, in: Lecture Notes in Mathematics, vol. 1985, Springer-Verlag, 2010. [7] C. Xenophontos, Optimal mesh design for the finite element approximation of reaction–diffusion problems, Internat. J. Numer. Methods Engrg. 53 (2002) 929–943. [8] P. Cnstantinou, C. Xenophontos, Finite element analysis of an exponentially graded mesh for singularly perturbed problems, Comput. Methods Appl. Math. 15 (2) (2015) 135–143. [9] P. Grisvard, Elliptic Problems in Nonsmooth Domains, Pitman Publishing, 1985. [10] R.B. Kellogg, M. Stynes, Corner singularities and boundary layers in a simple convection–diffusion problem, J. Differential Equations 213 (2005) 81–120. [11] T. Linß, Analysis of a Galerkin finite element method on a Bakhvalov–Shishkin mesh for a linear convection–diffusion problem, IMA J. Numer. Anal. 20 (2000) 621–632.
1540
C. Xenophontos et al. / Computers and Mathematics with Applications 72 (2016) 1532–1540
[12] Q. Lin, N. Yan, A. Zhou, A rectangle test for interpolated finite elements, in: Proc. Syst. Sci. Eng., Culture Publishing Co., Great Wall (HK), 1991, pp. 217–229. [13] M. Stynes, L. Tobiska, Using rectangular Qp elements in the SDFEM for a convection–diffusion problem with a boundary layer, Appl. Numer. Math. 58 (2008) 1789–1802. [14] T. Apel, Anisotropic Finite Elements: Local Estimates and Applications, in: Advances in Numerical Mathematics, B. G. Teubner, Stuttgart, 1999. [15] E. O’Riordan, M. Stynes, A uniformly convergent Galerkin method on a Shishkin mesh for a convection–diffusion problem, J. Math. Anal. Appl. 214 (1997) 36–54. [16] H.-G. Roos, Lj. Teofanov, Z. Uzelac, Graded meshes for higher order FEM, J. Comput. Math. 33 (2015) 1–16.