A defect–correction parameter-uniform numerical method for a singularly perturbed convection–diffusion problem in one dimension

A defect–correction parameter-uniform numerical method for a singularly perturbed convection-diffusion problem in one dimension∗ J.L. Gracia†

E. O’Riordan‡

Abstract Monotone second order parameter–robust numerical methods for singularly perturbed differential equations can be designed using the principles of defect–correction. However, the proofs of second order parameter–uniform convergence can be difficult, even in one dimension. In this paper, we examine a variant of the standard defectcorrection scheme. This variant is proposed in order to simplify the analysis of a defect–correction method in the case of a one dimensional convection–diffusion problem. Numerical results are presented to validate the theoretical results. Keywords: Singular perturbation, defect-correction, Shishkin mesh, second order. AMS subject classifications: 65L70, 65L12.

1

Introduction

Singularly perturbed differential equations are often characterized by the presence of a small parameter (denoted here by ε) multiplying the highest order derivatives appearing in the differential equation. Two of the most commonly studied classes of second order singularly perturbed problems are problems of reaction-diffusion type (no first derivatives present) and problems of convection-diffusion type. Classical finite difference or finite element methods constructed on uniform meshes are known to be unsuitable for solving singularly perturbed problems [8, 19]. Their inadequacy is most pronounced in the case of convection-diffusion problems. Large spurious oscillations are observed in the numerical approximations [8] when the standard central difference operator is employed on a uniform mesh. Parameter-uniform numerical methods [8, 17] are methods whose numerical approximations U N satisfy error bounds of the form kuε − U N k∞ ≤ Cp N −p ,

p > 0,

where uε is the solution of the continuous problem, k · k∞ is the maximum pointwise norm, (N +1) (independent of ε) is the number of mesh elements used in each coordinate direction ∗ This research was partially supported by the project MEC/FEDER MTM 2004-019051 and the grant EUROPA XXI of the Caja de Ahorros de la Inmaculada. † Departamento de Matem´ atica Aplicada, Universidad de Zaragoza, Spain. email: [email protected] ‡ School of Mathematical Sciences, Dublin City University, Dublin 9, Ireland. email: [email protected]

1

and C is a constant which is independent of both parameters N and ε. In other words, the numerical approximations U N converge to uε for all values of ε in the range 0 < ε ≤ 1. Parameter-uniform numerical methods can be constructed using standard finite difference operators and piecewise-uniform Shishkin meshes [8, 17, 20]. For singularly perturbed reaction-diffusion problems given by   Lε u ≡ εu00 − bu = f, x ∈ Ω = (0, 1), b(x) ≥ β > 0, ∀x ∈ [0, 1], (1.1)  u(0) = u0 , u(1) = u1 , the standard central difference scheme on an appropriate Shishkin piecewise-uniform mesh achieves an asymptotic parameter uniform error bound of the form [18] ku − U N k ≤ C(N −1 ln N )2 . Second order (up to logarithmic factors) convergence for singularly perturbed convectiondiffusion problems given by   Lε u ≡ εu00 + au0 − bu = f, x ∈ Ω = (0, 1), a ≥ α > 0, b ≥ 0, ∀x ∈ [0, 1], (1.2)  u(0) = u0 , u(1) = u1 , can be designed using central difference finite difference operators on a appropriately chosen Shishkin mesh [2]. However, in the case of a convection-diffusion problem the central difference operator is not monotone for ε sufficiently small. In the case of a uniform mesh the oscillations in the numerical approximations become unbounded as ε → 0, but on an appropriate piecewise-uniform Shishkin mesh the oscillations are small and bounded [8, 15]. To preserve monotonicity, upwind finite difference schemes are often used in the case of convection-diffusion problems, which results in only first order convergence being obtained [17, 20]. It is well-known that monotone second order numerical methods can be designed using the principles of defect–correction [7, 10]. The standard defect-correction approach consists of using the monotone first order upwind operator corrected by the second order centered difference operator on an appropriate Shishkin mesh. Proofs of parameter-uniform convergence in the case of the standard defect-correction approach exist in the literature [9, 14]. However, the proofs of higher order parameter–uniform convergence can be difficult, even in one dimension. In [21], the defect-correction technique is used to improve the approximations of the solutions of singularly perturbed equations of parabolic and elliptic types on a segment and on a strip, respectively, when there are not convective terms present. Defect correction techniques have also been used to increase the parameter-uniform order of convergence in time in the case of singularly perturbed reaction-diffusion and convection diffusion parabolic problems [11, 12, 13, 14]. In this paper, we examine a variant of the standard defect-correction approach to generate numerical approximations to the solution of problem (1.2). The variant is designed primarily to generate a simpler proof of second order uniform convergence in the case of a defect-correction method. Monotone second order finite difference schemes on Shishkin 2

meshes for convection-diffusion problems already exist [1, 3, 4, 6, 16]. The main interest in this paper lies not in establishing a monotone second order scheme for a convectiondiffusion problem, but rather in providing a simple proof of second order convergence of the numerical approximations generated from a defect-correction algorithm. The variant defect-correction approach adopted in this paper differs from the standard in several ways. Firstly, the hybrid difference operator only employs a first order upwind operator outside the boundary layer and when ε ≤ CN −1 . In all other cases, the uncorrected hybrid operator uses the second order central difference operator for a non uniform mesh D ± Uj =

hj+1 hj D+ Uj + D − Uj , hj + hj+1 hj + hj+1

where D+ , D− denote the classical forward and backward differences and hj , hj+1 the local step sizes of the mesh. This choice ensures that the error in the regular component of the uncorrected approximation is uniformly second order when ε ≥ CN −1 . Secondly, the subsequent analysis of the corrected approximations primarily involves the analysis of the regular component outside the boundary layer in the special case of ε ≤ CN −1 . Note that when the approximations are being corrected, a first order central difference operator is required at the transition point D0 Uj =

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

In [15], Kopteva examined the effect of using different central difference approximations for the first derivative term au0 in (1.2) on a piecewise-uniform Shishkin mesh. Three possible choices are D± U1 , 0.5(D+ + D− )U2 and D0 U3 . Note these three operators differ only at the transition point, where the Shishkin mesh is non-uniform. Using exact solutions of a discrete constant coefficient problem, Kopteva [15] showed that all three numerical solutions Ui (xj ), i = 1, 2, 3 oscillated in the coarse mesh region, but that in the coarse mesh region U1 and U2 are unbounded in the sense that ¶ µ µ ¶ 1 1 and |U2 (xi )| ≈ O |U1 (xi )| ≈ O for ε σ, εN (εN )2 while if the approximation D0 is used at the transition point then |U3 (xi ) − u(xi )| ≤ C(N −1 ln N )2 . In this paper, since the numerical methods are monotone, the numerical solutions will not oscillate. However, if the first order approximation is corrected using D± or 0.5(D+ + D− ) at the transition point the corrected solution is not ε-uniformly bounded. Moreover, when D0 is used in the correcting stage, the order of uniform convergence is increased from first order to essentially second order. See Remarks 8 and 14 for further details. In this paper C denotes a positive constant independent of the singular and discretization parameters. Moreover, we assume that the coefficients of (1.2) satisfy a, b, f ∈ C 4 (0, 1). Then, the solution can be decomposed into a sum of regular v and singular components w u = v + w, 3

satisfying [8] Lε v = f, v(0) = v ∗ (0), v(1) = u(1),

(1.3a)

∗

Lε w = 0, w(0) = u(0) − v (0), w(1) = 0,

(1.3b)

v ∗ (0)

and is suitably chosen so that the derivatives of the components satisfy the following parameter-explicit bounds |v (j) (x, ε)| ≤ C, 0 ≤ j ≤ 3, |w(x, ε)| ≤ C exp(−αx/ε),

2

|w

(j)

ε|v (4) (x, ε)| ≤ C,

(x, ε)| ≤ Cε

−j

,

1 ≤ j ≤ 4.

(1.3c) (1.3d)

The hybrid scheme

To approximate the solution of problem (1.2), we use a finite difference scheme defined on a piecewise-uniform Shishkin mesh (see for example [8, 17] for motivation for this choice of mesh). The domain is partitioned into two subintervals [0, σ] and [σ, 1], where the transition point σ is taken to be σ = min {1/2, σ0 ε ln N }, σ0 = 3/α. (2.4a) Assuming N is an even number, we divide each of the two subdomains [0, σ], [σ, 1] into N/2 ¯N uniform mesh intervals. Then, the set of mesh points Ω ε is given by ½ ih, if 0 ≤ i ≤ N/2, xi = (2.4b) σ + (i − N/2)H, if N/2 ≤ i ≤ N, where h = 2σ/N is the fine mesh step and H = 2(1 − σ)/N is the coarse mesh step. Any mesh step is denoted by 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 the following hybrid three point finite difference scheme LN Uj ≡ rj− Uj−1 + rjc Uj + rj+ Uj+1 = fj , U0 = u0 , UN = u1 , where LN

0 < j < N,

(2.4c)

 N  Lcd , if 1 ≤ j < N/2, −1 < ε, LN ≡ cd , if N/2 ≤ j < N and kakN  N −1 Lup , if N/2 ≤ j < N and kakN ≥ ε.

N The central difference operator LN cd and the upwind operator Lup are defined by 2 ± LN cd Uj ≡ εδ Uj + aj D Uj − bj Uj = fj , 2 + LN up Uj ≡ εδ Uj + aj D Uj − bj Uj = fj ,

where δ2U

j

1 = ~j

µ

¶ Uj − Uj−1 Uj+1 − Uj Uj − Uj−1 Uj+1 − Uj − , D + Uj = , D− Uj = , hj+1 hj hj+1 hj hj hj+1 D ± Uj = D+ Uj + D − Uj . hj + hj+1 hj + hj+1

Assume that N ≥ N0 where N0 is sufficiently large so that σ0 max{kak, 2} ln N0 < N0 . 4

(2.5)

Lemma 1. Assume (2.5), then the finite difference operator in (2.4c) satisfies a discrete comparison principle of the form: For any mesh function Z, if Z(0) ≥ 0, Z(1) ≥ 0 and LN Z ≤ 0 then Z ≥ 0, ∀xi . N N Proof. In the case where LN ≡ LN up the proof is standard. When L ≡ Lcd , note that we can rewrite the central difference operator in the form ¶ µ a(xj )hj+1 N Lcd Zj ≡ ε − δ 2 Zj + aj D+ Zj − bj Zj . 2

For mesh points in the layer xj < σ, (2.5) implies that µ ¶ µ ¶ a(xj )hj+1 kakσ0 ln N kakσ ≥ε 1− > 0. ε− ≥ε 1− 2 εN N For xj ≥ σ and kakN −1 < ε, we have that ε−

a(xj )hj+1 kak ≥ε− > 0. 2 N

In all cases where LN ≡ LN cd , we have that 2 + LN cd Zj ≡ µδ Zj + aj D Zj − bj Zj ,

µ > 0,

and a standard proof-by-contradiction argument can be applied. We can deduce the following truncation error bounds for the two different difference operators employed in LN : On an arbitrary mesh, we have that ³ ´ (3) (3) k(LN − L)u(x )k ≤ C ε~ ku k + h h ku k , j j j j+1 cd ³ ´ (3) (2) k(LN up − L)u(xj )k ≤ C ε~j ku k + hj+1 ku k , and on a uniform mesh with step size h ³ ´ 2 (4) 2 (3) k(LN − L)u(x )k ≤ C εh ku k + h ku k , j cd ³ ´ 2 (4) (2) k(LN up − L)u(xj )k ≤ C εh ku k + hku k . Note the increase in the order of the truncation error in the case of a uniform mesh. Let us examine the nodal error in the case of a uniform mesh xj+1 − xj = 1/N, ∀j. That is, consider the case where ε and N are such that 2σ0 ε ln N ≥ 1.

(2.6)

In this case the hybrid operator LN is the central difference operator at all the mesh points. Then, the truncation error is bounded as follows |LN (U − u)| ≤ 5

h2 . ε3

In the case of (2.6), this truncation error bound crudely yields the error bound |(U − u)| ≤ CN −2 (ln N )3 ,

if

σ = 0.5.

(2.7)

The remainder of the paper deals with the non-uniform mesh case. Consider a decomposition of the discrete solution U = V + W, where the discrete regular component V and the discrete singular component W satisfy LN V = LN U = f, LN W

= 0,

V (0) = v(0), V (1) = v(1),

(2.8)

W (0) = w(0), W (1) = w(1).

(2.9)

Lemma 2. If σ 6= 0.5 and N ≥ N0 , the error associated with the regular component satisfies ½ CN −1 (1 − xj ), if kakN −1 ≥ ε, |v(xj ) − Vj | ≤ CN −2 , if kakN −1 < ε. Proof. Note that since v is the regular component |LN (V − v)(xi )| ≤ CN −1 and thus |(V − v)(xi )| ≤ CN −1 (1 − xi ). In the case kakN −1 < ε and σ 6= 0.5, we will sharpen this error bound. In this case, the operator is the central difference on a non–uniform mesh. We have the following bounds for the truncation error ½ C(N −1 ε + N −2 ), if j = N/2, N |L (v(xj ) − Vj )| ≤ CN −2 , otherwise. Note that the truncation error is second order except at the transition point. To establish a bound of second order on the nodal error, we employ the following discrete barrier function ψj = CN −2 θ(zj ) + CN −2 (1 − xj ), where

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

The discrete derivatives of this function are    0,  0, if 0 ≤ z < σ, if 0 < z ≤ σ, −1 −1 D+ θ(z) = D− θ(z) = , if σ ≤ z < 1, , if σ < z ≤ 1,   1−σ 1−σ  0, if 0 < z < σ,     −hN , if z = σ, D± θ(z) = 2(1 − σ)     −1 , if σ < z < 1. 1−σ Then, the discrete barrier function satisfies   0, j 6= N/2, εδ 2 ψj ≤ C N −1 ε  − , j = N/2, 1−σ 6

D± ψj ≤ −CN −2 .

Therefore, |LN (v(xj ) − Vj )| ≤ −LN ψj . Applying the discrete comparison principle yields the desired bound |v(xj ) − Vj | ≤ ψj ≤ CN −2 .

Now, we analyze the error (W − w) associated with the singular component. From the bound |w| ≤ C exp(−αx/ε), we deduce that for xi ≥ σ |w(xi )| ≤ CN −3 . Motivated by this bound we consider the monotonically decaying discrete barrier function Ψ(xj ) which is the solution of the constant coefficient difference scheme (εδ 2 + αD+ )Ψ(xj ) = 0,

0 < j < N,

Ψ(0) = 1, Ψ(1) = 0.

(2.10)

This barrier function has the following properties (see [8]) Ψ > 0, D+ Ψ < 0, δ 2 Ψ > 0,

0 < xj < 1,

and using the inequality ¶ µ 2σ0 ln N −N/2 ≤ CN −σ0 , 1+ N

if

2σ0 ln N < N,

we deduce (as in pages 48, 49 of [8]) that |Ψ(xj )| ≤ CN −3 , xj ≥ σ − h. Lemma 3. If σ 6= 0.5 and N ≥ N0 , then the singular component W satisfies the bound |W (xj )| ≤ CN −3 , xj ≥ σ − h. Proof. Note that |w(0)|Ψ(0) ≥ W (0) and |w(0)|Ψ(1) ≥ W (1) = 0. Also, ¡ ¢ (εδ 2 + aD+ − b)Ψ(xj ) = (a − α)D+ − b Ψ(xj ) < 0, ¡ ¢ (εδ 2 + aD± − b)Ψ(xj ) = (a − α)D± + α(D± − D+ ) − b Ψ(xj ) < 0, since 2(D± − D+ ) = −hj+1 δ 2 . Hence, |w(0)|Ψ(xj ) ≥ |W (xj )|. Lemma 4. If σ 6= 0.5 and N ≥ N0 , the error associated with the singular component satisfies ½ CN −2 (ln N )3 , if xj < σ, |w(xj ) − Wj | ≤ CN −3 , if xj ≥ σ. 7

Proof. For xj ≥ σ, using the previous lemma and the bound on w in (1.3), we deduce that |w(xj ) − Wj | ≤ |w(xj )| + |Wj | ≤ CN −3 . For mesh points in the fine mesh region xj < σ the central difference operator leads to |LN (w − W )(xj )| ≤ C

h2 , ε3

(w − W )(0) = 0, |(w − W )(σ)| ≤ CN −3 .

Use the barrier function C(σ − xj )h2 ε−3 + CN −3 , to finish. Theorem 5. If N ≥ N0 , the error associated with the hybrid scheme is ½ CN −1 , if kakN −1 ≥ ε, |u(xj ) − Uj | ≤ −2 3 CN (ln N ) , if kakN −1 < ε.

(2.11)

Remark 6. Note that in the case of ε small relative to N −1 , the first order error bound is caused by the regular component.

3

The defect–correction scheme

To achieve almost second order uniform convergence, we only correct the solution when it is necessary using the defect correction principle. Note that when kakN −1 ≥ ε, then the use of the upwind operator (to maintain a monotone operator) in the coarse mesh has resulted in first order convergence in the regular component of the error. We will correct this defective numerical solution with a central difference operator (slightly modified at the transition point). Define LN cd∗ by ½ N Lcd Z(xj ), if j 6= N/2, N Lcd∗ Z(xj ) ≡ 2 0 εδ Z(xj ) + aj D Z(xj ) − bj Z(xj ), if j = N/2, where D0 Z(xj ) =

Z(xj+1 ) − Z(xj−1 ) . hj+1 + hj

It is worth noting that if D± is used instead of D0 at the transition point, then the resulting numerical method is not parameter-uniformly convergent. See Table 5 in §4. We correct the numerical solution U from (2.4) as follows LN ϑ = f − LN cd∗ U, in Ω, ϑ = 0, in ∂Ω,

(3.12a)

and compute the corrected solution from b = U + ϑ, if kakN −1 ≥ ε, U b = U, U otherwise. We decompose the correction as follows ϑ = ϑv + ϑw , 8

(3.12b)

where

LN ϑv = f − LN on ∂Ω, cd∗ V, in Ω, ϑv = 0, LN ϑw = −LN W, in Ω, ϑ = 0, on ∂Ω. w cd∗

(3.13)

In this way, we have that b k ≤ kv − (V + ϑv )k + kw − (W + ϑw )k. ku − U Lemma 7. Assume that N ≥ N0 and ε ≤ kakN −1 . The singular component satisfies kw − (W + ϑw )k ≤ CN −2 (ln N )3 . Proof. The triangular inequality and lemma 4 establish the bound kw − (W + ϑw )k ≤ kw − W k + kϑw k ≤ CN −2 (ln N )3 + kϑw k. To bound the second term in the above inequality we use the uniform stability of the hybrid scheme and (3.13) kϑw k ≤ k(LN )−1 kkLN ϑw k ≤ CkLN ϑw k = CkLN cd∗ W k. Note that if 1 ≤ j < N/2 N LN cd∗ Wj = Lcd Wj = 0,

and if N/2 ≤ j < N, taking into account that LN up Wj = 0, we obtain N N 0 + LN cd∗ Wj = (Lcd∗ − Lup )Wj = aj (D − D )Wj = aj

hj (D− Wj − D+ Wj ). hj + hj+1

Hence, using hj ≤ hj+1 , |LN cd∗ Wj | ≤ CN

j+1 X

|Wi |.

i=j−1

From lemma 3 we conclude that kϑw k ≤ CN −2 , and therefore the result follows. Remark 8. Consider the special case of constant a(x) = α, b ≡ 0 and ε ≤ αN −1 . The discrete singular component W generated by the hybrid scheme can be written in the form ½ W (σ) + (W (0) − W (σ))Z1 (xi ), if xi ≤ σ, W (xi ) = W (σ)Z2 (xi ), if xi ≥ σ, where

εδ 2 Z1 (xi ) + αD0 Z1 (xi ) = 0,

xi ∈ (0, σ),

Z1 (0) = 1, Z1 (σ) = 0,

and

εδ 2 Z

xi ∈ (σ, 1),

Z2 (σ) = 1, Z2 (1) = 0.

2 (xi )

+

αD+ Z2 (xi )

= 0,

Solving for Z1 exactly we get Z1 (xi ) = 1 −

1 − λi , 1 − λN/2

where 9

λ=

2−ρ , 2+ρ

ρ=

αh . ε

Recall that |W (σ)| ≤ CN −3 and note that α α 1 ζ D Z1 (σ) = − G(N ) = − , ε ε 1 − 3 lnNN 1 − ζ −

à ζ=

1 − 3 lnNN 1+

!N/2 ≤ CN −3 ,

3 lnNN

using the inequalities µ ¶ k ln N −N 1+ ≤ CN −k N

and

µ ¶ k ln N N 1− ≤ CN −k . N

Thus, εD− Z1 (σ) is independent of ε and µ +

|D W (σ)| ≤ CN Then, at the transition point, ¯ ¯ N ¯ ¯ ¯L ϑw (σ)¯ = ¯α ¯

−2

,

but

−

D W (σ) = O

G(N ) ε

¶ .

¯ ¯ h + − (D W (σ) − D W (σ))¯¯ ≤ CN −3 ln N. h+H

0 ± + − In the definition of the correcting operator LN cd∗ if D is replaced by D or by 0.5(D + D ) ∗ then at the transition point the corresponding correction ϑw satisfies ¶ µ G(N ) N ∗ . L ϑw (σ) = O ε

Either of these replacements for D0 result in LN ϑ∗w (σ) being unbounded as ε → 0 at the transition point, but |LN ϑ∗w (xi )| ≤ CN −2 , xi 6= σ. In the last section, we show the adverse effect that this term produces in the uniform errors when D± is used to correct the solution of the basic scheme. Before analyzing the error in the corrected regular component, we first estimate the discrete derivative of the error, scaled appropriately in the fine mesh region. We denote the error in the regular component and the truncation error in the regular component by e(xi ) = (V − v)(xi ),

τV (xi ) = LN (V − v)(xi ).

Lemma 9. Assume that N ≥ N0 and ε ≤ kakN −1 . The error in the regular component satisfies ½ CN −1 , if xi ≥ σ, |D+ e(xi )| ≤ CN −1 /ε, if xi < σ. Proof. We will use the bounds given in Lemma 2. From these bounds and using ε ≤ CN −1 we get that ε|D+ e(xi )| ≤ CN −1 (1 − xi ), xi ≥ σ. At the transition point, we have that 2ε(D+ − D− )e(σ) + aN/2 (h + H)D+ e(σ) − (h + H)bN/2 e(σ) = (h + H)τV (σ), 10

(3.14)

which using |(h + H)D+ e(σ)| ≤ 2(|e(σ)| + |e(σ + H)|) yields ε|D− e(σ)| ≤ CN −1 . For the mesh points in the layer xi < σ, we have that ε(D+ − D− )e(xi ) + ai hD± e(xi ) − hbi e(xi ) = hτV (xi ).

(3.15)

For any i ≤ N/2−1, we sum these equations from j = i to j = N/2−1, use |e(xi )| ≤ CN −1 , ε|D− e(σ)| ≤ CN −1 and N/2−1

X

N/2−3 ±

±

±

al D el = aN/2−1 D eN/2−1 + aN/2−2 D eN/2−2 +

l=i

X

(al D± el ± al+2 el+1 ) =

l=i

N/2−3 X ¡ ¢ = aN/2−1 eN/2 + aN/2−2 eN/2−1 + (al − al+2 )el+1 − ai+1 ei − ai ei−1 /(2h) = l=i

= O(N −1 /h), to get that all mesh points

ε|D+ e(xi )| ≤ CN −1 ,

0 ≤ xi < 1.

We improve these crude bounds following an argument developed in [5]. From (3.15) and D0 ≡ (D+ + D− )/2 on a uniform mesh, we have that for xi < σ (1 + ρi /2)D+ e(xi ) = (1 − ρi /2)D− e(xi ) +

h (bi e(xi ) + τV )(xi ), ε

ρi =

ai h . ε

In the layer region (0, σ), ρi ≤ CN −1 ln N, and therefore, 0 < 1 − ρi /2 < 1, 1 − ρi /2 1 − ρ/2 1 < < , ρ = αh/ε. 1 + ρi /2 1 + ρ/2 1+ρ Thus for points in the layer xi < σ |D+ e(xi )| ≤

1 ρ |D− e(xi )| + C N −1 . 1+ρ 1+ρ

Using this expression recursively, we have +

−

−N/2

|D e(σ − h)| = |D e(σ)| ≤ C(1 + ρ)(1 + ρ)

N/2−2 X ρ −1 N (1 + ρ)−i |D e(0)| + C 1+ρ +

i=0

≤ C(1 + ρ) ≤ CN

−1

−N/2

+

|D e(0)| + CN

−1

+

|D e(0)| + CN −1 .

From (3.14), we have that |D+ e(σ)| ≤

ε h+H CN −1 |D− e(σ)| + CN −1 ≤ ε|D+ e(0)| + CN −1 ≤ CN −1 . ε + αH ε + αH ε + αH 11

Finally, we use the same argument for the equations xi = σ + H to xj < 1 to obtain |D+ e(xj )| ≤ C|D+ e(σ)| + CN −1 ≤ CN −1 ,

σ < xj < 1.

Now we estimate the error ec (xi ) = (V + ϑv − v)(xi ) in the corrected regular solution. Lemma 10. Assume that N ≥ N0 and that σ 6= 0.5. The truncation error associated with the regular component of the corrected solution satisfies kLN ec k ≤ CN −1 ln N. Proof. First, we calculate an expression for the truncation error N N N N N cd LN ec = LN e + LN ϑv = LN e + f − LN cd∗ V = (L − Lcd∗ )e − Lcd∗ v + f = (L − Lcd∗ )e + τV , (3.16) cd N where τV (xi ) = (f − Lcd∗ v)(xi ). Note that   0, if xi < σ, ai hi (LN − LN + − + 0 cd∗ )e(xi ) = (D − D )e(xi ), if xi ≥ σ.  ai (D − D )e(xi ) = hi + hi+1

If xi < σ, the mesh is uniform and trivially |τVcd (xi )| ≤ CN −2 . If xi > σ, by using that the mesh is uniform and lemma 9, we have |LN ec (xi )| ≤ |τVcd (xi )| +

ai (|D+ e(xi )| + |D+ e(xi−1 )|) ≤ CN −1 . 2

Finally, if xi = σ, using lemma 9 and h/(h + H) = O(ε ln N ), we get that |LN ec (xi )| ≤ |τVcd (σ)| +

aN/2 h (|D+ e(σ)| + |D+ e(σ − h)|) ≤ CN −1 ln N. h+H

Lemma 11. Assume that N ≥ N0 and σ 6= 0.5. The error in the regular component associated with the corrected solution satisfies |eci (xi )| ≤ CN −1 ln N (1 − xi ) ½ + c

|D e (xi )| ≤

CN −1 ln N, if xi ≥ σ, CN −1 ln N/ε, if xi < σ.

Proof. The first inequality follows from the previous lemma. Using the same ideas as in the proof of lemma 9, we deduce that ε|D+ ec (xi )| ≤ CN −1 ,

0 ≤ xi < 1,

|D− ec (σ)| ≤ CN −1 |D+ ec (0)| + CN −1 .

Then, at the transition point |D+ ec (σ)| ≤

ε h+H h |D− ec (σ)|+CN −1 +aN/2 |(D+ −D− )e(σ)| ≤ CN −1 ln N. ε + αH ε + αH ε + αH 12

We finish the proof in the same way as in the proof lemma 9 since the new term satisfies aj (D+ − D− )e(xj ) = O(N −1 ),

xj > σ.

The previous lemma yields an error bound of the regular component for the corrected solution which is not an improvement on the error bound of the uncorrected solution. In the next lemma, this first order error bound is sharpened to indicate an increased order of convergence in the corrected solution. The key ingredient in the proof is summarized in the following observation. On a uniform mesh, if ai D+ ei − bi ei = τˆi , i < N

eN = 0,

then, we establish below that ¯ ¯ −1 ¯ ¯NX ¯ ¯ |ei | ≤ CN −1 ¯ τˆk ¯ ≤ max |ˆ τk |. ¯ ¯ k k=i

Hence a first order truncation error of the form maxk |ˆ τk | ≤ CN −1 can produce a second order error bound if ¯N −1 ¯ ¯X ¯ ¯ ¯ τˆk ¯ ≤ CN −1 . ¯ ¯ ¯ k=i

Lemma 12. Assume that N ≥ N0 and σ 6= 0.5. The regular component associated with the corrected solution satisfies the error bound kv − (V + ϑv )k ≤ CN −2 ln N. Proof. Firstly, we establish the bound in the interval [σ + H, 1]. From (3.16), we have for l = i, . . . , N − 1 ε al (D+ ecl − D+ ecl−1 ) + al D+ ecl − bl ecl = τlcd + (D+ − D− )el . H 2 Summing all these equations from l = i to l = N − 1, results in the equation N −1 X

(al D+ ecl

l=i

−

bl ecl )

N −1 N −1 X X ε al + + c + c − = − (D eN −1 − D ei−1 ) + (D − D )el + τlcd . H 2 l=i

l=i

We bound each term of this expression separately. Note that, by lemma 11 and εN ≤ C, ε + c D ei−1 = O(N −1 ln N ) i > N/2, H

|

N −1 X l=i

τlcd | ≤

N −1 X

|τlcd | ≤ CN −1 .

l=i

Moreover, N −1 X l=i

al (D+ − D− )el = aN −1 D+ eN −1 +

N −2 X

(al − al+1 )D+ el − ai D− ei = O(N −1 ).

l=i

13

Thus we have the equation al D+ ecl − bl ecl = gl ,

where |

j X

gl | ≤ CN −1 ln N ;

i ≤ j ≤ N − 1.

(3.17)

l=i

From (3.17), we have ecl

µ ¶µ ¶ Hgl Hbl −1 c = el+1 − 1+ . al al

Hence eci

ÃN −2 µ Ãk−1 µ ¶ ! ¶ ! N −1 Y X Y Hbk −1 c Hbl −1 Hgk−1 = eN −1 − . 1+ 1+ ak al ak−1 k=i

Let βk−1 =

1 ak−1

N −2 X k=i

k=i+1

k−1 Yµ

1+

l=i

βk gk =

NX −i−2

Hbl al

l=i

¶−1 , then

à i+k ! N −2 X X gl = O(N −1 ln N ). gl (βi+k − βi+k+1 ) + βN −1 l=i

l=i

k=0

From this expression and |ecN −1 | ≤ CN −2 ln N , we have that |eci | ≤ CN −2 ln N,

i > N/2,

with C a positive constant independent of N . This is easily extended to |eci | ≤ CN −2 ln N, i ≥ N/2 using lemma 11. Finally, in [0, σ], it holds that ½

|LN ec (xi )| ≤ CN −2 , xi ∈ (0, σ), ec (0) = 0, |ec (σ)| ≤ CN −2 ln N.

Then, the discrete minimum principle on [0, σ], yields |ec (xi )| ≤ CN −2 ln N,

xi ∈ (0, σ].

An immediate consequence of theorem 5, lemmas 7, 12 and the bound (2.7) is the following parameter uniform error bound. Theorem 13. Assume that N ≥ N0 . The error associated with the defect–correction scheme satisfies, b k ≤ CN −2 ln3 N, ku − U b is the numerical approximation given in where u is the solution of the problem (1.2) and U (3.12b). 14

4

Numerical experiments

In this section we show the numerical results obtained with several schemes: upwind, hybrid and some variants of schemes of defect–correction type. We consider the following constant coefficient problem ½ 00 εu + u0 = −ex (1 + ε), x ∈ (0, 1), (4.18) u(0) = 1, u(1) = 1 − e, whose exact solution is known. For any value of N , the maximum pointwise errors EεN and the ε–uniform errors are calculated using b kΩN , E N = max EεN , EεN = ku − U (4.19) ε

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

EεN , Eε2N

pN = log2

EN . E 2N

(4.20)

In all of the following tables of numerical results the mesh employed is the Shishkin mesh defined in (2.4b). In the first row we display the maximum errors EεN and in the second N and row the orders of convergence pN ε . In the last row we display the uniform-errors E the corresponding uniform orders of convergence pN . In table 1, the standard upwind scheme is seen to be essentially first order for all values of the parameter ε. We can observe first order of convergence for large values of ε and that it degrades down to (N −1 ln N ) for ε small relative to N . In table 2, the classical defect-correction method of simple upwinding corrected with the central difference operator,  N  Lup U = f, in Ω, U0 = u0 , UN = u1 , LN ϑ = f − LN in Ω, ϑ0 = ϑN = 0, cd∗ U,  bup U = U + ϑ, is seen to be second order up to logarithmic factors. As in table 1, the logarithmic factor is not observed for large values of ε. In table 3 the hybrid finite difference scheme (2.4c) is examined. Note that the orders are second order for large values of ε and first order for small values of ε (relative to N −1 ). Note the absence of any logarithmic factors in the uniform orders of convergence. The results in this table support the theoretical bound established in Theorem 5. In table 4 we show the results obtained with the numerical scheme (3.12b). Now, the hybrid scheme (2.4c) is corrected with the central difference scheme when ε is small. Note that the orders are now second order up to logarithmic factors. The results in this table are similar to those presented in table 2 for the standard defect-correction method.

15

Remark 14. The hybrid defect-correction scheme used in table 4 can be modified by replacing the operator D0 with D± at the transition point  N L U = fj , in Ω, U0 = u0 , UN = u1 ,    N  L ϑ = f − LN in Ω, ϑ0 = ϑN = 0, cd U, ( (4.21) b = U + ϑ, if kakN −1 ≥ ε, U     b = U, U otherwise. In Figure 1, the numerical solution and the correction ϑ generated by applying the numerical method (4.21) to problem (4.18) are displayed for N = 32 and ε = 10−6 . In table 5 the maximum errors and the orders of convergence are displayed for this method applied to problem (4.18). In this table it is seen that the errors continue to increase as ε → 0. This is not a parameter-uniform numerical method. 1 0

0.5

−0.2 −0.4

0 −0.6

−0.5

−0.8 −1

−1 −1.2

−1.5

−1.4 −1.6

−2 −1.8

−2.5

0

0.2

0.4

0.6

0.8

−2

1

0

0.2

0.4

0.6

0.8

1

Figure 1: Corrected numerical solution and the correction ϑ generated by the numerical method (4.21) applied to problem (4.18) with N = 32 and ε = 10−6 . Throughout the theoretical sections of this paper, we have used σ0 = 3/α. However, numerical experiments show that it is sufficient to take σ0 = 2/α. Concretely, for the test problem (4.18), in table 6 we present the uniform errors and their corresponding orders of convergence associated with the numerical scheme defined in §3. We take σ0 = 0.5, 1, 2, 3, 4 and the parameters ε and N vary in the same range of values than in the previous tables. We observe that if σ0 ≤ 1 the method does not achieve second order of convergence. On the other hand, if σ0 ≥ 2, the method has uniform rate of convergence of second order. Finally, we demonstrate numerically that the correction of the scheme (3.12b) occurs mainly in the regular component. We define the following auxiliary problem ½ 00 εu + u0 = −ex (1 + ε), x ∈ (0, 1), (4.22) u(0) = 0, u(1) = 1 − e. Note that the value of the solution of (4.22) at x = 0 coincides with the solution of the reduced problem (ε = 0) associated with problem (4.18) ½ 0 u = −ex , x ∈ (0, 1), (4.23) u(1) = 1 − e. 16

Table 1: Upwind scheme on the mesh (2.4a,b): Maximum pointwise errors EεN , maximum N uniform errors E N , orders of convergence pN ε and the order of ε–uniform convergence p . ε 20 2−1 −2

2

2−3 −4

2

−5

2

−6

2

2−7 −8

2

−9

2

2−10 2−11 2

−12

2

−13

. . . .

2−29 2−30 EN pN

N=32 1.544E-3 1.016 1.311E-3 0.925 1.022E-2 0.940 2.925E-2 0.906 6.353E-2 0.870 7.599E-2 0.510 6.981E-2 0.478 6.599E-2 0.450 6.463E-2 0.449 6.392E-2 0.449 6.356E-2 0.449 6.338E-2 0.449 6.329E-2 0.449 6.324E-2 0.450 . . . . 6.320E-2 0.450 6.320E-2 0.450

N=64 7.638E-4 1.008 6.903E-4 0.962 5.331E-3 0.970 1.560E-2 0.958 3.475E-2 0.928 5.335E-2 0.637 5.014E-2 0.607 4.831E-2 0.604 4.733E-2 0.604 4.681E-2 0.604 4.655E-2 0.604 4.641E-2 0.604 4.635E-2 0.604 4.631E-2 0.604 . . . . 4.628E-2 0.604 4.628E-2 0.604

N=128 3.798E-4 1.004 3.543E-4 0.982 2.721E-3 0.985 8.030E-3 0.978 1.826E-2 0.959 3.431E-2 0.783 3.291E-2 0.713 3.179E-2 0.709 3.115E-2 0.707 3.080E-2 0.706 3.063E-2 0.705 3.053E-2 0.705 3.049E-2 0.705 3.047E-2 0.705 . . . . 3.044E-2 0.705 3.044E-2 0.705

N=256 1.893E-4 1.002 1.794E-4 0.991 1.375E-3 0.992 4.078E-3 0.989 9.398E-3 0.980 1.994E-2 0.963 2.007E-2 0.766 1.945E-2 0.764 1.908E-2 0.763 1.889E-2 0.763 1.878E-2 0.763 1.873E-2 0.762 1.870E-2 0.762 1.869E-2 0.762 . . . . 1.868E-2 0.762 1.868E-2 0.762

N=512 9.454E-5 1.001 9.026E-5 0.995 6.911E-4 0.996 2.054E-3 0.994 4.765E-3 0.990 1.023E-2 0.980 1.180E-2 0.807 1.146E-2 0.804 1.125E-2 0.802 1.113E-2 0.802 1.107E-2 0.801 1.104E-2 0.801 1.103E-2 0.801 1.102E-2 0.801 . . . . 1.101E-2 0.801 1.101E-2 0.801

N=1024 4.724E-5 1.001 4.527E-5 0.998 3.464E-4 0.998 1.031E-3 0.997 2.400E-3 0.995 5.185E-3 0.990 6.745E-3 0.835 6.563E-3 0.832 6.448E-3 0.831 6.386E-3 0.830 6.353E-3 0.830 6.336E-3 0.830 6.328E-3 0.830 6.323E-3 0.830 . . . . 6.319E-3 0.830 6.319E-3 0.830

N=2048 2.361E-5 1.000 2.267E-5 0.999 1.734E-4 0.999 5.166E-4 0.999 1.204E-3 0.997 2.610E-3 0.995 3.782E-3 0.855 3.687E-3 0.852 3.625E-3 0.852 3.591E-3 0.851 3.573E-3 0.851 3.564E-3 0.851 3.559E-3 0.851 3.557E-3 0.851 . . . . 3.554E-3 0.851 3.554E-3 0.851

N=4096 1.180E-5 1.000 1.134E-5 0.999 8.678E-5 1.000 2.585E-4 0.999 6.031E-4 0.999 1.310E-3 0.998 2.091E-3 0.871 2.043E-3 0.868 2.009E-3 0.867 1.990E-3 0.867 1.981E-3 0.867 1.976E-3 0.867 1.973E-3 0.867 1.972E-3 0.867 . . . . 1.971E-3 0.867 1.971E-3 0.867

N=8192 5.901E-6

6.320E-2 0.450

4.628E-2 0.604

3.044E-2 0.705

1.868E-2 0.762

1.101E-2 0.801

6.319E-3 0.830

3.554E-3 0.851

1.971E-3 0.867

1.081E-3

5.674E-6 4.340E-5 1.293E-4 3.018E-4 6.559E-4 1.143E-3 1.119E-3 1.101E-3 1.091E-3 1.086E-3 1.083E-3 1.082E-3 1.081E-3 . . . . 1.081E-3 1.081E-3

In Figure 2 we show the correction ϑ of problems (4.18) and (4.22) associated with the numerical scheme (3.12b) for N = 64 and ε = 10−6 . The values are quite similar (see table 7) and they support the theoretical analysis developed in the previous sections.

References [1] V.B. Andr´eiev and I.A. Savin, The uniform convergence with respect to a small parameter of A.A. Samarskii’s monotone scheme and its modification, Comp. Maths Math. Phys., 35 (1995) 581–591. [2] V.B. Andr´eiev, N.V. Kopteva, A study of difference schemes with the first derivative approximated by a central difference ratio, Comp. Maths. Math. Phys., 36 (8) (1996) 1065-1078. [3] V.B. Andr´eiev, N.V. Kopteva, On the convergence, uniform with respect to a small parameter, of monotone three-point finite-difference approximations, Differential Equations, 34 (7) (1998) 921-928. 17

Table 2: Standard defect correction on the mesh (2.4a,b): Maximum pointwise errors EεN , maximum uniform errors E N , orders of convergence pN ε and the order of ε–uniform convergence pN . ε 20 2−1 −2

2

−3

2

2−4 −5

2

−6

2

−7

2

2−8 −9

2 2

−10

2−11 2−12 2

−13

. . . .

2−29 2−30 EN pN

N=32 8.359E-5 1.994 1.724E-4 1.987 5.782E-4 1.949 2.445E-3 1.910 7.493E-3 1.851 1.263E-2 1.473 1.453E-2 1.519 1.568E-2 1.539 1.632E-2 1.547 1.665E-2 1.551 1.682E-2 1.553 1.691E-2 1.553 1.695E-2 1.554 1.697E-2 1.554 . . . . 1.699E-2 1.554 1.699E-2 1.554

N=64 2.099E-5 1.997 4.349E-5 1.994 1.497E-4 1.973 6.507E-4 1.954 2.077E-3 1.912 4.550E-3 1.518 5.069E-3 1.550 5.399E-3 1.561 5.584E-3 1.567 5.684E-3 1.571 5.734E-3 1.573 5.760E-3 1.574 5.773E-3 1.574 5.779E-3 1.574 . . . . 5.786E-3 1.574 5.786E-3 1.574

N=128 5.261E-6 1.998 1.092E-5 1.997 3.813E-5 1.987 1.679E-4 1.976 5.520E-4 1.954 1.589E-3 1.682 1.732E-3 1.609 1.830E-3 1.617 1.884E-3 1.621 1.913E-3 1.622 1.927E-3 1.623 1.935E-3 1.624 1.939E-3 1.624 1.941E-3 1.624 . . . . 1.943E-3 1.624 1.943E-3 1.624

N=256 1.317E-6 1.999 2.736E-6 1.998 9.622E-6 1.993 4.269E-5 1.988 1.425E-4 1.976 4.951E-4 1.950 5.677E-4 1.661 5.963E-4 1.668 6.126E-4 1.671 6.213E-4 1.672 6.257E-4 1.672 6.279E-4 1.672 6.290E-4 1.672 6.295E-4 1.672 . . . . 6.301E-4 1.672 6.301E-4 1.672

N=512 3.293E-7 2.000 6.847E-7 1.999 2.417E-6 1.997 1.076E-5 1.994 3.621E-5 1.988 1.282E-4 1.974 1.795E-4 1.701 1.876E-4 1.707 1.925E-4 1.710 1.950E-4 1.711 1.963E-4 1.711 1.970E-4 1.712 1.973E-4 1.712 1.975E-4 1.712 . . . . 1.977E-4 1.712 1.977E-4 1.712

N=1024 8.235E-8 2.000 1.713E-7 2.000 6.056E-7 1.998 2.702E-6 1.997 9.128E-6 1.994 3.262E-5 1.987 5.520E-5 1.733 5.747E-5 1.738 5.884E-5 1.740 5.957E-5 1.741 5.996E-5 1.741 6.015E-5 1.742 6.025E-5 1.742 6.030E-5 1.742 . . . . 6.035E-5 1.742 6.035E-5 1.742

N=2048 2.059E-8 2.000 4.283E-8 2.000 1.516E-7 1.999 6.769E-7 1.998 2.292E-6 1.997 8.231E-6 1.993 1.660E-5 1.758 1.722E-5 1.763 1.761E-5 1.764 1.782E-5 1.765 1.793E-5 1.765 1.799E-5 1.766 1.802E-5 1.766 1.803E-5 1.766 . . . . 1.804E-5 1.766 1.804E-5 1.766

N=4096 5.148E-9 2.000 1.071E-8 2.000 3.792E-8 2.000 1.694E-7 1.999 5.742E-7 1.998 2.067E-6 1.997 4.909E-6 1.778 5.076E-6 1.782 5.184E-6 1.784 5.243E-6 1.784 5.274E-6 1.785 5.290E-6 1.785 5.298E-6 1.785 5.302E-6 1.785 . . . . 5.306E-6 1.785 5.306E-6 1.785

N=8192 1.287E-9

1.699E-2 1.554

5.786E-3 1.574

1.943E-3 1.624

6.301E-4 1.672

1.977E-4 1.712

6.035E-5 1.742

1.804E-5 1.766

5.306E-6 1.785

1.540E-6

2.677E-9 9.483E-9 4.237E-8 1.437E-7 5.180E-7 1.431E-6 1.476E-6 1.506E-6 1.522E-6 1.531E-6 1.535E-6 1.538E-6 1.539E-6 . . . . 1.540E-6 1.540E-6

[4] V.B. Andr´eiev, N.V. Kopteva, Uniform with respect to a small parameter convergence of difference schemes for a convection-diffusion problem, in: J.J.H. Miller, G.I. Shishkin & L. Vulkov, Eds. Analytical and Computational Methods for Convection-Dominated and Singularly Perturbed Problems, Nova Science Publishers, New York, 2000, 133-139. [5] S. Bellew, E. O’Riordan, A parameter robust numerical method for a system of two singularly perturbed convection–diffusion equations, Applied Numerical Mathematics 51 (2004) 171-186. [6] 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. [7] V. Ervin and W. Layton, An analysis of a defect-correction method for a model convection-diffusion equation, SIAM J. Numer. Anal., 26 (1) (1989) 169-179. 18

Table 3: Hybrid scheme (2.4c), uncorrected, on the mesh (2.4a,b): Maximum pointwise errors EεN , maximum uniform errors E N , orders of convergence pN ε and the order of ε– uniform convergence pN . ε 20 2−1 −2

2

2−3 −4

2

−5

2

−6

2

2−7 −8

2

−9

2

2−10 2−11 −12

2

−13

2

. . . .

2−29 2−30 EN pN

N=32 6.093E-5 1.999 1.522E-4 2.000 5.467E-4 2.002 2.068E-3 2.007 8.043E-3 2.030 3.273E-2 2.790 4.366E-2 1.169 4.970E-2 1.100 5.293E-2 1.071 5.467E-2 1.059 5.555E-2 1.053 5.599E-2 1.050 5.621E-2 1.048 5.632E-2 1.047 . . . . 5.644E-2 1.047 5.644E-2 1.047

N=64 1.524E-5 2.000 3.804E-5 2.000 1.365E-4 2.000 5.144E-4 2.001 1.970E-3 2.007 4.731E-3 1.565 1.941E-2 3.592 2.319E-2 1.075 2.520E-2 1.046 2.624E-2 1.033 2.678E-2 1.027 2.705E-2 1.024 2.719E-2 1.022 2.725E-2 1.022 . . . . 2.732E-2 1.021 2.732E-2 1.021

N=128 3.811E-6 2.000 9.509E-6 2.000 3.412E-5 2.000 1.285E-4 2.000 4.899E-4 2.002 1.599E-3 1.729 1.610E-3 1.620 1.100E-2 4.386 1.220E-2 1.028 1.283E-2 1.016 1.314E-2 1.010 1.330E-2 1.008 1.338E-2 1.006 1.342E-2 1.006 . . . . 1.346E-2 1.005 1.346E-2 1.005

N=256 9.528E-7 2.000 2.377E-6 2.000 8.531E-6 2.000 3.212E-5 2.000 1.223E-4 2.000 4.825E-4 2.002 5.236E-4 1.665 5.262E-4 1.667 5.986E-3 5.171 6.342E-3 1.010 6.524E-3 1.005 6.615E-3 1.002 6.661E-3 1.001 6.685E-3 1.001 . . . . 6.708E-3 1.000 6.708E-3 1.000

N=512 2.382E-7 2.000 5.943E-7 2.000 2.133E-6 2.000 8.030E-6 2.000 3.057E-5 2.000 1.205E-4 2.000 1.651E-4 1.699 1.657E-4 1.700 1.662E-4 1.701 3.150E-3 5.943 3.251E-3 1.004 3.302E-3 1.002 3.328E-3 1.001 3.341E-3 1.000 . . . . 3.354E-3 1.000 3.354E-3 1.000

N=1024 5.955E-8 2.000 1.486E-7 2.000 5.332E-7 2.000 2.007E-6 2.000 7.643E-6 2.000 3.011E-5 2.000 5.083E-5 1.727 5.101E-5 1.728 5.113E-5 1.728 5.120E-5 1.728 1.621E-3 6.712 1.649E-3 1.002 1.663E-3 1.001 1.670E-3 1.000 . . . . 1.677E-3 1.000 1.677E-3 1.000

N=2048 1.489E-8 2.000 3.715E-8 2.000 1.333E-7 2.000 5.019E-7 2.000 1.911E-6 2.000 7.526E-6 2.000 1.536E-5 1.750 1.540E-5 1.751 1.543E-5 1.751 1.545E-5 1.751 1.546E-5 1.752 8.233E-4 7.486 8.311E-4 1.001 8.350E-4 1.000 . . . . 8.389E-4 1.000 8.389E-4 1.000

N=4096 3.723E-9 1.997 9.289E-9 1.998 3.332E-8 2.000 1.255E-7 2.000 4.777E-7 2.000 1.881E-6 2.000 4.564E-6 1.770 4.576E-6 1.771 4.584E-6 1.771 4.589E-6 1.771 4.592E-6 1.771 4.593E-6 1.771 4.152E-4 8.269 4.173E-4 1.001 . . . . 4.195E-4 1.000 4.195E-4 1.000

N=8192 9.326E-10

5.644E-2 1.047

2.732E-2 1.021

1.346E-2 1.005

6.708E-3 1.000

3.354E-3 1.000

1.677E-3 1.000

8.389E-4 1.000

4.195E-4 1.000

2.097E-4

2.325E-9 8.332E-9 3.137E-8 1.194E-7 4.703E-7 1.338E-6 1.341E-6 1.343E-6 1.345E-6 1.345E-6 1.346E-6 1.346E-6 2.086E-4 . . . . 2.097E-4 2.097E-4

[8] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan and G.I. Shishkin, Robust computational techniques for boundary layers, Chapman and Hall/CRC, Boca Raton, FL., 2000. [9] A. Frohner, T. Linß and H.-G. Roos, Defect correction on Shishkin-type meshes, Num. Algor. 26 (2001) 281-299. [10] P.W. Hemker, An accurate method without directional bias for the numerical solution of a 2-D elliptic singular perturbation problem, in Theory and Applications of Singular Perturbations, W. Eckhaus and E.M. de Jager, eds., Lecture Notes in Mathematics 942, Springer-Verlag, Berlin, 1982. [11] P.W. Hemker, G.I. Shishkin and L.P. Shishkina, The use of defect correction for the solution of parabolic singular perturbation problems, ZAMM, 77 (1) (1997) 59-74. [12] P.W. Hemker, G.I. Shishkin and L.P. Shishkina, ε-uniform schemes with high-order time-accuracy for parabolic singular perturbation problems, IMA J. Numer. Anal., 20 (1) (2000) 99-121. 19

Table 4: Hybrid scheme (2.4c), corrected, on the mesh (2.4a,b): Maximum pointwise errors EεN , maximum uniform errors E N , orders of convergence pN ε and the order of ε–uniform convergence pN . ε 20 2−1 −2

2

2−3 −4

2

−5

2

−6

2

2−7 −8

2

−9

2

2−10 2−11 −12

2

−13

2

. . . .

2−29 2−30 EN pN

N=32 6.093E-5 1.999 1.522E-4 2.000 5.467E-4 2.002 2.068E-3 2.007 8.043E-3 2.030 1.443E-2 1.609 1.534E-2 1.565 1.586E-2 1.574 1.613E-2 1.578 1.627E-2 1.579 1.634E-2 1.580 1.637E-2 1.581 1.639E-2 1.581 1.640E-2 1.581 . . . . 1.641E-2 1.581 1.641E-2 1.581

N=64 1.524E-5 2.000 3.804E-5 2.000 1.365E-4 2.000 5.144E-4 2.001 1.970E-3 2.007 4.731E-3 1.565 5.184E-3 1.687 5.327E-3 1.617 5.404E-3 1.620 5.444E-3 1.621 5.464E-3 1.622 5.473E-3 1.622 5.478E-3 1.622 5.481E-3 1.622 . . . . 5.483E-3 1.622 5.483E-3 1.622

N=128 3.811E-6 2.000 9.509E-6 2.000 3.412E-5 2.000 1.285E-4 2.000 4.899E-4 2.002 1.599E-3 1.729 1.610E-3 1.620 1.736E-3 1.722 1.758E-3 1.651 1.769E-3 1.652 1.775E-3 1.653 1.778E-3 1.653 1.780E-3 1.653 1.780E-3 1.653 . . . . 1.781E-3 1.653 1.781E-3 1.653

N=256 9.528E-7 2.000 2.377E-6 2.000 8.531E-6 2.000 3.212E-5 2.000 1.223E-4 2.000 4.825E-4 2.002 5.236E-4 1.665 5.262E-4 1.667 5.600E-4 1.753 5.630E-4 1.687 5.646E-4 1.687 5.654E-4 1.687 5.658E-4 1.688 5.660E-4 1.688 . . . . 5.662E-4 1.688 5.662E-4 1.688

N=512 2.382E-7 2.000 5.943E-7 2.000 2.133E-6 2.000 8.030E-6 2.000 3.057E-5 2.000 1.205E-4 2.000 1.651E-4 1.699 1.657E-4 1.700 1.662E-4 1.701 1.749E-4 1.772 1.753E-4 1.715 1.755E-4 1.715 1.756E-4 1.715 1.757E-4 1.715 . . . . 1.758E-4 1.715 1.758E-4 1.715

N=1024 5.955E-8 2.000 1.486E-7 2.000 5.332E-7 2.000 2.007E-6 2.000 7.643E-6 2.000 3.011E-5 2.000 5.083E-5 1.727 5.101E-5 1.728 5.113E-5 1.728 5.120E-5 1.728 5.340E-5 1.788 5.347E-5 1.739 5.350E-5 1.739 5.351E-5 1.739 . . . . 5.353E-5 1.739 5.353E-5 1.739

N=2048 1.489E-8 2.000 3.715E-8 2.000 1.333E-7 2.000 5.019E-7 2.000 1.911E-6 2.000 7.526E-6 2.000 1.536E-5 1.750 1.540E-5 1.751 1.543E-5 1.751 1.545E-5 1.751 1.546E-5 1.752 1.602E-5 1.802 1.603E-5 1.760 1.603E-5 1.760 . . . . 1.604E-5 1.760 1.604E-5 1.760

N=4096 3.723E-9 1.997 9.289E-9 1.998 3.332E-8 2.000 1.255E-7 2.000 4.777E-7 2.000 1.881E-6 2.000 4.564E-6 1.770 4.576E-6 1.771 4.584E-6 1.771 4.589E-6 1.771 4.592E-6 1.771 4.593E-6 1.771 4.733E-6 1.814 4.735E-6 1.777 . . . . 4.736E-6 1.777 4.736E-6 1.777

N=8192 9.326E-10

1.641E-2 1.581

5.483E-3 1.622

1.781E-3 1.653

5.662E-4 1.688

1.758E-4 1.715

5.353E-5 1.739

1.604E-5 1.760

4.736E-6 1.777

1.381E-6

2.325E-9 8.332E-9 3.137E-8 1.194E-7 4.703E-7 1.338E-6 1.341E-6 1.343E-6 1.345E-6 1.345E-6 1.346E-6 1.346E-6 1.381E-6 . . . . 1.381E-6 1.381E-6

[13] P.W. Hemker, G.I. Shishkin and L.P. Shishkina, High-order time-accurate schemes for parabolic singular perturbation problems with convection, Russ. J. Numer. Anal. Math. Modelling, 17 (1) (2002) 1-24. [14] P.W. Hemker, G.I. Shishkin and L.P. Shishkina, Novel defect correction high-order, in space and time, accurate schemes for parabolic singularly perturbed convection-diffusion problems, Computational Methods in Applied Mathematics, 3 (3) (2003) 387-404. [15] N.V. Kopteva, A study of certain difference schemes for convection-diffusion equations, Diplom Thesis, Moscow State University, 1993. [16] N.V. Kopteva, Uniform convergence with respect to a small parameter of a four-point scheme for the one-dimensional stationary convection-diffusion equation, Differential Equations, 32 (7) (1996) 951-957. [17] J.J.H. Miller, E. O’Riordan and G.I. Shishkin, Fitted numerical methods for singular perturbation problems, World-Scientific, (Singapore), 1996. 20

Table 5: Scheme (4.21) on the mesh (2.4a,b): Maximum pointwise errors EεN and orders of convergence pN ε . ε 20 2−2 2

−4

2−6 2

−8

−10

2

−12

2

2−14 −16

2

−18

2

2−20 2−22 −24

2

−26

2

2−28 −30

2

N=32 6.093E-5 1.999 5.467E-4 2.002 8.043E-3 2.030 1.500E-2 1.560 1.623E-2 1.600 1.769E-2 1.680 2.236E-2 1.940 4.251E-2 2.615 1.368E-1 3.529 5.320E-1 3.910 2.115E+0 3.951 8.445E+0 3.958 3.377E+1 3.959 1.351E+2 3.960 5.402E+2 3.960 2.161E+3 3.960

N=64 1.524E-5 2.000 1.365E-4 2.000 1.970E-3 2.007 5.087E-3 1.660 5.352E-3 1.625 5.522E-3 1.641 5.828E-3 1.697 6.940E-3 1.893 1.185E-2 2.472 3.538E-2 3.419 1.368E-1 3.942 5.436E-1 4.006 2.171E+0 4.013 8.680E+0 4.015 3.472E+1 4.015 1.389E+2 4.015

N=128 3.811E-6 2.000 3.412E-5 2.000 4.899E-4 2.002 1.610E-3 1.620 1.736E-3 1.650 1.770E-3 1.658 1.798E-3 1.668 1.868E-3 1.710 2.136E-3 1.864 3.309E-3 2.345 8.896E-3 3.259 3.384E-2 3.915 1.345E-1 4.019 5.370E-1 4.027 2.147E+0 4.029 8.587E+0 4.029

N=256 9.528E-7 2.000 8.531E-6 2.000 1.223E-4 2.000 5.236E-4 1.665 5.530E-4 1.734 5.611E-4 1.690 5.658E-4 1.691 5.708E-4 1.698 5.867E-4 1.729 6.514E-4 1.849 9.294E-4 2.242 2.243E-3 3.092 8.291E-3 3.854 3.294E-2 4.015 1.315E-1 4.024 5.259E-1 4.026

N=512 2.382E-7 2.000 2.133E-6 2.000 3.057E-5 2.000 1.651E-4 1.699 1.662E-4 1.701 1.740E-4 1.716 1.752E-4 1.717 1.759E-4 1.718 1.770E-4 1.724 1.808E-4 1.748 1.965E-4 1.843 2.631E-4 2.166 5.735E-4 2.940 2.037E-3 3.775 8.087E-3 4.005 3.229E-2 4.016

N=1024 5.955E-8 2.000 5.332E-7 2.000 7.643E-6 2.000 5.083E-5 1.727 5.113E-5 1.728 5.296E-5 1.776 5.328E-5 1.741 5.347E-5 1.740 5.358E-5 1.741 5.383E-5 1.746 5.479E-5 1.766 5.862E-5 1.842 7.474E-5 2.112 1.489E-4 2.810 5.038E-4 3.687 1.996E-3 3.993

N=2048 1.489E-8 2.000 1.333E-7 2.000 1.911E-6 2.000 1.536E-5 1.750 1.543E-5 1.751 1.546E-5 1.752 1.594E-5 1.761 1.600E-5 1.761 1.603E-5 1.761 1.605E-5 1.761 1.611E-5 1.765 1.635E-5 1.781 1.729E-5 1.845 2.123E-5 2.073 3.910E-5 2.700 1.254E-4 3.597

N=4096 3.723E-9 1.997 3.332E-8 2.000 4.777E-7 2.000 4.564E-6 1.770 4.584E-6 1.771 4.592E-6 1.771 4.705E-6 1.806 4.721E-6 1.779 4.731E-6 1.778 4.736E-6 1.778 4.740E-6 1.779 4.755E-6 1.782 4.813E-6 1.796 5.046E-6 1.849 6.016E-6 2.044 1.036E-5 2.609

N=8192 9.326E-10 8.332E-9 1.194E-7 1.338E-6 1.343E-6 1.345E-6 1.346E-6 1.376E-6 1.379E-6 1.381E-6 1.382E-6 1.383E-6 1.386E-6 1.401E-6 1.459E-6 1.698E-6

[18] 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. [19] H.-G. Roos, M. Stynes and L. Tobiska, Numerical methods for singularly perturbed differential equations. Convection–diffusion and flow problems, Springer–Verlag, New York, 1996. [20] G.I. Shishkin, Discrete approximation of singularly perturbed elliptic and parabolic equations, Russian Academy of Sciences, Ural section, Ekaterinburg, 1992. (In Russian) [21] G.I. Shishkin, Method of improving the accuracy of the approximate solutions to singularly perturbed equations by defect correction, Russ. J. Numer. Anal. Math. Modelling, 11 (6) (1996) 539-557.

21

Table 6: Hybrid scheme (2.4c), corrected, on the mesh (2.4a,b): Uniform errors E N and uniform orders of convergence pN for problem (4.18) for several values of σ0 . ε σ0 = 0.5 σ0 = 1 σ0 = 2 σ0 = 3 σ0 = 4

N=32 1.668E-1 0.504 3.766E-2 1.111 9.239E-3 1.676 1.641E-2 1.581 2.776E-2 1.579

N=64 1.176E-1 0.490 1.743E-2 1.068 2.891E-3 1.678 5.483E-3 1.622 9.291E-3 1.625

N=128 8.377E-2 0.486 8.313E-3 1.040 9.036E-4 1.695 1.781E-3 1.653 3.013E-3 1.634

N=256 5.983E-2 0.486 4.043E-3 1.023 2.791E-4 1.717 5.662E-4 1.688 9.707E-4 1.680

N=512 4.271E-2 0.488 1.990E-3 1.012 8.491E-5 1.737 1.758E-4 1.715 3.030E-4 1.706

0.03

0.03

0.025

0.025

0.02

0.02

0.015

0.015

0.01

0.01

0.005

0.005

N=1024 3.045E-2 0.490 9.865E-4 1.007 2.548E-5 1.756 5.353E-5 1.739 9.287E-5 1.733

N=2048 2.168E-2 0.493 4.909E-4 1.004 7.546E-6 1.772 1.604E-5 1.760 2.793E-5 1.755

N=4096 1.541E-2 0.494 2.448E-4 1.002 2.209E-6 1.787 4.736E-6 1.777 8.275E-6 1.774

N=8192 1.094E-2 1.223E-4 6.399E-7 1.381E-6 2.420E-6

0

0 0

0.2

0.4

0.6

0.8

0

1

0.2

0.4

0.6

0.8

1

Figure 2: Values of the correction ϑ for problems (4.18) and (4.22) for N = 64 and ε = 10−6 .

Table 7: Values of the correction ϑ for problems (4.18) and (4.22) for N = 64 and ε = 10−6 . Mesh points x0 x1 x2 x3 x4 x5 x6 x7 x8 x9 ... ...

Problem (4.18) 0 8.332E-3 1.394E-2 1.773E-2 2.027E-2 2.199E-2 2.315E-2 2.393E-2 2.445E-2 2.480E-2 ... ...

Problem (4.22) 0 8.334E-3 1.395E-2 1.773E-2 2.028E-2 2.200E-2 2.315E-2 2.393E-2 2.446E-2 2.481E-2 ... ...

Mesh points ... ... x55 x56 x57 x58 x59 x60 x61 x62 x63 x64

22

Problem (4.18) ... ... 1.009E-2 9.106E-3 8.088E-3 7.038E-3 5.955E-3 4.837E-3 3.684E-3 2.494E-3 1.267E-3 0

Problem (4.22) ... ... 1.009E-2 9.106E-3 8.088E-3 7.038E-3 5.955E-3 4.837E-3 3.684E-3 2.494E-3 1.267E-3 0