A parameter robust second order numerical method for a singularly perturbed two-parameter problem

A parameter robust higher order numerical method for a singularly perturbed two–parameter problem∗ J.L. Gracia†

E. O’Riordan‡

M.L. Pickett§

Abstract In this paper a second order monotone numerical method is constructed for a singularly perturbed ordinary differential equation with two small parameters affecting the convection and diffusion terms. The monotone operator is combined with a piecewise-uniform Shishkin mesh. An asymptotic error bound in the maximum norm is established theoretically whose error constants are show to be independent of both singular perturbation parameters. Numerical results are presented which support the theoretical results. Keywords: Singular perturbation, uniform convergence, Shishkin mesh, high order. AMS subject classifications: 65N12, 65N30, 65N06

1

Introduction

Differential equations with a small parameter 0 < ε ≤ 1 multiplying the highest order derivatives are called singularly perturbed differential equations. Typically, the solutions of such equations have steep gradients (proportional to ε−p , p > 0) in narrow layer regions of the domain. Classical numerical methods are inappropriate for singularly perturbed problems [5]. For any numerical method applied to a singularly perturbed problem, it is informative to establish an asymptotic pointwise error bound of the form ku − U N kΩN ≤ Cp N −p , p > 0, where U N are the numerical approximations, N is the number of mesh elements used in each coordinate direction, u is the solution of the continuous problem and kvkD = maxx∈D |v(x)| is the maximum pointwise norm. A numerical method is said to be parameter-uniform if the error constant Cp is independent of any perturbation parameters and the mesh parameter N . A common approach to constructing a parameter-uniform numerical method is to combine a standard monotone finite difference operator with an appropriate piecewise-uniform mesh (Shishkin mesh [5, 11]) where half the mesh points are concentrated in the layer regions. In the case of the singularly perturbed reaction-diffusion problem εu00 − bu = f (x), b(x) ≥ β > 0,

u(0), u(1) given,

∗ This research was partially supported by the project MEC/FEDER MTM 2004-019051, the grant EUROPA XXI of the Caja de Ahorros de la Inmaculada and by the National Center for Plasma Science and Technology Ireland. † Departamento de Matem´ atica Aplicada, Universidad de Zaragoza, Spain. email: [email protected] ‡ School of Mathematical Sciences, Dublin City University, Glasnevin, Dublin 9, Ireland. email: [email protected] § School of Mathematical Sciences, Dublin City University, Glasnevin, Dublin 9, Ireland. email: [email protected]

1

the standard central difference operator on an appropriate p Shishkin mesh (transition points between coarse and fine mesh are taken to be σ = min{0.25, 2 ε/β ln N }) produces a second order (up to logarithmic factor) parameter-uniform error bound [7] of the form ku − U N kΩN ≤ C(N −1 ln N )2 . In [6] the authors analyze a monotone finite difference scheme of fourth order (up to logarithmic factor) using a classical higher order compact finite difference operator on a Shishkin mesh. In the case of the singularly perturbed convection-diffusion problem εu00 + au0 = f (x), a(x) ≥ α > 0,

u(0), u(1) given,

the standard upwind difference operator on an appropriate Shishkin mesh (transition point of the form σ = min{0.5, αε ln N }) produces a first order (up to logarithmic factor) parameter-uniform error bound [5] ku − U N kΩN ≤ CN −1 ln N. For the above convection-diffusion problem, Stynes and Roos [13] established that a numerical method composed of the central difference operator in the layer region (0, σ) combined with the mid-point scheme [1] outside the layer region [σ, 1) on a Shishkin mesh with σ = min{0.5, 2ε α ln N } is a monotone numerical method; and when σ < 0.5, it satisfies a parameter-uniform error bound of the form ½ CN −1 (ε + N −1 ), if xi ∈ [σ, 1], ku − U kΩN ≤ C(N −1 ln N )2 , if xi ∈ [0, σ). This is essentially second order in the layer and first order outside the layer. Other high order parameter-uniform methods for the convection-diffusion problem, include a weighted monotone finite difference scheme [3], the non-monotone centered difference scheme [2] and the HODIE [4] on appropriate Shishkin meshes. In this paper, we construct a monotone numerical method for the two parameter problem εu00 + µau0 − bu = f (x), a(x) ≥ α > 0, b(x) ≥ β > 0,

u(0), u(1) given,

where 0 < ε ≤ 1, 0 ≤ µ ≤ 1. This problem encompasses both the reaction-diffusion problem when µ = 0 and the convection-diffusion problem when µ = 1. The nature of the continuous solution of this two parameter problem was examined by O’Malley [8], where the ratio of µ2 to ε was identified as significant. Parameter-uniform methods on a uniform mesh were constructed in [12]. In [9] and [10], the standard upwind finite difference operator on two different choices of Shishkin mesh was shown to be parameter-uniform of first order. In this paper, we construct a second order parameter-uniform numerical method on the Shishkin mesh given in [10] for this two parameter problem. The finite difference operator is a combination of the central difference, midpoint and standard upwinded difference operators. These three finite difference operators are monotone in various subdomains of the parameter space P = {(ε, µ, N ) | 0 < ε ≤ 1, 0 ≤ µ ≤ 1, N ≥ N0 }. When analysing the errors in the regular component in the coarse mesh region, the following interplay between the requirements of monotonicity and high order truncation error are worth noting. The standard upwinded operator is always monotone and has a second order truncation error when both ε and µ are relatively small so that ε, µ ≤ CN −1 . The central difference operator is monotone if ε is relatively large and µ is relatively small so that N ε ≥ C1 µ and it has second order truncation error away from the transition points. The midpoint scheme is monotone for all ε and for µ relatively large so that µN ≥ C2 . For ε ≤ CN −1 , the midpoint scheme has second order truncation error at all the mesh points. All of these ingredients are combined in this paper to produce a numerical scheme which is monotone and is uniformly second order (up to logarithmic factors) for all values of the parameters ε and µ. 2

2

The continuous problem

Consider the singularly perturbed boundary value problem ½

Lε,µ u ≡ εu00 + µau0 − bu = f, u(0) = u0 , u(1) = u1 ,

x ∈ Ω = (0, 1),

(2.1)

where 0 < ε, µ ≤ 1 are two small parameters. The coefficients a, b and f are sufficiently smooth and a(x) ≥ α > 0,

b(x) ≥ β > 0,

∀x ∈ Ω,

b γ = min . Ω a

The proof of the following comparison principle for the differential operator is standard. Lemma 1 (Comparison Principle). Let v ∈ C 2 (Ω). If v(0) ≥ 0, v(1) ≥ 0 and Lε,µ v(x) ≤ 0, ∀ x ∈ Ω, then v(x) ≥ 0, ∀ x ∈ Ω. An immediate consequence of this comparison principle is the following parameter uniform bound on the solution u. Lemma 2. If u is the solution of the boundary value problem (2.1), then ·

¸ 1 ||u|| ≤ max{|u0 |, |u1 |} + ||f || , β

∀ x ∈ Ω,

where k · k denotes the pointwise maximum norm. The derivatives of the solution of this problem satisfy the following parameter-explicit bounds. Theorem 3. Assuming that a, b, f ∈ C 2 (Ω), the derivatives of the solution u of (2.1) satisfy the following bounds ||u

(k)

|| ≤

||u(3) || (4)

||u

||

C √ k ( ε)

≤

C √ ( ε)3

≤

C √ 4 ( ε)

Ã

µ 1+

Ã

µ 1+

Ã

µ 1+

µ √ ε µ √ ε µ √ ε

¶k ! max{||u||, ||f ||}, ¶3 ! ¶4 !

k = 1, 2,

1 max{||u||, ||f ||} + ||f 0 ||, ε max{||u||, ||f ||} +

1 0 1 ||f || + kf 00 k, 2 ε ε

where the constant C depends only on ka(i) k, kb(i) k, for i = 0, 1, 2 and is independent of ε and µ. Proof. Obtain the bounds on the first two derivatives using the argument in Lemma 2.2 of [9]. The bounds on the third and fourth derivatives come from differentiating the differential equation in (2.1) and using the inequality n X

ti ≤ (n + 1)(1 + tn ),

i=0

3

t ≥ 0.

3

Decomposition of the continuous solution

The solution of the boundary value problem (2.1) can be decomposed into the following sum u = v + wL + wR ,

(3.2a)

where Lv = f, LwL = 0,

v(0), v(1) suitably chosen, wL (0) = u(0) − v(0) − wR (0), wL (1) = 0,

LwR = 0, wR (0) suitably chosen, wR (1) = u(1) − v(1).

(3.2b) (3.2c) (3.2d)

The function v is called the regular component and wL and wR are the left and right, respectively, singular components of the solution u. A major part of this section will involve establishing that the first three derivatives of the regular component are bounded independently of the perturbation parameters, which is a natural extension of the decomposition given in [10] for the solution of the time-dependent version of problem (2.1). We first consider the case of αµ2 ≤ γε. We have the following secondary decomposition of the regular component v √ √ √ v(x; ε, µ) = v0 (x) + εv1 (x; ε, µ) + ( ε)2 v2 (x; ε, µ) + ( ε)3 v3 (x; ε, µ), (3.3a) where √ µ µ εv000 + √ av00 , bv2 = εv100 + √ av10 , (3.3b) ε ε √ µ Lε,µ v3 = − εv200 − √ av20 , v3 (0; ε, µ) = v3 (1; ε, µ) = 0. (3.3c) ε √ √ We see that v(0; ε, µ) = v0 (0) + εv1 (0; ε, µ) + εv2 (0; ε, µ) and v(1; ε, µ) = v0 (1) + εv1 (1; ε, µ) + εv2 (1; ε, µ). Assuming sufficient smoothness on the coefficients (a ∈ C 6 (Ω), b, f ∈ C 8 (Ω)) and noting that αµ2 ≤ γε, we see that v0 and it’s derivatives up to order eight, v1 and it’s derivatives up to sixth order and v2 and its derivatives up to order four are bounded independently of ε and µ. We proceed to analyse the final term v3 (x; ε, µ). Using the minimum principle for Lε,µ and a suitable barrier function we obtain (see [9] Lemma 2.1), −bv0 = f,

bv1 =

√

© ª 1 kv3 k ≤ max |v3 (0)|, |v3 (1)| + (kv200 k + kv20 k). β Applying the bounds on v2 we therefore have kv3 k ≤ C. Using the differential equation (3.3c) and the Mean Value Theorem on an interval of width noting that αµ2 ≤ γε, we obtain (see [9] Lemma 2.2), © ª C C (k) kv3 k ≤ √ k max ||v3 ||, ||v20 ||, ||v200 || ≤ √ k , ( ε) ( ε)

√

ε and

k = 1, 2.

Differentiating the equation (3.3c) and using the above bounds we also obtain C (k) kv3 k ≤ √ k , ( ε)

k = 3, 4.

Substituting all of these bounds for v0 (x; µ), v1 (x; ε, µ), v2 (x; ε, µ) and v3 (x; ε, µ) into the equation for v(x; ε, µ) gives us µ ¶ 1 (k) kv k ≤ C 1 + √ k−3 , for 0 ≤ k ≤ 4, when αµ2 ≤ γε. (3.4) ( ε) 4

When αµ2 ≥ γε we consider the following alternative secondary decomposition of the regular component v(x; ε, µ) = v0 (x; µ) + εv1 (x; µ) + ε2 v2 (x; µ) + ε3 v3 (x; ε, µ), (3.5a) where Lµ v 0 Lµ v 1 Lµ v 2 Lε,µ v3 (x; ε, µ)

= f (x),

v0 (1; µ) chosen in (3.7),

(3.5b)

= −v000 (x; µ), = −v100 (x; µ), = −v200 (x; µ),

v1 (1; µ) chosen in (3.8), v2 (1; µ) = 0, v3 (0; ε, µ) = v3 (1; ε, µ) = 0,

(3.5c) (3.5d) (3.5e)

with Lµ z ≡ µaz 0 −bz. We see that v(0; ε, µ) = v0 (0; µ)+εv1 (0; µ)+ε2 v2 (0, µ). The following lemmas establish that when v0 (1; µ) and v1 (1; µ) are chosen correctly the first three derivatives of v0 (x; µ) and the first derivative of v1 (x; µ) are bounded independent of µ. Lemma 4. If v0 satisfies the first order differential equation (3.5b) then there exists a value for v0 (1; µ) such that the following bounds hold µ ¶ 1 (i) kv0 k ≤ C 1 + i−3 , for 0 ≤ i ≤ 7. µ Proof. We start by noting that since a > 0 and b > 0 we can establish the following ¯ ¯ ¯ ¯ ¯ If Lµ z ¯ ≤0 and z(1) ≥ 0, then z ¯¯ ≥ 0, [0,1)

(3.6)

[0,1]

using a simple proof by contradiction argument. We further decompose v0 (x; µ) as follows v0 (x; µ) = s0 (x) + µs1 (x) + µ2 s2 (x) + µ3 s3 (x; µ),

(3.7a)

where f as0 (x) s0 (x) = − , s1 (x) = 0 , b b Lµ s3 (x; µ) = −as02 (x),

as01 (x) , b s3 (1; µ) = 0.

s2 (x) =

(3.7b) (3.7c)

We see that v0 (1; µ) = s0 (1)+µs1 (1)+µ2 s2 (1) and assuming sufficient smoothness of the coefficients, we have (i) (i) (i) ks0 k ≤ C, ks1 k ≤ C and ks2 k ≤ C for 0 ≤ i ≤ 3. Using (3.6) we deduce that ks3 k ≤ (1/β)kas02 k ≤ C and from (3.7c) we obtain (i)

ks3 k ≤

C , µi

for

1 ≤ i ≤ 3. (i)

We use these bounds for s0 (x), s1 (x), s2 (x) and s3 (x; µ) to obtain ||v0 || ≤ C for 0 ≤ i ≤ 3 and then (3.5b) to obtain the required result on the higher derivatives of v0 . Lemma 5. If v1 satisfies the first order differential equation (3.5c) then there exists a value for v1 (1; µ) such that the following bounds hold ¶ µ 1 (i) kv1 k ≤ C 1 + i−1 , for 0 ≤ i ≤ 5. µ 5

Proof. We decompose v1 (x; µ) as follows v1 (x; µ) = ρ0 (x) + µρ1 (x) + µ2 ρ2 (x; µ)

(3.8a)

where v000 , b Lµ ρ2 (x; µ) = −aρ01 (x),

aρ00 (x) , b ρ2 (1; µ) = 0.

ρ0 (x) =

ρ1 (x) =

(3.8b) (3.8c)

We see that v1 (1; µ) = ρ0 (1) + µρ1 (1) and assuming sufficient smoothness of the coefficients, we have ¶ µ 1 C (i) (i) kρ0 k ≤ C 1 + i−1 and kρ1 k ≤ i , for 0 ≤ i ≤ 2. µ µ Using (3.6) and (3.8c) we can also obtain (i)

kρ2 k ≤

C , µi+1

for

0 ≤ i ≤ 2. (i)

We use these bounds for ρ0 (x), ρ1 (x), and ρ2 (x; µ) and their derivatives to obtain ||v1 || ≤ C(1 + µ1−i ), for i = 0, 1, 2. The required result for 0 ≤ i ≤ 5 follows by differentiating the differential equation (3.5c) for v1 . Lemma 6. If v2 satisfies the first order differential equation (3.5d) then the following bounds hold (i)

kv2 k ≤

C , µi+1

for

0 ≤ i ≤ 4.

Proof. The proof follows using (3.6), the differential equation (3.5d) and the bounds in Lemma 5. Lemma 7. If v3 satisfies the first order differential equation (3.5e) then the following bounds hold ³ µ ´i µ 1 ¶ (i) kv3 k ≤ C , for 0 ≤ i ≤ 4. ε µ3 Proof. Using the minimum principle for Lε,µ and a suitable barrier function we obtain (see [9] Lemma 2.1), © ª 1 kv3 k ≤ max |v3 (0)|, |v3 (1)| + kv200 k. β Applying the bounds in Lemma 6 we therefore have kv3 k ≤

C . µ3

(3.9)

Using the differential equation (3.5e) and the Mean Value Theorem on an interval of width obtain (see [9] Lemma 2.2), (k)

kv3 k ≤

³ µ ´k ´ © ª C ³ √ k 1+ √ max ||v3 ||, ||v200 || , ( ε) ε

k = 1, 2.

Simplifying this expression using the previous Lemma and (3.9) (k)

kv3 k

≤

³ µ ´k ´ 1 C ³ √ k 1+ √ , µ3 ( ε) ε 6

k = 1, 2.

√

ε we

Differentiating the differential equation for v3 and applying these bounds gives us (3)

kv3 k ≤

C ε3

and

(4)

kv3 k ≤

Cµ . ε4

(3.10)

Substituting all of these bounds for √ v0 (x; µ), v1 (x; ε, µ), v2 (x; ε, µ) and v3 (x; ε, µ) into the equation for v(x; ε, µ) and noting that µ ≥ C ε gives us à µ ¶(3−i) ! ε (i) kv k ≤ C 1 + , for 0 ≤ i ≤ 4. (3.11) µ In order to fully specify the decomposition given in (3.2), it is required to specify the values of v(0), v(1) and wR (0). In the case of αµ2 ≤ γε, v(0), v(1) are specified in (3.3) and wR (0) = 0. In the case of αµ2 ≥ γε v(0), v(1) are specified in (3.5), (3.7) and (3.8) and wR (0) is specified below in (3.12). Lemma 8. When µ2 ≥

γε α ,

the value wR (0) can be specified so that wR satisfies the following bounds (i)

||wR || ≤

C µi

0 ≤ i ≤ 3.

Proof. Consider the following decomposition wR (x; ε, µ) = w0 (x; µ) + εw1 (x; µ) + ε2 w2 (x; µ) + ε3 w3 (x; ε, µ)

(3.12a)

where v(1) = v0 (1) + εv1 (1) is given in (3.7) and Lµ w0 εLµ w1 εLµ w2 εLε,µ w3

= = = =

0, w0 (1; µ) = u(1) − v(1), (Lµ − Lε,µ )w0 , w1 (1; µ) = 0, (Lµ − Lε,µ )w1 , w2 (1; µ) = 0, (Lµ − Lε,µ )w2 , w3 (0; ε, µ) = w3 (1; ε, µ) = 0.

(3.12b) (3.12c) (3.12d) (3.12e)

We start by analysing w0 (x). Using (3.6) applied to w0 (x) = w0 (1)ψ(x), where Lµ ψ = 0, ψ(1) = 1, we obtain the following bounds C (i) ||w0 || ≤ i , 0 ≤ i ≤ 5. (3.13) µ Using this method again for w1 (x) and w2 (x) we obtain (i)

||w1 || ≤

C , µi+2

0 ≤ i ≤ 4,

(i)

and

||w2 || ≤

C , µi+4

Finally we consider w3 , we can apply Lemma 2 to obtain ||w3 || ≤

C . µ6

From Theorem 3 we have the following bounds for 1 ≤ i ≤ 2 ³ µ ´i ´ 1 C ³ (i) ||w3 || ≤ √ i 1 + √ . µ6 ( ε) ε Finally we obtain

³ µ ´3 ´ 1 1 1 C ³ (3) + . ||w3 || ≤ √ 3 1 + √ 6 µ ε µ7 ( ε) ε

The required bounds follow using (3.12) and the inequality µ2 ≥ 7

γε α .

0 ≤ i ≤ 3.

(3.14)

Lemma 9. [10] The singular components wL and wR satisfy the following bounds |wL (x)| ≤ Ce−θ1 x , where

(

√ γα √ , ε αµ ε ,

θ1 =

µ2 ≤ µ2 ≥

|wR (x)| ≤ Ce−θ2 (1−x) , (

γε α , γε α ,

θ2 =

√ γα √ , 2 ε γ 2µ ,

µ2 ≤ µ2 ≥

γε α , γε α .

Note that the derivatives of the singular components wL and wR satisfy the bounds given in Theorem 3.

4

Discrete Problem

To approximate the solution of problem (2.1), we employ a finite difference scheme defined on a piecewise uniform Shishkin mesh [5]. This mesh is defined as follows. Let N be a positive integer and a multiple of 4. Divide the interval Ω = [0, 1] into three subintervals [0, σ1 ], [σ1 , 1 − σ2 ] and [1 − σ2 , 1], where the transition parameters are given by   o o n n √ √ 4 ε 4 ε γε 2  min 1 , √  min 1 , √ ln N , if µ ≤ , ln N , if µ2 ≤ γε α α , n 4 γα o n 4 γα o σ1 = σ = 2 γε 4µ γε 1 1 4ε 2 2  min ,  min , ln N , if µ ≥ α , if µ ≥ α . 4 µα ln N , 4 γ Place N/4 + 1, N/2 + 1 and N/4 + 1 mesh points respectively in [0, σ1 ], [σ1 , 1 − σ2 ] and [1 − σ2 , 1]. Denote the step sizes in each subinterval by H1 = 4σ1 /N , H2 = 2(1 − σ1 − σ2 )/N and H3 = 4σ2 /N , respectively. The mesh points are given by  if 0 ≤ i ≤ N/4,  iH1 , σ1 + (i − N/4)H2 , if N/4 ≤ i ≤ 3N/4, (4.16) xi =  1 − σ2 + (i − 3N/4)H3 , if 3N/4 ≤ i ≤ N. We set hj+1 = xj+1 − xj for j = 0, . . . , N − 1, and ~j = (hj + hj+1 )/2 for j = 1, . . . , N − 1. On this mesh, we define a finite difference scheme that uses the central difference, upwind and midpoint schemes. This scheme is given by LN Uj ≡ rj− Uj−1 + rjc Uj + rj+ Uj+1 = QN (fj ), where

LN LN LN LN LN LN

≡ LN cd , ≡ LN cd , ≡ LN mp , ≡ LN up , ≡ LN cd , ≡ LN mp ,

if if if if if if

1 ≤ j < N/4, N/4 < j < 3N/4, N/4 < j < 3N/4, N/4 < j < 3N/4, 3N/4 < j ≤ N − 1, 3N/4 < j ≤ N − 1,

and and and and and

µH2 kak < 2ε, µH2 kak ≥ 2ε, kbkH2 < 2µα, µH2 kak ≥ 2ε, kbkH2 ≥ 2µα, µH3 kak < 2ε, µH3 kak ≥ 2ε.

At the left transition point σ1 we use LN ≡ LN if xj = σ1 = 0.25, cd , if xj = σ1 < 0.25, and kbkH2 < 2µα,

LN ≡ LN mp ,

LN ≡ LN up ,

otherwise.

At the right transition point 1 − σ2 we use LN ≡ LN cd ,

if xj = 1 − σ2 = 0.75,

µH3 kak < 2ε,

LN ≡ LN mp ,

if xj = 1 − σ2 = 0.75,

µH3 kak ≥ 2ε,

LN ≡ LN mp ,

if xj = 1 − σ2 > 0.75, and kbkH3 < 2µα, LN ≡ LN up , 8

otherwise.

(4.17)

Associated with each of these finite difference operators, we have the following finite difference scheme 2 0 LN cd Uj ≡ εδ Uj + µaj D Uj − bj Uj = fj , N 2 Lup Uj ≡ εδ Uj + µaj D+ Uj − bj Uj = fj , 2 LN aj D+ Uj − bj Uj = f¯j , mp Uj ≡ εδ Uj + µ¯ with z¯j = (zj + zj+1 )/2, and µ ¶ Uj − Uj−1 1 Uj+1 − Uj 2 δ Uj = − , ~j hj+1 hj

D+ Uj =

Uj+1 − Uj , hj+1

D0 Uj =

Uj+1 − Uj−1 . hj + hj+1

Thus the discrete problem is: LN U = QN (fj ), xi ∈ ΩN , U (0) = u(0), U (1) = u(1).

(4.18a) (4.18b)

Note that the elements in the system matrix LN are as follows ε µaj + ε µaj c − , r = + , r = −rj+ − rj− − bj , if LN ≡ LN cd , hj ~j 2~j j hj+1 ~j 2~j j µaj ε ε , r+ = + , rc = −rj+ − rj− − bj , if LN ≡ LN rj− = up , hj ~j j hj+1 ~j hj+1 j µ¯ aj bj+1 c ε ε , r+ = + − , rj = −rj+ − rj− − ¯bj , if LN ≡ LN rj− = mp , hj ~j j hj+1 ~j hj+1 2 rj− =

and N QN (fj ) = fj , if LN ≡ LN cd or Lup , QN (fj ) = f¯j , if LN ≡ LN mp .

In the left layer region (0, σ1 ), note that µkakH1 2µkakσ1 8kak ln N = ≤ . 2ε εN αN In the right-hand layer region (1 − σ2 , 1), note that µkakH3 2µkakσ2 8kak ln N = ≤ , 2ε εN αN kbkH3 2kbkσ2 8kbk ln N = ≤ , 2αµ αµN αγN

for

αµ2 ≤ γε.

for

αµ2 ≥ γε.

To guarantee a monotone difference operator LN , we impose the following mild assumption on the minimum number of mesh points ¾ ½ kak kbk −1 , . (4.19) N (ln N ) > 8 max α αγ Lemma 10. Assume (4.19). The finite difference scheme (4.17) satisfies a discrete minimum principle of the form: For any mesh function Z: If Z(0) ≥ 0, Z(1) ≥ 0, and LN Z ≤ 0,

then Z(xi ) ≥ 0, ∀xi .

Also, for any mesh function Z kZk ≤

1 N kL Zk + max{|Z(0)|, |Z(1)|}. β 9

(4.20)

Proof. We will simply check that the choice of the different discretizations used in the definition of the difference scheme (4.17) allow us to establish the following inequalities on the entries in the system matrix rj− > 0, rj+ > 0, rj− + rjc + rj+ < 0. (4.21) From these sign patterns it follows that the matrix LN is the negative of an M -matrix. In the case of the upwind operator LN up , note that the conditions (4.21) are automatically satisfied. In the case of the central difference operator LN cd , the conditions (4.21) are satisfied if ~j rj− =

ε µkak − > 0. hj 2

For the mesh points outside the left layer region xi > σ1 , the sign pattern of rj− > 0 is guaranteed as LN cd is only used when µH2 kak < 2ε or µH3 kak < 2ε. Using the assumption (4.19), one can check that rj− > 0 for xi ≤ σ1 . In the case of the mid-point operator LN mp , the conditions (4.21) are satisfied if b(xj+1 )hj+1 hj+1 rj+ > µ¯ aj − > 0. 2 The condition kbkH2 < 2µα (when the mid-point operator is used outside the layers) ensures that rj+ > 0 for xj ∈ [σ1 , 1 − σ2 ). In the right-hand layer region [1 − σ2 , 1), using (4.19), we can establish that rj+ > 0 for xj ∈ [1 − σ2 , 1). We can deduce the following truncation error bounds for two of the three different difference operators employed in LN : On an arbitrary mesh, we have that (3) k + µ~i kak ku(2) k, k(LN cd − L)uk ≤ ε~i ku (3) k + µhi+1 kak ku(2) k, k(LN up − L)uk ≤ ε~i ku

and on a uniform mesh with step size h 2 (4) k(LN k + µh2 kak ku(3) k, cd − L)uk ≤ εh ku 2 (4) k(LN k + µhkak ku(2) k. up − L)uk ≤ εh ku

Note the increase in the order of the truncation error in the case of a uniform mesh. Let us finally look at the mid-point scheme µ 00 ¶ 00 ¯j − LN u = ε u (xj+1 ) − u (xj ) + u00 (xj ) − δ 2 u(xj ) + LN (U − u) = f mp mp 2 +0.5µ(a(xj )u0ε (xj ) − a(xj )D+ uj + a(xj+1 )u0ε (xj+1 ) − a(xj+1 )D+ uj ). Hence, we have the following bound on a non-uniform mesh (3) k(LN k + C(kak,ka0 k) µh2i+1 (ku(3) k + ku(2) k), mp − L)uk ≤ ε~i ku

and on a uniform mesh there is no increase in the order as (3) k(LN k + C(kak,ka0 k) µh2 (ku(3) k + ku(2) k), mp − L)uk ≤ εhku

with C(kak,ka0 k) a positive constant that depends on kak and ka0 k. Note that the mid-point scheme has a higher order of truncation error in the convective (µau0 ) term than the centered difference operator on a non-uniform mesh. This is why the mid-point scheme is chosen at the transition points. 10

5

Error analysis

The solution of the discrete problem (4.18) can be decomposed into the following sum U = V + WL + WR ,

(5.22a)

where LN V = LN U, LN WL = 0, LN WR = 0,

V (0) = v(0), V (1) = v(1), WL (0) = wL (0), WL (1) = 0, WR (0) = wR (0), WR (1) = wR (1).

(5.22b) (5.22c) (5.22d)

Lemma 11. Assume (4.19). The regular component of the error satisfies the following error bound kv − V k ≤ CN −2 ,

(5.23)

where v is the solution of (3.2b) and V is the solution of (5.22b). Proof. If the mesh is uniform (16σ1 σ2 = 1) , then |LN (v − V )(xj )| ≤ CN −2 , and the result follows from the uniform stability given in (4.20). In the case when the mesh is not uniform (16σ1 σ2 < 1), away from the transition points the mesh is uniform and in all cases one can check that |LN (v − V )(xj )| ≤ CN −2 , xj 6= σ1 , 1 − σ2 . If the mesh is nonuniform, then at the transition points we either employ the mid-point scheme or if 2µα ≤ kbkhj+1 we use the upwind scheme. Either way, we have the following truncation error bound ½ CN −2 , if j 6= N/4, 3N/4, N |L (v − V )(xj )| ≤ CN −1 (ε + N −1 ), otherwise. Define the barrier function where

ψj = CN −2 (θ(zj ) + 1),

 1, if 0 ≤ z ≤ σ1 ,    z − σ1  1− , if σ1 ≤ z ≤ 1 − σ2 , θ(zj ) = 2(1 − σ1 − σ2 )  1−z    , if 1 − σ2 ≤ z ≤ 1. 2σ2

Noting that 1/σ2 ≥ 4, we have that ½ 0, j 6= N/4, 3N/4, 2 εδ ψj = O(−N −1 ε), otherwise,

D0 ψj ≤ 0,

D+ ψj ≤ 0.

Combining this with (4.20) completes the proof. Motivated by the pointwise bounds on the continuous singular components wL , wR , consider the following barrier functions  j  N    Y  Y (1 + θR hi )−1 , 0 ≤ j < N, (1 + θL hi )−1 , 1 ≤ j ≤ N, ΨL,j = ΨR,j =    i=1  i=j+1 1, j = 0, 1, j = N, 11

where

 √ γα γε   √ , , if µ2 ≤ α 2 ε θL = γε   µα , if µ2 ≥ , 2ε α First we prove the following technical result.

 √ γα   √ , 2 ε θR = γ   , 2µ

γε , α γε if µ2 ≥ . α if µ2 ≤

Lemma 12. The barrier functions ΨL,j , ΨR,j satisfy the inequalities LN ΨL,j ≤ 0,

LN ΨR,j ≤ 0.

Proof. We begin by applying the operator to the left-hand barrier function ΨL,j µ ¶ 1 LN ΨL,j = ΨL,j (1 + θL hj )rj− + rjc + rj+ = 1 + θL hj+1 Ã Ã !! hj+1 rj+ − + − c = ΨL,j rj + rj + rj + θL hj rj − . 1 + θL hj+1 Now we analyze each of the different discretizations used in the definition of the operator LN . First, in the case of the central difference operator we have 2 LN cd ΨL,j ≤ ΨL,j+1 (2εθL − µaj θL − bj ). √ √ If µ2 ≤ γε/α, let θL = γα/(2 ε) to prove that ³ γα ´ ³ γα ´ − bj − µaj θL ≤ ΨL,j+1 − bj ≤ 0. LN cd ΨL,j ≤ ΨL,j+1 2 2

If µ2 ≥ γε/α, use θL = µα/2ε to get µ LN cd ΨL,j ≤ ΨL,j+1

µ2 α (α − aj ) − bj 2ε

¶ ≤ 0.

Now consider the upwind scheme. Then we obtain 2 LN up ΨL,j ≤ ΨL,j+1 (2εθL − µaj θL − bj ),

and we argue as for LN cd . Finally for the midpoint scheme we can deduce that 2 LN aj θL − bj /2), mp ΨL,j ≤ ΨL,j+1 (2εθL − µ¯

and we finish as for LN cd . Applying the operator to the right hand barrier function we obtain à à N

L ΨR,j = ΨR,j

rj−

+

rjc

+

rj+

+ θR

hj+1 rj+

hj rj− − 1 + θR hj

!! .

If we use the central difference or the upwind scheme, we have ¡ 2 ¢ LN ΨR,j ≤ ΨR,j 2εθR + µaj θR − bj . Now it is not necessary to consider two cases depending on the relation between ε and µ. In both cases, we obtain µ ¶ ´ ³ γα a γ bj j N + − bj ≤ ΨR,j (γaj − bj ) ≤ ΨR,j aj γ − ≤ 0, L ΨR,j ≤ ΨR,j 2 2 aj 12

since γ = min (b/a). If we use the midpoint scheme, then µ µ ¶ µ ¶¶ ¡ 2 ¢ aj bj aj+1 bj+1 N ¯ L ΨR,j ≤ ΨR,j 2εθR + µ¯ aj θR − bj ≤ ΨR,j γ− + γ− ≤ 0, 2 aj 2 aj+1 and this completes the proof. An immediate consequence of the previous lemma and the minimum principle is the following result, which implies that the discrete singular components WL , WR are computationally small (≤ CN −2 ) outside their respective layer regions. Corollary 13. The barrier functions ΨL,j and ΨR,j and the singular components WL and WR satisfy |WL,j | ≤ |WL,0 | ΨL,j , |WR,j | ≤ |WR,N | ΨR,j . Moreover, it holds

ΨL,j ≤ CN −2 , ΨR,j ≤ CN −2 ,

if if

N/4 ≤ j ≤ N, 0 ≤ j ≤ 3N/4.

(5.24)

Proof. The bounds of (5.24) are obtained in a standard way. For j ≥ N/4, the component WL satisfies µ ¶−N/4 4σ1 ΨL,j ≤ ΨL,N/4 = 1 + θL = ((1 + 8N −1 ln N )−N/8 )2 ≤ CN −2 , N where we have used the inequality ln(1 + t) > t(1 − t/2) to deduce (1 + 8N −1 ln N )−N/8 ≤ 8N −1 . We use the same argument to bound the barrier function ΨR,j . Lemma 14. Assume (4.19). The left singular component of the error satisfies the following estimate ½ CN −2 ln3 N, if γε ≤ αµ2 , (5.25) kwL − WL k ≤ −1 2 C(N ln N ) , if γε ≥ αµ2 , where wL is the solution of (3.2c) and WL is the solution of (5.22c). (k)

2 −k/2 Proof. and √ First consider the uniform mesh case of σ1 = 1/4. If µ ≤ γε/α, using kwL k ≤ Cε 1/ ε ≤ C ln N , we deduce that (4)

(3)

|LN (wL − WL )(xj )| ≤ CN −2 (εkwL k + µkwL k) ≤ CN −2 /ε ≤ C(N −1 ln N )2 . (k)

If µ2 ≥ γε/α, using kwL k ≤ C(µ/ε)k and µ/ε ≤ C ln N , we deduce that (4)

(3)

|LN (wL − WL )(xj )| ≤ CN −2 (εkwL k + µkwL k) ≤ CN −2 µ4 /ε3 ≤ CµN −2 ln3 N. To complete the proof in the case of a non-uniform mesh, we split the argument into two cases depending of the localization of the mesh point. In the first case xj ∈ [σ1 , 1). The triangular inequality, lemma 9 and (5.24) yield |(wL − WL )(xj )| ≤ |wL (xj )| + |WL (xj )| ≤ C(e−θ1 xj + N −2 ) ≤ C(e−θ1 σ1 + N −2 ). Using the fact that θ1 σ1 = 4 ln N , we conclude that |(wL − WL )(xj )| ≤ CN −2 ,

if xj ∈ [σ1 , 1).

If xj ∈ (0, σ1 ), we calculate a bound on the truncation error of the form |LN (wL − WL )(xj )|

= ≤

(4)

(3)

2 2 |LN cd (wL − WL )(xj )| ≤ C(εH1 kwL k + µH1 kwL k) ≤ (4) (3) −2 2 2 CN (εσ1 kwL k + µσ1 kwL k).

13

Note that this is a second order truncation error bound because the scheme is the √ central difference (k) operator on a uniform mesh. If µ2 ≤ γε/α, from kwL k ≤ Cε−k/2 and σ1 = O( ε ln N ), it follows that µ ¶ µ −1 2 √ 1 + |LN (w − W )(x )| ≤ C(N ln N ) ≤ C(N −1 ln N )2 . L L j cd ε (k)

If µ2 ≥ γε/α, from kwL k ≤ C( µε )k and σ1 = O(ε/µ ln N ), we conclude that |LN cd (wL − WL )(xj )| ≤ C

µ2 −1 (N ln N )2 . ε

We use this bound to obtain an appropriate bound on the error in the layer region (0, σ1 ). Consider the barrier function ³ µ´ ψj = C N −2 + (N −1 ln N )2 (σ1 − xi ) . ε Taking into account that δ 2 ψj = 0 and D0 ψj < 0, it follows that LN cd ψj

µ ¶ µ2 −1 2 = − bj ψj + C (N ln N ) ≤ −|LN cd (wL − WL )(xj )|. ε

Therefore the discrete minimum principle on [0, σ1 ], yields ³ µ´ |(wL − WL )(xj )| ≤ ψj ≤ C N −2 + (N −1 ln N )2 σ1 ≤ CN −2 ln3 N. ε

Lemma 15. Assume (4.19). The right singular component of the error satisfies the following estimate kwR − WR k ≤ C(N −1 ln N )2 , (5.26) where wR is the solution of (3.2d) and WR is the solution of (5.22d). Proof. The proof is similar to the proof given for the left boundary layer component. If xj ∈ (0, 1 − σ2 ], again we have |(wR − WR )(xj )| ≤ |wR (xj )| + |WR (xj )| ≤ C(e−θ2 σ2 + N −2 ) ≤ CN −2 . If xj ∈ (1 − σ2 , 1), the truncation error associated with the central differences satisfies (4)

(3)

(4)

(3)

2 2 −2 2 |LN σ2 (εkwR k + µkwR k), cd (wR − WR )(xj )| ≤ C(εH3 kwR k + µH3 kwR k) ≤ CN

and the truncation error associated with the midpoint scheme satisfies (3)

(3)

(2)

2 |LN mp (wR − WR )(xj )| ≤ C(εH3 kwR k + µH3 (kwR k + kwR k)) ≤ (3) (2) (3) (2) ≤ C(µH32 kwR k + µH32 kwR k) ≤ CN −2 µσ22 (kwR k + kwR k),

since the midpoint scheme is activated when µH3 kak ≥ 2ε. We begin the proof with the uniform mesh case, where σ2 = 1/4. If µ2 ≤ γε/α, then using √ (k) 1/ ε ≤ C ln N and the bounds kwR k ≤ Cε−k/2 , we have that (4)

(3)

−2 |LN (εkwR k + µkwR k) ≤ CN −2 /ε ≤ C(N −1 ln N )2 , cd (wR − WR )(xj )| ≤ CN (3) (2) N −2 |Lmp (wR − WR )(xj )| ≤ CN µ(kwR k + kwR k) ≤ CN −2 µ/ε3/2 ≤ C(N −1 ln N )2 .

14

(k)

(4)

0 If µ2 ≥ γε/α, using the facts that kwR k ≤ Cµ−k , k ≤ 3 , εwR = (−µwR +bwr )00 and 1/µ ≤ C ln N , we deduce that (4)

(3)

−2 |LN (εkwR k + µkwR k) ≤ CN −2 /µ2 ≤ C(N −1 ln N )2 , cd (wR − WR )(xj )| ≤ CN (3) (2) N −2 |Lmp (wR − WR )(xj )| ≤ CN µ(kwR k + kwR k) ≤ CN −2 /µ2 ≤ C(N −1 ln N )2 .

Now consider the non-uniform mesh case of σ2 < 1/4. We begin with the central difference scheme. √ (k) If µ2 ≤ γε/α, from kwR k ≤ Cε−k/2 and σ2 = O( ε ln N ), it follows that µ ¶ µ N −1 2 |Lcd (wR − WR )(xj )| ≤ C(N ln N ) 1 + √ ≤ C(N −1 ln N )2 . ε If µ2 ≥ γε/α, from kw(k) k ≤ Cµ−k and σ2 = O(µ ln N ), we have that µ ¶ ε N −1 2 |Lcd (wR − WR )(xj )| ≤ C(N ln N ) 1 + 2 ≤ C(N −1 ln N )2 . µ To finish the proof we consider the midpoint scheme. The midpoint scheme is not used in the γε 2 case µ2 ≤ γε α under the assumption (4.19). Hence, we consider only the case of µ ≥ α . From (k) (4) 0 kwR k ≤ Cµ−k , k ≤ 3 , εwR = (−µwR + bwr )00 and σ2 = O(µ ln N ) we get that (3)

(2)

−2 |LN µσ22 (kwR k + kwR k) ≤ C(N −1 ln N )2 . mp (wR − WR )(xj )| ≤ CN

Thus, in all cases,

|LN (wR − WR )(xj )| ≤ C(N −1 ln N )2 .

Complete the proof using (4.20). Theorem 16. Assume that N is sufficiently large so that (4.19) is satisfied. The maximum pointwise error satisfies the following parameter-uniform error bound ½ CN −2 ln3 N, if γε ≤ αµ2 , (5.27) ku − U k ≤ −1 2 C(N ln N ) , if γε ≥ αµ2 , where u is the solution of the continuous problem (2.1) and U is the solution of the discrete problem (4.18). Proof. The result follows from the triangular inequality, Lemmas 11, 14 and 15.

6

Numerical experiments

In this final section, we apply the numerical method (4.18) to the following constant coefficients problem εu00 (x) + µu0 (x) − u(x) = −x, x ∈ Ω = (0, 1), (6.28) whose exact solution and boundary conditions are given by u(x) = (x + µ) + where

(1 + µ)em2 /(2ε) + 1 − µ −xm1 /(2ε) 1 + µ + (1 − µ)e−m1 /(2ε) (1−x)m2 /(2ε) e − e , 1 − e−m1 /ε 1 − e−m1 /ε m1,2 = µ2 ±

p

µ2 + 4ε.

Figures 1–2 display solutions of this problem (6.28) for several values of the singular perturbation parameters ε and µ. 15

1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

(a) ε = 10

0.3

−2

0.4

0.5

0.6

0.7

and µ = 10

0.8

0.9

0

1

−1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

(b) ε = 10−2 and µ = 10−7

Figure 1: Exact solutions of problem (6.28) for different values of ε and µ . 1

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6

0.6

0.5

0.5

0.4

0.4

0.3

0.3

0.2

0.2

0.1

0.1

0

0

0.1

0.2

(a) ε = 10

0.3

−3

0.4

0.5

0.6

and µ = 10

0.7

0.8

0.9

0

1

−1

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

(b) ε = 10−3 and µ = 10−7

Figure 2: Exact solutions of problem (6.28) for different values of ε and µ . N For any value of N , the maximum pointwise errors Eε,µ , the ε–uniform errors EµN and the N (ε, µ)–uniform errors E are calculated using N Eε,µ = ku − U kΩN ,

N EµN = max Eε,µ , ε

N E N = max Eε,µ ε,µ

(6.29)

where u is the exact solution of (6.28) and U is the numerical solution of the finite difference scheme (4.18). From these values the local orders of convergence pN ε,µ , the local order of ε-uniform convergence N N pµ and the local order of (ε, µ)–uniform convergence p are calculated using pN ε,µ = log2

N Eε,µ , 2N Eε,µ

pN µ = log2

EµN , Eµ2N

pN = log2

EN . E 2N

(6.30)

In Tables 1 we display the (ε, µ)–uniform errors E N and the corresponding orders pN for N = 26 , 28 , . . . , 215 and N max max Eε,µ , µ∈A ε∈C1

with A = {20 , 2−1 . . . , 2−20 }, Bλ = λ × {20 , 2−1 , . . . , 2−60 } and Cλ = {ε ∈ Bλ |2−20 µ2 ≤ ε}. We consider this set of values in order to elude the precision error of the computer (all calculations N are in double precision) and to ensure that for each value of µ the Eε,µ stabilize (smaller values of ε 16

N do not change Eε,µ ) for all values of N . For example, we observed that for µ = 1 it is necessary to consider values of ε ≤ 2−17 and for µ = 2−20 it is necessary to consider ε ≤ 2−60 before stabilization with respect to ε is observed. If ε is too small (e.g., ε = 2−60 , µ = 1) then rounding errors pollute the calculation of the order of parameter-uniform convergence. N Also, a surface plot of the maximum pointwise errors Eε,µ is displayed in Figure 3 for ε = −26 −28 −46 −13 −14 −25 N 2 ,2 ,...,2 , µ = 2 ,2 ,...,2 and N = 64. Observe that as ε → 0 the errors Eε,µ begin increasing, then decrease and finally stabilize at a constant value.

0.045

Maximum errors

0.04 0.035 0.03 0.025 0.02 0.015 12 10 8 6

−5

x 10

4 2

2

4

Mu

6

8

10

12

14

−9

x 10

Epsilon

N Figure 3: Surface plot of the maximum pointwise errors Eε,µ for N = 64 within the range 0 < µ < −4 −9 1.2 × 10 and 0 < ε < 1.4 × 10

Table 1: The (ε, µ)–uniform errors E N and the (ε, µ)–uniform orders pN for various values of N computed over the parameter range µ = 2−p , ε = 2−q and 2−20 ≤ µ ≤ 1, 2−20 µ2 ≤ ε ≤ 1. EN pN

N=64

N=128

N=256

N=512

N=1024

N=2048

N=4096

N=8192

N= 16384

N= 32768

4.577E-2 1.180

2.021E-2 1.370

7.817E-3 1.523

2.719E-3 1.639

8.728E-4 1.723

2.643E-4 1.159

1.184E-4 1.809

3.380E-5 1.831

9.499E-6 1.865

2.607E-6

The results are in line with the theoretical error bound given in Theorem 16, although we observe that the (ε, µ)–uniform orders of convergence are not always increasing. To understand this anomalous behaviour, we display in Table 2 for µ = 2−10 the ε–uniform errors for ε ≤ µ2 . N In Table 2, the maximum error Eε,µ=2 −10 lies along two distinct diagonals. The corresponding N orders of convergence pε,µ=2−10 for these errors are bigger than two due to the fact that the midpoint scheme is being replaced by the central difference scheme. Below the diagonals the operator is the midpoint (the orders of convergence are decreasing because there are related to the bound N −1 ε if ε/µ2 ≤ CN −1 ) and central difference scheme above the diagonal (the orders of convergence are seen to be increasing). Nevertheless, we observe a jump in the ε–uniform orders of convergence for N = 2048. The reason for this jump is that we can only consider a discrete set of values for the parameters ε and µ. In Table 3 we show the results when we compute over a more detailed range of ε, ε = {0.6, 0.7, 0.8, 0.9, 1} × {2−21 , . . . , 2−30 } and µ = 2−10 . The jump in the ε–uniform orders has been moved to N = 8192 and damped for this larger set of parameter values. In Table 4 we display the 17

(ε, µ)–uniform errors when the singular perturbation parameters take the values µ∈A

and

[

ε∈

Cλ .

(6.31)

λ=0.6,0.7,0.8,0.9,1

where A = {20 , 2−1 . . . , 2−20 }, Bλ = λ × {20 , 2−1 , . . . , 2−60 } and Cλ = {ε ∈ Bλ |2−20 µ2 ≤ ε}. Again we observe that if we take a larger set of parameter values, the computed orders of convergence are closer to the order of convergence established theoretically in this paper. N N Table 2: The maximum pointwise errors Eε,µ=2 −10 and the orders of convergence pε,µ=2−10 for several values of ε.

ε

N=64

N=128

N=256

N=512

N=1024

N=2048

N=4096

N=8192

N= 16384

N= 32768

2−21

4.546E-2 2.510 4.538E-2 1.177 4.336E-2 1.325 4.120E-2 1.462 3.966E-2 1.565 3.873E-2 1.632 3.822E-2 1.671 3.795E-2 1.691 3.782E-2 1.702 3.775E-2 1.708

7.981E-3 1.586 2.007E-2 2.391 1.730E-2 1.156 1.496E-2 1.280 1.340E-2 1.401 1.250E-2 1.493 1.201E-2 1.552 1.175E-2 1.586 1.162E-2 1.604 1.156E-2 1.613

2.659E-3 1.646 3.827E-3 1.626 7.764E-3 2.271 6.158E-3 1.189 5.075E-3 1.328 4.441E-3 1.454 4.096E-3 1.546 3.916E-3 1.604 3.824E-3 1.637 3.778E-3 1.654

8.497E-4 1.692 1.240E-3 1.688 1.609E-3 1.686 2.701E-3 2.197 2.021E-3 1.221 1.621E-3 1.367 1.402E-3 1.492 1.288E-3 1.582 1.230E-3 1.636 1.200E-3 1.667

2.630E-4 1.723 3.847E-4 1.723 5.001E-4 1.722 5.889E-4 1.721 8.669E-4 2.144 6.284E-4 1.259 4.984E-4 1.406 4.304E-4 1.530 3.956E-4 1.615 3.780E-4 1.667

7.967E-5 1.748 1.165E-4 1.748 1.516E-4 1.748 1.786E-4 1.748 1.961E-4 1.748 2.626E-4 1.153 1.880E-4 1.297 1.491E-4 1.444 1.292E-4 1.564 1.191E-4 1.645

2.371E-5 1.769 3.469E-5 1.769 4.514E-5 1.769 5.318E-5 1.769 5.839E-5 1.769 1.181E-4 2.712 7.652E-5 1.183 5.480E-5 1.331 4.368E-5 1.477 3.806E-5 1.594

6.958E-6 1.786 1.018E-5 1.786 1.325E-5 1.786 1.561E-5 1.786 1.714E-5 1.786 1.802E-5 1.786 3.370E-5 2.651 2.179E-5 1.211 1.569E-5 1.361 1.261E-5 1.506

2.017E-6 1.801 2.952E-6 1.801 3.841E-6 1.801 4.525E-6 1.801 4.969E-6 1.801 5.226E-6 1.801 5.365E-6 1.801 9.414E-6 2.593 6.111E-6 1.236 4.440E-6 1.388

5.789E-7

4.546E-2

2.007E-2

7.764E-3

2.701E-3

8.669E-4

2.626E-4

1.181E-4

3.370E-5

9.414E-6

1.180

1.370

1.523

1.639

1.723

1.153

1.809

1.840

1.859

2−22 2−23 2−24 2−25 2−26 2−27 2−28 2−29 2−30 E N−10 2 pN−10 2

8.471E-7 1.102E-6 1.299E-6 1.426E-6 1.500E-6 1.540E-6 1.560E-6 2.595E-6 1.696E-6 2.595E-6

Table 3: The (ε, µ)–uniform errors E N and the (ε, µ)–uniform orders pN for various values of N computed over the parameter range µ = 2−10 , ε = {0.6, 0.7, 0.8, 0.9, 1.0} × {2−21 , ..., 2−30 }. E N−10 2 pN−10 2

N=64

N=128

N=256

N=512

N=1024

N=2048

N=4096

N=8192

N= 16384

N= 32768

4.580E-2

2.074E-2

8.253E-3

3.154E-3

1.122E-3

3.731E-4

1.181E-4

3.370E-5

1.072E-5

2.952E-6

1.114

1.329

1.388

1.491

1.588

1.660

1.809

1.652

1.861

Table 4: The (ε, µ)–uniform errors E N and the (ε, µ)–uniform orders pN for various values of N computed over the parameter range (6.31) for problem (6.28). EN pN

N=64

N=128

N=256

N=512

N=1024

N=2048

N=4096

N=8192

N= 16384

N= 32768

4.611E-2 1.143

2.088E-2 1.329

8.309E-3 1.388

3.175E-3 1.492

1.129E-3 1.588

3.756E-4 1.666

1.184E-4 1.809

3.380E-5 1.642

1.083E-5 1.820

3.067E-6

As in [5] we have estimated the order of convergence corresponding to the theoretical error bounds (N −1 ln N )2 and N −2 ln3 N . These values are given in Table 5. On comparing these rates with the computed rates given in Table 4, the computed rates appear to be in close agreement with the theoretical rates established in Theorem 16. 18

Table 5: Orders of convergence corresponding to different theoretical error bounds. (N −1 ln N )2 N −2 ln3 N

N=64

N=128

N=256

N=512

N=1024

N=2048

N=4096

N=8192

N= 16384

1.526 1.277

1.593 1.383

1.644 1.461

1.683 1.522

1.715 1.570

1.741 1.609

1.762 1.642

1.780 1.669

1.796 1.693

The second test problem that we consider is the variable coefficients problem ½ εu00 + µ(1 + x)u0 − u = (1 + x)2 , x ∈ (0, 1), u(0) = u(1) = 0.

(6.32)

To estimate the errors we use the following variant of the double mesh principle. Compute N 2N ˜2j Dε,µ = kUjN − U k∞,ΩN ,

N DN = max Dε,µ , ε,µ

˜ 2N }2N is the solution obtained on a mesh containing the mesh points of the original where {U j=0 j Shishkin mesh and its midpoints xj+1/2 = (xj+1 + xj )/2, j = 0, . . . , N − 1. We note that we use the ˜ 2N 2N same difference operators to calculate {UjN }N j=0 and {Uj }j=0 . To estimate the errors we compare the values of both computed solutions on the coarser mesh ΩN without using polynomial interpolation. From these differences we estimate the (ε, µ)–uniform errors using DN and the corresponding (ε, µ)–uniform orders of convergence pN are estimated using q N = log2

DN . D2N

In Table 6 we show the estimatedS(ε, µ)–uniform errors when the singular perturbation parameters take the values µ ∈ A and ε ∈ λ=0.6,0.7,0.8,0.9,1 Cλ . We observe the same behaviour as in the previous test problem. N Table 6: The (ε, µ)–uniform differences Dε,µ and the computed orders of (ε, µ)–uniform convergence N q for problem (6.32). DN qN

N=64

N=128

N=256

N=512

N=1024

N=2048

N=4096

N=8192

N= 16384

N= 32768

1.196E-1 1.303

4.848E-2 1.343

1.911E-2 1.405

7.215E-3 1.499

2.552E-3 1.592

8.467E-4 1.667

2.667E-4 1.798

7.670E-5 1.665

2.418E-5 1.852

6.698E-6

References [1] L. R. Abrahamsson, H. B. Keller and H. O. Kreiss, Difference approximations for singularly perturbations of systems of ordinary differential equations, Numer. Math., 22 (1974) 367-391. [2] V. B. Andr´eiev, N. V. Kopteva, A study of difference schemes with the first derivate approximated by a central difference ratio, Comp. Maths. Math. Phys., 36(8) (1996) 1065-1078. [3] V. B. Andr´eiev, N. V. Kopteva, Uniform with respect to a small parameter convergence of difference schemes for a convection-diffusion problem, in Analytical and Computational Methods for Convection-Dominated and Singularly Perturbed Problems (J. J. H. Miller, G. I. Shishkin & L. Vulkov, eds.), Nova Science Publishers, New York, 2000, 133-139. [4] C. Clavero, J .L. Gracia and F. Lisbona, High order methods on Shishkin meshes for singular perturbation problems of convection-diffusion type, Numer. Algor. 22 (1999) 73-97. 19

[5] P. A. Farrell, A. F. Hegarty, J. J. H Miller, E. O’Riordan and G. I. Shishkin , Robust computational techniques for boundary layers, Applied Mathematics 16 (2000). [6] J. L. Gracia, F. Lisbona and C. Clavero, High order ε-uniform methods for singularly perturbed reaction-diffusion problems, Lecture Notes in Computer Science, 1988 (2001) 350–358. [7] J. J. H. Miller, E. O’Riordan, G. I. Shishkin and L. P. Shishkina, Fitted mesh methods for problems with parabolic boundary layers, Mathematical Proceedings of the Royal Irish Academy, 98A(2) 1998 173-190. [8] R. E. O’Malley Jr, Two-parameter singular perturbation problems for second order eqautions, J. Math. Mech. 16 (1967) 1143-1164. [9] E. O’Riordan, M. L. Pickett and G. I. Shishkin, Singularly perturbed problems modeling reactionconvection-diffusion processes Computational Methods in Applied Mathematics, 3(3) (2003) 424-442. [10] E. O’Riordan, M. L. Pickett and G. I. Shishkin, Parameter-uniform finite difference schemes for singularly perturbed parabolic diffusion-convection-reaction problems , Dublin City University, School of Mathematical Sciences preprint, 2004, MS-04-08. [11] G. I. Shishkin, Discrete approximation of singularly perturbed elliptic and parabolic equations, Russian Academy of Sciences, Ural Section, Ekaterinburg, 1992. [12] G. I. Shishkin and V. A. Titov, A difference scheme for a differential equation with two small parameters at the derivatives, Chisl. Metody Mekh. Sploshn. Sredy, 7(2) (1976) 145-155, in Russian. [13] M. Stynes and H.-G. Roos, The midpoint upwind scheme, Appl. Num. Math., 23 (1997) 361-374.

20