- Email: [email protected]
Richardson extrapolation for a convection–diffusion problem using a Shishkin mesh
Applied Numerical Mathematics 45 (2003) 315–329 www.elsevier.com/locate/apnum
Richardson extrapolation for a convection–diffusion problem using a Shishkin mesh Maria Caridad Natividad a,∗ , Martin Stynes b a University of the Philippines, Diliman, Quezon City, Philippines 1101 b National University of Ireland, Cork, Ireland
Abstract We consider a convection–diffusion two-point boundary value problem on a piecewise-uniform Shishkin mesh, and show that when simple upwinding is used, a version of Richardson extrapolation improves the accuracy of the computed solution (measured in the discrete L∞ norm) from O(N −1 ln N) to O(N −2 ln2 N), where N + 1 mesh points are used. Numerical tests are provided to support these theoretical results. 2002 IMACS. Published by Elsevier Science B.V. All rights reserved.
1. Introduction Consider the linear convection–diffusion problem Lε uε (x) := −εuε (x) + a(x)uε (x) + b(x)uε (x) = f (x) (Pε ) uε (1) = u1 , uε (0) = u0 ,
for x ∈ (0, 1),
where a(·), b(·), f (·) ∈ C 2 [0, 1], a(·) > 2α > 0, 0 < ε 1 and u0 and u1 are given. Here α is some fixed positive constant. The solution uε of (Pε ) lies in C 4 [0, 1] and typically has a boundary layer at x = 1. Thus the problem (Pε ) is singularly perturbed. Many numerical methods (see, e.g., [2,9–11]) have been used to solve (Pε ). Most of these methods yield solutions that are first-order accurate (measured in the discrete L∞ norm); a few higher-order methods are described in [9], but these have complicated constructions and it is not known how to extend them to elliptic convection–diffusion problems in two or more variables. An exception is the defect correction method of [2], but its analysis seems fairly demanding. The aim of the present paper * Corresponding author. The work of M.C. Natividad was supported by the Center of Excellence grant awarded to the
Department of Mathematics by the Philippine Commission on Higher Education. E-mail addresses: [email protected] (M.C. Natividad), [email protected] (M. Stynes). 0168-9274/02/$30.00 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. doi:10.1016/S0168-9274(02)00212-X
316
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
is to construct and analyse a simple postprocessing technique that on a Shishkin mesh converts the almost-first-order convergence of the standard upwinding method into almost-second-order convergence. It seems feasible that this technique can be extended to elliptic convection–diffusion problems posed in more than one dimension. Richardson extrapolation is a well-known classical postprocessing procedure where two computed solutions approximating a particular quantity are averaged to provide a third and better approximation. It has been used [3,4] to improve the accuracy of computed solutions to non-singularly perturbed two-point boundary problems. Its application to the numerical solution of singularly perturbed reaction–diffusion problems using a Bakhvalov mesh was examined in [12]. In this paper we shall consider Richardson extrapolation applied to upwinding on a Shishkin mesh. This piecewise-uniform mesh is much simpler to construct than the exponentially graded Bakhvalov mesh, but (unlike the Bakhvalov mesh) has a certain nonuniformity in the local truncation error that complicates its analysis. We shall show how extrapolation reduces the nodal errors from O(N −1 ln N) to O(N −2 ln2 N), where N + 1 points are used in the mesh. The paper is organized as follows. Section 2 recalls pertinent properties of the solution uε of (Pε ), describes the Shishkin mesh and our upwind discretization of (Pε ), and gives some technical results that will be used later. The heart of our paper is Section 3, which shows how to construct a Richardson extrapolant of the computed solution and proves that this extrapolated solution is almost-second-order convergent in the discrete L∞ norm. Finally, Section 4 gives a numerical example that confirms the theoretical results. Throughout our analysis, C (sometimes subscripted) is a generic positive constant that is independent of ε, N , and the mesh. Note that an unsubscripted C may take different values in different places, but a subscripted C is a fixed constant that does not change throughout the paper. We say that aε,N = O(bε,N ) if for some C and some positive N0 independent of ε, we have |aε,N | C|bε,N | for all N N0 . Analogously, aε,N = o(bε,N ) means that given any C > 0, there exists N0 = N0 (C) such that |aε,N | C|bε,N | for all N N0 . Remark 1.1. In the analysis below we shall frequently assume that ε N −1 . If ε N −1 , then in practice the problem (Pε ) is not difficult to solve computationally, so the assumption is not restrictive. It can be replaced by the hypothesis that ε C0 N −1 for some fixed constant C0 without altering the conclusions of the paper.
2. Background material We begin with a basic stability result for Lε . Lemma 2.1. Let i 5 be a positive integer. Suppose that Lε ψ = φ on (0, 1), where φ ∈ C i [0, 1] satisfies (k) φ (x C 1 + ε −k−1 e−α(1−x)/ε for 0 k i and all x ∈ (0, 1) (2.1) and some constant C > 0. Suppose also that |ψ(0)| C and |ψ(1)| C. Then there exists a constant C such that (k) ψ (x) C 1 + ε −k e−α(1−x)/ε for 0 k i + 1 and all x ∈ [0, 1]; (2.2) ε
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
furthermore, (i+2) ψ (x) Cε −1 1 + ε −i−1 e−α(1−x)/ε ε
for all x ∈ [0, 1].
317
(2.3)
If instead of (2.1) one has the stronger assumption that (k) φ (x) Cε −k−1 e−α(1−x)/ε for 0 k i and all x ∈ (0, 1),
(2.4)
then (2.2) can be replaced by the conclusion (k) ψ (x) Cε −k e−α(1−x)/ε for 0 k i + 1 and all x ∈ [0, 1].
(2.5)
ε
Proof. Inequality (2.2) is proved in [5, Lemma 2.3]. Then (2.3) follows when Lε ψ = φ is differentiated i times, on solving for ψε(i+2) (x) and invoking (2.2). One obtains (2.5) easily from an inspection of the argument leading to (2.2). ✷ Linß [6] decomposes the solution uε of (Pε ) under minimal regularity hypotheses on the data a(·), b(·), f (·). His argument produces a splitting uε = vε + wε , where (k) v (x) C 1 + ε 3−k for all x ∈ [0, 1] and 0 k 4, (2.6a) ε (k) −k −α(1−x)/ε w (x) Cε e for all x ∈ [0, 1] and 0 k 4, (2.6b) ε and Lε vε = f, vε (0) = u0 , Lε wε = 0, wε (0) = 0,
v (1) = u1 − wε (1), ε wε (1) C.
(2.7a) (2.7b)
The functions vε and wε are called the smooth and layer parts of the solution uε . To solve (Pε ) numerically we shall use a piecewise-uniform Shishkin mesh [1,7,9] containing N + 1 points, where N 4 is an even positive integer. This mesh has a transition point 1 − τ , where 1 2ε τ = min , ln N . 2 α In fact we assume that 2ε ln N, (2.8) τ= α as otherwise N is very large compared to 1/ε, which means that the problem is easily solved numerically. Divide the subinterval [0, 1 − τ ] into N/2 equal subdivisions of width H , where 2(1 − τ ) . N Similarly divide the subinterval [1 − τ, 1] into N/2 equal subdivisions of width h, where H=
(2.9)
4ε ln N 2τ = . (2.10) N αN Then the Shishkin mesh is SN,τ = {xi }N i=0 , where x0 = 0, xN = 1, and the mesh widths hi := xi − xi−1 satisfy hi = H for i = 1, . . . , N/2 and hi = h for i = 1 + N/2, . . . , N . The notation SN,τ indicates that this piecewise-uniform mesh is entirely determined by the two user-chosen parameters N and τ . Define the standard first-order and second-order finite difference operators h=
318
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
UεN (xi ) − UεN (xi−1 ) , h i N 2 Uε (xi+1 ) − UεN (xi ) UεN (xi ) − UεN (xi−1 ) 2 N δ Uε (xi ) = − . hi+1 + hi hi+1 hi D − UεN (xi ) =
The upwind difference operator D − UεN (xi ) is used to approximate u (xi ) in (Pε ) since it is more stable than a standard central difference operator [1,7,9]. Now we form an upwind discretization on SN,τ of the problem (Pε ): Find the mesh function UεN such that N N N Pε Lε Uε (xi ) := −εδ 2 UεN (xi ) + a(xi )D − UεN (xi ) + b(xi )UεN (xi ) = f (xi ) for i = 1, . . . , N − 1, with UεN (0) = u0 , UεN (1) = u1 . The finite difference operator LN ε satisfies the following discrete maximum principle. N Lemma 2.2 (Discrete Maximum Principle). Let {φi }N i=0 and {ψi }i=0 be mesh functions. Assume that N N φ0 ψ0 and φN ψN . Then Lε φi Lε ψi for i = 1, . . . , N − 1 implies that φi ψi for all i.
Proof. See [7, pp. 60–61].
✷
If the conditions of Lemma 2.2 are satisfied, then we say that {φi } is a barrier function for {ψi }. The Lemma also holds true when applied to {xi : xi ∈ [1 − τ, 1]} mutatis mutandis. In particular, this Lemma implies that (PεN ) has a unique solution UεN . The error between the true solution uε and the computed solution UεN obeys the following error bound: Theorem 2.1. Let uε be the solution to (Pε ) and UεN the solution to (PεN ) computed on SN,τ . Then there exists a constant C such that N U (xi ) − uε (xi ) CN −1 ln N for all xi ∈ [0, 1]. ε
Proof. See [8].
✷
We close this Section with two technical Lemmas that will be needed in Section 3. For i = 0, 1, . . . , N , define on SN,τ the mesh functions i i
αhj αhj , Si = , 1+ 1+ Si = ε 2ε j =1 j =1 where by convention S0 = S0 = 1. Lemma 2.3. There exists a positive constant C1 such that for i = 1, . . . , N − 1, C1 C1 Si and LN Si . LN ε Si ε Si ε + αhi 2ε + αhi Moreover, there exists a positive constant C2 such that for i = N/2 + 1, . . . , N − 1, we have −1 LN ε Si C2 ε Si
and
−1 LN ε Si C2 ε Si .
(2.11)
(2.12)
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
319
Proof. Set Si = ij =1 (1 + βhj /ε), where β can be α or α/2. A calculation yields −2β 2 hi Si Si C Si
= + a(x )β + b(xi ) Si −2β 2 + a(xi )β S LN i ε i hi+1 + hi ε + βhi ε + βhi ε + βhi for some C, since a(·) > 2α 2β. This proves (2.11), and (2.12) then follows since h = O(ε) by (2.10). ✷ Lemma 2.4. Let SN,τ = {xi }N i=0 with hi = xi − xi−1 for all i. Then there exists a constant C such that N
−1
N
αhj −1 CN −1 , 1+ 2ε j =N/2+1
(2.13)
and N
αhj −1 CN −2 . 1+ ε j =N/2+1
(2.14)
Proof. First, N
2 ln N −N/2 αhj −1 = 1+ CN −1 , 1+ 2ε N j =N/2+1 by [11, Lemma 3.1]. Inequality (2.14) is proved similarly. Next, for 0 x 1, it is easy to verify that ln(1 + x) x − x 2 /4. Hence N
2 ln N ln2 N N αhj −1 2 ln N −N/2 − ln N + , = − ln 1 + 1+ = ln 1 + ln 2ε N 2 N 2N j =N/2+1 and the first inequality in (2.13) follows.
✷
3. Extrapolation of UεN The solution UεN (xi ) is computed on the mesh SN,τ . To improve its accuracy by extrapolation, we also solve (PεN ) on the Shishkin mesh S2N,τ . This mesh has 2N subintervals and the same transition point 1 − τ as SN,τ . The two meshes are nested, that is, SN,τ = {xi } ⊂ S2N,τ = {x˜i }. Clearly S2N,τ is obtained from SN,τ by bisecting each mesh interval, so on S2N,τ one has x˜i − x˜i−1 = H /2 for x˜i ∈ (0, 1 − τ ] and x˜i − x˜i−1 = h/2 for x˜i ∈ (1 − τ, 1]. ε2N . We denote the solution of the discrete problem (PεN ) on S2N,τ by U An inspection of the proof of Theorem 2.1 reveals that the ln N factor in this error estimate comes from the formula τ = (2ε ln N)/α of (2.8). Thus it is more illuminating to write the conclusion of this Theorem as UεN (xi ) − uε (xi ) = C3 N −1 ln N + RN (xi ) = C3 N −1 (ατ/2ε) + RN (xi )
for all xi ∈ SN,τ ,
(3.1)
320
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
ε2N is solved where C3 is some fixed constant and the remainder term RN is o(N −1 ln N). Recalling that U using the same transition point 1 − τ , it follows that 2N (x˜i ) for all x˜i ∈ S2N,τ , ε2N (x˜i ) − uε (x˜i ) = C3 (2N)−1 (ατ/2ε) + R U
(3.2)
2N is o(N −1 ln N). Eliminating the O(N −1 ) terms from the right-hand of (3.1) and where the remainder R (3.2) gives 2N ε (xi ) − UεN (xi ) = RN (xi ) − 2R 2N (xi ) = o N −1 ln N for all xi ∈ SN,τ . (3.3) uε (xi ) − 2U We shall therefore use the simple extrapolation formula ε2N (xi ) − UεN (xi ) 2U
for xi ∈ SN,τ ,
(3.4)
ε2N (xi ). which we expect will yield an approximation of uε (xi ) more accurate than either UεN (xi ) or U N N Decompose the solution Uε of (Pε ) on SN,τ into smooth and layer mesh functions analogously to (2.7) by setting UεN = VεN + WεN , where N LN ε Vε = f, N LN ε Wε
= 0,
VεN (0) = vε (0), WεN (0)
= wε (0),
VεN (1) = vε (1), WεN (1)
= wε (1).
(3.5a) (3.5b)
ε2N + W ε2N . Then UεN − uε = (VεN − vε ) + (WεN − wε ) and U ε2N − uε = ε2N = V Similarly one can set U 2N 2N ε − wε ). ε − vε ) + (W (V 3.1. Extrapolation of VεN Consider first the error VεN (xi ) − vε (xi ) in the smooth part of uε . Set ζ(x) = a(x)vε (x)/2 for all x ∈ [0, 1]. Lemma 3.1 (Bound on truncation error of vε ). Assume that ε N −1 . Then for all xi ∈ (0, 1), 2 N LN ε Vε − vε (xi ) = ζ(xi )hi + O H . Proof. Recall from (2.6a) that |vε(3)(·)| C. It is straightforward to deduce the result from this bound, some Taylor expansions (cf. [7, p. 21]), and ε N −1 H . ✷ Following the classical approach of Keller [3,4], define the function E(x) to be the solution of the boundary value problem Lε E(x) = ζ(x) for all x ∈ (0, 1), E(0) = E(1) = 0.
(3.6a) (3.6b)
Now a ∈ C 2 [0, 1], and vε ∈ C 4 [0, 1] has its derivatives bounded as in (2.6a), so ζ ∈ C 2 [0, 1] with |ζ(x)| C and |ζ (x)| C for all x ∈ [0, 1]. By [6, Theorem 1], one can decompose E as E = η + ϑ , where η and ϑ are the smooth and layer parts of E, and for all x ∈ (0, 1), (k) η (x) C 1 + ε 2−k for 0 k 3, (3.7a) (k) −k −α(1−x)/ε ϑ (x) Cε e for 0 k 3, (3.7b) with
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
Lε η(x) = ζ(x) η(0) = ϑ(0) = 0
and and
Lε ϑ(x) = 0 for all x ∈ (0, 1), η(1) = −ϑ(1).
321
(3.8a) (3.8b)
We are now ready to show how extrapolation works for the smooth part vε . Let us first consider the interval [0, 1 − τ ]. Lemma 3.2. Assume that ε N −1 . Then VεN (xi ) − vε (xi ) = H E(xi ) + O N −2 for all xi ∈ [0, 1 − τ ]. Proof. Fix xi ∈ (0, 1). By (3.7a), the truncation error formulas of [7, p. 21], and the inequality hj H for all j , we get N LN ε η(xi ) = Lε η(xi ) + Lε η(xi ) − Lε η(xi ) = Lε η(xi ) + O(H ).
Hence, recalling (3.8a),
2 H LN ε η(xi ) = H ζ(xi ) + O H .
Consequently Lemma 3.1 yields LN ε
VεN
− vε − H η (xi ) =
O H2
for xi ∈ (0, 1 − τ ], 2 for xi ∈ (1 − τ, 1), (h − H )ζ(xi ) + O H 2 for xi ∈ (0, 1 − τ ], O H = O(H ) for xi ∈ (1 − τ, 1).
Define the mesh function Zi = C3 N −2 (1 + xi ) + H Si
N
j =1
αhj 1+ 2ε
−1
(3.9)
for i = 0, . . . , N,
where the constant C3 will be specified later. Then N
αhj −1 N N −2 Lε Si 1+ for i = 1, . . . , N − 1. Lε Zi C3 2αN + H 2ε j =1 From Lemmas 2.3 and 2.4 and ε N −1 , for 0 < i N/2 it follows that −2 C3 αH 2 /2, LN ε Zi 2C3 αN
(3.10)
and for N/2 < i < N , LN ε Zi Now
C3 C2 H ε
−1
N
αhj −1 C3 C2 H. 1+ 2ε j =N/2+1
(3.11)
Z0 0 = VεN (0) − vε (0) − H η(0)
(3.12)
by (3.5a) and (3.8b). Also ZN C3 H V N (1) − vε (1) − H η(1)
(3.13)
ε
322
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
by (3.5a) and (3.7b), provided that C3 is sufficiently large (independently of ε and N ). From (3.9)–(3.13), one can choose the constant C3 so that Zi is a barrier function for ±[VεN (xi ) − vε (xi ) − H η(xi )], and Lemma 2.2 yields (3.14) Zi VεN (xi ) − vε (xi ) − H η(xi ) for all i. But for i = 1, . . . , N/2, Lemma 2.4 shows that N
αhj −1 −2 1+ C N −2 + H N −1 CN −2 . Zi C3 2N + H 2ε j =i+1
(3.15)
Thus for i = 1, . . . , N/2, N V (xi ) − vε (xi ) − H E(xi ) V N (xi ) − vε (xi ) − H η(xi ) + H ϑ(xi ) ε ε Zi + 2N −1 ϑ(xi ) CN −2 , where we used E = η + ϑ , (3.14), (3.15) and (3.7b).
✷
It is now a simple matter to show that the extrapolation formula (3.4) improves the accuracy of VεN in (0, 1 − τ ]. Lemma 3.3. Assume that ε N −1 . For all xi ∈ [0, 1 − τ ], 2N vε (xi ) − 2V ε (xi ) − VεN (xi ) CN −2 for some constant C. Proof. Let xi ∈ [0, 1 − τ ]. Since the subinterval meshwidths of S2N,τ are half those of SN,τ and the function E(·) depends only on τ , Lemma 3.2 implies that ε2N (xi ) − vε (xi ) = H E(xi ) + O N −2 . V 2 This relation and Lemma 3.2 yield 2N 2N ε (xi ) − VεN (xi ) = −2 V ε − vε (xi ) + VεN − vε (xi ) = O N −2 . ✷ vε (xi ) − 2V Next, consider the approximation of vε (xi ) in (1 − τ, 1). Lemma 3.4. Assume that ε N −1 . For all xi ∈ [1 − τ, 1], 2N vε (xi ) − 2V ε (xi ) − VεN (xi ) CN −2 ln N for some constant C. Proof. Define the function G(x) on [1 − τ, 1] by G(1 − τ ) = 1, G(1) = 0, and Lε G(x) = 0 for 1 − τ < x < 1. N We work initially on the mesh SN,τ . Define a discrete approximation GN ε of G by Gε (1 − τ ) = 1, N N N Gε (1) = 0, and Lε Gε (xi ) = 0 for N/2 < i < N . From a standard convergence analysis [8] for the upwind method, G(xi ) − GN (xi ) CN −1 ln N for N/2 i N. (3.16) ε
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
323
Set
Φ(xi ) = VεN (xi ) − vε (xi ) − H E(1 − τ ) GN ε (xi )
for N/2 i N.
−2
2 −2 ln N) Then Φ(xN ) = 0, Φ(xN/2 ) = O(N ) by Lemma 3.2, and LN ε Φ(xi ) = O(hi ) + O(H ) = O(N −1 −2 by Lemma 3.1 and ε N . Consequently one can use the barrier function C(1 + xi )N ln N to get Φ(xi ) CN −2 ln N for N/2 i N. (3.17)
Combining (3.16) and (3.17) and observing that |E(1 − τ )| C is implied by (3.7a) and (3.7b), we have shown that on the mesh SN,τ , for N/2 i N. (3.18) VεN (xi ) − vε (xi ) = H E(1 − τ ) G(xi ) + O N −2 ln N Similarly, on the mesh S2N,τ one gets V˜ε2N (x˜i ) − vε (x˜i ) = (H /2)E(1 − τ ) G(x˜i ) + O N −2 ln N
for N i 2N.
(3.19)
But x˜2i = xi for N/2 i N . From (3.18) and (3.19), we have 2N ε (xi ) − VεN (xi ) = 2 vε − V ε2N (xi ) − vε − VεN (xi ) = O N −2 ln N vε (xi ) − 2V for N/2 i N .
✷
3.2. Extrapolation of WεN Now we consider the error WεN (xi ) − wε (xi ). The two subintervals [0, 1 − τ ] and [1 − τ, 1] will be analysed separately. First, consider [0, 1 − τ ]. Lemma 3.5. For all xi ∈ [0, 1 − τ ], 2N wε (xi ) − 2W ε (xi ) − WεN (xi ) CN −2 for some constant C. Proof. Let xi ∈ [0, 1 − τ ]. With only slight modifications, the argument leading to [10, (3.2)] shows that ε2N )(xi )| CN −2 . Hence |(wε − WεN )(xi )| CN −2 . Similarly |(wε − W 2N wε (xi ) − 2W ε (xi ) − WεN (xi ) = 2 wε − W ε2N (xi ) − wε − WεN (xi ) CN −2 . ✷ To analyse the effect of extrapolation on (1 − τ, 1), define the function F on [1 − τ, 1] by the boundary value problem 2 Lε F (x) = − εa(x)wε (x) α F (1 − τ ) = F (1) = 0.
for x ∈ (1 − τ, 1),
(3.20a) (3.20b)
Note that F depends on τ but is otherwise independent of N . Lemma 3.6. For xi ∈ [1 − τ, 1], we have WεN (xi ) − wε (xi ) = N −1 ln N F (xi ) + O N −2 ln2 N .
(3.21)
324
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
Proof. Recalling (2.6b) and applying Lemma 2.1 to (3.20), for 0 k 3 and x ∈ (0, 1) we have (k) F (x) Cε −k e−α(1−x)/ε . (3.22) Let xi ∈ (1 − τ, 1). By a Taylor expansion of LN ε F (xi ) about x = xi , and an appeal to (2.10) and (3.22), we get −1 −α(1−x )/ε −1 i+1 N ln N . (3.23) LN ε F (xi ) = Lε F (xi ) + O ε e Recalling (2.7b) and (3.5b), more Taylor expansions show that, for some x˜1 ∈ (xi , xi+1 ) and x˜2 , x˜3 ∈ (xi−1 , xi ), N ε (4) 1 1 wε (x˜1 ) + wε(4) (x˜2 ) h2 − a(xi )wε (xi )h + a(xi )wε (x˜3 )h2 LN ε Wε − wε (xi ) = − 4! 2 3! −1 −α(1−x )/ε −2 2 2 −1 i+1 = − εN ln N a(xi )wε (xi ) + O ε e N ln N , (3.24) α by (2.6b) and (2.10). From (3.23), (3.20) and (3.24) one sees that −1 −1 ln N Lε F (xi ) + O ε −1 e−α(1−xi+1 )/ε N −2 ln2 N N ln N LN ε F (xi ) = N −1 −α(1−x )/ε −2 2 N i+1 N ln N . = LN ε Wε − wε (xi ) + O ε e That is, for all xi ∈ (1 − τ, 1), N N L W − wε − N −1 ln N F (xi ) C4 ε −1 e−α(1−xi+1 )/ε N −2 ln2 N, ε ε
(3.25)
for some constant C4 . For i = N/2, . . . , N , consider the mesh function N −2 2
αhj −1 −2 1+ , Γi := C5 N (1 + xi ) + N ln N Si ε j =1 where the constant C5 will be chosen later. Invoking (2.12), we see that N −1 −2 2
αhj −1 N . 1+ Lε Γi C2 C5 ε N ln N ε j =i+1
(3.26)
But for each i N/2, e
−α(1−xi+1 )/ε
= e
N
αhi+1 /ε
e−αhj /ε
j =i+1
e
4N −1 ln N
N
j =i+1
αhj 1+ ε
−1
N
αhj −1 e . 1+ ε j =i+1
If C5 e4 C4 /C2 , then by (3.25)–(3.27), LN Γi LN W N − wε − N −1 ln N F (xi ) ε
ε
ε
4
(3.27)
for xi ∈ (1 − τ, 1).
From (3.5b) and (3.20b), ΓN 0 = |WεN (1) − wε (1) − (N −1 ln N)F (1)|. The proof of Lemma 3.5 shows that |(wε − WεN )(1 − τ )| C6 N −2 for some constant C6 ; hence, recalling (3.20b), ΓN/2 C5 N −2 |(WεN − wε − (N −1 ln N)F )(1 − τ )| provided that C5 C6 .
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
325
Choose C5 = max{e4 C4 /C2 , C6 }. Then the above calculations establish that Γi is a barrier function for ±[WεN (xi ) − wε (xi ) − (N −1 ln N)F (xi )] when 1 − τ xi 1. Applying Lemma 2.2 on [1 − τ, 1] yields N W (xi ) − wε (xi ) − N −1 ln N F (xi ) Γi C5 N −2 ln2 N ε for i = N/2, . . . , N , as desired.
✷
Now we demonstrate how extrapolation improves the accuracy of WεN (xi ) on [1 − τ, 1]. Lemma 3.7. For all xi ∈ [1 − τ, 1], 2N wε (xi ) − 2W ε (xi ) − WεN (xi ) CN −2 ln2 N for some constant C. Proof. Let xi ∈ [1 − τ, 1]. An inspection of the proof of Lemma 3.6 shows that the factor ln N in (3.21) comes from the definition of h in (2.10). Consequently, rewriting (3.21) so as to make explicit its dependence on the choice of N and τ , one has WεN (xi ) − wε (xi ) = N −1 (τ α/2ε)F (xi ) + O N −2 (τ α/2ε)2 . ε2N is computed on the mesh S2N,τ which has the same transition point 1 − τ as SN,τ , it follows that As W ε2N (xi ) − wε (xi ) = (2N)−1 (τ α/2ε)F (xi ) + O N −2 (τ α/2ε)2 . W Combining these two formulas, we obtain 2N ε (xi ) − WεN (xi ) = 2 wε − W ε2N (xi ) − wε − WεN (xi ) = O N −2 ln2 N . wε (xi ) − 2W
✷
We come finally to our main result. Theorem 3.1 (error after extrapolation). Assume that ε N −1 . There exists a constant C such that 2N uε (xi ) − 2U ε (xi ) − UεN (xi ) CN −2 ln2 N for all xi ∈ SN,τ . Proof. For each i, 2N ε (xi ) − UεN (xi ) uε (xi ) − 2U 2N 2N ε (xi ) − VεN (xi ) + wε (xi ) − 2W ε (xi ) − WεN (xi ) . = vε (xi ) − 2V Combining the results of Lemmas 3.3, 3.4, 3.5 and 3.7, we are done.
✷
Remark 3.1. One could continue the process analysed in this paper to obtain still more accuracy by repeated extrapolation, as is done in classical Richardson extrapolation. To be precise, if a(·), b(·), f (·) ∈ C m [0, 1] for some m 2, then one can extrapolate m − 1 times. This means: choose a Shishkin mesh SN,τ as in Section 2 but with τ = mεα −1 ln N , then, setting UεN,0 = UεN , compute UεN,j (xi ) =
ε2N,j −1 (xi ) − UεN,j −1 (xi ) 2j U 2j −1 − 1
for j = 1, . . . , m − 1 and i = 0, 1, . . . , N.
326
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
ε2N,j −1 is computed on a mesh with the same transition point as UεN,j −1 but having twice as (Here U many mesh points.) The above choice of τ yields N −m in Lemma 3.5. The final result, which generalizes Theorem 3.1, is that uε (xi ) − U N,j (xi ) CN −m lnm N for all xi ∈ SN,τ . ε
Table 1 Maximum nodal errors before and after extrapolation ε
N
64
128
256
512
1024
10−1
before extrapolation rate after extrapolation rate
2.05e–02 1.21 7.91e–04 2.44
1.07e–02 1.21 2.12e–04 2.41
5.44e–03 1.19 5.51e–05 2.38
2.75e–03 1.17 1.40e–05 2.34
1.38e–03 1.16 3.54e–06 2.31
10−2
before extrapolation rate after extrapolation rate
1.55e–01 0.82 2.63e–02 1.15
9.95e–02 0.91 1.41e–02 1.97
5.99e–02 1.08 4.70e–03 2.16
3.22e–02 1.11 1.35e–03 2.21
1.67e–02 1.13 3.69e–04 2.25
10−3
before extrapolation rate after extrapolation rate
1.50e–01 0.82 2.58e–02 1.14
9.63e–02 0.74 1.39e–02 1.66
6.38e–02 0.84 5.49e–03 1.72
3.92e–02 0.92 2.04e–03 1.85
2.80e–02 0.95 6.90e–04 1.91
10−4
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.89e–04 1.91
10−5
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.89e–04 1.91
10−6
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.89e–04 1.91
10−7
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.89e–04 1.98
10−8
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.87e–04 1.83
10−9
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.49e–03 1.70
3.91e–02 0.92 2.07e–03 2.16
2.27e–02 0.96 5.81e–04 2.00
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
327
We do not present the general theory for this result as the extra notation involved would obscure the central ideas that we have described here.
4. Numerical example Consider the following singularly perturbed convection–diffusion problem in one dimension: −εuε (x) + 1 + x(1 − x) uε (x) = f (x) for x ∈ (0, 1), uε (0) = uε (1) = 0, where f (x) is chosen in such a way that the real solution uε is 1 − e−(1−x)/ε πx . − cos −1/ε 1−e 2 This function displays typical boundary-layer behaviour. Our theory assumes that ε N −1 , but this inequality is not always satisfied in the numerical experiments. We take α = 0.199. Table 1 shows the maximum nodal errors between the real and computed solutions, before and after extrapolation, for a range of values of ε and N . The rates of convergence r are computed in the standard way, assuming that convergence of the form N −r lnr N is attained; to be precise, if Eε,N is the maximum nodal error for given values of ε and N , then the associated rate in Table 1 is (ln Eε,N − ln Eε,2N )/ln(2 ln N/ln 2N ). Figs. 1 and 2 give these errors in graphical form, together with benchmark curves for certain rates of convergence to show that the computed errors decrease at approximately these rates. uε (x) =
Fig. 1. Comparison of nodal errors before and after extrapolation for ε = 10−2 .
328
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
Fig. 2. Comparison of nodal errors before and after extrapolation for ε = 10−8 .
The results show that extrapolation decreases the nodal errors significantly and increases the rate of convergence of the numerical method from O(N −1 ln N) to O(N −2 ln2 N), which supports the theoretical bounds of Theorems 2.1 and 3.1. Remark 4.1. Richardson extrapolation on Shishkin meshes is easily implemented. To compute Uε2N we must solve a linear system with 2N equations, which costs approximately four times as much as solving the original linear system for UεN . Thus the rate of convergence of the simple upwind method is raised from almost-first-order to almost-second-order at the cost of five times as much computation; this tradeoff between cost and accuracy is clearly worthwhile.
Acknowledgement We are grateful to two unknown referees who read the original version of the paper carefully and made constructive suggestions for its improvement.
References [1] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman & Hall/CRC, Boca Raton, FL, 2000. [2] A. Fröhner, H-G. Roos, The ε-uniform convergence of a defect correction method on a Shishkin mesh, Appl. Numer. Math. 37 (2001) 79–94.
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
329
[3] H.B. Keller, Accurate difference methods for linear ordinary differential systems subject to linear constraints, SIAM J. Numer. Anal. 6 (1969) 8–30. [4] H.B. Keller, Numerical Methods for Two-Point Boundary Value Problems, Dover, New York, 1992. [5] R.B. Kellogg, A. Tsan, Analysis of some difference approximations for a singular perturbation problem without turning points, Math. Comp. 32 (1978) 1025–1039. [6] T. Linß, Solution decompositions for linear convection–diffusion problems, Z. Anal. Anwendungen 21 (2002) 209–214. [7] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimension, World Scientific, Singapore, 1996. [8] H.-G. Roos, T. Linß, Sufficient conditions for uniform convergence on layer-adapted grids, Computing 63 (1999) 27–45. [9] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations: Convection– Diffusion and Flow Problems, Springer, Berlin, 1996. [10] M. Stynes, H.-G. Roos, The midpoint upwind scheme, Appl. Numer. Math. 23 (1997) 361–374. [11] M. Stynes, L. Tobiska, A finite difference analysis of a streamline diffusion method on a Shishkin mesh, Numer. Algorithms 18 (1998) 337–360. [12] R. Vulanovi´c, D. Herceg, N. Petrovi´c, On the extrapolation of a singularly perturbed boundary value problem, Computing 36 (1986) 69–79.
Richardson extrapolation for a convection–diffusion problem using a Shishkin mesh Maria Caridad Natividad a,∗ , Martin Stynes b a University of the Philippines, Diliman, Quezon City, Philippines 1101 b National University of Ireland, Cork, Ireland
Abstract We consider a convection–diffusion two-point boundary value problem on a piecewise-uniform Shishkin mesh, and show that when simple upwinding is used, a version of Richardson extrapolation improves the accuracy of the computed solution (measured in the discrete L∞ norm) from O(N −1 ln N) to O(N −2 ln2 N), where N + 1 mesh points are used. Numerical tests are provided to support these theoretical results. 2002 IMACS. Published by Elsevier Science B.V. All rights reserved.
1. Introduction Consider the linear convection–diffusion problem Lε uε (x) := −εuε (x) + a(x)uε (x) + b(x)uε (x) = f (x) (Pε ) uε (1) = u1 , uε (0) = u0 ,
for x ∈ (0, 1),
where a(·), b(·), f (·) ∈ C 2 [0, 1], a(·) > 2α > 0, 0 < ε 1 and u0 and u1 are given. Here α is some fixed positive constant. The solution uε of (Pε ) lies in C 4 [0, 1] and typically has a boundary layer at x = 1. Thus the problem (Pε ) is singularly perturbed. Many numerical methods (see, e.g., [2,9–11]) have been used to solve (Pε ). Most of these methods yield solutions that are first-order accurate (measured in the discrete L∞ norm); a few higher-order methods are described in [9], but these have complicated constructions and it is not known how to extend them to elliptic convection–diffusion problems in two or more variables. An exception is the defect correction method of [2], but its analysis seems fairly demanding. The aim of the present paper * Corresponding author. The work of M.C. Natividad was supported by the Center of Excellence grant awarded to the
Department of Mathematics by the Philippine Commission on Higher Education. E-mail addresses: [email protected] (M.C. Natividad), [email protected] (M. Stynes). 0168-9274/02/$30.00 2002 IMACS. Published by Elsevier Science B.V. All rights reserved. doi:10.1016/S0168-9274(02)00212-X
316
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
is to construct and analyse a simple postprocessing technique that on a Shishkin mesh converts the almost-first-order convergence of the standard upwinding method into almost-second-order convergence. It seems feasible that this technique can be extended to elliptic convection–diffusion problems posed in more than one dimension. Richardson extrapolation is a well-known classical postprocessing procedure where two computed solutions approximating a particular quantity are averaged to provide a third and better approximation. It has been used [3,4] to improve the accuracy of computed solutions to non-singularly perturbed two-point boundary problems. Its application to the numerical solution of singularly perturbed reaction–diffusion problems using a Bakhvalov mesh was examined in [12]. In this paper we shall consider Richardson extrapolation applied to upwinding on a Shishkin mesh. This piecewise-uniform mesh is much simpler to construct than the exponentially graded Bakhvalov mesh, but (unlike the Bakhvalov mesh) has a certain nonuniformity in the local truncation error that complicates its analysis. We shall show how extrapolation reduces the nodal errors from O(N −1 ln N) to O(N −2 ln2 N), where N + 1 points are used in the mesh. The paper is organized as follows. Section 2 recalls pertinent properties of the solution uε of (Pε ), describes the Shishkin mesh and our upwind discretization of (Pε ), and gives some technical results that will be used later. The heart of our paper is Section 3, which shows how to construct a Richardson extrapolant of the computed solution and proves that this extrapolated solution is almost-second-order convergent in the discrete L∞ norm. Finally, Section 4 gives a numerical example that confirms the theoretical results. Throughout our analysis, C (sometimes subscripted) is a generic positive constant that is independent of ε, N , and the mesh. Note that an unsubscripted C may take different values in different places, but a subscripted C is a fixed constant that does not change throughout the paper. We say that aε,N = O(bε,N ) if for some C and some positive N0 independent of ε, we have |aε,N | C|bε,N | for all N N0 . Analogously, aε,N = o(bε,N ) means that given any C > 0, there exists N0 = N0 (C) such that |aε,N | C|bε,N | for all N N0 . Remark 1.1. In the analysis below we shall frequently assume that ε N −1 . If ε N −1 , then in practice the problem (Pε ) is not difficult to solve computationally, so the assumption is not restrictive. It can be replaced by the hypothesis that ε C0 N −1 for some fixed constant C0 without altering the conclusions of the paper.
2. Background material We begin with a basic stability result for Lε . Lemma 2.1. Let i 5 be a positive integer. Suppose that Lε ψ = φ on (0, 1), where φ ∈ C i [0, 1] satisfies (k) φ (x C 1 + ε −k−1 e−α(1−x)/ε for 0 k i and all x ∈ (0, 1) (2.1) and some constant C > 0. Suppose also that |ψ(0)| C and |ψ(1)| C. Then there exists a constant C such that (k) ψ (x) C 1 + ε −k e−α(1−x)/ε for 0 k i + 1 and all x ∈ [0, 1]; (2.2) ε
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
furthermore, (i+2) ψ (x) Cε −1 1 + ε −i−1 e−α(1−x)/ε ε
for all x ∈ [0, 1].
317
(2.3)
If instead of (2.1) one has the stronger assumption that (k) φ (x) Cε −k−1 e−α(1−x)/ε for 0 k i and all x ∈ (0, 1),
(2.4)
then (2.2) can be replaced by the conclusion (k) ψ (x) Cε −k e−α(1−x)/ε for 0 k i + 1 and all x ∈ [0, 1].
(2.5)
ε
Proof. Inequality (2.2) is proved in [5, Lemma 2.3]. Then (2.3) follows when Lε ψ = φ is differentiated i times, on solving for ψε(i+2) (x) and invoking (2.2). One obtains (2.5) easily from an inspection of the argument leading to (2.2). ✷ Linß [6] decomposes the solution uε of (Pε ) under minimal regularity hypotheses on the data a(·), b(·), f (·). His argument produces a splitting uε = vε + wε , where (k) v (x) C 1 + ε 3−k for all x ∈ [0, 1] and 0 k 4, (2.6a) ε (k) −k −α(1−x)/ε w (x) Cε e for all x ∈ [0, 1] and 0 k 4, (2.6b) ε and Lε vε = f, vε (0) = u0 , Lε wε = 0, wε (0) = 0,
v (1) = u1 − wε (1), ε wε (1) C.
(2.7a) (2.7b)
The functions vε and wε are called the smooth and layer parts of the solution uε . To solve (Pε ) numerically we shall use a piecewise-uniform Shishkin mesh [1,7,9] containing N + 1 points, where N 4 is an even positive integer. This mesh has a transition point 1 − τ , where 1 2ε τ = min , ln N . 2 α In fact we assume that 2ε ln N, (2.8) τ= α as otherwise N is very large compared to 1/ε, which means that the problem is easily solved numerically. Divide the subinterval [0, 1 − τ ] into N/2 equal subdivisions of width H , where 2(1 − τ ) . N Similarly divide the subinterval [1 − τ, 1] into N/2 equal subdivisions of width h, where H=
(2.9)
4ε ln N 2τ = . (2.10) N αN Then the Shishkin mesh is SN,τ = {xi }N i=0 , where x0 = 0, xN = 1, and the mesh widths hi := xi − xi−1 satisfy hi = H for i = 1, . . . , N/2 and hi = h for i = 1 + N/2, . . . , N . The notation SN,τ indicates that this piecewise-uniform mesh is entirely determined by the two user-chosen parameters N and τ . Define the standard first-order and second-order finite difference operators h=
318
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
UεN (xi ) − UεN (xi−1 ) , h i N 2 Uε (xi+1 ) − UεN (xi ) UεN (xi ) − UεN (xi−1 ) 2 N δ Uε (xi ) = − . hi+1 + hi hi+1 hi D − UεN (xi ) =
The upwind difference operator D − UεN (xi ) is used to approximate u (xi ) in (Pε ) since it is more stable than a standard central difference operator [1,7,9]. Now we form an upwind discretization on SN,τ of the problem (Pε ): Find the mesh function UεN such that N N N Pε Lε Uε (xi ) := −εδ 2 UεN (xi ) + a(xi )D − UεN (xi ) + b(xi )UεN (xi ) = f (xi ) for i = 1, . . . , N − 1, with UεN (0) = u0 , UεN (1) = u1 . The finite difference operator LN ε satisfies the following discrete maximum principle. N Lemma 2.2 (Discrete Maximum Principle). Let {φi }N i=0 and {ψi }i=0 be mesh functions. Assume that N N φ0 ψ0 and φN ψN . Then Lε φi Lε ψi for i = 1, . . . , N − 1 implies that φi ψi for all i.
Proof. See [7, pp. 60–61].
✷
If the conditions of Lemma 2.2 are satisfied, then we say that {φi } is a barrier function for {ψi }. The Lemma also holds true when applied to {xi : xi ∈ [1 − τ, 1]} mutatis mutandis. In particular, this Lemma implies that (PεN ) has a unique solution UεN . The error between the true solution uε and the computed solution UεN obeys the following error bound: Theorem 2.1. Let uε be the solution to (Pε ) and UεN the solution to (PεN ) computed on SN,τ . Then there exists a constant C such that N U (xi ) − uε (xi ) CN −1 ln N for all xi ∈ [0, 1]. ε
Proof. See [8].
✷
We close this Section with two technical Lemmas that will be needed in Section 3. For i = 0, 1, . . . , N , define on SN,τ the mesh functions i i
αhj αhj , Si = , 1+ 1+ Si = ε 2ε j =1 j =1 where by convention S0 = S0 = 1. Lemma 2.3. There exists a positive constant C1 such that for i = 1, . . . , N − 1, C1 C1 Si and LN Si . LN ε Si ε Si ε + αhi 2ε + αhi Moreover, there exists a positive constant C2 such that for i = N/2 + 1, . . . , N − 1, we have −1 LN ε Si C2 ε Si
and
−1 LN ε Si C2 ε Si .
(2.11)
(2.12)
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
319
Proof. Set Si = ij =1 (1 + βhj /ε), where β can be α or α/2. A calculation yields −2β 2 hi Si Si C Si
= + a(x )β + b(xi ) Si −2β 2 + a(xi )β S LN i ε i hi+1 + hi ε + βhi ε + βhi ε + βhi for some C, since a(·) > 2α 2β. This proves (2.11), and (2.12) then follows since h = O(ε) by (2.10). ✷ Lemma 2.4. Let SN,τ = {xi }N i=0 with hi = xi − xi−1 for all i. Then there exists a constant C such that N
−1
N
αhj −1 CN −1 , 1+ 2ε j =N/2+1
(2.13)
and N
αhj −1 CN −2 . 1+ ε j =N/2+1
(2.14)
Proof. First, N
2 ln N −N/2 αhj −1 = 1+ CN −1 , 1+ 2ε N j =N/2+1 by [11, Lemma 3.1]. Inequality (2.14) is proved similarly. Next, for 0 x 1, it is easy to verify that ln(1 + x) x − x 2 /4. Hence N
2 ln N ln2 N N αhj −1 2 ln N −N/2 − ln N + , = − ln 1 + 1+ = ln 1 + ln 2ε N 2 N 2N j =N/2+1 and the first inequality in (2.13) follows.
✷
3. Extrapolation of UεN The solution UεN (xi ) is computed on the mesh SN,τ . To improve its accuracy by extrapolation, we also solve (PεN ) on the Shishkin mesh S2N,τ . This mesh has 2N subintervals and the same transition point 1 − τ as SN,τ . The two meshes are nested, that is, SN,τ = {xi } ⊂ S2N,τ = {x˜i }. Clearly S2N,τ is obtained from SN,τ by bisecting each mesh interval, so on S2N,τ one has x˜i − x˜i−1 = H /2 for x˜i ∈ (0, 1 − τ ] and x˜i − x˜i−1 = h/2 for x˜i ∈ (1 − τ, 1]. ε2N . We denote the solution of the discrete problem (PεN ) on S2N,τ by U An inspection of the proof of Theorem 2.1 reveals that the ln N factor in this error estimate comes from the formula τ = (2ε ln N)/α of (2.8). Thus it is more illuminating to write the conclusion of this Theorem as UεN (xi ) − uε (xi ) = C3 N −1 ln N + RN (xi ) = C3 N −1 (ατ/2ε) + RN (xi )
for all xi ∈ SN,τ ,
(3.1)
320
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
ε2N is solved where C3 is some fixed constant and the remainder term RN is o(N −1 ln N). Recalling that U using the same transition point 1 − τ , it follows that 2N (x˜i ) for all x˜i ∈ S2N,τ , ε2N (x˜i ) − uε (x˜i ) = C3 (2N)−1 (ατ/2ε) + R U
(3.2)
2N is o(N −1 ln N). Eliminating the O(N −1 ) terms from the right-hand of (3.1) and where the remainder R (3.2) gives 2N ε (xi ) − UεN (xi ) = RN (xi ) − 2R 2N (xi ) = o N −1 ln N for all xi ∈ SN,τ . (3.3) uε (xi ) − 2U We shall therefore use the simple extrapolation formula ε2N (xi ) − UεN (xi ) 2U
for xi ∈ SN,τ ,
(3.4)
ε2N (xi ). which we expect will yield an approximation of uε (xi ) more accurate than either UεN (xi ) or U N N Decompose the solution Uε of (Pε ) on SN,τ into smooth and layer mesh functions analogously to (2.7) by setting UεN = VεN + WεN , where N LN ε Vε = f, N LN ε Wε
= 0,
VεN (0) = vε (0), WεN (0)
= wε (0),
VεN (1) = vε (1), WεN (1)
= wε (1).
(3.5a) (3.5b)
ε2N + W ε2N . Then UεN − uε = (VεN − vε ) + (WεN − wε ) and U ε2N − uε = ε2N = V Similarly one can set U 2N 2N ε − wε ). ε − vε ) + (W (V 3.1. Extrapolation of VεN Consider first the error VεN (xi ) − vε (xi ) in the smooth part of uε . Set ζ(x) = a(x)vε (x)/2 for all x ∈ [0, 1]. Lemma 3.1 (Bound on truncation error of vε ). Assume that ε N −1 . Then for all xi ∈ (0, 1), 2 N LN ε Vε − vε (xi ) = ζ(xi )hi + O H . Proof. Recall from (2.6a) that |vε(3)(·)| C. It is straightforward to deduce the result from this bound, some Taylor expansions (cf. [7, p. 21]), and ε N −1 H . ✷ Following the classical approach of Keller [3,4], define the function E(x) to be the solution of the boundary value problem Lε E(x) = ζ(x) for all x ∈ (0, 1), E(0) = E(1) = 0.
(3.6a) (3.6b)
Now a ∈ C 2 [0, 1], and vε ∈ C 4 [0, 1] has its derivatives bounded as in (2.6a), so ζ ∈ C 2 [0, 1] with |ζ(x)| C and |ζ (x)| C for all x ∈ [0, 1]. By [6, Theorem 1], one can decompose E as E = η + ϑ , where η and ϑ are the smooth and layer parts of E, and for all x ∈ (0, 1), (k) η (x) C 1 + ε 2−k for 0 k 3, (3.7a) (k) −k −α(1−x)/ε ϑ (x) Cε e for 0 k 3, (3.7b) with
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
Lε η(x) = ζ(x) η(0) = ϑ(0) = 0
and and
Lε ϑ(x) = 0 for all x ∈ (0, 1), η(1) = −ϑ(1).
321
(3.8a) (3.8b)
We are now ready to show how extrapolation works for the smooth part vε . Let us first consider the interval [0, 1 − τ ]. Lemma 3.2. Assume that ε N −1 . Then VεN (xi ) − vε (xi ) = H E(xi ) + O N −2 for all xi ∈ [0, 1 − τ ]. Proof. Fix xi ∈ (0, 1). By (3.7a), the truncation error formulas of [7, p. 21], and the inequality hj H for all j , we get N LN ε η(xi ) = Lε η(xi ) + Lε η(xi ) − Lε η(xi ) = Lε η(xi ) + O(H ).
Hence, recalling (3.8a),
2 H LN ε η(xi ) = H ζ(xi ) + O H .
Consequently Lemma 3.1 yields LN ε
VεN
− vε − H η (xi ) =
O H2
for xi ∈ (0, 1 − τ ], 2 for xi ∈ (1 − τ, 1), (h − H )ζ(xi ) + O H 2 for xi ∈ (0, 1 − τ ], O H = O(H ) for xi ∈ (1 − τ, 1).
Define the mesh function Zi = C3 N −2 (1 + xi ) + H Si
N
j =1
αhj 1+ 2ε
−1
(3.9)
for i = 0, . . . , N,
where the constant C3 will be specified later. Then N
αhj −1 N N −2 Lε Si 1+ for i = 1, . . . , N − 1. Lε Zi C3 2αN + H 2ε j =1 From Lemmas 2.3 and 2.4 and ε N −1 , for 0 < i N/2 it follows that −2 C3 αH 2 /2, LN ε Zi 2C3 αN
(3.10)
and for N/2 < i < N , LN ε Zi Now
C3 C2 H ε
−1
N
αhj −1 C3 C2 H. 1+ 2ε j =N/2+1
(3.11)
Z0 0 = VεN (0) − vε (0) − H η(0)
(3.12)
by (3.5a) and (3.8b). Also ZN C3 H V N (1) − vε (1) − H η(1)
(3.13)
ε
322
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
by (3.5a) and (3.7b), provided that C3 is sufficiently large (independently of ε and N ). From (3.9)–(3.13), one can choose the constant C3 so that Zi is a barrier function for ±[VεN (xi ) − vε (xi ) − H η(xi )], and Lemma 2.2 yields (3.14) Zi VεN (xi ) − vε (xi ) − H η(xi ) for all i. But for i = 1, . . . , N/2, Lemma 2.4 shows that N
αhj −1 −2 1+ C N −2 + H N −1 CN −2 . Zi C3 2N + H 2ε j =i+1
(3.15)
Thus for i = 1, . . . , N/2, N V (xi ) − vε (xi ) − H E(xi ) V N (xi ) − vε (xi ) − H η(xi ) + H ϑ(xi ) ε ε Zi + 2N −1 ϑ(xi ) CN −2 , where we used E = η + ϑ , (3.14), (3.15) and (3.7b).
✷
It is now a simple matter to show that the extrapolation formula (3.4) improves the accuracy of VεN in (0, 1 − τ ]. Lemma 3.3. Assume that ε N −1 . For all xi ∈ [0, 1 − τ ], 2N vε (xi ) − 2V ε (xi ) − VεN (xi ) CN −2 for some constant C. Proof. Let xi ∈ [0, 1 − τ ]. Since the subinterval meshwidths of S2N,τ are half those of SN,τ and the function E(·) depends only on τ , Lemma 3.2 implies that ε2N (xi ) − vε (xi ) = H E(xi ) + O N −2 . V 2 This relation and Lemma 3.2 yield 2N 2N ε (xi ) − VεN (xi ) = −2 V ε − vε (xi ) + VεN − vε (xi ) = O N −2 . ✷ vε (xi ) − 2V Next, consider the approximation of vε (xi ) in (1 − τ, 1). Lemma 3.4. Assume that ε N −1 . For all xi ∈ [1 − τ, 1], 2N vε (xi ) − 2V ε (xi ) − VεN (xi ) CN −2 ln N for some constant C. Proof. Define the function G(x) on [1 − τ, 1] by G(1 − τ ) = 1, G(1) = 0, and Lε G(x) = 0 for 1 − τ < x < 1. N We work initially on the mesh SN,τ . Define a discrete approximation GN ε of G by Gε (1 − τ ) = 1, N N N Gε (1) = 0, and Lε Gε (xi ) = 0 for N/2 < i < N . From a standard convergence analysis [8] for the upwind method, G(xi ) − GN (xi ) CN −1 ln N for N/2 i N. (3.16) ε
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
323
Set
Φ(xi ) = VεN (xi ) − vε (xi ) − H E(1 − τ ) GN ε (xi )
for N/2 i N.
−2
2 −2 ln N) Then Φ(xN ) = 0, Φ(xN/2 ) = O(N ) by Lemma 3.2, and LN ε Φ(xi ) = O(hi ) + O(H ) = O(N −1 −2 by Lemma 3.1 and ε N . Consequently one can use the barrier function C(1 + xi )N ln N to get Φ(xi ) CN −2 ln N for N/2 i N. (3.17)
Combining (3.16) and (3.17) and observing that |E(1 − τ )| C is implied by (3.7a) and (3.7b), we have shown that on the mesh SN,τ , for N/2 i N. (3.18) VεN (xi ) − vε (xi ) = H E(1 − τ ) G(xi ) + O N −2 ln N Similarly, on the mesh S2N,τ one gets V˜ε2N (x˜i ) − vε (x˜i ) = (H /2)E(1 − τ ) G(x˜i ) + O N −2 ln N
for N i 2N.
(3.19)
But x˜2i = xi for N/2 i N . From (3.18) and (3.19), we have 2N ε (xi ) − VεN (xi ) = 2 vε − V ε2N (xi ) − vε − VεN (xi ) = O N −2 ln N vε (xi ) − 2V for N/2 i N .
✷
3.2. Extrapolation of WεN Now we consider the error WεN (xi ) − wε (xi ). The two subintervals [0, 1 − τ ] and [1 − τ, 1] will be analysed separately. First, consider [0, 1 − τ ]. Lemma 3.5. For all xi ∈ [0, 1 − τ ], 2N wε (xi ) − 2W ε (xi ) − WεN (xi ) CN −2 for some constant C. Proof. Let xi ∈ [0, 1 − τ ]. With only slight modifications, the argument leading to [10, (3.2)] shows that ε2N )(xi )| CN −2 . Hence |(wε − WεN )(xi )| CN −2 . Similarly |(wε − W 2N wε (xi ) − 2W ε (xi ) − WεN (xi ) = 2 wε − W ε2N (xi ) − wε − WεN (xi ) CN −2 . ✷ To analyse the effect of extrapolation on (1 − τ, 1), define the function F on [1 − τ, 1] by the boundary value problem 2 Lε F (x) = − εa(x)wε (x) α F (1 − τ ) = F (1) = 0.
for x ∈ (1 − τ, 1),
(3.20a) (3.20b)
Note that F depends on τ but is otherwise independent of N . Lemma 3.6. For xi ∈ [1 − τ, 1], we have WεN (xi ) − wε (xi ) = N −1 ln N F (xi ) + O N −2 ln2 N .
(3.21)
324
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
Proof. Recalling (2.6b) and applying Lemma 2.1 to (3.20), for 0 k 3 and x ∈ (0, 1) we have (k) F (x) Cε −k e−α(1−x)/ε . (3.22) Let xi ∈ (1 − τ, 1). By a Taylor expansion of LN ε F (xi ) about x = xi , and an appeal to (2.10) and (3.22), we get −1 −α(1−x )/ε −1 i+1 N ln N . (3.23) LN ε F (xi ) = Lε F (xi ) + O ε e Recalling (2.7b) and (3.5b), more Taylor expansions show that, for some x˜1 ∈ (xi , xi+1 ) and x˜2 , x˜3 ∈ (xi−1 , xi ), N ε (4) 1 1 wε (x˜1 ) + wε(4) (x˜2 ) h2 − a(xi )wε (xi )h + a(xi )wε (x˜3 )h2 LN ε Wε − wε (xi ) = − 4! 2 3! −1 −α(1−x )/ε −2 2 2 −1 i+1 = − εN ln N a(xi )wε (xi ) + O ε e N ln N , (3.24) α by (2.6b) and (2.10). From (3.23), (3.20) and (3.24) one sees that −1 −1 ln N Lε F (xi ) + O ε −1 e−α(1−xi+1 )/ε N −2 ln2 N N ln N LN ε F (xi ) = N −1 −α(1−x )/ε −2 2 N i+1 N ln N . = LN ε Wε − wε (xi ) + O ε e That is, for all xi ∈ (1 − τ, 1), N N L W − wε − N −1 ln N F (xi ) C4 ε −1 e−α(1−xi+1 )/ε N −2 ln2 N, ε ε
(3.25)
for some constant C4 . For i = N/2, . . . , N , consider the mesh function N −2 2
αhj −1 −2 1+ , Γi := C5 N (1 + xi ) + N ln N Si ε j =1 where the constant C5 will be chosen later. Invoking (2.12), we see that N −1 −2 2
αhj −1 N . 1+ Lε Γi C2 C5 ε N ln N ε j =i+1
(3.26)
But for each i N/2, e
−α(1−xi+1 )/ε
= e
N
αhi+1 /ε
e−αhj /ε
j =i+1
e
4N −1 ln N
N
j =i+1
αhj 1+ ε
−1
N
αhj −1 e . 1+ ε j =i+1
If C5 e4 C4 /C2 , then by (3.25)–(3.27), LN Γi LN W N − wε − N −1 ln N F (xi ) ε
ε
ε
4
(3.27)
for xi ∈ (1 − τ, 1).
From (3.5b) and (3.20b), ΓN 0 = |WεN (1) − wε (1) − (N −1 ln N)F (1)|. The proof of Lemma 3.5 shows that |(wε − WεN )(1 − τ )| C6 N −2 for some constant C6 ; hence, recalling (3.20b), ΓN/2 C5 N −2 |(WεN − wε − (N −1 ln N)F )(1 − τ )| provided that C5 C6 .
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
325
Choose C5 = max{e4 C4 /C2 , C6 }. Then the above calculations establish that Γi is a barrier function for ±[WεN (xi ) − wε (xi ) − (N −1 ln N)F (xi )] when 1 − τ xi 1. Applying Lemma 2.2 on [1 − τ, 1] yields N W (xi ) − wε (xi ) − N −1 ln N F (xi ) Γi C5 N −2 ln2 N ε for i = N/2, . . . , N , as desired.
✷
Now we demonstrate how extrapolation improves the accuracy of WεN (xi ) on [1 − τ, 1]. Lemma 3.7. For all xi ∈ [1 − τ, 1], 2N wε (xi ) − 2W ε (xi ) − WεN (xi ) CN −2 ln2 N for some constant C. Proof. Let xi ∈ [1 − τ, 1]. An inspection of the proof of Lemma 3.6 shows that the factor ln N in (3.21) comes from the definition of h in (2.10). Consequently, rewriting (3.21) so as to make explicit its dependence on the choice of N and τ , one has WεN (xi ) − wε (xi ) = N −1 (τ α/2ε)F (xi ) + O N −2 (τ α/2ε)2 . ε2N is computed on the mesh S2N,τ which has the same transition point 1 − τ as SN,τ , it follows that As W ε2N (xi ) − wε (xi ) = (2N)−1 (τ α/2ε)F (xi ) + O N −2 (τ α/2ε)2 . W Combining these two formulas, we obtain 2N ε (xi ) − WεN (xi ) = 2 wε − W ε2N (xi ) − wε − WεN (xi ) = O N −2 ln2 N . wε (xi ) − 2W
✷
We come finally to our main result. Theorem 3.1 (error after extrapolation). Assume that ε N −1 . There exists a constant C such that 2N uε (xi ) − 2U ε (xi ) − UεN (xi ) CN −2 ln2 N for all xi ∈ SN,τ . Proof. For each i, 2N ε (xi ) − UεN (xi ) uε (xi ) − 2U 2N 2N ε (xi ) − VεN (xi ) + wε (xi ) − 2W ε (xi ) − WεN (xi ) . = vε (xi ) − 2V Combining the results of Lemmas 3.3, 3.4, 3.5 and 3.7, we are done.
✷
Remark 3.1. One could continue the process analysed in this paper to obtain still more accuracy by repeated extrapolation, as is done in classical Richardson extrapolation. To be precise, if a(·), b(·), f (·) ∈ C m [0, 1] for some m 2, then one can extrapolate m − 1 times. This means: choose a Shishkin mesh SN,τ as in Section 2 but with τ = mεα −1 ln N , then, setting UεN,0 = UεN , compute UεN,j (xi ) =
ε2N,j −1 (xi ) − UεN,j −1 (xi ) 2j U 2j −1 − 1
for j = 1, . . . , m − 1 and i = 0, 1, . . . , N.
326
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
ε2N,j −1 is computed on a mesh with the same transition point as UεN,j −1 but having twice as (Here U many mesh points.) The above choice of τ yields N −m in Lemma 3.5. The final result, which generalizes Theorem 3.1, is that uε (xi ) − U N,j (xi ) CN −m lnm N for all xi ∈ SN,τ . ε
Table 1 Maximum nodal errors before and after extrapolation ε
N
64
128
256
512
1024
10−1
before extrapolation rate after extrapolation rate
2.05e–02 1.21 7.91e–04 2.44
1.07e–02 1.21 2.12e–04 2.41
5.44e–03 1.19 5.51e–05 2.38
2.75e–03 1.17 1.40e–05 2.34
1.38e–03 1.16 3.54e–06 2.31
10−2
before extrapolation rate after extrapolation rate
1.55e–01 0.82 2.63e–02 1.15
9.95e–02 0.91 1.41e–02 1.97
5.99e–02 1.08 4.70e–03 2.16
3.22e–02 1.11 1.35e–03 2.21
1.67e–02 1.13 3.69e–04 2.25
10−3
before extrapolation rate after extrapolation rate
1.50e–01 0.82 2.58e–02 1.14
9.63e–02 0.74 1.39e–02 1.66
6.38e–02 0.84 5.49e–03 1.72
3.92e–02 0.92 2.04e–03 1.85
2.80e–02 0.95 6.90e–04 1.91
10−4
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.89e–04 1.91
10−5
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.89e–04 1.91
10−6
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.89e–04 1.91
10−7
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.89e–04 1.98
10−8
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.48e–03 1.72
3.91e–02 0.92 2.04e–03 1.85
2.27e–02 0.95 6.87e–04 1.83
10−9
before extrapolation rate after extrapolation rate
1.49e–01 0.82 2.57e–02 1.14
9.60e–02 0.74 1.39e–02 1.66
6.36e–02 0.84 5.49e–03 1.70
3.91e–02 0.92 2.07e–03 2.16
2.27e–02 0.96 5.81e–04 2.00
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
327
We do not present the general theory for this result as the extra notation involved would obscure the central ideas that we have described here.
4. Numerical example Consider the following singularly perturbed convection–diffusion problem in one dimension: −εuε (x) + 1 + x(1 − x) uε (x) = f (x) for x ∈ (0, 1), uε (0) = uε (1) = 0, where f (x) is chosen in such a way that the real solution uε is 1 − e−(1−x)/ε πx . − cos −1/ε 1−e 2 This function displays typical boundary-layer behaviour. Our theory assumes that ε N −1 , but this inequality is not always satisfied in the numerical experiments. We take α = 0.199. Table 1 shows the maximum nodal errors between the real and computed solutions, before and after extrapolation, for a range of values of ε and N . The rates of convergence r are computed in the standard way, assuming that convergence of the form N −r lnr N is attained; to be precise, if Eε,N is the maximum nodal error for given values of ε and N , then the associated rate in Table 1 is (ln Eε,N − ln Eε,2N )/ln(2 ln N/ln 2N ). Figs. 1 and 2 give these errors in graphical form, together with benchmark curves for certain rates of convergence to show that the computed errors decrease at approximately these rates. uε (x) =
Fig. 1. Comparison of nodal errors before and after extrapolation for ε = 10−2 .
328
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
Fig. 2. Comparison of nodal errors before and after extrapolation for ε = 10−8 .
The results show that extrapolation decreases the nodal errors significantly and increases the rate of convergence of the numerical method from O(N −1 ln N) to O(N −2 ln2 N), which supports the theoretical bounds of Theorems 2.1 and 3.1. Remark 4.1. Richardson extrapolation on Shishkin meshes is easily implemented. To compute Uε2N we must solve a linear system with 2N equations, which costs approximately four times as much as solving the original linear system for UεN . Thus the rate of convergence of the simple upwind method is raised from almost-first-order to almost-second-order at the cost of five times as much computation; this tradeoff between cost and accuracy is clearly worthwhile.
Acknowledgement We are grateful to two unknown referees who read the original version of the paper carefully and made constructive suggestions for its improvement.
References [1] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, Chapman & Hall/CRC, Boca Raton, FL, 2000. [2] A. Fröhner, H-G. Roos, The ε-uniform convergence of a defect correction method on a Shishkin mesh, Appl. Numer. Math. 37 (2001) 79–94.
M.C. Natividad, M. Stynes / Applied Numerical Mathematics 45 (2003) 315–329
329
[3] H.B. Keller, Accurate difference methods for linear ordinary differential systems subject to linear constraints, SIAM J. Numer. Anal. 6 (1969) 8–30. [4] H.B. Keller, Numerical Methods for Two-Point Boundary Value Problems, Dover, New York, 1992. [5] R.B. Kellogg, A. Tsan, Analysis of some difference approximations for a singular perturbation problem without turning points, Math. Comp. 32 (1978) 1025–1039. [6] T. Linß, Solution decompositions for linear convection–diffusion problems, Z. Anal. Anwendungen 21 (2002) 209–214. [7] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimension, World Scientific, Singapore, 1996. [8] H.-G. Roos, T. Linß, Sufficient conditions for uniform convergence on layer-adapted grids, Computing 63 (1999) 27–45. [9] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations: Convection– Diffusion and Flow Problems, Springer, Berlin, 1996. [10] M. Stynes, H.-G. Roos, The midpoint upwind scheme, Appl. Numer. Math. 23 (1997) 361–374. [11] M. Stynes, L. Tobiska, A finite difference analysis of a streamline diffusion method on a Shishkin mesh, Numer. Algorithms 18 (1998) 337–360. [12] R. Vulanovi´c, D. Herceg, N. Petrovi´c, On the extrapolation of a singularly perturbed boundary value problem, Computing 36 (1986) 69–79.