Numerical approximation of solution derivatives in the case of singularly perturbed time dependent reaction–diffusion problems

Accepted Manuscript Numerical approximation of solution derivatives in the case of singularly perturbed time dependent reaction–diffusion problems J.L. Gracia, E. O’Riordan PII: DOI: Reference:

S0377-0427(14)00268-4 http://dx.doi.org/10.1016/j.cam.2014.05.023 CAM 9692

To appear in:

Journal of Computational and Applied Mathematics

Received date: 1 March 2013 Revised date: 13 December 2013 Please cite this article as: J.L. Gracia, E. O’Riordan, Numerical approximation of solution derivatives in the case of singularly perturbed time dependent reaction–diffusion problems, Journal of Computational and Applied Mathematics (2014), http://dx.doi.org/10.1016/j.cam.2014.05.023 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Numerical approximation of solution derivatives in the case of singularly perturbed time dependent reaction–diffusion problems J.L. Gracia∗ and E. O’Riordan† May 30, 2014

Abstract Numerical approximations to the solution of a linear singularly perturbed parabolic problem are generated using a classical finite difference operator on a piecewise-uniform Shishkin mesh. First order convergence of these numerical approximations in an appropriately weighted C 1 -norm is established. Numerical results are given to illustrate the theoretical error bounds. AMS Classification: 65L11,65L12 Key Words: Singular perturbation, approximation of derivatives, reaction– diffusion problem, Shishkin mesh.

1

Introduction

Singularly perturbed differential equations are typically characterized by the presence of a small positive parameter (denoted here by ε) multiplying some or all of the highest derivatives present in the differential equation. These equations arise in numerous applications [9, 10], once the physical variables have been reformulated into dimensionless variables. Fluid dynamics, chemical kinetics, combustion and control theory are sample areas from science where singularly perturbed problems naturally appear. ∗ Corresponding author: IUMA - Department of Applied Mathematics, University of Zaragoza, Spain. [email protected] † School of Mathematical Sciences, Dublin City University, Ireland. [email protected]

1

It has been well established [2] that classical numerical approaches to solving these problems have many deficiencies, especially if one is interested in pointwise global accuracy or in estimating physically meaningful quantities such as the flux. Nodal accuracy is a necessary, but not a sufficient condition for global accuracy. In addition, the level of nodal accuracy may depend crucially on the location of the mesh points. In the context of singularly perturbed problems, high levels of accuracy can be obtained if there are no mesh points located within the thin layer regions. This effect is caused by the fact that within the layer regions the magnitude of the first derivative of the solution can be significantly larger compared to the size of the solution derivatives outside the layers. Moreover, the size of any layer region shrinks to zero as the singular perturbation parameter tends to zero. Hence, unless the mesh is layer-adapted, for sufficiently small values of the perturbation parameter there will be no mesh points within the layer. For this reason, we are interested in measuring the pointwise accuracy of numerical approximations at all points within the domain of the continuous problem, both within and outside the layers. Hence, our interest is not only in parameter-uniform nodal convergence, but in parameter-uniform global convergence [2]. Since the seminal papers of Bakhvalov [1] and Il’in [5], parameter-uniform numerical methods have been developed for several classes of singularly perturbed problems. These parameter-uniform methods generate numerical approximations, whose pointwise accuracy is retained independently of the singular perturbation parameter. In [12] Shishkin established that a uniform mesh is not an adequate discretization of the domain, if one wants to design a parameter-uniform numerical method for a class of singularly perturbed time dependent reaction-diffusion problems. For the approximations to have a pointwise accuracy independent of the perturbation parameter, a non-uniform (layer-adapted) mesh is a necessary component for any parameter-uniform numerical method, when dealing with a class of problems of the form − εuxx + bu + ut = f, u = g,

(x, t) ∈ G := (0, 1) × (0, 1], ¯ \ G. (x, t) ∈ G

(1a) (1b)

Moreover, in [12] Shishkin created his now well-known piecewise-uniform mesh to generate such parameter-uniform approximations to the solutions of problems from this class of time dependent reaction-diffusion problems. Shishkin [14] and others have subsequently continued to show that piecewiseuniform meshes can be used to create parameter-uniform numerical methods 2

for an extensive class of singularly perturbed problems. In this paper, we return to the same class of singularly perturbed reaction-diffusion problems as in [12], to examine other attributes of the Shishkin mesh. For the problem class (1), it was established in [7] that for a numerical solution U N generated using a standard finite difference operator and an appropriate piecewise-uniform Shishkin mesh, one has a global error bound of the form ¯ ∥G ≤ C(N −1 ln N )2 + CM −1 , ∥u − U (2) ¯ denotes a piecewise bilinear interpolant over the domain G ¯ of the where U N,M discrete solution U , N, M denotes the number of mesh points used in the space and time directions, respectively; and we define the global pointwise norm ∥u∥G := ess sup(x,t)∈G |u(x, t)|

Throughout this paper, C denotes a generic constant that is independent of the singular perturbation parameter ε and of all discretization parameters. An additional feature of Shishkin meshes is that accurate approximations to the scaled flux can be easily generated. Numerical evidence for this was given in [2, 7]. In this paper a proof of the convergence of the scaled discrete derivatives is given. As the derivatives of the solution are only large in the boundary layers, it is only necessary to scale the derivatives in the boundary layer regions. For this reason, in this paper, we will derive an error estimate in the following weighted C 1 -norm ∥v∥1,κ,G := ∥κ(x)vx ∥G + ∥vt ∥G + ∥v∥G ,

(3)

with κ(x) :=

{ √ ε, if

x ∈ (0,

1, otherwise.

√

ε β

ln 1ε ) ∪ (1 −

√

ε β

ln 1ε , 1),

Shishkin et al. [13], [14, pp.353-354] have highlighted the difficulties in selecting a weighted-C 1 norm to appropriately measure the relative error in numerical approximations to the solutions of singularly perturbed problems. In particular, observe that the weight κ(x) in the √ above norm is excessively √ small in the region where C ε << α|x − p| ≤ βε ln(1/ε), p = 0, 1. Nevertheless, outside this subregion, the norm ∥ · ∥1,κ,G is an appropriate norm to use in the context of problem (1). This paper is structured as follows: In §2 the continuous problem is stated and parameter-explicit bounds on the derivatives of the solution are 3

given. In §3, the numerical method is described and in §4 scaled nodal approximations of the space and time derivatives, respectively, are deduced. The main result of the paper, which establishes the error estimate ¯ − u∥1,κ,G ≤ CN −1 ln N, ∥U is given in §5. Some numerical results are presented in the final section of the paper.

2

Continuous problem

Consider the following class of singularly perturbed parabolic problems Lε u := −εuxx + b(x, t)u + c(x, t)ut = f (x, t), in G := Ω × (0, T ], (4a) u = 0,

on ΓB ∪ ΓL ∪ ΓR ,

b(x, t) ≥ ∥ct ∥G + β > 0,

0 < ε ≤ 1;

c(x, t) ≥ c0 > 0,

(4b)

(4c)

where 0 < ε ≤ 1, Ω := (0, 1), ΓB := {(x, 0) | 0 ≤ x ≤ 1}, ΓL := {(0, t) | 0 ≤ ¯ t ≤ T }, ΓR := {(1, t) | 0 ≤ t ≤ T } and Γ := G\G. Since the problem is linear, there is no loss in generality in assuming zero boundary/initial conditions. The constraint (4c) on the coefficient b(x, t) can be absorbed into the time variable by using the change of variable u = veγt , where γ > 0 is sufficiently large. We assume that the data of the problem satisfy sufficient regularity and compatibility conditions so that the solution of problem (4) is such that u ∈ C 6+γ (G)1 and for the analysis presented below to be applicable. It is well–known that the differential operator associated with (4) satisfies a comparison principle. From this, one can establish the stability estimate ∥u∥ ≤

t ∥f ∥G . c0

√ Layers of width O( ε ln(1/ε)) occur on both sides of the domain and we have the following bounds on the solution and its components. 1

The space C n+γ (D) is the set of all functions, whose derivatives of order n are H¨ older continuous of degree γ > 0. That is, C n+γ (D) := {z :

∂ i+j z ∈ C γ (D), 0 ≤ i + 2j ≤ n}. ∂xi ∂tj

4

Theorem 1. The solution of (4) can be written in the form u = v+wL +wR , where the regular component v ∈ C 6+γ (G) satisfies Lε v = f, in G,

v = u, on ΓB ,

and v = v ∗ can be specified on the boundary ΓR ∪ ΓL , so that

∂ k+m v

k m ≤ C(1 + ε2−k/2 ), 0 ≤ k + 2m ≤ 6. ∂x ∂t G

(5a)

(5b)

The singular components wL , wR satisfy the problems

Lε wL = 0, in G, wL = 0, on ΓB ∪ ΓR , wL = u − v on ΓL ,

Lε wR = 0, in G, wR = 0, on ΓB ∪ ΓL , wR = u − v on ΓR ,

(6a) (6b)

and for all points (x, t) ∈ G and 0 ≤ k + 2m ≤ 6

√ ∂ k+m w β R −k/2 − ε (1−x) , e k m (x, t) ≤ Cε ∂x ∂t √ ∂ k+m w L − βε x . k m (x, t) ≤ Cε−k/2 e ∂x ∂t

(6c) (6d)

Proof. Apply the arguments given in [7] and extend to the case where u ∈ C 6+γ (G). Sufficient compatibility conditions so that u ∈ C 6+γ (G) are given in [11, Appendix]. The regular component v is written in the form v = v0∗ + εv1∗ + ε2 v2∗ , where v0∗ is the solution of the reduced problem L∗0 v0∗ := b∗ (x, t)v0∗ + c∗ (x, t)(v0∗ )t = f ∗ (x, t), (x, t) ∈ G∗ ,

v ∗ (x, 0) = 0, x ∈ (−d, 1 + d);

defined on an extended domain G∗ := (−d, 1 + d) × (0, T ], d > 0. The correcting terms v1∗ , v2∗ are defined as the solutions of L∗0 v1∗ = (v0∗ )xx , (x, t) ∈ G∗ , v1∗ (x, 0) = 0, x ∈ (−d, 1 + d); ¯ ∗ \ G∗ . L∗ε v2∗ = (v1∗ )xx , (x, t) ∈ G∗ , v2∗ (x, t) = 0, (x, t) ∈ G The introduction of the extended domain avoids the imposition of artificial compatibility conditions in this decomposition of the solution into regular and singular components (see, for example [14], or [8, Chapter 12] for a discussion of this extension of the domain in the context of singularly perturbed problems). To deduce the pointwise bounds on the derivatives of the components v, wL , wR , classical a priori bounds from [6] are combined with the stretched variables x/ε, (1 − x)/ε, as described in [7], to complete the proof. 5

3

Numerical scheme

Consider a uniform mesh in time ω ¯ M = {kτ, 0 ≤ k ≤ M, τ = T /M } ¯ N [2] in space on which numerical and a piecewise–uniform Shishkin mesh Ω approximations of the solution of problem (4) are generated. The Shishkin mesh for the reaction-diffusion problem is fitted to the two boundary layers by splitting the space domain as follows: [0, σ] ∪ [σ, 1 − σ] ∪ [1 − σ, 1].

(7a)

The mesh points are distributed in the ratio N/4 : N/2 : N/4 across these three subintervals. The transition point σ is taken to be √ ε σ := min{0.25, 2 ln N }. (7b) β We confine our attention to the case where σ < 0.25. The classical case, when σ = 0.25, is dealt with using the same arguments as here and the fact that ε−1/2 ≤ C ln N in this case. ¯ N = {xi } is given by The grid in the space variable Ω  if 0 ≤ i ≤ N/4,  ih, σ + (i − N/4)H, if N/4 ≤ i ≤ 3N/4, xi = (8)  (1 − σ) + (i − 3N/4)h, if 3N/4 ≤ i ≤ N,

where the step sizes are h := 4σ/N and H := 2(1 − 2σ)/N . We denote the local step sizes by hj := xj − xj−1 for j = 1, . . . , N , and we define the following sets of mesh points ¯ N,M := Ω ¯N × ω ¯ N,M ∩ G, G ¯ M , GN,M := G

¯ N,M \GN,M . ΓN,M := G

We use a classical finite difference operator on this mesh to produce the following numerical method: { N,M L Ui,j := (−εδx2 + b(xi , tj )I + c(xi , tj )Dt− )Ui,j = f (xi , tj ), (xi , tj ) ∈ GN,M , U0,j = 0, UN,j = 0, j = 0, . . . , M ; Ui,0 = 0, i = 0, . . . , N, (9) where Ui,j := U (xi , tj ) and Dt− Ui,j :=

Ui,j − Ui,j−1 Ui+1,j − Ui,j Ui,j − Ui−1,j , Dx+ Ui,j := , Dx− Ui,j := , τ hi+1 hi 1 hi + hi+1 δx2 Ui,j := (Dx+ Ui,j − Dx− Ui,j ), ~i := . ~i 2 6

We also introduce the following notation Ωi := (xi−1 , xi ),

Ri,j := Ωi × (tj−1 , tj ).

Throughout the analysis in this paper we assume that (10)

C1 N ≤ M ≤ C2 N.

4

Nodal approximations of derivatives

We denote the nodal error by e(xi , tj ) := U (xi , tj ) − u(xi , tj ), and the associated truncation error by T (xi , tj ) := LN,M e(xi , tj ). We have the following truncation error bounds |(Dt− u − ut )(xi , tj )| ≤ CN −1 ∥utt (xi , t)∥t∈(ti−1 ,ti ) , (11a)

|(δx2 u − uxx )(xi , tj )| ≤ min{~i |uxxx (x, tj )|Ωi ∪Ωi+1 , |uxx (x, tj )|Ωi ∪Ωi+1 },(11b) and if hi = hi+1 , then the following bound holds |(δx2 u−uxx )(xi , tj )| ≤ min{h2i |uxxx (x, tj )|Ωi ∪Ωi+1 , |uxx (x, tj )|Ωi ∪Ωi+1 }, (11c) which will be used throughout the paper. Using the bounds on the components of the solution u given in Theorem 1 we easily deduce from (11a) that |(Dt− u − ut )(xi , tj )| ≤ CN −1 , (xi , tj ) ∈ [0, 1] × [τ, T ]. (12) We follow the approach outlined in [3, Appendix A.2] and [4, §8.2] to deduce nodal approximations of the time derivatives. Lemma 1. Assume (10). Then, the following bound holds, |(Dt− U − ut )(xi , tj )| ≤ CN −1 ln N,

xi ∈ [0, 1], tj > 0,

where U is the solution generated by the numerical method (9) and u is the solution of the continuous problem (4). Proof. Given the bound (12), we only need to bound the quantity |Dt− (U − u)|. The crude truncation error bound ∥LN,M (U − u)∥GN,M ≤ CN −1 , can be deduced from the bounds on the components given in Theorem 1 and the truncation error bounds (11). Using a barrier function we obtain the nodal error bound |(U − u)(xi , tj )| ≤ Ctj N −1 , 7

xi ∈ [0, 1].

From this, one deduces that at the first time level |Dt− (U − u)(xi , τ )| ≤ CN −1 ,

xi ∈ [0, 1].

Along the side boundaries we note that Dt− (U − u)(xi , tj ) = 0,

¯ N,M , (xi , tj ) ∈ (ΓL ∪ ΓR ) ∩ G

tj ≥ τ.

At the interior points (xi , tj ) ∈ GN,M ∩ {tj ≥ 2τ }, we will estimate |LN,M (Dt− (U − u)(xi , tj )|, by first reversing the order of the operators LN,M and Dt− . To this end, we use the identity Dt− (Pi,j Qi,j ) ≡ Pi,j Dt− Qi,j + Qi,j−1 Dt− Pi,j ,

(13)

and we define a minor modification to the operator LN,M denoted by ˜ N,M Zi,j := (LN,M + (c(xi , tj−1 ) − c(xi , tj ))D− + (D− c(xi , tj ))I)Zi,j , L t t to obtain ˜ N,M D− (U − u)| = |D− (LN,M (U − u)) − D− (b(U − u)) + bD− (U − u) |L t t t t ≤ |Dt− (LN,M (U − u))| + CN −1

≤ ε|Dt− (δx2 u − uxx )| + C|Dt− (Dt− u − ut )| + CN −1 . Using the solution decomposition in Theorem 1 and the truncation error bounds (11) we deduce that ˜ N,M D− (U − u)| ≤ CN −1 . |L t Use the discrete comparison principle to complete the proof. To bound the quantity Dx− (U − u) at the mesh points, we bound a ˆ N,M (Dx− (U −u))∥, where the finite difference operator quantity of the form ∥L N,M ˆ L is monotone and is defined below in (14). Then we use a discrete comparison principle to bound Dx− (U − u). We define the discrete error flux to be − Ui,j := Dx− e(xi , tj ), if 0 < xi ≤ 1.

We identify a discrete problem associated with the error flux defined over the region ¯ N,M := G ¯ N,M ∩ {[H, 1] × [0, T ]}. G H

GN,M := GN,M ∩ {(H, 1) × (0, T ]}; H 8

In addition, we introduce a new finite difference operator δˆx2 defined by ) ~i 1 ( hi+1 + δˆx2 Zi,j := Dx − Dx− Zi,j , ~i hi ~i−1

which has the property that

δˆx2 Dx− Zi,j ≡ Dx− δx2 Zi,j .

Using the identities (13) and Dx− (LN,M e(xi , tj )) = Dx− T (xi , tj ), we see that ¯ N,M , the quantity U − satisfies for all mesh points within the region G i,j H ˆ N,M U − = Dx− T (xi , tj ) − e(xi−1 , tj )Dx− b(xi , tj ), (xi , tj ) ∈ GN,M ,(14) L i,j H where for the internal points (xi , tj ) ∈ GN,M , H

ˆ N,M Zi,j := −εδˆx2 Zi,j + b(xi , tj )Zi,j + c(xi , tj )D− Zi,j , L t

ˆ N,M Zi,j := Zi,j for (xi , tj ) ∈ G ¯ N,M \ GN,M . and L H H

Remark 1. The discrete operator δˆx2 on the piecewise uniform mesh is defined as follows  2 δ , if xi ̸= 1 − σ, 1 − σ + h, σ, σ + H,   x      1( )   2(1 − 2σ)Dx+ − Dx− , if xi = σ,    h       1( + ) Dx − 2(1 − 2σ)Dx− , if xi = σ + H, δˆx2 ≡ H      )  1(   if xi = 1 − σ, 4σDx+ − Dx− ,   H         1 (D+ − 4σD− ), if xi = 1 − σ + h. x h x

When bounding the term Dx− T (xi , tj ) we will make use of the following truncation error bounds ∫ ∫ tj 1 tj |Dx− (ut − Dt− u)(xi , tj )| = Dx− utt (xi , s)dsdr τ t=tj−1 s=r If hi−1 = hi = hi+1 , then

≤ Cτ ∥uttx ∥Ri,j .

|Dx− (uxx − δx2 u)(xi , tj )| ≤ Ch2i ∥uxxxxx (x, tj )∥x∈(xi−2 ,xi+1 ) . 9

(15a)

(15b)

In the next result, we bound the scaled discrete error flux, using a weight√ ing factor of ε across the entire domain. Theorem 2. Assume (10). Then, the following bound holds √ ε|Dx− (U − u)(xi , tj )| ≤ CN −1 (ln N ), 0 < xi ≤ 1, tj ≥ 0;

(16)

where U is the solution generated by the numerical method (9) and u is the solution of the continuous problem (4). Proof. Using the fact that utt (0, t) = 0, we have for some η ∈ (0, xi ), ∫ ∫ tj xi tj − |(ut − Dt u)(xi , tj )| = uttx (η, s)dsdt τ t=tj−1 s=t ∫ tj ∫ x i tj xi τ ∥uttx ∥(0,xi )×(tj−1 ,tj ) ds dt ≤ C √ . ≤ τ tj−1 s=t ε Hence, one can deduce the nodal error bound, xi 1 − xi |(U − u)(xi , tj )| ≤ C((N −1 ln N )2 + τ min{ √ , √ }), ε ε which allows one prove that the boundary error fluxes are bounded by √ max ε|Dx− (U − u)(p, tj )| ≤ CN −1 ln N. (17) p=x1 ,xN

Using the truncation error bounds (15), and the following estimate of the second order central difference 2 |δxx w(xi , tj )| ≤

max

x∈(xi−1 ,xi+1 )

10

|wxx (x, tj )|,

(18)

we can deduce that ˆ N,M U − | ≤ CN −1 (1 + ∥uttx ∥R ) |L i,j i,j  √ −1 ln N, εN ) if xi = x1 , 1, (   √   β √ √  − ε xi  ), if xi ∈ (x1 , σ), N√−1 ln N ( ε + ( ε)−1 e    −1  ε N √ if xi = σ, 1 − σ + h, ln N + ε ln N , +C −1  ε+N , if xi = 1 − σ, σ + H,    −1 ,  N if xi ∈ (σ + H, 1 − σ),   √  β  √ √  −1 − ε (1−xi ) N ln N ( ε + ( ε)−1 e ), if xi ∈ (1 − σ + h, 1).  √ −1 ( εN ) ln N, if xi = x1 , 1,   √   β √ √  − ε xi  N√−1 ln N ( ε + ( ε)−1 e ), if xi ∈ (x1 , σ),    −1  ε N √ if xi = σ, 1 − σ + h, ln N + ε ln N , ≤ CN −1 + C −1  ε+N , if xi = 1 − σ, σ + H,    −1 ,  N if xi ∈ (σ + H, 1 − σ),   √   √ −1 − βε (1−xi ) √  −1 ), if xi ∈ (1 − σ + h, 1). N ln N ( ε + ( ε) e

Except for the four mesh points σ, σ + H, 1 − σ, 1 − σ + h, this bound is √ of the form O((N −1 ln N )/ ε) within the layers and it is O(N −1 ) outside the layers. We construct a suitable barrier function to reflect this difference between the fine and coarse mesh regions. Consider the following four mesh functions: { xi { xi , if xi ≤ σ σ+H , if xi ≤ σ + H σ ; R2 (xi ) := ; R1 (xi ) := 1, if σ < xi ≤ 1 1, if σ + H < xi ≤ 1 { { 1, if xi ≤ 1 − σ + h, 1, if xi ≤ 1 − σ R3 (xi ) := ; R (x ) := 4 i 1−xi 1−xi , if 1 − σ + h < xi ≤ 1, , if 1 − σ < x ≤ 1 i σ−h σ

11

which satisfy the inequalities { βN ε hσ = (4 ln N )2 , if xi = σ ; ˆ N,M R1 (xi ) ≥ L 0, if xi ̸= σ  2ε(1−2σ)   H(H+σ) ≥ εN 2 , if xi = σ + H N,M ˆ L R2 (xi ) ≥ 0, if xi < σ + H ;   β, if xi > σ + H  4ε  H ≥ 2εN, if xi = 1 − σ ˆ N,M R3 (xi ) ≥ 0, if xi > 1 − σ ; L  β, if xi < 1 − σ { βN ε h(σ−h) ≥ (4 ln N )2 , if xi = 1 − σ + h ˆ N,M R4 (xi ) ≥ L 0, if xi ̸= 1 − σ + h.

From the previous lower bounds, we form the barrier function √ √ CN −1 ln N ε(R1 (xi ) + R4 (xi )) + CN −1 (R2 (xi ) + R3 (xi )) + C( εN )−1 ln N, to complete the proof. We can sharpen this error bound, by removing the scaling factor of from outside the computational layer regions (σ, 1 − σ]. Theorem 3. Assume (10). Then, for all tj ≥ 0, √ ε|(Dx− (U − u))(xi , tj )| ≤ CN −1 , if xi ∈ (0, σ] ∪ (1 − σ, 1], − −1 |(Dx (U − u))(xi , tj )| ≤ CN (ln N ), if σ < xi ≤ 1 − σ,

√ ε

(19)

where U is the solution generated by the numerical method (9) and u is the solution of the continuous problem (4). Proof. Define the local mesh Peclet numbers as follows: √ √ β β ξ1 := H, ρ1 := h. ε ε Let us consider the following two mesh functions  (1 + 0.5ρ1 )−i , if xi ≤ σ,    (1 + 0.5ρ1 )−N/4 (1 + ξ1 )N/4−i , if σ ≤ xi ≤ 0.5, P1 (xi ) =  (1 + 0.5ρ1 )−N/4 (1 + ξ1 )i−3N/4 , if 0.5 ≤ xi < 1 − σ,   (1 + 0.5ρ1 )i−N , if xi ≥ 1 − σ.   1, if xi ≤ σ, 0, if σ < xi ≤ 1 − σ, Q1 (xi ) =  1, if xi < 1 − σ. 12

We give some properties for the functions P1 and Q1 . By symmetry, we only need to consider the points where xi < xN/2 . At the transition point σ, we have that N −1 ≤ P1 (σ) ≤ 4N −1 , and for the remaining mesh points that  √ √ ε(1 + 0.5ρ1 )Dx+ P1 = −0.5√βP1  √ − εDx P1 = −0.5 βP1 if xi < σ,  ε(1 + 0.5ρ1 )δx2 P1 = 0.25βP1  √ √ ε(1 + ξ1 )Dx+ P1 = −√βP1  √ − εDx P1 = − βP1 if σ < xi < xN/2 .  ε(1 + ξ1 )δx2 P1 = βP1

These expressions allow us establish the following lower bounds  1, if xi = x1 , 1,  √   β  − x   P1 (xi ) ≥ e √ ε i if xi ∈ (x1 , σ),    β   P (x ) ≥ e− ε (1−xi ) if x ∈ (1 − σ + h, 1), 1 i i ˆ N,M P1 (xi ) ≥ C L 1 , if x = σ, 1 − σ + h, −  i ln N  √ √   −N εP1 ≥ − ε, if xi = 1 − σ, σ + H,    −1 ,  N if xi ∈ (σ + H, 1 − σ) and xi ̸= xN/2 ,    √ − ε(1 + ξ)−N/4 , if xi = xN/2 .  if xi ̸= σ, σ + H, 1 − σ, 1 − σ + h,  0, Nε ˆ N,M Q1 (xi ) ≥ C , if xi = σ, 1 − σ + h, L  hN ε − H , if xi = σ + H, 1 − σ. Hence, we can construct a majoring function for the mesh function U − − |Ui,j |

≤ CN

Noting that

−1

ln N

4 ∑

Rj (xi ) + C

j=1

) CN −2 N −1 ln N ( √ P1 (xi ) + √ Q1 (xi ) . ε ε

N −1 N −1 N −1 √ |P1 (σ + H)| ≤ √ ≤ CN −1 , ε ε 1+ξ

we deduce the following estimates for the discrete error flux { √ − ε|Ui,j | ≤ CN −1 ln N, if xi ∈ (0, σ] ∪ [1 − σ, 1], − |Ui,j | ≤ CN −1 ln N, otherwise, which completes the proof. 13

√ In the next result, we show that it is not necessary to scale by ε in the √ √ 1 ε 1 outer region ( βε ln , 1 − ln ), even if the mesh points lie within the ε β ε fine mesh. Theorem 4. Assume (10). Then, for all tj ≥ 0, √ √ ε 1 ε 1 − −1 |Dx (U − u)(xi , tj )| ≤ CN ln N, xi ∈ ( ln , 1 − ln ), β ε β ε

(20)

where U is the solution generated by the numerical method (9) and u is the solution of the continuous problem (4). Proof. From the bounds given in Theorem 3, we only need to consider the case where N −2 ≤ ε, as then ln 1ε ≤ 2 ln N . Moreover, we only need to consider the mesh points in the region √ √ √ √ ε 1 ε ε ε 1 ln , 2 ln N ) ∪ (1 − 2 ln N, 1 − ln ). ( β ε β β β ε √ √ √ ε 1 Let us examine the subregion ( βε ln 1ε , 2 βε ln N ). Since xi > β ln ε , √ there exists some θ > 1 such that xi ≥ θ βε ln 1ε . For N sufficiently large, we note that

(1 + ρ)−1 ≤ e−ρ/θ , ρ ≤ θ ln θ,

for θ > 1.

Hence, looking at the barrier function P1 , it follows in this case that √ P1 (xi ) ≤ (1 + 0.5ρ1 )−i ≤ e−0.5iρ1 /θ ≤ ε. Using this inequality we establish that, in this particular case, the following estimate for the discrete error flux |Dx− (U − u)(xi , tj )| ≤ CN −1 ln N.

5

Global accuracy in weighted C 1 -norm

In this section, we examine the global accuracy, in the weighted C 1 -norm ∥ · ∥1,χ,G defined in (3), of the bilinear interpolant ¯ (x, t) := U

N∑ −1,M

U (xi , tj )ϕi (x)ψj (t),

i,j=1

14

¯ (x, t) ∈ G,

where ϕi (x), ψj (t) are piecewise linear basis functions in space and time, defined by the nodal values of ϕi (xk ) = δi,k = ψi (tk ). Note the following bound on the bilinear interpolant g¯ of a function g (see e.g. [15, Lemma 4.1]) ∫ xi 2 ∥g − g¯∥Rij ≤ C min{hi ∥gxx ∥Rij , max |gx (s, t)|ds} t∈[tj−1 ,tj ] xi−1 ∫ tj

+C min{τ 2 ∥gtt ∥Rij ,

max

x∈[xi−1 ,xi ] tj−1

|gt (x, s)|ds}. (21)

Theorem 5. Assume (10). Then, ¯ − u∥1,κ,G ≤ CN −1 ln N, ∥U

(22)

¯ is the bilinear interpolant of the nuwhere u is the solution of (4) and U merical solution on the piecewise-uniform mesh defined by (7). Proof. Using the decomposition u = v +wL +wR and splitting the argument to inside and outside the computational layer regions, we have ∥u − u ¯∥G ≤ C(N −1 ln N )2 . Hence, the following global error estimate follows: ¯ ∥G ≤ CN −1 ln N. ∥u − U Note that ¯ −u (U ¯)t (x, t) = ¯ −u (U ¯)x (x, t) =

N −1 ∑ i=1 M ∑ j=1

Dt− (U − u)(xi , tj )ϕi (x),

t ∈ (tj−1 , tj ],

Dx− (U − u)(xi , tj )ψj (t),

x ∈ (xi−1 , xi ].

Using Lemma 1 for the time derivatives and Theorems 3 and 4 for the space derivatives, we have that ¯ −u ∥U ¯∥1,κ,G ≤ CN −1 ln N. The interpolation error ∥u − u ¯∥1,κ,G is next to be bounded. First, we note that for t ∈ (tj − 1, tj ], (¯ u−u)t (x, t) =

N −1 ∑ i=1

(

N −1 ∑ ) Dt− u(xi , tj )−ut (xi , t) ϕi (x)+ ut (xi , t)ϕi (x)−ut (x, t). i=1

15

Consequently, in the rectangular cell Rij , we have that ∥(u − u ¯)t ∥Rij ≤ Cτ ∥utt ∥Rij + min{hi ∥uxt ∥Rij , ∥ut ∥Rij }. By using this estimate, the decomposition u = v + wL + wR and splitting the argument to inside and outside the layer region, we obtain that ∥(u − u ¯)t ∥G ≤ CN −1 . Secondly, for x ∈ (xi−1 , xi ], we have (¯ u−u)x (x, t) =

M ∑ (

)

Dx− u(xi , tj )−ux (x, tj )

j=1

ψj (t)+

M ∑

ux (x, tj )ψj (t)−ux (x, t).

j=1

Therefore, in the rectangular cell Rij , we obtain ∥(u − u ¯)x ∥Rij ≤ min{hi ∥uxx ∥Rij , ∥ux ∥Rij } + min{τ ∥uxt ∥Rij , ∥ux ∥Rij }. (23) To deduce bounds on the term ∥(u−¯ u)x ∥Rij , we employ the decomposition of the solution u = v + wL + wR . The regular component satisfies the following estimate ∥(v − v¯)x ∥Rij ≤ CN −1 . For the layer component wL√, we again split the argument to inside and outside the layer region (0, βε ln 1ε ] × (0, T ] and deal with the two cases

of ε ≤ N −2 and ε ≥ N −2 . We observe the following: If ε ≤ N −2 then 2 ln N ≤ ln 1ε and in this case we obtain √ √ ε 1 −1 ∥(wL )x ∥Rij ≤ C ε ≤ CN , xi ≥ ln , β ε √ ε 1 ∥(wL )x ∥Rij ≤ Cε−1/2 N −2 , σ ≤ xi < ln , β ε hi ∥(wL )xx ∥Rij + τ ∥(wL )xt ∥Rij

≤ Cε−1/2 N −1 ln N,

In the second case, where ε ≥ N −2 , then hi ∥(wL )xx ∥Rij + τ ∥(wL )xt ∥Rij

√ ≤ C(hi + τ ε),

hi ∥(wL )xx ∥Rij + τ ∥(wL )xt ∥Rij

≤ Cε−1/2 N −1 ln N,

xi < σ. √

ε 1 ln , β ε √ ε 1 xi ≤ ln . β ε

xi >

Similar estimates are obtained for the right boundary layer component wR . Combining all these bounds completes the proof. 16

6

Numerical experiments

Consider the following test problem: − εuxx + (1 + x + t)u + ut = 43 x3 t2 (1 − x)3 ,

u(0, t) = u(1, t) = t ;

u(x, 0) = (4x(1 − x)) .

3

3

(24a) (24b)

The computed solution is displayed in Figure 1 for ε = 2−20 and N = M = 64, where boundary layers are visible near x = 0 and x = 1. To estimate

1 0.8 0.6 0.4 0.2 0 0 0.5 1

0

0.2

0.4

0.6

0.8

1

Time Variable

Space variable

Figure 1: Computed solution generated by the numerical scheme (9) for ε = 2−20 and N = M = 64 the errors in the weighted C 1 -norm ∥ · ∥1,κ,G of the numerical scheme (9) for any fixed value of the singular perturbation parameter ε, we use a variant of the double mesh principle (see [2]): Given a numerical approximation U N,M ¯ N,M , we also generate the numerical solution on a generated over a mesh G 2048,2048 fine Shishkin mesh U with N, M < 2048, and the global two–mesh differences are computed: ¯ N,M − U ¯ 2048,2048 ∥ ¯ N,M ¯ 2048,2048 , EεN,M := ∥U 1,κ,G ∪G

¯ denotes the bilinear interpolant of the numerical solution. From where U these values, we compute the approximate order of convergence using QN,M := log2 (EεN,M /Eε2N,2M ). ε The uniform global errors in the norm ∥ · ∥1,χ,G and their orders of convergence are estimated as follows E N,M := max EεN,M , ε∈S

QN,M := log2 (E N,M /E 2N,2M ), 17

with S = {20 , 2−1 , 2−2 , . . . , 2−30 }. In Table 1, we present the computed global EεN,M and uniform global E N,M errors, measured in the global norm ∥ · ∥1,κ,G , for N = M = 2j , j = 4, 5, 6, 7, 8, 9 with their corresponding orders of convergence. The numerical results are in agreement with the bounds given in Theorem 5.

Acknowledgments This research was partially supported by the project MEC/FEDER MTM 2010-16917 and the Diputaci´on General de Arag´on.

References [1] N.S. Bakhvalov, On the optimization of methods for boundary-value problems with boundary layers, J. Numer. Meth. Math. Phys., 9 (1969), no. 4, 841–859 (in Russian). [2] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan and G.I. Shishkin, Robust computational techniques for boundary layers, Chapman & Hall, 2000. [3] P.W. Hemker, G.I. Shishkin and L.P. Shishkina, ε-uniform schemes with high–order time–accuracy for parabolic singular perturbation problems, IMA J. Num. Anal., 20 (2000), no. 1, 99–121. [4] P.W. Hemker, G.I. Shishkin and L.P. Shishkina, High-order timeaccuracy schemes for parabolic singular perturbation problems with convection, Russian J. Numer. Anal. Math. Modelling, 17 (2002), no. 1, 1-24. [5] A.M. Il’in, Differencing scheme for a differential equation with a small parameter affecting the highest derivative. Math. Notes, 6 (1969) no. 2, 596–602. [6] O.A. Ladyzhenskaya, V.A. Solonnikov and N.N. Ural’tseva, Linear and quasilinear equations of parabolic type, Transactions of Mathematical Monographs, 23, American Mathematical Society, 1968. [7] J.J.H. Miller, E. O’Riordan, G.I. Shishkin and L.P. Shishkina, Fitted mesh methods for problems with parabolic boundary layers. Math. Proc. Royal Irish Acad., 98A (1998), no. 2, 173–190. 18

[8] J.J.H. Miller, E. O’Riordan and G.I. Shishkin, Fitted numerical methods for singular perturbation problems, Revised Edition, World-Scientific, Singapore, 2012. [9] K.W. Morton, Numerical solution of convection diffusion problems, Chapman and Hall, London, 1995. [10] R.E. O’Malley Jr, Singular perturbation methods for ordinary differential equations. Springer–Verlag, New York, 1991. [11] E. O’Riordan and G. I. Shishkin, Solution decompositions for singularly perturbed elliptic problems, Dublin City University School of Mathematical Sciences Preprint MS-06-01, (2006). [12] G.I. Shishkin, Approximation of solutions of singularly perturbed boundary value problems with a parabolic boundary layer, USSR Comput. Maths. Math. Phys., 29 (1989), no. 4, 1–10. [13] G.I. Shishkin, Discrete approximations of solutions and derivatives for a singularly perturbed parabolic convection-diffusion equation, J. Comput. Appl. Math., 166, (2004), no. 1, 247–266. [14] G.I. Shishkin and L.P. Shishkina, Discrete approximation of singular perturbation problems, Chapman and Hall/CRC Press, 2008. [15] M. Stynes and E. O’Riordan, A uniformly convergent Galerkin method on a Shishkin mesh for a convection-diffusion problem. J. Math. Anal. Appl. 214, (1997), no. 1, 36-54.

19

Table 1: Computed global errors in ∥·∥1,κ,G estimated by EεN,M and uniform global errors E N,M with their corresponding computed orders of convergence QN,M , QN,M for the test problem (24). ε N=M=16 1.774E+001 0.183 8.758E+000 0.322 3.919E+000 0.469 1.759E+000 0.648 8.937E-001 0.939 1.117E+000 0.904 1.434E+000 0.847 1.932E+000 0.841 2.528E+000 0.797 3.108E+000 0.664 3.093E+000 0.429 3.073E+000 0.432 3.051E+000 0.436 3.051E+000 0.436 3.051E+000 0.436 3.051E+000 0.436 3.051E+000 0.436 3.051E+000 0.436 3.052E+000 0.436 3.052E+000 0.436 3.052E+000 0.436

N=M=32 1.562E+001 0.306 7.007E+000 0.459 2.832E+000 0.599 1.123E+000 0.729 4.663E-001 0.874 5.970E-001 0.970 7.971E-001 0.963 1.078E+000 0.949 1.455E+000 0.919 1.962E+000 0.885 2.297E+000 0.662 2.278E+000 0.617 2.255E+000 0.623 2.255E+000 0.623 2.255E+000 0.623 2.255E+000 0.623 2.256E+000 0.623 2.256E+000 0.623 2.256E+000 0.623 2.256E+000 0.623 2.256E+000 0.623

N=M=64 1.263E+001 0.468 5.098E+000 0.614 1.870E+000 0.730 6.773E-001 0.840 2.544E-001 0.938 3.048E-001 1.017 4.090E-001 1.024 5.583E-001 1.025 7.695E-001 1.021 1.062E+000 1.007 1.452E+000 0.974 1.486E+000 0.754 1.464E+000 0.769 1.464E+000 0.769 1.464E+000 0.769 1.464E+000 0.769 1.464E+000 0.769 1.464E+000 0.769 1.464E+000 0.769 1.464E+000 0.769 1.464E+000 0.769

N=M=128 9.134E+000 0.659 3.332E+000 0.779 1.127E+000 0.871 3.785E-001 0.950 1.328E-001 1.024 1.506E-001 1.082 2.011E-001 1.085 2.744E-001 1.092 3.792E-001 1.094 5.285E-001 1.094 7.391E-001 1.089 8.808E-001 0.880 8.593E-001 0.913 8.593E-001 0.913 8.593E-001 0.913 8.593E-001 0.913 8.593E-001 0.913 8.593E-001 0.913 8.593E-001 0.913 8.593E-001 0.913 8.593E-001 0.913

N=M=256 5.784E+000 0.903 1.941E+000 0.993 6.166E-001 1.060 1.959E-001 1.119 6.530E-002 1.173 7.117E-002 1.213 9.481E-002 1.213 1.288E-001 1.216 1.777E-001 1.220 2.475E-001 1.222 3.473E-001 1.224 4.785E-001 1.190 4.563E-001 1.089 4.563E-001 1.089 4.563E-001 1.089 4.563E-001 1.089 4.563E-001 1.089 4.563E-001 1.089 4.563E-001 1.089 4.563E-001 1.089 4.563E-001 1.089

N=M=512 3.094E+000

. . .

. . .

. . .

. . .

. . .

. . .

. . .

ε = 2−30

3.052E+000 0.436 1.774E+001 0.183

2.256E+000 0.623 1.562E+001 0.306

1.464E+000 0.769 1.263E+001 0.468

8.593E-001 0.913 9.134E+000 0.659

4.563E-001 1.089 5.784E+000 0.903

2.145E-001

ε = 2−0 ε = 2−1 ε=2

−2

ε=2

−3

ε = 2−4 ε=2

−5

ε=2

−6

ε = 2−7 ε = 2−8 ε=2 ε=2

−9

−10

ε = 2−11 ε=2

−12

ε=2

−13

ε = 2−14 ε = 2−15 ε=2

−16

ε=2

−17

ε = 2−18 ε=2

−19

ε=2

−20

N,M

E QN,M

20

9.757E-001 2.958E-001 9.021E-002 2.895E-002 3.070E-002 4.091E-002 5.542E-002 7.624E-002 1.061E-001 1.487E-001 2.098E-001 2.145E-001 2.145E-001 2.145E-001 2.145E-001 2.145E-001 2.145E-001 2.145E-001 2.145E-001 2.145E-001

3.094E+000