Exponentially graded mesh for a singularly perturbed problem with two small parameters

Accepted Manuscript Exponentially graded mesh for a singularly perturbed problem with two small parameters

Helena Zarin

PII: DOI: Reference:

S0168-9274(17)30139-3 http://dx.doi.org/10.1016/j.apnum.2017.06.003 APNUM 3221

To appear in:

Applied Numerical Mathematics

Received date: Accepted date:

28 November 2016 3 June 2017

Please cite this article in press as: H. Zarin, Exponentially graded mesh for a singularly perturbed problem with two small parameters, Appl. Numer. Math. (2017), http://dx.doi.org/10.1016/j.apnum.2017.06.003

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Exponentially graded mesh for a singularly perturbed problem with two small parameters Helena Zarin1 Department of Mathematics and Informatics, Faculty of Sciences, University of Novi Sad, Trg Dositeja Obradovi´ca 4, 21000 Novi Sad, Serbia

Abstract A one–dimensional singularly perturbed boundary value problem with two small perturbation parameters is numerically solved on an exponentially graded mesh. Using an h−version of the standard Galerkin method with higher order polynomials, we prove a robust convergence in the corresponding energy norm. Numerical experiments support theoretical findings. Keywords: singularly perturbed problem, two small parameters, Galerkin finite element method, exponentially graded mesh 2010 MSC: 65L11, 65L60, 65L70 1. Introduction We consider the following singularly perturbed boundary value problem ⎧ ⎪ ⎨ −ε1 u (x) + ε2 b(x)u (x) + c(x)u(x) = f (x), x ∈ Ω := (0, 1), ⎪ ⎩

u(0) = 0,

(1)

u(1) = 0,

with two small perturbation parameters 0 < ε1 , ε2  1. Let the data functions b, c, f ¯ = [0, 1] and be sufficiently smooth on Ω b(x) ≥ β > 0,

c(x) ≥ γ > 0,

c(x) −

ε2  b (x) ≥  > 0, 2

x ∈ Ω,

(2)

where β, γ,  are constants. There are two limiting cases for the perturbation parameter ε2 : • for ε2 = 0, the problem (1)–(2) is the well–known reaction–diffusion problem 1/2

1/2

whose solution has two boundary layers of the widths O(ε1 | ln ε1 |); Email address: [email protected] (Helena Zarin) Research of this author is supported by the Ministry of Education, Science and Technological Development of the Republic of Serbia under grant 174030. 1

Preprint submitted to Applied Numerical Mathematics

June 7, 2017

• when ε2 = 1, (1)–(2) reduces to the convection–diffusion with a solution exhibiting a layer at x = 1 of the width O(ε1 | ln ε1 |). For both problem classes there is a vast literature devoted to their numerical solving, see for instance books [5, 11] and references therein. Here we are interested in the reaction–convection–diffusion problem (1)–(2) and the purely two–parameter case 0 < ε1 , ε2  1. O’Malley in [8] examines the asymptotic structure of the solutions to (1) and establishes a relation between ε1 and ε2 that directly influences the layer structure in the vicinity of x = 0 and x = 1. In [4, 6, 7, 9, 13], the authors consider various robust finite difference schemes on piecewise uniform meshes of Shishkin type. Bakhvalov–type meshes for higher–order finite p+1/2

difference methods appear in [15, 16], for a special ε2 = ε1

, p > 0. The problem

(1)–(2) is analyzed in the context of finite elements as well. For example, a robust error estimate in an energy norm for the standard Galerkin method can be found in [1] for a Bakhvalov–type mesh, and in [2] for recursively graded meshes. The streamline–diffusion finite element method on the Shishkin mesh for (1)–(2) has been presented in [12] where authors prove (almost) second order of convergence in the maximum norm. In this paper we study numerical analysis of (1)–(2) on exponentially graded mesh (eXp–mesh). Its construction has been presented in [17], where it was applied in a combination with the Galerkin finite element method for the reaction–diffusion problem. It has been further used for solving reaction–diffusion systems, [19], and for two–dimensional convection–diffusion problems, [18]. The main goal of this paper is twofolds. First we generalize the results from [3], where separate error bounds have been derived for reaction–diffusion and convection– diffusion case. Furthermore, for an h−version of the standard Galerkin method with higher order polynomials, we prove a robust convergence of the optimal order in an energy norm. Unlike the Bakhvalov–type mesh that also leads to optimal convergence, [1], here we do not use Cl´emant quasi–interpolant. The rest of the paper is organized as follows. In Section 2 we review properties of the solution to (1)–(2), and we present the finite element method and the exponentially graded mesh, here adapted to resolve two boundary layers. Section 3 is devoted to the error analysis split into interpolation error and discretization error. The main 2

assertion is Theorem 3 on parameter–uniform convergence of optimal order for the Galerkin method and the eXp–mesh applied to (1)–(2). In the last section we present numerical experiments that illustrate theoretical results, and their comparison with selected layer–adapted meshes. Notation. With C we denote a generic positive constant independent of the perturbation parameters ε1 , ε2 and the number of degrees of freedom N . For a set D ⊆ Ω, Pk (D) is the set of polynomials defined on D, of the highest degree k ≥ 1. Moreover, on D we use the standard notation for Banach spaces Lq (D), Sobolev spaces H q (D), norms  · Lq (D) ,  · H q (D) and seminorms | · |H q (D) . The scalar product in L2 (D) is denoted with (·, ·)D ; we write (·, ·) when D = Ω. 2. The problem discretization Following [12], the solution u to the problem (1)–(2) can be decomposed as u = S + E0 + E1 , with ()

()

|S () (x)| ≤ C, |E0 (x)| ≤ Cμ0 e−pμ0 x , |E1 (x)| ≤ Cμ1 e−pμ1 (1−x) ,

(3)

for x ∈ Ω and 0 ≤  ≤ q, where q depends on the smoothness of the data functions. Here μ0 , μ1 define the widths of boundary layers. They are given with μ0 = − max λ0 (x),

μ1 = min λ1 (x),

x∈[0,1]

x∈[0,1]

where λ0 (x), λ1 (x) represent the solutions of the characteristic equation associated to the differential operator in (1). Since |λ0 | < λ1 , the layer at x = 1 is stronger than the layer at x = 0. Moreover, from [12, 14] we have the following properties μ0 ≤ μ1 , ε2 μ0 ≤ C,

1/2

max{μ−1 0 , ε1 μ1 } ≤ C(ε2 + ε1 ), 1/2

ε2 (ε1 μ1 )−1/2 ≤ Cε2 ,

(4) (5)

which will be frequently used in the error analysis. 2.1. Finite element method The standard weak formulation for the starting problem is: find u ∈ H01 (Ω) such that a(u, v) = (f, v),

for all v ∈ H01 (Ω), 3

with H01 (Ω) = {w ∈ H 1 (Ω) : w(0) = w(1) = 0} and the bilinear form a(w, v) = ε1 (w , v  ) + ε2 (bw , v) + (cw, v),

w, v ∈ H01 (Ω).

On an arbitrary mesh 0 = x0 < x1 < · · · < xN = 1, N ∈ N, let k ∈ N be the polynomial degree and VkN = {v ∈ L2 (Ω) : v|[xi−1 ,xi ] ∈ Pk ([xi−1 , xi ]), i = 1, 2, . . . , N } the finite element space that consists of kth degree piecewise polynomial functions defined on the general mesh. The discrete problem that defines the standard Galerkin finite element method reads: find uN ∈ VkN such that a(uN , v N ) = (f, v N ),

for all v N ∈ VkN .

(6)

Due to the last assumption in (2), the bilinear form a(·, ·) is coercive in the energy norm v2E = ε1 |v|2H 1 (Ω) + v2L2 (Ω) .

(7)

Consequently, the discrete problem (6) has a unique solution. 2.2. Discretization mesh Let N ∈ N, N ≥ 8, be divisible by 4 and τ ≥ 1 a user–chosen parameter. We need two mesh generating functions ϕ0 and ϕ1 in order to resolve two boundary layers. These functions are given with ϕ0 (t) = − ln (1 − 4Θ0 t) , ϕ1 (t) = − ln (1 − 4Θ1 (1 − t)) . The constants Θ0 , Θ1 from [3, 17], are here adapted to depend on both perturbation parameters ε1 , ε2 , i.e. Θj = 1 − exp



−pμj τ (k + 1)

 ∈ (0, 1),

The mesh points xi , i = 0, 1, . . . , N , are now defined ⎧ τ ⎪ ⎪ (k + 1)ϕ0 (ti ), ⎪ ⎪ 2pμ0 ⎪ ⎪ ⎪ ⎪   ⎨ x3N/4+1 − xN/4−1 N +1 , i− xi = xN/4−1 + ⎪ N/2 + 2 4 ⎪ ⎪ ⎪ ⎪ ⎪ τ ⎪ ⎪ ⎩1 − (k + 1)ϕ1 (ti ), 2pμ1 4

j = 0, 1.

as i = 0, 1, . . . ,

N − 1, 4

i=

N N 3N , + 1, . . . , , 4 4 4

i=

3N + 1, . . . , N. 4

(8)

If we introduce Ω0 = (0, xN/4−1 ),

Ωc = (xN/4−1 , x3N/4+1 ),

Ω1 = (x3N/4+1 , 1),

(9)

¯ 0 = [0, xN/4−1 ] and Ω ¯ 1 = [x3N/4+1 , 1] the mesh points are gradually then on the sets Ω ¯ c = [xN/4−1 , x3N/4+1 ] they are equidistant. Notice that Ωc is distributed, while on Ω well defined since ϕ0 (ϕ1 ) is monotonically increasing (decreasing) function and xN/4−1 <

τ τ 1 (k + 1)ϕ0 (1/4) = = 1 − (k + 1)ϕ1 (3/4) < x3N/4+1 . 2pμ0 2 2pμ1

Similarly to [3, 18], in the sequel we assume τ (k + 1) ln(N − 4) ≤ 1, pμj

j = 0, 1,

(10)

which does not represent a restriction in practical implementation, due to 1  μ0 in the singularly perturbed case. Under this assumption, the boundary layer functions exp(−pμ0 x) and exp(−pμ1 (1 − x)) from (3) are sufficiently small at the points xN/4−1 and x3N/4+1 , respectively. This stems from   −pμj −1 −1 + 4N −1 ∈ [4N −1 , 5N −1 ], 1 − Θj + 4N Θj = (1 − 4N ) exp τ (k + 1)

(11)

for j = 0, 1, and  max e−pμ0 xN/4−1 , e−pμ1 (1−x3N/4+1 ) τ (k+1)/2

≤ CN −τ (k+1)/2 . ≤ 1 − Θ0 + 4N −1 Θ0

(12)

Notice that from (11) we also have max |ϕj | =

4Θj 4Θj ≤ < N, −1 1 − Θj + 4N Θj 4N −1

j = 0, 1,

(13)

which directly influences the properties of the mesh step sizes hi = xi − xi−1 , i = 1, 2, . . . , N , cf. [10]. For the subsets of indices I0 = {1, 2, . . . , N/4 − 1}, Ic = {N/4, N/4 + 1 . . . , 3N/4 + 1}, I1 = {3N/4 + 2, 3N/4 + 3, . . . , N }, from the definition (8) of the mesh points and from (13), we directly obtain ⎧ τ ⎪ (k + 1)N −1 max |ϕj | ≤ Cμ−1 i ∈ Ij , j = 0, 1, ⎨ hi ≤ j , 2pμj ⎪ ⎩ h ≤ 2N −1 , i ∈ Ic . i 5

(14)

Introducing mesh characterizing functions ψj = e−ϕj ,

j = 0, 1,

we can derive other properties of element widths from the layer region. For example, for i ∈ I0 we also have τ −1 (k + 1)N −1 max (−ψ0 (t)eϕ0 (t) ) ≤ Cμ−1 max |ψ0 |eϕ0 (ti ) 0 N t∈[ti−1 ,ti ] 2pμ0   2pμ0 xi −1 −1 = Cμ0 N exp , τ (k + 1)

hi ≤

(15)

while for i ∈ I1 we can similarly get τ −1 (k + 1)N −1 max −ψ1 (t)eϕ1 (t) ≤ Cμ−1 max |ψ1 |eϕ1 (ti−1 ) 1 N t∈[ti−1 ,ti ] 2pμ1   2pμ1 (1 − xi−1 ) −1 −1 . (16) = Cμ1 N exp τ (k + 1)

hi ≤

Notice that max |ψj | = 4Θj < 4. The previous inequalities (15)–(16) are further used to estimate interpolation error. 3. Error analysis Let uI ∈ VkN be the standard nodal interpolant of the solution function u. Then we can split the error u − uN in the following way u − uN = η + χ,

η = u − uI ,

χ = uI − uN ∈ VkN .

In the sequel, we separately analyze energy norms of the interpolation error η and the discretization error χ on the exponentially graded mesh (8). 3.1. Interpolation error Following the solution decomposition u = S + E0 + E1 from (3), for the interpolation function uI ∈ VkN of u, we write uI = S I + E0I + E1I . The standard interpolation result |g − g I |H q (xi−1 ,xi ) ≤ Chk+1−q |g|H k+1 (xi−1 ,xi ) , i

(17)

for k ∈ N, 0 ≤ q ≤ k + 1 and a function g ∈ H k+1 (xi−1 , xi ), applied to ηS := S − S I gives



−(k+1−q) |ηS |H q (Ω) ≤ C μ0 + N −(k+1−q) . 6

(18)

Here we have used |S|H k+1 (Ω) ≤ C from (3), the bounds (14) for the mesh steps and μ0 ≤ μ1 . Now we derive H q −estimates for η0 := E0 − E0I , separately on Ω0 and Ωc ∪ Ω1 . Let (xi−1 , xi ) ⊂ Ω0 be an arbitrary element. The local interpolation error (17) and the properties (3) of the solution component E0 , together yield |η0 |2H q (xi−1 ,xi )

≤

2(k+1−q) Chi |E0 |2H k+1 (xi−1 ,xi ) 2(k+1−q)

= Chi

2(k+1−q) 2(k+1) Chi μ0

≤

2k+1

μ0

e−2pμ0 xi−1 − e

2(k+1−q) 2k+1 −2pμ0 xi pμ0 hi μ0 e e

= Chi



xi

e−2pμ0 x dx

xi−1

 −2pμ0 xi

sinh(pμ0 hi ).

(19)

Since μ0 hi ≤ C, we have epμ0 hi ≤ C and also  sinh(pμ0 hi ) ≤ Cμ0 hi = C(ϕ0 (ti ) − ϕ0 (ti−1 )) = C 



ti

=C ti−1

(−ψ0 (t))eϕ0 (t) dt 

≤ CN −1 exp

2pμ0 xi τ (k + 1)

≤ C exp



ti ti−1

ϕ0 (t)dt

2pμ0 xi τ (k + 1)



ti ti−1

(−ψ0 (t))dt

.

(20)

Now, from (19), (20) and (15) we get |η0 |2H q (xi−1 ,xi )

≤

2(k+1−q) 2k+1 −1 −2pμ0 xi Chi μ0 N e

 exp

2pμ0 xi τ (k + 1)



≤ Cμ2q−1 N −2(k−q)−3 · Λ 0 with

 Λ = exp

2pμ0 xi τ (k + 1)

2(k+1−q)

 exp

2pμ0 xi τ (k + 1)

 e−2pμ0 xi .

If we choose the parameter τ to be larger than (2k + 3)/(k + 1), we deduce Λ ≤ 1 since 4p(k + 1 − q)μ0 xi 2pμ0 xi 2pμ0 xi + − 2pμ0 xi = (−2q + 2k + 3 − τ (k + 1)) ≤ 0. τ (k + 1) τ (k + 1) τ (k + 1) Taking sum over all intervals from Ω0 we obtain  |η0 |H q (Ω0 ) =



1/2 |η0 |2H q (xi−1 ,xi )

1/2 ≤ C μ2q−1 N −2(k−q)−3 N 0

i∈I0 q−1/2

≤ Cμ0

N −(k+1−q) .

(21)

7

On Ωc ∪ Ω1 we estimate |η0 |H q (Ωc ∪Ω1 ) using the triangle inequality. First   xi  2q 2 2 |E0 |H q (Ωc ∪Ω1 ) = |E0 |H q (xi−1 ,xi ) ≤ Cμ0 e−2pμ0 x dx i∈Ic ∪I1

≤ Cμ2q−1 0

i∈Ic ∪I1



e−2pμ0 xi−1 − e

 −2pμ0 xi

xi−1

e−2pμ0 xN/4−1 ≤ Cμ2q−1 0

i∈Ic ∪I1

Cμ2q−1 N −(2k+3) . 0

≤

(22)

Here we have used (3), (12) and previously established assumption on τ . The inverse inequality I |E0I |H q (xi−1 ,xi ) ≤ Ch−q i E0 L2 (xi−1 ,xi ) ,

together with the stability E0I L∞ (xi−1 ,xi ) ≤ CE0 L∞ (xi−1 ,xi ) , (14), (3) and (12), imply |E0I |2H q (Ωc ∪Ω1 ) =



|E0I |2H q (xi−1 ,xi ) ≤ C

i∈Ic ∪I1

≤C





h−2q E0I 2L2 (xi−1 ,xi ) i

i∈Ic ∪I1

h−2q+1 E0 2L∞ (xi−1 ,xi ) ≤ C(N 2q−1 + μ2q−1 )N −τ (k+1)+1 1 i

i∈Ic ∪I1

≤ C(N 2q−1 + μ2q−1 )N −2(k+1) . 1

(23)

From (22) and (23) we conclude q−1/2

|η0 |H q (Ωc ∪Ω1 ) ≤ Cμ0

q−1/2

N −(k+3/2) + CN q−k−3/2 + Cμ1

N −(k+1) .

(24)

Analogously, for the interpolation error η1 := E1 − E1I it holds q−1/2

N −(k+1−q) ,

q−1/2

N −(k+3/2) + CN q−k−3/2 + Cμ0

|η1 |H q (Ω1 ) ≤ Cμ1 |η1 |H q (Ω0 ∪Ωc ) ≤ Cμ1

(25) q−1/2

N −(k+1) .

(26)

Collecting (18), (21) and (24)–(26) for q = 0, 1, and with (4), we can summarize the result on the energy norm of the interpolation error on the discretization mesh (8). Theorem 1. Let u be the solution of the reaction–convection–diffusion problem (1)– (2) and uI ∈ VkN its interpolant constructed on the mesh (8). Under the assumption (10) and for τ≥

2k + 3 , k+1

the interpolation error satisfies 1/2

u − uI E ≤ C(ε2 + ε1 )1/2 N −k + CN −(k+1) . 8

3.2. Discretization error We now turn to the estimation of χE , with χ = uI − uN ∈ VkN . The bilinear form a(·, ·) is coercive with respect to the energy norm (7). Since the finite element method (6) posesses the Galerkin orthogonality property a(u − uN , v N ) = 0,

for all v N ∈ VkN ,

then we can start the analysis with min{1, }χ2E ≤ a(χ, χ) = −a(η, χ) ≤ ε1 |(η  , χ )| + ε2 |(bη  , χ)| + |(cη, χ)|.

(27)

The first and the third term in (27) can be easily estimated with the Cauchy– Schwarz inequality ε1 |(η  , χ )| + |(cη, χ)| ≤ ε1 |η|H 1 (Ω) |χ|H 1 (Ω) + CηL2 (Ω) χL2 (Ω) ≤ CηE χE ,

(28)

but the second term needs to be handled in more details. Recalling the notation from the previous subsection, let η = ηS + η0 + η1 . Then we can write (bη  , χ) = (bηS , χ) + (bη0 , χ)Ω0 ∪Ωc − (η0 , b χ)Ω1 − (η0 , bχ )Ω1 − (η1 , b χ) − (η1 , bχ ).

(29)

Notice that (η0 , b χ)Ω1 and (η1 , b χ) can be bounded as in (28). The remaining terms, multiplied with ε2 , will be considered separately. A straightforward consequence of (18) is ε2 |(bηS , χ)| ≤ Cε2 |ηS |H 1 (Ω) χL2 (Ω) ≤ Cε2 N −k χE .

(30)

Similarly, from (21), (24) and ε2 μ0 ≤ C we get ε2 |(bη0 , χ)Ω0 ∪Ωc | ≤ Cε2 |η0 |H 1 (Ω0 ∪Ωc ) χL2 (Ω0 ∪Ωc ) 1/2

1/2

≤ Cε2 (μ0 N −k + N −k−1/2 )χE ≤ Cε2 N −k χE . In order to estimate ε2 (η0 , bχ )Ω1 we shall use η0 L∞ (Ω1 ) ≤ CE0 L∞ (Ω1 ) ≤ Ce−pμ0 xN/4−1 ≤ CN −τ (k+1)/2 ,

9

(31)

that follows from the L∞ −stability of the interpolation operator, (3) and (12). Now −1/2

ε2 |(η0 , bχ )Ω1 | ≤ Cε2 η0 L2 (Ω1 ) |χ|H 1 (Ω1 ) ≤ Cε2 ε1 −1/2

≤ Cε2 ε1

η0 L2 (Ω1 ) χE

η0 L∞ (Ω1 ) (meas(Ω1 ))1/2 χE .

Each subinterval from Ω1 has the length hi ≤ Cμ−1 1 . Thus meas(Ω1 ) ≤ C(max hi )N ≤ Cμ−1 1 N, i∈Ic

and consequently −1/2

ε2 |(η0 , bχ )Ω1 | ≤ Cε2 ε1

1/2

1/2 N −τ (k+1)/2 (μ−1 χE ≤ Cε2 N −(k+1) χE , 1 N)

(32)

due to (5) and the assumption on τ from Theorem 1. Finally, it remains to estimate the last term ε2 (η1 , bχ ). On the set Ω1 , the Cauchy– Schwarz inequality, the L2 −bound from (25) and (5) yield −1/2

ε2 |(η1 , bχ )Ω1 | ≤ Cε2 η1 L2 (Ω1 ) |χ|H 1 (Ω1 ) ≤ Cε2 ε1

η1 L2 (Ω1 ) χE

1/2

≤ Cε2 (ε1 μ1 )−1/2 N −(k+1) χE ≤ Cε2 N −(k+1) χE .

(33)

On the equidistant part of the mesh, for |χ|H 1 (Ωc ) we use the inverse inequality and hi = O(N −1 ), i ∈ Ic . Therefore ε2 |(η1 , bχ )Ωc | ≤ Cε2 η1 L2 (Ωc ) |χ|H 1 (Ωc ) ≤ Cε2 N −(k+3/2) N χL2 (Ωc ) ≤ Cε2 N −(k+1/2) χE .

(34)

In order to use the inverse inequality on Ω0 as well, we need the lower bounds of the mesh step sizes. Since Θj , j = 0, 1, can be bounded from below with a constant 1 − exp(−p/(τ (k + 1))) in the singularly perturbed regime, we get −1 −1 −1 hi ≥ Cμ−1 min ϕ0 = Cμ−1 Θ0 ≥ Cμ−1 . 0 N 0 N 0 N

Hence −1/2

ε2 |(η1 , bχ )Ω0 | ≤ Cε2 η1 L2 (Ω0 ) |χ|H 1 (Ω0 ) ≤ Cε2 (μ0 −1/2

≤ Cμ0

N −k χE ,

N −(k+1) )(μ0 N )χL2 (Ω0 ) (35)

invoking ε2 μ0 ≤ C and (26). 10

Collecting (29)–(35), from (4) we deduce −1/2

ε2 |(bη  , χ)| ≤ C(ε2 ηE + μ0

N −k )χE .

At the end, combining the previous inequality with (28) and after canceling with χE in (27), we come to the final estimate of the discretization error. Theorem 2. Let u be the solution of the reaction–convection–diffusion problem (1)– (2), uI ∈ VkN its interpolant constructed on the mesh (8), and uN the finite element approximation of u from (6). Under the assumption (10) and for τ ≥ (2k+3)/(k+1), for the discretization error it holds 1/2

uI − uN E ≤ C(ε2 + ε1 )1/2 N −k + Cu − uI E . The main result on robustness of the Galerkin finite element method (6) on the exponentially graded mesh (8) for the two–parameter singularly perturbed problem (1)–(2) follows directly from the triangle inequality, Theorem 1 and Theorem 2. Theorem 3. Let u be the solution of the reaction–convection–diffusion problem (1)– (2) and uN its finite element approximation from (6) constructed on the discretization mesh (8). Under the assumptions of Theorem 1, the energy norm of the error satisfies 1/2

u − uN E ≤ C(ε2 + ε1 )1/2 N −k + CN −(k+1) . The previous estimate represents a generalization of [3] where the energy norm error O(N −k ) has been proved for the same method and the corresponding eXp–mesh, separately for reaction–diffusion and convection–diffusion problems. 4. Numerical tests In this section we present the results of numerical experiments for the Galerkin method (6) applied to the eXp–mesh (8). We take the test problem from [12] −ε1 u (x) + ε2 u (x) + u(x) = cos πx,

x ∈ Ω,

u(0) = u(1) = 0,

with the known exact solution. For fixed perturbation parameters ε1 , ε2 , polynomial degree k and the number of mesh subintervals N , we compute the energy norm of the error eN := u − uN E and the rate of convergence pN with the standard formula pN =

ln eN − ln e2N . ln 2 11

All integrals are calculated using Gauss–Legendre 5−point quadrature. The results are also compared with energy norm errors for (6) applied to the standard Shishkin mesh from [12], the Bakhvalov mesh from [1] and the recursive Duran mesh from [2]. We first distinguish two different choices of perturbation parameters. Taking ε1 = 10−10 and ε2 = 10−3 , the problem solution has a stronger layer at x = 1 than the layer at x = 0, which illustrates the case when ε1  ε22  1, μ0 = O(ε−1 2 ) and μ1 = O(ε2 ε−1 1 ). For polynomial degrees k = 1, 2, 3, the results on the Shishkin, Bakhvalov and the eXp–mesh are presented in Table 1, showing the smallest errors and the predicted rate of convergence for the latter mesh. Table 2 contains similar results but for ε1 = 10−5 and ε2 = 10−8 , that corresponds to the case when ε22  −1/2

ε1  1, μ0 = O(ε1

−1/2

) and μ1 = O(ε1

), where the boundary layers are similar to

the reaction–diffusion case. We also calculate energy norm errors eN R for a fixed k and different number of mesh elements N , now considered on a range of perturbation parameters such that N −5 −10 eN ; ε2 = 10−3 , . . . , 10−10 }. R = max{u − u E : ε1 = 10 , . . . , 10

The results on Shishkin, Bakhvalov and eXp–mesh are depicted in Figure 1. In Table 3 and Table 4 we compare energy norm errors for the Duran mesh and the eXp–mesh. Here we first choose the mesh parameter h ∈ (0, 1) from the Duran mesh, then we calculate the number of mesh subintervals M0 on [0, 1/2] and M1 on [1/2, 1] in the Duran mesh, and at the end we choose N for the eXp–mesh such that it is divisible by 4 and close to M0 + M1 . For the polynomial degree k = 1 and the proposed choices of perturbation parameters, the eXp–mesh produces slightly smaller errors when ε1 = 10−10 and ε2 = 10−3 (Table 3), while it is moderately outperformed by the Duran mesh for ε1 = 10−5 and ε2 = 10−8 (Table 4). The presented numerical results confirm robustness of the standard Galerkin method on the exponentially graded mesh for the reaction–convection–diffusion problem (1)–(2). As already indicated in [3], the eXp–mesh does not require an explicit definition of transition points, which is one of the key properties of layer–adapted Shishkin– and Bakhvalov–type meshes, cf. [5]. The advantage of the optimal eXp– mesh over the non–optimal standard Shishkin mesh is evident, both in theory and in our experiments. Though the Bakhvalov mesh for (1)–(2) is optimal too, [1], the 12

Table 1: Error and rate of convergence, ε1 = 10−10 , ε2 = 10−3 , k = 1, 2, 3, for Shishkin (S–), Bakhvalov (B–) and exponentially graded (eXp–) mesh.

N

S–mesh

B–mesh

eXp–mesh

k=1

u − uN E

pN

u − uN E

pN

u − uN E

pN

26

1.376(−2)

0.786

3.485(−3)

1.010

2.275(−3)

1.009

27

7.983(−3)

0.823

1.731(−3)

1.002

1.130(−3)

1.002

28

4.513(−3)

0.838

8.639(−4)

1.001

5.642(−4)

1.001

29

2.525(−3)

0.851

4.318(−4)

1.000

2.820(−4)

1.000

210

1.400(−3)

0.863

2.158(−4)

1.000

1.410(−4)

1.000

211

7.693(−4)

−

1.079(−4)

−

7.049(−5)

−

k=2

u − uN E

pN

u − uN E

pN

u − uN E

pN

26

9.524(−3)

1.162

9.417(−4)

1.969

2.079(−4)

1.999

27

4.256(−3)

1.417

2.406(−4)

1.993

5.198(−5)

2.006

28

1.594(−3)

1.581

6.044(−5)

2.003

1.294(−5)

2.021

29

5.326(−4)

1.669

1.508(−5)

2.015

3.189(−6)

2.046

210

1.674(−4)

1.718

3.731(−6)

2.036

7.723(−7)

2.043

211

5.091(−5)

−

9.097(−7)

−

1.874(−7)

−

k=3

u − uN E

pN

u − uN E

pN

u − uN E

pN

26

7.189(−4)

2.123

1.531(−5)

2.997

1.112(−5)

2.999

27

1.650(−4)

2.343

1.918(−6)

2.999

1.392(−6)

3.000

28

3.252(−5)

2.463

2.399(−7)

3.000

1.740(−7)

3.000

29

5.894(−6)

2.536

2.999(−8)

3.000

2.175(−8)

3.000

210

1.017(−6)

2.585

3.749(−9)

2.999

2.720(−9)

2.997

211

1.694(−7)

−

4.689(−10)

−

3.407(−10)

−

error analysis on the eXp–mesh relies on more standard arguments and avoids using a Cl´emant quasi–interpolant. The Duran mesh from [2], as a recursively graded mesh without transition points, in our experiments produces smaller energy norm errors in comparison to the eXp–mesh. Nevertheless, for higher order polynomials we expect the errors and the order of convergence for the Duran mesh to deteriorate due to a logarithmic dependence of the mesh parameter h.

13

Table 2: Error and rate of convergence, ε1 = 10−5 , ε2 = 10−8 , k = 1, 2, 3, for Shishkin (S–), Bakhvalov (B–) and exponentially graded (eXp–) mesh.

N

S–mesh

B–mesh

eXp–mesh

k=1

u − uN E

pN

u − uN E

pN

u − uN E

pN

26

3.244(−2)

0.741

8.728(−3)

0.999

5.686(−3)

1.002

27

1.941(−2)

0.793

4.367(−3)

1.000

2.839(−3)

1.000

28

1.119(−2)

0.825

2.183(−3)

1.003

1.419(−3)

1.000

29

6.320(−3)

0.846

1.090(−3)

1.003

7.096(−4)

0.998

210

3.515(−3)

0.862

5.435(−4)

1.002

3.552(−4)

0.994

211

1.934(−3)

−

2.713(−4)

−

1.784(−4)

−

k=2

u − uN E

pN

u − uN E

pN

u − uN E

pN

26

6.005(−3)

1.424

4.704(−4)

1.997

3.605(−4)

1.999

27

2.238(−3)

1.564

1.178(−4)

1.999

9.020(−5)

2.000

28

7.570(−4)

1.642

2.948(−5)

1.999

2.256(−5)

2.000

29

2.425(−4)

1.690

7.376(−6)

1.999

5.639(−6)

2.000

210

7.514(−5)

1.723

1.845(−6)

2.000

1.410(−6)

1.990

211

2.276(−5)

−

4.613(−7)

−

3.549(−7)

−

k=3

u − uN E

pN

u − uN E

pN

u − uN E

pN

26

1.713(−3)

2.083

3.800(−5)

2.992

2.789(−5)

2.996

27

4.043(−4)

2.320

4.777(−6)

2.992

3.497(−6)

2.999

28

8.095(−5)

2.454

6.004(−7)

2.981

4.374(−7)

2.999

29

1.477(−5)

2.532

7.606(−8)

2.970

5.470(−8)

2.984

210

2.553(−6)

2.584

9.710(−9)

2.975

6.914(−9)

2.391

211

4.259(−7)

−

1.235(−9)

−

1.318(−9)

−

14

eRN 0.1   

  

0.001



  

  

105





  



  

  

 



  



  



  



  

107

  



  



Smesh k1





Bmesh k1

 



eXpmesh k1





Smesh k2



Bmesh k2

  



eXpmesh k2



Smesh k3



100

Bmesh k3



109



eXpmesh k3

N

1000

Figure 1: Energy norm eN R with k = 1, 2, 3, for Shishkin (S–), Bakhvalov (B–) and exponentially graded (eXp–) mesh.

Table 3: Error for ε1 = 10−10 , ε2 = 10−3 , k = 1, on Duran (D–) and exponentially graded (eXp–) mesh.

D–mesh

eXp–mesh

h

M0

M1

u − uN E

N

u − uN E

2−1

19

41

3.866(−2)

60

2.010(−2)

2−2

36

77

1.596(−2)

116

1.031(−2)

2−3

72

150

6.855(−3)

224

5.329(−3)

2−4

150

302

3.135(−3)

452

2.639(−3)

2−5

316

615

1.495(−3)

932

1.280(−3)

2−6

671

1265

7.303(−4)

1936

6.161(−4)

References [1] M. Brdar, H. Zarin, A singularly perturbed problem with two parameters on a Bakhvalov-type mesh, J. Comput. Appl. Math. 292 (2016), 307–319. [2] M. Brdar, H. Zarin, On graded meshes for a two–parameter singularly perturbed problem, Appl. Math. Comput. 282 (2016), 97–107. [3] P. Constantinou, C. Xenophontos, Finite Element Analysis of an Exponentially Graded Mesh for Singularly Perturbed Problems, Comput. Methods Appl. Math. 15(2) (2015), 135–143.

15

Table 4: Error for ε1 = 10−5 , ε2 = 10−8 , k = 1, on Duran (D–) and exponentially graded (eXp–) mesh.

D–mesh

eXp–mesh

h

M0

M1

u − uN E

N

u − uN E

2−1

16

16

8.229(−2)

32

9.460(−2)

2−2

30

30

3.673(−2)

60

5.015(−2)

2−3

62

62

1.671(−2)

124

2.422(−2)

2−4

131

131

7.830(−3)

262

2.050(−2)

2−5

279

279

3.768(−3)

560

5.361(−3)

2−6

596

596

1.845(−3)

1192

2.521(−3)

[4] 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. [5] T. Linss, Layer–Adapted Meshes for Reaction–Convection–Diffusion Problems, Springer–Verlag Berlin Heidelberg, 2010. [6] 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. [7] S. Natesan, J. L. Gracia, C. Clavero, Singularly perturbed boundary–value problems with two small parameters – a defect correction approach, Proceedings of the International Conference on Boundary and Interior Layers – Computational and Asymptotic Methods (BAIL 2004, ONERA-Centre de Toulouse, France), 1–6. [8] R. E. O’Malley Jr, Introduction to singular perturbations, Academic Press, New York, 1974. [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] H.-G. Roos, T. Linss, Sufficient conditions for uniform convergence on layer-adapted grids, Computing 63 (1999), 27–45.

16

[11] H.-G. Roos, M. Stynes, L. Tobiska, Robust Numerical Methods for Singularly Perturbed Differential Equations, Springer Series in Computational Mathematics 24, Springer–Verlag Berlin Heidelberg, 2008. [12] 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. [13] K. Surla, Z. Uzelac, Lj. Teofanov, The discrete minimum principle for quadratic spline discretization of a singularly perturbed problem, Math. Comput. Simulat. 79(8) (2009), 2490–2505. [14] Lj. Teofanov, H.-G. Roos, An elliptic singularly perturbed problem with two parameters I: Solution decomposition, J. Comput. Appl. Math. 206 (2007), 1082–1097. [15] R. Vulanovi´c, A higher–order scheme for quasilinear boundary value problems with two small parameters, Computing 67(4) (2001), 287–303. [16] R. Vulanovi´c, Special meshes and higher–order schemes for singularly perturbed boundary values problems, Novi Sad Journal of Mathematics 31(1) (2001), 1–7. [17] C. Xenophontos, Optimal mesh design for the finite element approximation of reaction– diffusion problems, Int. J. Num. Meth. Eng. 53 (2002), 929–943. [18] C. Xenophontos, S. Franz, L. Ludwig, Finite element approximation of convection– diffusion problems using an exponentially graded mesh, Computers and Mathematics with Applications 72 (2016), 1532–1540. [19] C. Xenophontos, L. Oberbroeckling, A numerical study on the finite element solution of singularly perturbed systems of reaction–diffusion problems, Appl. Math. Comp. 187 (2007), 1351–1367.

17