Hybrid difference scheme for singularly perturbed reaction-convection-diffusion problem with boundary and interior layers

Applied Mathematics and Computation 314 (2017) 237–256

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Hybrid difference scheme for singularly perturbed reaction-convection-diffusion problem with boundary and interior layersR T. Prabha a,∗, M. Chandru a,b, V. Shanthi a a b

Department of Mathematics, National Institute of Technology, Tiruchirappalli, India Department of Mathematics, Vignan’s University, Guntur-522216, Andhra Pradesh, India

a r t i c l e

i n f o

MSC: 65L10 CR G1.7 Keywords: Singular perturbation problem Boundary and interior layers Reaction-convection-diffusion Two parameter problem Hybrid difference scheme

a b s t r a c t A singularly perturbed second order ordinary differential equation having two small parameters with a discontinuous source term is considered. The presence of two parameters gives rise to boundary layers of different widths and the discontinuous source term generates interior layers on both sides of the discontinuous point. Theoretical bounds are derived. The problem is solved numerically with finite difference methods on a Shishkin mesh. The discretization combines a five point second order scheme at the interior layer together with the standard central, mid-point and upwind difference scheme for other regions. This combination is used in order to obtain almost second order convergence for the considered problem. Parameter uniform error bounds for the numerical approximation are established. Numerical results are presented to illustrate the convergence of the numerical approximations. © 2017 Elsevier Inc. All rights reserved.

1. Introduction Singularly perturbed differential equations arise in many areas of applied mathematics such as transport phenomena in chemistry and biology [1], chemical reactor theory [2–4], lubrication theory [5], elasticity [6], theory of plates and shells [7], fluid dynamics [8], D-C motors [9], etc. Certain types of problem arise in the models of chemical reactors and the governing equations, which determine chemical species concentrations and fluid temperature, based on conservation laws involving chemical reactions, diffusion, advection and external sources [10]. The differential equation depends on a small positive parameter (ε ) multiplying the highest derivative term. When the parameter tends to zero (ε → 0) the problem has a limiting solution, which is the solution of the reduced problem [11] and the regions of non-uniform convergence lie near the boundary, known as boundary layers. These problems have steep gradients in the narrow layer regions of the domain in consideration. This causes severe hurdles in the computations for classical numerical methods. In order to capture the layers, a large number of special purpose methods have been developed by the researchers to provide accurate numerical solutions, which cover second order equations with single parameter for smooth [11–13] and non-smooth data [14–17]. Two parameter Singularly Perturbed Problems (SPPs) are not considered in broad area and robust numerical methods for solving two parameter problem are not much identified when compared to one parameter SPPs. In two parameter R ∗

The second and third authors wish to thank Department of Science and Technology, New Delhi for financial support of project SR/FTP/MS-039/2012. Corresponding author. E-mail addresses: [email protected] (T. Prabha), [email protected] (M. Chandru).

http://dx.doi.org/10.1016/j.amc.2017.06.029 0 096-30 03/© 2017 Elsevier Inc. All rights reserved.

238

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

(reaction-convection-diffusion) problem both the diffusion and convection terms are multiplied by small parameters. Depending on the size of the parameters the solution of the problem may absorb exponential layers at the boundary points of the domain. Several authors have mapped convection-diffusion problem with smooth data. For instance, Shishkin and Titov [18] applied an exponential fitted difference scheme on an equidistribution mesh. O’Riordan et al. [19] proved a first order parameter uniform method with standard upwind finite difference operator. Linß [20] has suggested a streamline diffusion finite element method for reaction-diffusion-convection type problem. A robust a posteriori error estimate in the maximum norm is derived. An adaptive moving mesh method for linear 1-D parabolic reaction-convection-diffusion initialboundary value problems with two small parameters is examined [21]. Higher order finite difference scheme on a piecewise uniform mesh of Shishkin and Bakhavalov type is constructed for solving quasi-linear boundary value problems with small parameters [4]. Gracia et al. [22] have discussed a parameter uniform second order method on Shishkin mesh. A novel fitted operator method is constructed by Patidar [23] to obtain a parameter uniform first order convergence. But not many results and analysis are known about parameter uniform numerical methods to solve two parameter problems with non-smooth data. The non-smoothness would cause interior layers with different scales. Analysis and robust numerical methods to this type of problem are very challenging. Shanthi et al. [24] considered a simple fitted mesh method to solve reactionconvection-diffusion problem with non-smoothness occurring on the source term and derived first order convergence. For a similar type of problem the authors have used a hybrid difference scheme with average technique on an adaptive mesh to improve the order of convergence [25]. Bhakvalov and Shishkin meshes are studied by Mohapatra [26] for the convection-reaction-diffusion problem with discontinuity in the convection coefficient. But the results are not improved. Clavero et al. [27] have considered a one-dimensional two parameter singularly perturbed parabolic problem with the source term having discontinuity of first kind on the degeneration line. This type of problems are found in simple models of diffraction. The authors also have presented some collected works on SPPs of PDE with discontinuous data and degenerating convective terms. The problem which we have considered falls under the nature of this type. In this paper, we are interested to improve the order of accuracy for the more general form of linear one-dimensional second order SPP with a discontinuity in source term. With the motivation from Gracia et al. [22] and Shanthi et al. [24], we consider a singularly perturbed reaction-convection-diffusion problem in one dimension with a discontinuous source term of the form

Lu(x ) ≡ ε u (x ) + μa(x )u (x ) − b(x )u(x ) = f (x ), u ( 0 ) = u0 ,

∀ x ∈ (− ∪ + ),

u ( 1 ) = u1 .

(1) (2)

The notation  = [0, 1], − = (0, d ) and + = (d, 1 ) are introduced for convenience. The coefficients a(x) and b(x) are sufficiently smooth functions in  and f(x) is sufficiently smooth in (− ∪ + ) ∪ {0, 1}. Also, f(x) and its derivatives have a discontinuity at d ∈  = (0, 1 ), |[f](d)| ≤ C, 0 < ε < <1, 0 ≤ μ ≤ 1, a(x) ≥ α > 0 and b(x) ≥ β > 0. Under these assumptions, the SPP (1) and (2) has a solution u(x ) ∈ C 0 () ∩ C 1 () ∩ C 2 (− ∪ + ), when μ = 1 the problem is the well known convection-diffusion problem [16] and when μ = 0, we get the reaction-diffusion problem [14,17]. In the present study the following cases are dealt with Case (i): Case (ii):

√

√

αμ ≤ ρε and

√



√

αμ ≥ ρε, where ρ = min 



b( x ) . a (x )

Throughout this study C denotes a generic positive constant independent of nodal points, mesh size (N) and the perturbation parameters ε , μ. All functions in the continuous maximum norm are denoted by

w D = max |w(x )|, x∈D

where D is a bounded closed interval [r, s]. The discrete maximum norm is defined as

W DN = max |W (xi )|, 0≤i≤N

N

N

where D is an arbitrary mesh on D. If the domain is evident, it may be simply denoted as ||.|| since the domain D is dropped. The structure of the paper is as follows. In Section 2, we establish an existence theorem for (1) and (2), comparison principle, stability result and some priori estimates on the solution and its derivatives. Section 3 presents a hybrid finite difference scheme to solve the discrete problem, which generates robust numerical approximation to the solution. A decomposition of the discrete solution is introduced and truncation error analysis is estimated in Section 4. This analysis gears the main theoretical results presented. ε − μ uniform convergence in the maximum norm of the approximations are generated by the numerical method in Section 5. Numerical examples are provided in Section 6 to illustrate the applicability of the method with maximum pointwise errors and rate of convergence in the form of tables and graphs. The major finding of the paper is presented in conclusion.

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

239

2. A priori bounds on the solution and its derivatives Some analytic properties like, the minimum principle, uniform stability estimate and bounds for the derivatives of the solution (1) and (2), for its regular and singular components are studied in this section. Section 2 is commenced with the following existence Lemma. Lemma 1. The SPP (1) and (2) has a solution u(x) which belongs to the class C 0 () ∩ C 1 () ∩ C 2 (− ∪ + ). Proof. By using the method of construction presented by Shanthi et al. [24], a solution of the SPP (1) and (2) can be obtained.  It is convenient to introduce the notation for jump discontinuity at d for any function w(x ) as [w](d ) = w(d+ ) − w(d− ). The operator L of (1) satisfies the following minimum principle on . Lemma 2. Let us suppose that a function u(x ) ∈ C 0 () ∩ C 2 (− ∪ + ) satisfies u(0 ) ≥ 0, u(1 ) ≥ 0, Lu(x ) ≤ 0, + ) and [u ](d) ≤ 0. Then u(x ) ≥ 0, ∀ x ∈ .

∀ x ∈ (− ∪

Proof. Let x∗ be any point at which u(x) attains its minimum value in . If u(x∗ ) ≥ 0, then the result is obvious. Suppose that u(x∗ ) < 0, with the assumptions considered on the boundary value, we have either x∗ ∈ (− ∪ + ) or x∗ = d. If x∗ ∈ (− ∪ + ) then u (x∗ ) = 0, u (x∗ ) ≥ 0 and hence

Lu(x∗ ) = ε u (x∗ ) + μa(x∗ )u (x∗ ) − b(x∗ )u(x∗ ) > 0, which is a contradiction. The only possibility remaining is that x∗ = d. We have that u (d ) = 0 and that there exists a small positive number δ such that u (x) ≥ 0 on [d, d + δ ]. Though u  (x) has a discontinuity at d, it is a continuous function on the x interval [d, d + δ ]. By the fundamental theorem on calculus we have u (x ) = d+0 u (s )ds and it follows that u (x˜) ≥ 0 for some x˜ ∈ [d, d + δ ]. Hence, Lu(x˜) > 0, which is a contradiction to our assumption. Therefore, u(x) ≥ 0 for all x ∈ .



An immediate consequence of the minimum principle is the following stability result. Lemma 3 and 4 can be derived by applying the methods given in [16]. Lemma 3. Let u(x) be a solution of (1) and (2), then

  1

u  ≤ max |u(0 )|, |u(1 )|, f \{d} . β



Lemma 4. Let u(x) be the solution of the problem (1) and (2), where |u(0)| ≤ C, |u(1)| ≤ C and the solution satisfies the following bounds



C

u(k) \{d} ≤ √ k

ε



μ 1+ √ ε

k

,

0 ≤ k ≤ 4.



Before going into the details about decomposition of u(x) into regular and singular components we consider the following observations. Let F(x) be a smooth function in (− ∪ + ), F(x) and its derivatives have a jump discontinuity at d ∈ . Find u(x ) ∈ C 1 () ∩ C 2 (− ∪ + ), such that



Lu(x ) = F (x ), x ∈ (− ∪ + ), u(0 ) = p, u(1 ) = q.

(3)

It can be proved that the problem (3) has a unique solution [16]. Let



F ∗ (k ) ( x ) =

F ( k ) ( x ), x ∈ ( 0, d ), F (k ) (d− ), at x = d,

and further let u∗l (x ) be the solution of



Lu∗l (x ) = F ∗ (x ), x ∈ (0, d ), u∗l (0 ) = p, u∗l (d ) = u(d ).

Similarly, one can define u∗r (x ) on the interval [d, 1]. Now



u (x ) =

u∗l (x ), x ∈ [0, d ), u∗l (d ) = u∗r (d ), u∗r (x ), x ∈ (d, 1].

It can be verified that the solution u(x) of the BVP (1) and (2), can be decomposed as u(x ) = v(x ) + wl (x ) + wr (x ), where

Lv ( x ) = f ( x ),

x ∈ (− ∪ + ),

v(0 ) = u(0 ), v(1 ) = u(1 ), v(d− ) and v(d+ ) are chosen.

240

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

and

Lwl (x ) = 0,

x ∈ (− ∪ + ),

wl ( 0 ) = u ( 0 ) − v ( 0 ), Lwr (x ) = 0,

wl ( 1 ) = 0,

x ∈ (− ∪ + ),

and wr (0 ) = 0,

wr ( 1 ) = u ( 1 ) − v ( 1 ),

[wr ](d ) = −[v ](d ) − [wl ](d ).

[wr ](d ) = −[v](d ) − [wl ](d ),

(4)

Hence, v(x ), wl (x ) and wr (x ) are discontinuous at x = d, but by (4) their sum is in C1 (). Note that

 v− ( x ), v+ ( x ),  w− r ( x ), and wr (x ) = w+ r ( x ),

x ∈ − , x ∈ + ,

v (x ) =



wl ( x ) =

w− ( x ), l w+ l

x ∈ − ,

( x ),

x ∈ + ,

x ∈ − , x ∈ + .

√

√

αμ ≤ ρε √ 2 √ 3 Let v(x ) = v0 (x ) + ε v1 (x ) + ε v2 (x ) + ε v3 (x ), where v0 (x ), v1 (x ), v2 (x ) and v3 (x ) to be the solutions of the

Consider the case (i):

√

following problems:

−b(x )v0 (x ) = f (x ),

x ∈ (− ∪ + ), √ μ b(x )v1 (x ) = √ a(x )v0 (x ) + ε v0 (x ), x ∈ (− ∪ + ), ε √ μ b(x )v2 (x ) = √ a(x )v1 (x ) + ε v1 (x ), x ∈ (− ∪ + ), ε Choose v3 (x ) ∈ C 0 () ∩ C 1 () ∩ C 2 (− ∪ + ), such that



√ Lv3 (x ) = −√με a(x )v2 (x ) − ε v2 (x ),

x ∈ (− ∪ + ),

v3 (0 ) = v3 (1 ) = 0, v3 (d− ) and v3 (d+ ) are chosen. Let

θ1 =

 √ρα 2

√

ε,

αμ , 2ε

if if

√

√

√

√

αμ ≤ ρε, αμ ≥ ρε,

 √ρα and

θ2 =

√ 2 ε

, if

,

if

ρ

2μ

√

√

√

√

αμ ≤ ρε, αμ ≥ ρε.

The following lemmas can be proved using the methods described in [22,24]. Lemma 5. When

√

√

αμ ≤ ρε, the regular component v(x ) satisfies the following bounds



1

(k )



v \{d} ≤ C 1 + √ k−3 , 0 ≤ k ≤ 4. ( ε) Lemma 6. When (k )

wl

√



√

αμ ≤ ρε the singular components wl (x ) and wr (x ) satisfy the bounds



C Ce−θ1 x , x ∈ − ,

\{d} ≤ √ k −θ1 (x−d ) Ce , x ∈ + , ε



C

wr(k) \{d} ≤ √

ε

0 ≤ k ≤ 4,

k

Ce−θ2 (d−x ) , x ∈ − , Ce−θ2 (1−x ) , x ∈ + ,

0 ≤ k ≤ 4.



√ √ Consider the case (ii): αμ ≥ ρε 2 Let v(x ) = v0 (x ) + εv1 (x ) + ε v2 (x ) + ε 3 v3 (x ), where v0 (x ), v1 (x ), v2 (x ) and v3 (x ) be the solution of the problems:

Lμ v0 ( x ) = f ( x ),

x ∈ (− ∪ + ),

v0 (d− ), v0 (1 ) are chosen, v1 (d− ), v1 (1 ) are chosen, v2 (d− ), v2 (1 ) are chosen,

Lμ v1 (x ) = −v0 (x ), x ∈ (− ∪ + ), Lμ v2 (x ) = −v1 (x ), x ∈ (− ∪ + ),

where Lμ z = μaz − bz. Choose v3 (x ) ∈ C 0 () ∩ C 1 () ∩ C 2 (− ∪ + ), such that

(5)

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256



Lv3 (x ) = −v2 (x ),

x ∈ (− ∪ + ),

(6)

v3 (0 ) = v3 (1 ) = 0, v3 (d− ) and v3 (d+ ) are chosen. Lemma 7. When

√

√

αμ ≥ ρε, the regular component v(x ) satisfies the following bounds



(k )

v \{d} ≤ C 1 + Lemma 8. When

241

√

 ε 3−k

,

μ

0 ≤ k ≤ 4.



(7)

√

αμ ≥ ρε the singular components wl (x ) and wr (x ) satisfy the bounds

wl(k) \{d} ≤ C (k )

wr \{d} ≤

 μ k 

Ce−θ1 x , x ∈ − , −θ1 (x−d ) Ce , x ∈ + ,

ε

0 ≤ k ≤ 4,



Ce−θ2 (d−x ) , x ∈ − , 0 ≤ k ≤ 4. μk Ce−θ2 (1−x) , x ∈ + , C



The unique solution of the considered problem (1) and (2) is given by

 − v (x ) + w−l (x ) + w−r (x ), x ∈ − , (d− ) + w−r (d− ) = v+ (d+ ) + w+l (d+ ) + w+r (d+ ) at x = d, u(x ) = v− (d− ) + w− l + + v (x ) + wl (x ) + w+r (x ), x ∈ + .

3. Discrete problem In this section, an appropriate piecewise uniform Shishkin mesh for the Boundary Value Problem (BVP) (1) and (2) is introduced and hybrid difference scheme is used on this mesh to obtain the numerical solution. On  a piecewise uniform mesh of N mesh intervals is constructed as follows. For simplicity, let N be a positive even integer. The domain  is divided into six subintervals as

 = [0, τ1 ] ∪ [τ1 , d − τ2 ] ∪ [d − τ2 , d] ∪ [d, d + τ3 ] ∪ [d + τ3 , 1 − τ4 ] ∪ [1 − τ4 , 1], for some τ 1 , τ 2 , τ 3 and τ 4 . The interior points of the mesh are denoted by

    N N N := xi : 1 ≤ i ≤ − 1 ∪ xi : + 1 ≤ i ≤ N − 1 . 2

2

N

Here xN/2 = d is chosen and  = {xi }N ∪ {d}. It is fitted to the problem by choosing τ 1 , τ 2 , τ 3 and τ 4 to be the 0 following functions of N, ε and μ

    ⎧ d 2 d 2 ⎪ ⎪ τ = min , ln N , τ = min , ln N , ⎨ 1 2 4 θ1 4 θ2     1−d 2 1−d 2 ⎪ ⎪ , ln N , τ4 = min , ln N . ⎩τ3 = min 4 θ1 4 θ2

(8)

where θ 1 and θ 2 are defined in the previous section. On the subintervals [0, τ1 ], [d − τ2 , d], [d, d + τ3 ] and [1 − τ4 , 1] a uniform mesh with N/8 + 1 mesh intervals is used, whereas [τ1 , d − τ2 ] and [d + τ3 , 1 − τ4 ] have a uniform mesh with N/4 + 1 mesh intervals. The step sizes of the subintervals are denoted by H1 = 8τ1 /N, H2 = 4(d − τ1 − τ2 )/N, H3 = 8τ2 /N, H4 = 8τ3 /N, H5 = 4(1 − d − τ3 − τ4 )/N and H6 = 8τ4 /N. Let Hi = xi − xi−1 , i = 1, 2, . . . , N be the mesh step. N

On the piecewise uniform mesh  , the operator of the BVP (1) and (2) is discretized as the combinations of mid-point, upwind and central difference scheme as given below:

LNmU (xi ) ≡ εδ 2U (xi ) + μa(xi )D+U (xi ) − b(xi )U (xi ) = f (xi ), LNu U (xi ) ≡ εδ 2U (xi ) + μa(xi )D+U (xi ) − b(xi )U (xi ) = f (xi ), LNc U (xi ) ≡ εδ 2U (xi ) + μa(xi )D0U (xi ) − b(xi )U (xi ) = f (xi ), where δ 2U (xi ) =

(D+ − D− )U (xi ) Hi

, D+ U ( xi ) =

U (xi+1 ) − U (xi ) U (xi ) − U (xi−1 ) U (xi+1 ) − U (xi−1 ) , D− U ( xi ) = , D0 U ( xi ) = and Hi+1 Hi Hi + Hi+1

z(xi ) + z(xi+1 ) . 2 At the point of discontinuity xN/2 = d, we use five point second order difference scheme, the operator LtN U (xN/2 ) is defined by z ( xi ) =

242

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

LtNU (xN/2 ) ≡

−U (xN/2+2 ) + 4U (xN/2+1 ) − 3U (xN/2 ) U (xN/2−2 ) − 4U (xN/2−1 ) + 3U (xN/2 ) − = 0, 2H 2H

(9)

−3U (x ) + 4U (x + H ) − U (x + 2H ) 3U (x ) − 4U (x − H ) + U (x − 2H ) where u (x ) ≈ for forward difference operator, u (x ) ≈ for 2H 2H backward difference operator and H = min{H3 , H4 }. √ √ Further, when αμ ≤ ρε

 ⎧ ( 0 , τ1 ) , ( d − τ2 , d ) , ⎪ N ⎪ L U ( x ) , if x ∈ ⎪c i i ⎪ ⎨ (d, d + τ3 ), (1 − τ4 , 1 ), ( τ1 , d − τ2 ) , LN U ( xi ) ≡ ⎪ LNc U (xi ), if xi ∈ with 2μ a Hk < ε , k = 2, 5, ⎪ ⎪ ⎪ ( d + τ3 , 1 − τ4 ) , ⎩N Lt U (xi ), if xi = d.

For

√

√

αμ ≥ ρε

⎧  ( 0 , τ1 ) , ⎪ N ⎪ L U ( x ) , if x ∈ ⎪ i i c ⎪ ⎪ ⎪ (d, d + τ3 ), ⎪ ⎪ ( τ1 , d − τ2 ) , ⎪ ⎪ ⎪ LNmU (xi ), if xi ∈ with 2 b Hk < μα , k = 2, 5, ⎪ ⎪ ⎨ ( d + τ3 , 1 − τ4 ) ,  ( τ1 , d − τ2 ) , LN U ( xi ) ≡ ⎪ LNu U (xi ), if xi ∈ with 2 b Hk ≥ μα , k = 2, 5, ⎪ ⎪ ⎪ ( d + τ3 , 1 − τ4 ) , ⎪  ⎪ ⎪ ( d − τ2 , d ) , ⎪ ⎪ ⎪ LNmU (xi ), if xi ∈ ⎪ ⎪ ( 1 − τ4 , 1 ) , ⎪ ⎩N Lt U (xi ), if xi = d.

At the transition points τ 1 and d + τ3 the scheme is given by

⎧ ⎧ ⎪ ⎨τ1 = d , ⎪ ⎪ ⎪ 4 N ⎪ Lc U (xi ), if xi = ⎪ ⎪ 5d ⎩ ⎪ ⎪ d + τ3 = , ⎨ 4 ⎧ d LN U ( xi ) ≡ ⎨τ1 < , 2 b H2 < μα , ⎪ ⎪ 4 ⎪ LNmU (xi ), if xi = ⎪ ⎪ 5d ⎩ ⎪ ⎪ d + τ3 < , 2 b H5 < μα , ⎪ ⎪ 4 ⎩N Lu U (xi ), otherwise.

At the transition points (d − τ2 ) and (1 − τ4 ) the scheme is given by

LN U ( xi ) ≡

⎧ 3d ⎪ LNc U (xi ), if xi = d − τ2 = , for 2μ a H3 < ε , ⎪ ⎪ 4 ⎪ ⎪ ⎪ ⎪ d ⎪ ⎪ and 1 − τ4 = 1 − , 2μ a H6 < ε , ⎪ ⎪ 4 ⎪ ⎪ ⎪ 3d ⎪ N ⎪ L U (x ), if xi = d − τ2 = , for 2μ a H3 ≥ ε , ⎪ ⎨m i 4 ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ LNmU (xi ), ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩N Lu U ( xi ),

d , 4 3d if xi = d − τ2 > , 4 d and 1 − τ4 > 1 − , 4 otherwise. and 1 − τ4 = 1 −

2μ a H6 ≥ ε , for 2 b H3 < μα , 2 b H6 < μα ,

Now the discrete scheme is given by

LN U ( xi ) = Q N f ( xi ), U ( 0 ) = u ( 0 ),

for i = 1, 2, . . . , N − 1,

U ( 1 ) = u ( 1 ).

(10) (11)

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

where

243

⎧ f (x ), if LN ≡ LNc or LNu , ⎪ ⎨ i Q N f (xi ) = f (xi ), if LN ≡ LNm , ⎪ ⎩ N N 0,

if L ≡ Lt .

The matrix associated with the above Eqs. (10) and (11) is not an M-matrix. We transform the equation to establish the monotonicity property by estimating U (xN/2−2 ) and U (xN/2+2 ) for LtN U (xi ). √ √ When αμ ≤ ρε , from the operator LN c , we get

















2H3 2ε U (xN/2−2 ) = H3 f (xN/2−1 ) + + H3 b(xN/2−1 ) × U (xN/2−1 ) − 2ε − H3 μa(xN/2−1 ) H3 U (xN/2+2 ) =

2H4 2ε H4 f (xN/2+1 ) + + H4 b(xN/2+1 ) × U (xN/2+1 ) − 2ε + H4 μa(xN/2+1 ) H4









2ε + H3 μa(xN/2−1 ) U (xN/2 ) , 2H3 2ε − H4 μa(xN/2+1 ) U (xN/2 ) . 2H4

Inserting the above expressions for U (xN/2−2 ) and U (xN/2+2 ) in LtN U (xN/2 ), we obtain



LNT1 U (xN/2 ) : =



2ε − H4 μa(xN/2+1 ) 2ε + H3 μa(xN/2−1 ) −6+ U (xN/2 ) 2ε + H4 μa(xN/2+1 ) 2ε − H3 μa(xN/2−1 )

 + = If

√

√



−4ε − 2H42 b(xN/2+1 ) + 4 U (xN/2+1 ) + 2ε + H4 μa(xN/2+1 )



2H32 f (xN/2−1 ) 2H42 f (xN/2+1 ) + . 2ε − H3 μa(xN/2−1 ) 2ε + H4 μa(xN/2+1 )

H3





H3 f (xN/2−1 ) +

ε



× U (xN/2−1 ) −

2ε + H3 μa(xN/2−1 ) H3 b(xN/2−1 ) + H3 2





× U (xN/2+1 ) −



 2ε H4

+ H4 b(xN/2+1 )







2ε − H4 μa(xN/2+1 ) U (xN/2 ) . 2H4

Inserting the above expressions in LtN U (xN/2 ), we obtain



LNT2 U (xN/2 ) : =

 

= Let us define

LNT U



(xN/2 ) ≡



2ε − H4 μa(xN/2+1 ) 2ε + H3 μa(xN/2−1 ) H32 b(xN/2 ) −6+ − U (xN/2 ) 2ε + H4 μa(xN/2+1 ) ε 2ε

+ +



2ε + H3 μa(xN/2−1 ) H3 b(xN/2−1 ) − U (xN/2 ) , H3 2

2H4 U (xN/2+2 ) = H4 f (xN/2+1 ) + 2ε + H4 μa(xN/2+1 )



(12) (13)

αμ ≥ ρε from the operators LNm and LNc , we get U (xN/2−2 ) =



−4ε − 2H32 b(xN/2−1 ) + 4 U (xN/2−1 ) 2ε − H3 μa(xN/2−1 )



−4ε − 2H42 b(xN/2+1 ) + 4 U (xN/2+1 ) 2ε + H4 μa(xN/2+1 ) −2ε − H3 μa(xN/2−1 )

ε

H32 f (xN/2−1 )

ε

+

LNT1 U (xN/2 ), if LNT2 U (xN/2 ), if



H 2 b(xN/2−1 ) − 3 + 4 U (xN/2−1 ) 2ε

2H42 f (xN/2+1 ) . 2ε + μH4 a(xN/2+1 ) √

√

√

√

(14)

(15)

αμ ≤ ρε, αμ ≥ ρε.

Now, the matrix row associated with LN U (xN/2 ) transforms (10) and (11) to M-Matrix. The discrete problem is given by T

LN∗ U (xi ) = Q∗N f (xi ), U ( 0 ) = u ( 0 ),

for i = 1, 2, . . . , N − 1,

U ( 1 ) = u ( 1 ).

(16)

(17)

244

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

where

 LN∗ U (xi ) ≡

and

LNU (xi ), for i = N/2, LNT U (xi ), for i = N/2,

⎧ √ 2H32 f (xN/2−1 ) 2H42 f (xN/2+1 ) √ ⎪ ⎪ + , if αμ ≤ ρε , i = N/2, ⎪ ⎨ 2ε − H3 μa(xN/2−1 ) 2ε + H4 μa(xN/2+1 ) Q∗N f (xi ) = H32 f (xN/2−1 ) √ 2H42 f (xN/2+1 ) √ ⎪ + , if αμ ≥ ρε , i = N/2, ⎪ ε 2ε + μH4 a(xN/2+1 ) ⎪ ⎩ Q N f (xi ), if i = N/2. N

On the piecewise-uniform mesh  , the elements in the system matrix LN ∗ are given by

LN∗ U (xi ) ≡ ri−U (xi−1 ) + ricU (xi ) + ri+U (xi+1 ) = Q∗N f (xi ),

(18)

where

ε

ri− =

Hi H i

ε

ri− =

Hi H i

ε

ri− =

Hi H i

−

μa ( x i ) 2H i

,

ri+ =

,

ri+ =

ε

ri+ =

,

ε Hi+1 H i

ε Hi+1 H i

+ +

Hi+1 H i

μa ( x i ) Hi+1

μa ( x i ) Hi+1

+

μa ( x i ) 2H i

,

ric = −ri+ − ri− − b(xi ), if LN ≡ LNc ,

,

ric = −ri+ − ri− − b(xi ), if LN ≡ LNu ,

−

b( x i ) c , ri = −ri+ − ri− − b(xi ), if LN ≡ LNm , 2

⎧ −4ε − 2H 2 b(xN/2−1 ) −4ε − 2H 2 b(xN/2+1 ) ⎪ ri− = + 4, ri+ = + 4, ⎪ ⎪ 2ε − H μa(xN/2−1 )  2ε + H μa(xN/2+1 ) ⎨ ⎪ ⎪ ⎪ ⎩



H2 H2 = − − 2b(xN/2−1 ) + , 2ε − μa(xN/2−1 )H 2ε + μa(xN/2+1 )H N N if L ≡ LT1 and H = H3 = H4 ,

ric

−ri−

ri+

⎧ −2ε − H μa(xN/2−1 ) H 2 b(xN/2−1 ) −4ε − 2H 2 b(xN/2+1 ) ⎪ ⎪ ri− = − + 4, ri+ = + 4, ⎪ ⎨ 2ε 2ε + H μa(xN/2+1 )  2 ε  H 2H 2 ric = −ri− − ri+ − b(xN/2+1 ) + , if LN ≡ LNT2 and ⎪ ⎪ ε 2ε + a(xN/2+1 )μH ⎪ ⎩ H = max{H3 , H4 }. In (0, τ 1 ) and (d, d + τ3 ), note that

2μ a H1 /ε = 16μ a τ1 /(ε N ) ≤ 64 a ln N/(α N ), For

√

2μ a H4 /ε = 16μ a τ3 /(ε N ) ≤ 64 a ln N/(α N ).

(19)

√

αμ ≤ ρε in (d − τ2 , d ) and (1 − τ4 , 1 ), note that 2μ a H3 /ε = 16μ a τ2 /(ε N ) ≤ 64 a ln N/(α N ),

For

√

2μ a H6 /ε = 16μ a τ4 /(ε N ) ≤ 64 a ln N/(α N ).

(20)

√

αμ ≥ ρε in (d − τ2 , d ) and (1 − τ4 , 1 ), note that 2 b H3 /(αμ ) = 16 b τ2 /(αμN ) ≤ 64 b ln N/(αρ N ), 2 b H6 /(αμ ) = 16 b τ4 /(αμN ) ≤ 64 b ln N/(αρ N ).

At xi = d, note that

4 b H 2 < ε ,

2μ a H < ε ,

√

√

αμ ≤ ρε, √ √ 2 b H < μα , 2μ a H4 < ε , 2μ a H3 < αμ2 /ρ , if αμ ≥ ρε. if

(21)

(22)

To guarantee that the operator LN ∗ is monotone, it is imposed that



N (ln N )−1 > 64 max



a b , . α αρ

(23)

The following lemma shows that the finite difference operator LN ∗ has properties analogous to those of the differential operator L.

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

245

N N Lemma 9. Suppose that a mesh function U(xi ) satisfies U(0) ≥ 0, U(1) ≥ 0, LN ∗ U (xi ) ≤ 0, ∀ xi ∈  and LT U (xN/2 ) ≤ 0. Then N

U(xi ) ≥ 0, ∀ xi ∈  . Proof. Different discretization used in the definition of the scheme (18) are discussed. It is important to show

ri+ > 0,

ri− > 0,

ri+ + ric + ri− < 0,

(24)

which form the coefficient matrix of (18). The methods followed by Gracia et al. [22] are applied to prove the smooth region N

 \{d}. At the point of discontinuity xN/2 = d, a hybrid difference scheme LNT U (xN/2 ) is used. Here the inequalities ri+ > 0, 1 + ri− > 0 are guaranteed since 4 b H2 < ε , 2μ a H < ε using (22). In the next case when LN T2 U (xN/2 ) scheme is applied ri > 0, − ri > 0 is guaranteed, because 2 b H < μα , 2μ a H4 < ε and 2μ a H3 < α μ2 /ρ from (22). These inequalities can be verified using (23). ri+ + ric + ri− < 0 is true for both the cases (i) and (ii). All the operators defined in (18) guarantee M-matrix. Therefore, LN ∗ U (xi ) satisfies discrete minimum principle.



N

Lemma 10. If U(xi ) is the solution of (16) and (17), then |U (xi )| ≤ C, ∀ xi ∈  . Proof. Define the mesh function

ψ ( xi ) = ϕ ( xi ) ± U ( xi ), where ϕ (xi ) = max{|U (0 )|, |U (1 )|, β1 Q∗N f N }. Clearly ψ (xi ) ∈ C 0 (), ψ (0) ≥ 0, ψ (1) ≥ 0. Now for each xi ∈ N

LNc ψ (xi ) = −b(xi )

ϕ (xi ) ± LNc U (xi ) ≤ 0,

similarly

LNm ψ (xi ) ≤ 0 and LNu ψ (xi ) ≤ 0. At the point xN/2 = d

LNT1 ψ (xN/2 ) = LNT1 ϕ (xN/2 ) ± LNT1 U (xN/2 ) ≤ 0, in a similar way

LNT2 ψ (xN/2 ) ≤ 0. Applying Lemma 9, it follows that N

ψ ( xi ) ≥ 0, ∀ xi ∈  . This leads to the required result N

|U (xi )| ≤ C, ∀ xi ∈  .  4. Truncation error analysis The solution of the discrete problems (16) and (17) can be decomposed as U (xi ) = V (xi ) + Wl (xi ) + Wr (xi ). Let us denote N

the nodal error at each mesh point xi ∈  by |e(xi )| = |U (xi ) − u(xi )|. To bound the nodal error |e(xi )|, the argument is divided into two main parts. Initially, we define mesh functions V − (xi ) and V + (xi ) which approximate V(xi ) respectively to the left and right sides of the point of discontinuity xi = d. Then, we construct mesh functions Wl− (xi ), Wl+ (xi ) and Wr− (xi ), Wr+ (xi ) to approximate respectively Wl (xi ) and Wr (xi ) on each side of xi = d. Using these mesh functions the nodal error |e(xi )| is bounded separately outside and inside the layers. Define the regular solution



V ( xi ) =

V − (xi ), 1 ≤ i ≤ N/2 − 1, V + (xi ), N/2 + 1 ≤ i ≤ N − 1.

Let V − (xi ) and V + (xi ) be respectively, the solutions of the following discrete problems

∀ xi ∈ (N ∩ − ), V (0 ) = v(0 ), V (d ) = v(d− ), LN∗ V − (xi ) = f (xi ), −

−

and

LN∗ V + (xi ) = f (xi ), V + (d ) = v(d+ ),

∀ xi ∈ (N ∩ + ), V ( 1 ) = v ( 1 ).

246

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

Define the singular solution



W (xi ) = Wl (xi ) + Wr (xi ) =

(Wl− + Wr− )(xi ), 1 ≤ i ≤ N/2 − 1, (Wl+ + Wr+ )(xi ), N/2 + 1 ≤ i ≤ N − 1.

Let Wl− (xi ), Wl+ (xi ), Wr− (xi ) and Wr+ (xi ) be respectively, the solutions of the following discrete problems

LN∗ Wl− (xi ) = 0 on 1 ≤ i ≤ N/2 − 1,

Wl− (0 ) = w− ( 0 ), l

LN∗ Wl+ (xi ) = 0 on N/2 + 1 ≤ i ≤ N − 1, LN∗ Wr− (xi ) = 0 on 1 ≤ i ≤ N/2 − 1, LN∗ Wr+

Wl+ (d ) = w+ ( d ), l

Wr− (0 ) = 0,

(xi ) = 0 on N/2 + 1 ≤ i ≤ N − 1,

Wl− (d ) = w− ( d ), l

Wr+

Wl+ (1 ) = 0,

Wr− (d ) = w− r ( d ),

(d ) = 0, Wr+ (1 ) = w+r (1 ).

Now, we set the discrete solution

⎧ − − − N − ⎨V (xi ) + Wl (xi ) + Wr (xi ), xi ∈ ( ∩  ), − + U (xi ) = V − (d ) + Wl (d ) + Wr− (d ) = V + (d ) + Wl (d ) + Wr+ (d ), ⎩ + V (xi ) + Wl+ (xi ) + Wr+ (xi ), xi ∈ (N ∩ + ).

Lemma 11. We have the following bounds on Wl− (xi ), Wl+ (xi ), Wr− (xi ) and Wr+ (xi )

| Wl− (xi ) |≤ C

i 

(1 + θ1 H j )−1 = ψli− , ψl0− = C,

j=1 i 

| Wl+ (xi ) |≤ C

(1 + θ1 H j )−1 = ψli+ ,

+ ψlN/ = C, 2

j=N/2+1 N/2 

| Wr− (xi ) |≤ C

(1 + θ2 H j )−1 = ψri− ,

− ψrN/ 2 = C,

(1 + θ2 H j )−1 = ψri+ ,

+ ψrN = C,

j=i+1 N 

| Wr+ (xi ) |≤ C

j=i+1

Proof. Define the barrier functions

φli− = ψli− ± Wl− (xi ) and φri− = ψri− ± Wr− (xi ), where

ψ = − li

and

ψ = − ri

i j=1

(1 + θ1 H j )−1 , 1 ≤ i ≤ N/2,

1,

i = 0,

N/2 1,

j=i+1

(1 + θ2 H j )−1 , 0 ≤ i < N/2, i = N/2,

− and θ 1 and θ 2 are defined in Section 2. To prove, LN ∗ ψli ≤ 0,

LN∗ ψli− =

 ψli− (1 + θ1 Hi )ri− + ric + 

=

ψ

− li

ri−

 +

ric

+

ri+

− θ1



1 1 + θ1 Hi+1 Hi+1 ri+

1 + θ1 Hi+1

− − LN ∗ ψri ≤ 0, apply discrete operator (18) on ψli to get

ri+

−

Hi ri−

.

As defined earlier, various discretization methods used in the operator LN ∗ are discussed here. Consider the central difference operator

LNc ψli− =

     H H ψli−+1 2ε θ12 i+1 − 1 + 2ε θ12 − μai θ1 − 2bi + μai θ1 (1 − θ1 Hi ) i+1 + bi θ1 Hi+1 ≤ 0,

and it reduces to

2H i



2H i



LNc ψli− ≤ ψli−+1 2εθ12 − μai θ1 − bi ≤ 0. Applying the value of θ 1 for both the cases we get

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

 LNc ψli− ≤

ψli−+1

≤

ψli−+1

and

√ ρα − 2 b i − μa i √ 2 2 ε

ρα

 ρα 2

 LNc ψli− ≤ ψli−+1

247





− 2bi ≤ 0,

μ2 α ( α − ai ) − 2 bi ≤ 0. 2ε

For the upwind scheme, we show that

LNu ψli− ≤ ψli−+1



 εθ12 − μai θ1 − bi ≤ 0,

for the mid-point scheme, we prove that

LNm ψli− ≤ ψli−+1



 εθ12 − μai θ1 − bi ≤ 0.

We complete the argument by following the methods applied for central difference scheme. All the above results show that − − − − N − LN ∗ ψli ≤ 0. Also, φl0 > 0, φlN/2 > 0 and L∗ φl (xi ) ≤ 0. Now, applying Lemma 10 in [22], we prove that φl (xi ) > 0 for 0 ≤ i ≤  i N/2 and hence Wl− (xi ) ≤ C j=1 (1 + θ1 H j )−1 . Now consider the right layer barrier function ψri− , operating the discrete operator (18) on ψri− we find

LN∗ ψri− =

   Hi ri− ψri− ri− + ric + ri+ − θ2 − Hi+1 ri+ . 1 + θ2 Hi

− Applying the central difference scheme to LN ∗ ψri gives

     ψri−  2ε θ22 Hi /H i − 1 + 2ε θ22 + μai θ2 − 2bi [1 + θ2 Hi ] − 2ε θ23 Hi 1 + θ2 Hi   ≤ ψri− 2εθ22 + μai θ2 − 2bi ≤ 0,

LNc ψri− =

for the upwind scheme





LNu ψri− ≤ ψri− 2εθ22 + μai θ2 − bi ≤ 0. Applying the value of θ 2 for both the cases



LNu

ψ ≤ ψ ( ρ ai − 2bi ) ≤ ψ − ri

− ri

− ri

2b ρ− i ai

≤ 0.

When the mid-point scheme is applied

LNm ψri− ≤ ψri−



 εθ22 + μai θ2 − bi ≤ 0.

− − − N − N − In both cases, LN m ψri ≤ 0. Combining all the above results would give L∗ ψri ≤ 0. Also, φr0 > 0, φrN/2 > 0 and L∗ φr (xi ) ≤ 0. N/2 − − Now, applying the Lemma 10 in [22], we prove that φr (xi ) ≥ 0 for 0 ≤ i ≤ N/2 and hence Wr (xi ) ≤ C j=i+1 (1 + θ2 H j )−1 .

Similarly Wli+ and Wri+ in the interval (N/2 + 1, N ) can be proved to obtain the required result.



Now, we examine the truncation error at the interior mesh points xi ∈ N \{d}. In fact, when the mesh is arbitrary we have

⎧ N (3 ) (2 ) ⎪ ⎨|(Lc − L )u(xi )| ≤ ε H i u + μH i a

u , |LN∗ e(xi )| ≡ |(LNu − L )u(xi )| ≤ ε H i u(3) + μHi+1 a

u(2) ,  (3 )  ⎪ ⎩|(LN − L )u(x )| ≤ ε H u(3) + C 2 (2 ) i i ( a , a ) μHi+1 u + u , m

where H i = (

Hi +Hi+1 2

) and on a uniform mesh with step size H, we have

⎧ N 2 (4 ) 2 (3 ) ⎨|(Lc − L )u(xi )| ≤ ε H u + μH a

u , N 2 ( 4) ( 2) N |L∗ e(xi )| ≡ |(Lu − L )u(xi )| ≤ ε H u + μH a

u ,   ⎩ N |(Lm − L )u(xi )| ≤ ε H u(3) + C( a , a ) μH 2 u(3) + u(2) ,

where C( a , a ) a positive constant depends on a and a . It can be noted that the mid-point scheme has higher order truncation error when compared with the central difference operator. Though the central difference operator and the upwind difference operator have first order truncation error, ε − μ ratio and the combination of the mid-point difference scheme at the transition point will provide a numerical scheme, that has almost second-order convergence for both the cases.

248

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

Lemma 12. At each mesh point xi ∈ N , the regular component of the truncation error satisfies the following estimate

V − v ≤ CN−2 . Proof. For both the cases

√

√

αμ ≤ ρε and

√

√

αμ ≥ ρε, when the mesh is uniform (τ1 = τ2 = d/4)

|LN∗ (V − − v− )(xi )| = |LN∗ V − (xi ) − Q∗N f (xi )|           d2 d 2 +    ≤ ε δ − 2 V (xi ) + μa(xi ) D − V (xi ) dx

dx

≤ C ε (xi+1 − xi )2 (v− )(4) + μa(xi )(xi+1 − xi ) (v− )(2)

|LN∗ (V − − v− )(xi )| ≤ CN−2 , for xi ∈ (N ∩ − ). Similarly when τ3 = τ4 = (1 − d )/4, we get

|LN∗ (V + − v+ )(xi )| ≤ CN−2 , for xi ∈ (N ∩ + ). In particular, d = 1/2 is a special case of uniform mesh (τ1 = τ2 = τ3 = τ4 = 1/8). Hence

|LN∗ (V − − v− )(xi )| ≤ CN−2 , for xi ∈ (N ∩ − ), and |LN∗ (V + − v+ )(xi )| ≤ CN −2 , for xi ∈ (N ∩ + ). When the mesh is non-uniform

|LN∗ (V − − v− )(xi )| ≤ CN−2 , for xi ∈ (N ∩ − )\{τ1 , d − τ2 }, |LN∗ (V + − v+ )(xi )| ≤ CN−2 , for xi ∈ (N ∩ + )\{d + τ3 , 1 − τ4 }. At the transition points, either mid-point difference scheme or upwind difference scheme is used. For both the schemes

|LN∗ (V − v )(xi )| ≤ CN−1 (ε + N−1 ). To improve the convergence at the transition points define the barrier function N

ψ (xi ) = CN−2t (xi ) + CN−2 ± (V − v )(xi ), for xi ∈  , where

⎧ 1, ⎪ ⎪ ⎪ ( x i − τ1 ) ⎪ ⎪ 1− , ⎪ ⎪ 2 ( d − τ1 − τ2 ) ⎪ ⎪ ⎪ ( d − x ) ⎪ i ⎪ ⎨1 − 2τ , 2 t ( xi ) = ( xi − d ) ⎪ 1− , ⎪ ⎪ 2 τ3 ⎪ ⎪ 1 − τ4 − x i ⎪ ⎪ ⎪ ⎪1 − 2(1 − (d + τ3 + τ4 )) , ⎪ ⎪ ⎪ ⎩ 1 − xi , τ4

if 0 ≤ xi ≤ τ1 , if

τ1 ≤ x i ≤ d − τ2 ,

if d − τ2 ≤ xi ≤ d, if d ≤ xi ≤ d + τ3 , if d + τ3 ≤ xi ≤ 1 − τ4 , if 1 − τ4 ≤ xi ≤ 1.

Now, ψ (0) ≥ 0, ψ (1) ≥ 0. To find LN ∗ ψ (xi ), we have



εδ 2 ψ (xi ) =

0,

for xi = τ1 ,

−O(N −1 ε ),

d − τ2 ,

d + τ3 and 1 − τ4 ,

otherwise,

D0 ψ (xi ) ≤ 0 and D+ ψ (xi ) ≤ 0, which shows LN ∗ ψ (xi ) ≤ 0. Using Lemma 9, we prove that ψ (xi ) ≥ 0. Combining the above results

V − v ≤ CN−2 , ∀ xi ∈ N .  Lemma 13. Assume (23). The left singular component of the error satisfies the following estimate



Wl − wl ≤

C (N −1 ln N )2 , if 3

CN −2 ln N,

if

√

√

αμ ≤ ρε, for xi ∈ N . √ αμ ≥ ρε,

√

Proof. Consider the uniform mesh (τ1 = τ2 = τ3 = τ4 = 1/8) when

√

|LN∗ (Wl− − w−l )(xi )| = |LNc Wl− − Lw−l |     ≤ CH 2 ε (w− )(4) (x ) + CH 2 μ(w− )(3) (x ) l

l

√

αμ ≤ ρε in the interval 1 ≤ i ≤ N/2 − 1, we have

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

249

|LN∗ (Wl− − w−l )(xi )| ≤ CN−2 (ε (w−l )(4) + μ (w−l )(3) ) ≤ CN−2 /ε ≤ C (N−1 ln N )2 . Similarly for the interval N/2 + 1 ≤ i ≤ N − 1,

|LN∗ (Wl− − w−l )(xi )| ≤ C (N−1 ln N )2 . When

√

√

αμ ≥ ρε for the same domains

|LN∗ (Wl− − w−l )(xi )| = |LNc Wl− − Lw−l |     ≤ CH 2 ε (w− )(4) (x ) + CH 2 μ(w− )(3) (x ) l

l

|LN∗ (Wl− − w−l )(xi )| ≤ CN−2 (ε (w−l )(4) + μ (w−l )(3) ) < C N −2 μ4 /ε 3 ≤ C μN −2 ln N, 3

in a similar way

|LN∗ (Wl− − w−l )(xi )| ≤ C μN−2 ln3 N. For non-uniform mesh, when xi ∈ (0, τ 1 ) the truncation error is

|LN∗ (Wl− − w−l )(xi )| ≤ |LNc (Wl− − w−l )(xi )| ≤ |ε ( w− ) + a(xi )μ(w− ) − b(xi )w−l − (εδ 2 + a(xi )μD0 − b(xi ))Wl− |  2 l − (4)  l 2  ≤ CH1 ε (w ) (x ) + CH1 μ(w− )(3) (x ) l

If

√

√

l

≤ CN −2 (ετ12 (w− )(4) + μτ12 (w−l )(3) ). l

(25)

αμ ≤ ρε, from (25)

 μ |LN∗ (Wl− − w−l )(xi )| ≤ C (N−1 ln N )2 1 + √ ≤ C (N −1 ln N )2 , ε

and similarly for xi ∈ (d, d + τ3 ), we obtain the same result

If

√

|LN∗ (Wl− − w−l )(xi )| ≤ C (N−1 ln N )2 . √

αμ ≥ ρε, from (25)

|LN∗ (Wl− − w−l )(xi )| ≤ C

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

Consider the barrier function [22]

 μ

(xi ) = C N−2 + (N−1 ln N )2 (τ1 − xi ) , ε

for xi ∈ (0, τ1 ).

Here, (0), (τ 1 ) are non-negative and LN c (xi ) < 0. Applying Lemma 10 from [22], it follows that (xi ) ≥ 0. Hence

 μ 3 |(Wl− − w−l )(xi )| ≤ (xi ) ≤ C N−2 + (N−1 ln N )2 τ1 ≤ CN −2 ln N. ε

Similarly we prove for xi ∈ (d, d + τ3 ) as

|(Wl− − w−l )(xi )| ≤ CN−2 ln3 N. When xi ∈ [τ 1 , d), we obtain

|LN∗ (Wl− − w−l )(xi )| ≤ |w−l (xi )| + |Wl− (xi )|  −θ x  −2 ≤C e



≤C e

1 i

−θ1 τ1

+N

+ N −2



|LN∗ (Wl− − w−l )(xi )| ≤ CN−2 .

Similarly, when xi ∈ [d + τ3 , 1 ),

|LN∗ (Wl− − w−l )(xi )| ≤ CN−2 . Combining all the above findings, we obtain the desired result.



Following the methods described in Lemma 13, the right singular components can be obtained. Lemma 14. Assume (23). The right singular component of the error satisfies the following estimate

Wr − wr ≤ C (N−1 ln N )2 , if

√

√

αμ ≤ ρε,

√

√

αμ ≥ ρε for xi ∈ N ,

where Wr and wr are right singular components of discrete and continuous decomposition.



250

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

Lemma 15. At the point of discontinuity xN/2 = d, the error e(xN/2 ) satisfies the following estimate

 2 3/2 N  L∗ (U − u )(xN/2 ) ≤ CH /ε ,2

C min H3 /μ3 , H42 μ3 /ε

Proof. Consider the case

√

√ √ if αμ ≤ ρε , √ √ , if αμ ≥ ρε .

 3

√

αμ ≤ ρε

  2 2  N   L∗ (U − u )(xN/2 ) = LN∗ U (xN/2 ) − 2H3 f (xN/2−1 ) − 2H4 f (xN/2+1 )   2ε − H3 μa(xN/2−1 ) 2ε + H4 μa(xN/2+1 )     2H32 f (xN/2−1 ) 2H42 f (xN/2+1 )  = LtNU (xN/2 ) − − 2ε − H3 μa(xN/2−1 ) 2ε + H4 μa(xN/2+1 )      + C (L − LNc )U (xN/2−1 ) + C (L − LNc )U (xN/2+1 )

for

√

≤ CH 2 /ε 3/2 , √

αμ ≥ ρε

  2  N  N H32 f (xN/2−1 ) 2 H f ( x ) N/ 2+1 4  L∗ (U − u )(xN/2 ) = L∗ U (xN/2 ) − −  ε 2ε + H4 μa(xN/2+1 )     H 2 f (xN/2−1 ) 2H42 f (xN/2+1 )  = LtNU (xN/2 ) − 3 − ε 2ε + H4 μa(xN/2+1 )      + C (L − LNc )U (xN/2−1 ) + C (L − LNc )U (xN/2+1 )   ≤ C min H32 /μ3 , H42 μ3 /ε 3 .

By using the Lemma 4 and the procedure adopted by Cen [28], we derive the desired result.



5. Error estimate This section presents the main contribution of the article, the theorem which conveys the ε –μ uniform convergence error estimate e = U − u . Theorem 1. Let u(xi ) be the solution of the continuous problem (1) and (2) and U(xi ) be the solution of the discrete problem (16) and (17). Then, for sufficiently large N, we have



C (N −1 ln N )2 , if

U − u ≤

CN −2 (ln N )3 ,

if

√

√

√

√

αμ ≤ ρε, αμ ≥ ρε.

Proof. From the results of Lemma 4, 12–14 and using the procedure adopted in [22] it follows that



|e ( xi )| ≤

2

√

3

√

CN −2 ln N, CN −2 ln N,

√

αμ ≤ ρε, N ∀ xi ∈  \{d}. √ αμ ≥ ρε,

To prove the error at xi = d: √ √ Case (i): αμ ≤ ρε , define the discrete barrier function φd∗ (xi ) by

εδ 2 φd∗ (xi ) + μα D+ φd∗ (xi ) − βφd∗ (xi ) = 0, ∀ xi ∈ N , φd ∗ ( 0 ) = 0 , φ d ∗ ( d ) = 1 , φ d ∗ ( 1 ) = 0 . From the minimum principle on the intervals [0, d] and [d, 1], we derive that 0 ≤ φd∗ ≤ 1 [19]. Also



LN∗ φd∗ (xi ) = (a(xi ) − α )μ ≤0

φd∗ (xi+1 ) − φd∗ (xi ) Hi+1



+ (β − b(xi ) )φd∗ (xi )

∀ x i ∈ N .

Define the ancillary continuous functions u1 , u2 by

ε u1 + μa(x )u1 − b(x )u1 = 0, u1 (0 ) = 0, u1 (d ) = 1, ε u2 + μa(x )u2 − b(x )u2 = 0, u2 (d ) = 1, u2 (1 ) = 0. It is to be observed that, the solutions of the above equations are

(26)

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

 u1 (x ) = eη (d−x )

sinh λx sinh λd

where

αμ η= and λ = 2ε Define





 and u2 (x ) = eη (d−x )

251



sinh λ(1 − x ) , sinh λ(1 − d )

 (αμ )2 + 4εβ . 2ε



u1 (x ), for x ∈ (0, d ),

u˜ (x ) =

u2 (x ), for x ∈ (d, 1 ).

Now

 LNT1 u1 (xN/2 ) = =



−8ε + 6H3 μa(xN/2−1 ) u1 (xN/2 ) + 2ε − H3 μa(xN/2−1 )





4ε − 2H32 b(xN/2−1 ) − 4H3 μa(xN/2−1 ) u1 (xN/2−1 ) 2ε − H3 μa(xN/2−1 )



−8ε + 6H3 μa(xN/2−1 ) 4ε − 2H32 b(xN/2 ) − 4H3 μa(xN/2−1 ) −ηH3 sinh λ(d − H3 ) + e 2ε − H3 μa(xN/2−1 ) 2ε − H3 μa(xN/2−1 ) sinh(λd )



= C1 + C2



e(η+λ)H3 1 − e−2λd



1 − e−2λ(d−H3 )







≥ C2∗ 1 − e−2λ(d−H3 ) , similarly



LNT1 u2 (xN/2 ) =

−8ε − 6H4 μa(xN/2+1 ) 4ε − 2H42 b(xN/2+1 ) + 4H3 μa(xN/2+1 ) −ηH4 sinh λ(1 − d − H4 ) + e 2ε + H4 μa(xN/2+1 ) 2ε + H4 μa(xN/2+1 ) sinh(λd )



= C3 + C4

e−(η+λ)H4 1 − e−2λ(1−d )



1 − e−2λ(1−d−H4 )





 N   LT u2 (xN/2 ) ≥ C4∗ 1 − e−2λ(1−d−H4 ) . 1

Hence

LNT1 u˜ (xN/2 ) = LNT1 u2 (xN/2 ) − LNT1 u1 (xN/2 )

   −C = − LNT1 u2 (xN/2 ) + LNT1 u1 (xN/2 ) ≤ √ . ε

Following the analysis given in [17] on the interval [0, d] and [d, 1], we obtain

|(φd∗ − u1 )(xi )| ≤ C (N−1 ln N )2 , for i ≤ N/2, |(φd∗ − u2 )(xi )| ≤ C (N−1 ln N )2 , for i ≥ N/2, and at the point of discontinuity i.e., for i = N/2

−C LNT1 φd∗ (xi ) ≤ √ .

ε

Define the mesh function

ω1± (xi ) = C N−2 ln2 N +

CH2

ε

φd ∗ ( x i ) ± e ( x i ) .

± ± N N ± It is seen clearly that ω1± (0 ) ≥ 0, ω1± (1 ) ≥ 0, LN ∗ ω1 (xi ) ≤ 0, ∀ xi ∈  and LT ω1 (xN/2 ) ≤ 0. Lemma 9 implies that ω1 (xi ) ≥ 0, 1

N

∀ xi ∈  . Hence

|e(xi )| ≤ CN−2 ln2 N. Case (ii):

√

√

αμ ≥ ρε, define the mesh function

 ( xi ) = ω ± ( xi ) ± e ( xi ), W 2 where

(27)

252

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

⎧ CH2 ⎪ ⎨C N−2 ln3 N + 33 (xi − d − τ2 ), for xi ∈ (d − τ2 , d], μ ω2± (xi ) = 2 3 ⎪ ⎩C N−2 ln3 N + C H4 μ (d + τ − x ), for x ∈ [d, d + τ ). 3 3 i i ε3 Applying mid-point difference operator over the domain (d − τ2 , d], we find

εδ 2 ω2± (xi ) = 0, D+ ω2± (xi ) > 0, LNm ω2± (xi ) = −b(xi )ω2± (xi ) +

CH32 a(xi )

μ

< 0.

Applying central difference operator for the domain [d, d + τ3 ), we get

εδ 2 ω2± (xi ) = 0, D0 ω2± (xi ) < 0, LNc ω2± (xi ) = −b(xi )ω2± (xi ) −

CH42 μ2 a(xi )

ε2

< 0,

and

ω2± (xN/2−1 ) < 0, ω2± (xN/2 ) > 0, ω2± (xN/2+1 ) < 0. Using (22), we have

LNT2 ω2 (xN/2 ) < 0.  ( d − τ ) ≥ 0, Now, W 2 domain, we get

 (d + τ ) ≥ 0 and LN W  (x ) ≤ 0 for x ∈ (d − τ , d + τ ). Applying Lemma 9 to W  (x ) in the same W 3 2 3 i i i ∗

⎧ 3 Cτ ⎪ ⎨ 3 2 2 , for xi ∈ (d − τ2 , d], μ N e ( xi ) ≤ 3 3 ⎪ C μ τ3 ⎩ , for xi ∈ [d, d + τ3 ). ε3 N2

In particular at xN/2 = d

|e(xN/2 )| ≤ CN−2 ln3 N.

(28)

From (26)–(28), the desired result is obtained.



6. Numerical examples To show the applicability and efficiency of the present method it has been implemented to the following test problems. Consider the singularly perturbed two parameter BVPs with discontinuous source term. Example 1.

ε u (x ) + 4μu (x ) − 2u(x ) = f (x ), x ∈ (− ∪ + ), u ( 0 ) = 0, u ( 1 ) = 0, where

 f (x ) =

0.7,

for

0 ≤ x ≤ 0.5,

−0.6, for

0.5 < x ≤ 1.

Example 2.

ε u (x ) + μ(1 + x )2 u (x ) − u(x ) = f (x ), x ∈ (− ∪ + ), u ( 0 ) = 0, u ( 1 ) = 0, where

 f (x ) =

2x + 1,

for 0 ≤ x ≤ 0.5,

−(3x + 4 ), for 0.5 < x ≤ 1.

To calculate the maximum pointwise errors and orders of convergence, we use the double mesh principle, which is followed by many researchers, one can refer [12] for further details. The error of double mesh difference is defined by

E N = max |U N (xi ) − U 2N (xi )|, xi ∈ 

N

and EεN,μ = max E N , ε,μ

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

253

Fig. 1. Plots of Numerical Solution and Error for ε = 10−6 , μ = 10−4 with N = 256 for Example 1.

Fig. 2. Plots of Numerical Solution and Error for ε = 10−6 , μ = 10−2 with N = 256 for Example 1.

Fig. 3. Loglog Plots of the Maximum Pointwise Error for the Different Values of ε and μ for Example 1.

where UN (xi ) and U2N (xi ) denote respectively, the numerical solutions obtained using N and 2N mesh intervals. In addition, the parameter–robust orders of convergence are calculated from



RNε,μ

= log2

EεN,μ Eε2,Nμ



.

The numerical solutions, errors and loglog plots are given in Figs. 1–3 for Example 1 and in Figs. 4–6 for Example 2. Tables 1–4 correspond respectively to Examples 1 and 2.

254

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

Fig. 4. Plot of Numerical Solution and Error for ε = 10−6 , μ = 10−4 with N = 256 for Example 2.

Fig. 5. Plot of Numerical Solution and Error for ε = 10−6 , μ = 10−2 with N = 256 for Example 2.

Fig. 6. Loglog Plots of the Maximum Pointwise Error for the Different Values of ε and μ for Example 2. Table 1 Maximum pointwise errors and orders of convergence for μ = 10−2 for Example 1.

ε / μ2 10−2 10−4 10−6 10−8 EεN,μ RNε,μ

Number of Mesh Points N 64

128

256

512

1024

1.15520E−02 1.15770E−02 1.15770E−02 1.15770E−02 1.15770E−02 1.320557955

4.63170E−03 4.63520E−03 4.63520E−03 4.63520E−03 4.63520E−03 1.598534285

1.53060E−03 1.53020E−03 1.53020E−03 1.53020E−03 1.53060E−03 1.582278679

5.11150E−04 5.11080E−04 5.11080E−04 5.11080E−04 5.11150E−04 1.656684640

1.62120E−04 1.61970E−04 1.61970E−04 1.61970E−04 1.62120E−04 1.679624714

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

255

Table 2 Maximum pointwise errors and orders of convergence for μ = 10−2 for Example 1. μ2 / ε

10−2 10−4 10−6 10−8 EεN,μ RNε,μ

Number of Mesh Points N 64

128

256

512

1024

5.86660E−04 3.97170E−06 4.64060E−10 4.64850E−14 5.86660E−04 1.543877732

2.01200E−04 1.15610E−06 1.34900E−10 1.35120E−14 2.01200E−04 2.239212774

4.26150E−05 2.80330E−07 3.28950E−11 3.29650E−15 4.26150E−05 2.132963040

9.71570E−06 6.86360E−08 8.07540E−12 8.08520E−16 9.71570E−06 2.118609988

2.23720E−06 1.69550E−08 1.99810E−12 2.0 0 050E−16 2.23720E−06 2.086611015

Table 3 Maximum pointwise errors and orders of convergence for μ = 10−2 for Example 2.

ε / μ2

Number of Mesh Points N 64

128

256

512

1024

10−2 10−4 10−6 10−8 EεN,μ RNε,μ

1.30750E−02 1.33240E−02 1.33270E−02 1.33270E−02 1.33270E−02 1.523689724

4.58240E−03 4.63450E−03 4.63500E−03 4.63500E−03 4.63500E−03 1.747363533

1.38050E−03 1.37930E−03 1.37930E−03 1.37930E−03 1.38050E−03 1.771713943

4.04300E−04 4.01740E−04 4.01710E−04 4.01710E−04 4.04300E−04 1.698752270

1.24550E−04 1.23920E−04 1.23910E−04 1.23910E−04 1.24550E−04 1.743605938

Table 4 Maximum pointwise errors and orders of convergence for μ = 10−2 for Example 2. μ2 / ε

10−2 10−4 10−6 10−8 EεN,μ RNε,μ

Number of Mesh Points N 64

128

256

512

1024

4.58740E−03 1.81560E−05 2.49400E−07 2.50330E−09 4.58740E−03 2.323349579

9.16570E−04 2.15680E−06 3.11410E−08 3.12910E−10 9.16570E−04 2.237623596

1.94350E−04 2.18580E−07 3.88210E−09 3.91130E−11 1.94350E−04 2.104789146

4.51830E−05 1.44950E−08 4.82630E−10 4.88880E−12 4.51830E−05 2.055494821

1.08690E−05 1.40180E−09 5.96710E−11 6.11040E−13 1.08690E−05 2.029896284

7. Conclusion A singularly perturbed second order ordinary differential equation having two parameters with a discontinuous source term was examined. Theoretical bounds on the derivatives, regular and singular components of the solution were derived. A hybrid monotone finite difference scheme was constructed on Shishkin mesh to obtain second-order convergence. Parameter uniform error bounds for the numerical approximation were established. Numerical results presented in this paper support the theoretical results. Further study on the problems (1) and (2) with condition |[a](d)| ≤ C is communicated for publication. Acknowledgment The authors wish to thank all the reviewers for their valuable comments and suggestions to improve the article. References [1] E. Böhl, Finite Modelle gewöhnlicher Randwertaufgaben, 51, Teubner, 1981. [2] J. Bigge, E. Bohl, Deformations of the bifurcation diagram due to discretization, Math. Comput. 45 (172) (1985) 393–403. [3] A.B. Vasil’Eva, Asymptotic methods in the theory of ordinary differential equations containing small parameters in front of the higher derivatives, USSR Comput. Math. Math. Phys. 3 (4) (1963) 823–863. [4] R. Vulanovic´ , A higher-order scheme for quasilinear boundary value problems with two small parameters, Computing 67 (4) (2001) 287–303. [5] R.C. DiPrima, Asymptotic methods for an infinitely long slider squeeze-film bearing, J. Lubr. Technol. 90 (1) (1968) 173–183. [6] L.S. Frank, Singular Perturbations in Elasticity Theory, 1, IOS Press, 1997. [7] C. Huan-wen, Some applications of the singular perturbation method to the bending problems of thin plates and shells, Appl. Math. Mech. 5 (4) (1984) 1449–1457. [8] H. Schlichting, K. Gersten, Boundary Layer Theory, Springer Science & Business Media, 2003. [9] E. Robert, et al., Introduction to Singular Perturbations, 14, Elsevier, 2012. [10] J. Chen, R.E. O’Malley, On the asymptotic solution of a two-parameter boundary value problem of chemical reactor theory, SIAM J. Appl. Math. 26 (4) (1974) 717–729. [11] E.P. Doolan, J.J.H. Miller, W.H.A. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, 1, Boole Press, 1980. [12] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Robust Computational Techniques for Boundary Layers, CRC Press, 20 0 0. [13] H.G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer, 1996.

256

T. Prabha et al. / Applied Mathematics and Computation 314 (2017) 237–256

[14] M. Chandru, T. Prabha, V. Shanthi, A hybrid difference scheme for a second-order singularly perturbed reaction-diffusion problem with non-smooth data, Int. J. Appl. Comput. Math. 1 (1) (2015) 87–100, doi:10.1007/s40819- 014- 0004- 8. [15] M. Chandru, V. Shanthi, Fitted mesh method for singularly perturbed robin type boundary value problem with discontinuous source term, Int. J. Appl. Comput. Math. 1 (3) (2015) 491–501, doi:10.1007/s40819-015- 0030- 1. [16] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Singularly perturbed convection-diffusion problems with boundary and weak interior layers, J. Comput. Appl. Math. 166 (1) (2004) 133–151, doi:10.1016/j.cam.2003.09.033. [17] P.A. Farrell, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Singularly perturbed differential equations with discontinuous source terms, in: Proceedings of Lozenetz, 20 0 0. [18] G.I. Shishkin, V.A. Titov, A difference scheme for a differential equation with two small parameters at the derivatives, Chisl. Metody Meh. Sploshn. Sredy 7 (2) (1976) 145–155. [19] E. O’Riordan, M.L. Pickett, G.I. Shishkin, Singularly perturbed problems modeling reaction-convection-diffusion processes, Comput. Methods Appl. Math. 3 (3) (2003) 424–442, doi:10.2478/cmam- 2003- 0028. [20] T. Linß, A posteriori error estimation for a singularly perturbed problem with two small parameters, Int. J. Numer. Anal. Model. 7 (3) (2010) 491–506. [21] P. Das, V. Mehrmann, Numerical solution of singularly perturbed convection-diffusion-reaction problems with two small parameters, BIT Numer. Math. (2015) 1–26, doi:10.1007/s10543- 015- 0559- 8. [22] J.L. Gracia, E. O’Riordan, M.L. Pickett, A parameter robust second order numerical method for a singularly perturbed two-parameter problem, Appl. Numer. Math. 56 (7) (2006) 962–980, doi:10.1016/j.apnum.2005.08.002. [23] K.C. Patidar, A robust fitted operator finite difference method for a two-parameter singular perturbation problem, J. Differ. Eq. Appl. 14 (12) (2008) 1197–1214. [24] V. Shanthi, N. Ramanujam, S. Natesan, Fitted mesh method for singularly perturbed reaction-convection-diffusion problems with boundary and interior layers, J. Appl. Math. Comput. 22 (1–2) (2006) 49–65, doi:10.1007/BF02896460. [25] M. Chandru, T. Prabha, V. Shanthi, A parameter robust higher order numerical method for singularly perturbed two parameter problems with non-smooth data, J. Comput. Appl. Math. 309 (2017) 11–27. [26] J. Mohapatra, Equidistribution grids for two-parameter convection-diffusion boundary value problems, J. Math. Model. 2 (1) (2014) 1–21. [27] C. Clavero, J.L. Gracia, G.I. Shishkin, L.P. Shishkina, An efficient numerical scheme for 1d parabolic singularly perturbed problems with an interior and boundary layers, J. Comput. Appl. Math. 318 (2017) 634–645, doi:10.1016/j.cam.2015.10.031. [28] Z. Cen, A hybrid difference scheme for a singularly perturbed convection-diffusion problem with discontinuous convection coefficient, Appl. Math. Comput. 169 (1) (2005) 689–699, doi:10.1016/j.amc.2004.08.051.