ε-Uniformly convergent fitted methods for the numerical solution of the problems arising from singularly perturbed general DDEs

Applied Mathematics and Computation 182 (2006) 119–139 www.elsevier.com/locate/amc

e-Uniformly convergent fitted methods for the numerical solution of the problems arising from singularly perturbed general DDEs Mohan K. Kadalbajoo a, Kailash C. Patidar

b,*

, Kapil K. Sharma

c

a

b

Department of Mathematics, Indian Institute of Technology Kanpur, Kanpur 208016, India Department of Mathematics and Applied Mathematics, University of the Western Cape, Private Bag X17, Bellville 7535, South Africa c MAB Universite´ Bordeaux 1, 351, cours de la Libr´ation, F-33405 Talence cedex, France

Abstract We consider some problems arising from singularly perturbed general differential difference equations. First we construct (in a new way) and analyze a ‘‘fitted operator finite difference method (FOFDM)’’ which is first order e-uniformly convergent. With the aim of having just one function evaluation at each step, attempts have been made to derive a higher order method via Shishkin mesh to which we refer as the ‘‘fitted mesh finite difference method (FMFDM)’’. This FMFDM is a direct method and e-uniformly convergent with the nodal error as Oðn2 ln2 nÞ which is an improvement over the existing direct methods (i.e., those which do not use any acceleration of convergence techniques, e.g., Richardson’s extrapolation or defect correction, etc.) for such problems on a mesh of Shishkin type that lead the error as Oðn1 ln nÞ where n denotes the total number of sub-intervals of [0, 1]. Comparative numerical results are presented in support of the theory.  2006 Elsevier Inc. All rights reserved. Keywords: Differential difference equations; Singular perturbations; Boundary value problems; Fitted operator methods; Fitted mesh methods; Shishkin mesh

1. Introduction Consider the singularly perturbed differential difference equation (SPDDE): ey 00 ðxÞ þ aðxÞy 0 ðxÞ þ aðxÞyðx  dÞ þ fðxÞyðxÞ þ bðxÞyðx þ gÞ ¼ f ðxÞ

on X ¼ ð0; 1Þ;

ð1:1Þ

under the interval and boundary conditions yðxÞ ¼ /ðxÞ

on  d 6 x 6 0;

yðxÞ ¼ cðxÞ

on 1 6 x 6 1 þ g;

*

ð1:2Þ

Corresponding author. E-mail addresses: [email protected] (M.K. Kadalbajoo), [email protected] (K.C. Patidar), [email protected] (K.K. Sharma). 0096-3003/$ - see front matter  2006 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2006.03.043

120

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

where a(x), a(x), f(x), b(x), f(x), /(x) and c(x) are sufficiently smooth functions, 0 < e  1 is the singular perturbation parameter and 0 < d = o(e), 0 < g = o(e) are the delay and the advance parameters (sometimes referred to as ‘negative shift’ and ‘positive shift’, respectively, as in [17,18]). Here it is worthwhile to mention that some of the analytical results (when d = g = 0) being proved in this paper, supplements the existing theory for FOFDMs and FMFDMs for singular perturbation problems (SPPs). The solution of (1.1) and (1.2) exhibits layer at the left or right end of the interval depending on whether a(x)  a(x)d + b(x)g > 0 or <0 on [0, 1]. In case, if a(x)  0 then one may have oscillatory solutions or two layers (one at each end) depending upon the cases whether a(x) + f(x) + b(x) is positive or negative, both of which are considered in [29]. The arguments for small delay problems are found throughout the literature on epidemics and population where these small shifts play an important role in the modeling of various real life phenomena [16]. For example, in the mathematical model for the determination of the expected first-exit time in the generation of action potential in nerve cells by random synaptic inputs in dendrites, the shifts are due to the jumps in the potential membrane which are very small. Numerous researchers have observed that often the shifts are very small and affect the solution significantly. By considering several numerical examples, Lange and Miura [18] have shown the effect of very small shifts (of the order of e) on the solution and pointed out that they drastically affect the solution and therefore cannot be neglected. Further studies of the effect of d and g on the layer behavior of the solution has been carried out by Kadalbajoo and Sharma [15]. For the occurrence and further motivation for solving such problems, the readers may refer to the introduction section in [14] and the references therein. Some of the other relevant references are [2,17,18,20,36,37]. In this paper, first we present some analytical results for the case when there is one boundary layer at the left end of the interval. The case when the layer occurs at right end, can be analyzed similarly. However, we do present some numerical results for the later case. On the other hand, we construct two types of methods: Fitted operator finite difference method (FOFDM) and fitted mesh finite difference method (FMFDM). The first aim of this paper is to show a new way of designing these FOFDMs. Like the other existing ways to get the similar schemes, its simplicity is explained by incorporating the method on some difficult problems than those considered earlier while still gaining e-uniform numerical results. The next step is about how to improve the order of convergence of the method? It is possible with the same class of methods but one requires more function evaluations at each step. See, e.g., [27] in which considering three evaluations of the function at each step, a fourth order method has been presented (here function means the right-hand side function f(x)). However, there are no other ways (to the best of our knowledge) to get a e-uniformly convergent method (in the sense of Definition 1.1) of high order with just one function evaluation per step via any FOFDM. To this end, it is natural to ask whether there are other class of methods which can give higher order convergence with only one function evaluation per step? The answer is yes. These are the carefully designed fitted mesh methods which use the same number of grid points (of course distributed differently). The FMFDM thus obtained is applied to the class of problems whose solution has a layer at left end. Since with some appropriate changes in the transition parameters of the mesh, the method is analogously applicable to the problems whose solution has a layer at right end, we omit that part. We also provide some e-uniform numerical results with this FMFDM for a problem when a(x)  0 in (1.1). As far as the SPDDEs are concerned, the idea of the FOFDM considered in this paper has successfully been implemented for the first time by the authors in [30] where they considered the problems having small delay. The FOFDM presented in this paper is obtained via so-called nonstandard finite difference method (NSFDM). The method is nonstandard due to the fact that the classical denominator h or h2 of the discrete first or second order derivative in the standard finite difference scheme is replaced by a nonnegative function / such that /ðzÞ ¼ z þ Oðz2 Þ or /ðzÞ ¼ z2 þ Oðz3 Þ as 0 < z ! 0:

ð1:3Þ

This denominator function replicates a number of significant properties of the governing differential equation. Interested readers may refer to a comprehensive list of works on NSFDMs (till early 2005) in a recent survey article [28] by the second author of this paper. On the other hand, a special section is devoted in the book chapter [21] towards the definition of these NSFDMs.

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

121

When it comes to fitted mesh methods, various options are possible for the appropriate mesh but to be as precise as possible, we restrict ourselves to the answer of the question of Stynes [39]: Miller has moved from [23] the question ‘‘what scheme should one use on a given mesh?’’ to [4] ‘‘what mesh should one use with a given scheme?’’. Various mesh selection strategies (either graded or piecewise uniform meshes) are available now-adays, however, the most successful and widely used ones are those of Bakhvalov-type and Shishkin-type meshes. Bakhvalov [1] was the first to use an a priori mesh to solve a singular perturbation problem. His mesh was generalized by Vulanovic´ [41]. A mesh of this kind is generated by a suitable smooth mesh generating function which in the layer corresponds to the inverse of the function describing the layer behavior of the continuous solution. Bakhvalov meshes have been successfully applied to many different types of singular perturbation problem (see, e.g., [42–46]). A new impetus to the a priori mesh approach was given by Shishkin [35], who discovered that simple piecewise equidistant meshes can also produce e-uniform results. A mesh of this kind uses a fine equidistant mesh inside each layer, and a coarse, again equidistant, mesh outside. The Shishkin meshes are discussed in detail in [31,38] as well as in [32] where, Roos surveyed results on layer-adapted meshes. It is shown in [46] that though Shishkin meshes are simpler than Bakhvalov meshes, the later produces much better numerical results. The reason is that Bakhvalov meshes are better adopted to the layer structure. Comparisons of these two meshes made in [19,33,46,47] have shown that the standard three-point difference schemes produce much better results on Bakhvalov meshes than on Shishkin meshes. However, it is also well known that a higher order scheme is much easier to construct on an equidistant mesh than on a fully nonequidistant one. The paper [48] shows that the real advantage of Shishkin meshes over Bakhvalov meshes is born out when applying more complicated higher order schemes. Since we are also working on higher order schemes, all of these comparisons suggest us therefore to use a mesh of Shishkin type for the problem under consideration. The rest of the paper is organized as follows. In Section 2, we derive some a priori results for the particular case of problem (1.1) when there is one boundary layer at left end. Section 3 is devoted to the construction and analysis of FOFDM. In Section 4, we design and analyze a FMFDM. The appropriate Shishkin mesh is also explained in the same section for the cases of one or two layers. Several numerical examples are given in Section 5. Discussion on the results as well as further research plans are indicated in Section 6. Notations and Terminology: Most of the notations and symbols we employ are fairly standard. However, to avoid notational complexities, sometimes we omit subscripted e, e.g., from ye and so on. Further, throughout the paper, M will denote a positive constant which may take different values in different equations (inequalities) but is always independent of e and the mesh size h. Definition 1.1 [24]. Consider a family of mathematical problems parameterized by e. Assume that each problem in the family has a unique solution denoted by ue, and that each ue is approximated by a sequence of n n numerical solutions fðme ; Xn Þg1 n , where me is defined on the mesh X :¼ fxi g0 and n is a discretization parameter. Then, the numerical solutions me are said to converge e-uniformly to the exact solution ue, if there exist a positive integer n0 and positive numbers M and p, where n0, M and p are all independent of n and e, such that for all n P n0 sup max jðue Þi  ðme Þi j 6 Mnp :

0
2. Some a priori estimates Taking the Taylor series expansions of the terms y(x  d) and y(x + g) in Eq. (1.1), we have yðx  dÞ  yðxÞ  dy 0 ðxÞ þ d2 y 00 ðxÞ=2

ð2:1Þ

yðx þ gÞ  yðxÞ þ gy 0 ðxÞ þ g2 y 00 ðxÞ=2:

ð2:2Þ

and

122

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

Using (2.1) and (2.2) in (1.1) and (1.2), we obtain C e ðxÞy 00 ðxÞ þ Ae ðxÞy 0 ðxÞ þ BðxÞyðxÞ  f ðxÞ; yð0Þ  /ð0Þ; yð1Þ  cð1Þ;

ð2:3Þ

where Ce(x)  Ce,d,g(x) = e + (a(x)d2 + b(x)g2)/2, Ae(x)  Ae,d,g(x) = a(x)  a(x)d + b(x)g and B(x) = (a(x) + f(x) + b(x)). Furthermore, the terms a(x), a(x), b(x), d and g are such that Ae(x) P A > 0 and B(x) 6 h < 0 where h is a positive constant and later on we will use the term Ce to denote the constant part of Ce(x) as and when a(x) and b(x) depend on x. (Here it is to be noted that since a(x) and b(x) are bounded and d, g are of the small order of e, we have C e ðxÞ ¼ OðeÞ.) Since (2.3) is an approximation version of (1.1), it is good to use different notation (say u(x)) for the solution of this approximate differential equation. Thus the problem (2.3) results into Lu  C e ðxÞu00 ðxÞ þ Ae ðxÞu0 ðxÞ þ BðxÞuðxÞ ¼ f ðxÞ; uð0Þ ¼ /ð0Þ ¼ /0 ;

uð1Þ ¼ cð1Þ ¼ c1 ;

ð2:4Þ

which differ from the original problem (1.1) and (1.2) by Oðd3 u000 ; g3 u000 Þ terms. Remark 2.1. It is important to note here that the case, when only delay is present, has already been treated in detail in our paper [30]. In a similar manner (as in [30]), one can study the case when only advance is present. There will merely be the changes in the location of the layers. However, the problem considered in this paper is more general than those considered there. Further, there the authors treated them via fitted operator methods only. We note that the operator L satisfies Lemma 2.2 (Continuous minimum principle). Assume that p(x) is any sufficiently smooth function satisfying p(0) P 0 and p(1) P 0. Then LpðxÞ 6 0 for all x 2 (0, 1) implies that p(x) P 0 for all x 2 [0, 1]. Proof. Let x* be such that p(x*) = minx2[0,1]p(x) and assume that p(x*) < 0. Clearly x* 62 {0, 1}, therefore p 0 (x*) = 0 and p00 (x*) P 0. Moreover, Lpðx Þ ¼ C e ðxÞp00 ðx Þ þ Ae ðx Þp0 ðx Þ þ Bðx Þpðx Þ P 0, which is a contradiction. It follows that p(x*) P 0 and thus p(x) P 0, "x 2 [0, 1]. h Furthermore, Lemma 2.3. Let u(x) be the solution of the problem (2.4), then kuk 6 h1 kf k þ maxðj/0 j; jc1 jÞ: Proof. We construct two barrier functions p± defined by p ðxÞ ¼ h1 kf k þ maxðj/0 j; jc1 jÞ  uðxÞ: Then p ð0Þ ¼ h1 kf k þ maxðj/0 j; jc1 jÞ  uð0Þ ¼ h1 kf k þ maxðj/0 j; jc1 jÞ  /0 P 0; p ð1Þ ¼ h1 kf k þ maxðj/0 j; jc1 jÞ  uð1Þ ¼ h1 kf k þ maxðj/0 j; jc1 jÞ  c1 P 0 and Lp ðxÞ ¼ C e ðxÞðp ðxÞÞ00 þ Ae ðxÞðp ðxÞÞ0 þ BðxÞp ðxÞ ¼ BðxÞðh1 kf k þ maxðj/0 j; jc1 jÞÞ  LuðxÞ ¼ BðxÞ½h1 kf k þ maxðj/0 j; jc1 jÞ  f ðxÞ 6 ðkf k  f ðxÞÞ þ BðxÞ maxðj/0 j; jc1 jÞ 6 0: Therefore, from Lemma 2.2, we obtain p±(x) P 0 for all x 2 [0, 1], which gives the required estimate.

h

Lemma 2.2 implies that the solution is unique and since the problem under consideration is linear, the existence of the solution is implied by its uniqueness. Furthermore, Lemma 2.3 gives the boundedness of the solution.

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

123

It is easy to prove that Lemma 2.4. The solution u(x) of the constant coefficients problem corresponding to (2.4) satisfies    Ax juðxÞj 6 M 1 þ exp  Ce and ðkÞ

ju ðxÞj 6 MðC e Þ

k



Ax exp  Ce

 8k P 1:

This result is useful in proving the following lemma: Lemma 2.5. The solution u(x) of (2.4) admits the decomposition uðxÞ :¼ ve;d;g ðxÞ þ we;d;g ðxÞ  ve ðxÞ þ we ðxÞ

ðsince both d; g depend on eÞ;

where the regular(smooth) component ve(x) (ur(x)) satisfies    Ax jve ðxÞj 6 M 1 þ exp  ; Ce    Ax 2k ðxÞj 6 M 1 þ ð C Þ exp  8k P 1 jvðkÞ e e Ce and the singular component we(x) (us(x)) satisfies   Ax k jwðkÞ ðxÞj 6 MðC Þ exp  8k P 0: e e Ce Proof. Plugging the three term asymptotic expansion for ve(x), viz., 2

ve ðxÞ ¼ v0 ðxÞ þ C e v1 ðxÞ þ ðC e Þ v2 ðxÞ

ð2:5Þ

into Eq. (2.4), we obtain the relations Ae ðxÞv00 ðxÞ þ BðxÞv0 ðxÞ ¼ f ðxÞ; v0 ð1Þ ¼ uð1Þ; Ae ðxÞv01 ðxÞ þ BðxÞv1 ðxÞ ¼ v000 ðxÞ; v1 ð1Þ ¼ 0; v002 ðxÞ ¼ 0;

v2 ð0Þ ¼ 0;

v2 ð1Þ ¼ 0;

Lve ¼ f ;

ve ð0Þ ¼ v0 ð0Þ þ C e v1 ð0Þ;

Lwe ¼ 0;

we ð0Þ ¼ uð0Þ  ve ð0Þ;

ve ð1Þ ¼ uð1Þ

ð2:6Þ

and we ð1Þ ¼ 0:

ð2:7Þ

Since v0 and v1 are independent of Ce, we obtain ðkÞ

jvi ðxÞj 6 M

8k P 0 and i ¼ 0; 1:

Lemma 2.4 then implies    Ax jv2 ðxÞj 6 M 1 þ exp  Ce and ðkÞ

jv2 ðxÞj 6 MðC e Þ

k

  Ax exp  Ce

ð2:8Þ

ð2:9Þ

8k P 1:

ð2:10Þ

Using (2.8)–(2.10) into (2.5) and the equations which are obtained via the successive differentiation of (2.5), we will get the desired bounds on the smooth component.

124

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

To obtain the bounds on the singular components we, consider two barrier functions g± defined by g ðxÞ ¼ ðjuð0Þj þ jve ð0ÞjÞ expðAx=C e Þ  we ðxÞ: Clearly, g±(0) and g±(1) are nonnegative. Further, Lwe ¼ 0 implies that Lg 6 0. Therefore, Lemma 2.2 implies that g± P 0 which gives   Ax jwe ðxÞj 6 M exp  ð2:11Þ ; where M ¼ ðjuð0Þj þ jve ð0ÞjÞ: Ce Now to find the bound on the derivatives of we(x), we introduce Rx Z x expðA ðtÞ=C e Þ dt  /ðxÞ ¼ R01 ðxÞ ¼ Ae ðsÞ ds: ; where A 0 expðA ðtÞ=C e Þ dt 0

ð2:12Þ

Clearly, /(0) = 0, /(1) = 1 and x 2 [0, 1] ) /(x) 6 1. Furthermore, L/ðxÞ 6 0. Therefore, from Lemma 2.2, we have /(x) P 0. Thus 0 6 /(x) 6 1. We write we(x) as we ðxÞ ¼ c1 /ðxÞ þ c2 ð1  /ðxÞÞ: Using (2.7) and the above values of /(x) at 0 and 1, we get c1 = 0, c2 = u(0)  ve(0). This implies that we ðxÞ ¼ ðuð0Þ  ve ð0ÞÞð1  /ðxÞÞ

ð2:13Þ

and therefore, jw0e ðxÞj 6 Mj/0 ðxÞj;

where M ¼ uð0Þ  ve ð0Þ:

ð2:14Þ

Eq. (2.12) then implies that jw0e ðxÞj 6 MðC e Þ

1

  Ax exp  : Ce

ð2:15Þ

Now using the relation w00e ðxÞ ¼ ðAe ðxÞw0e ðxÞ þ BðxÞwe ðxÞÞ=C e , the bounds on w00e can be obtained via (2.11) and (2.15). Finally, the successive differentiation in (2.7) and the use of the bounds on the derivatives obtained at each earlier steps gives the desired bounds on wðkÞ h e ðxÞ. This completes the proof. Remark 2.6. The first consequence of the above lemmas is the fact that  Z Z      xj  xj expðAt=C e Þ     0 juj  uj1 j ¼  u ðtÞ dt 6 M  1þ dt 6 Mðxj  xj1 Þ;    xj1  xj1 Ce meaning that the difference between the values of the exact solution at two successive grid points does not grow unboundedly. 3. Fitted operator finite difference method Many numerical methods constructed on uniform mesh have been mentioned in the extensive survey article [11]. Out of these, the most widely used ones are those mentioned in [3]. In these types of methods, one considers first the standard finite difference method in which a fitting parameter is introduced into the difference scheme. This parameter is then obtained by requiring that the proposed scheme satisfies the necessary conditions for e-uniform convergence. Apart from some other works, in [34] it has been observed that this idea does not always give the e-uniform methods. Another approach used in [12,13,27] rely on the idea of replacing the singular perturbation parameter e by the fitting factor. Then a difference scheme is obtained via a suitable standard method. Using some inverse monotonicity properties, this fitting factor is obtained then by requiring that the truncation error for the boundary layer function(s) should be equal to zero in the case of constant coefficient problems. Unlike the above two approaches, here we present a novel approach which does not require any of the above requirements (necessary condition for uniform convergence or inverse monotonicity of the matrix of

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

125

the system) and is fairly easy to construct. The construction of this new FOFDM, to which we may sometimes refer (in this paper) as nonstandard finite difference method (NSFDM), depends in many cases on the knowledge of the corresponding exact finite difference methods. Therefore, in this section, first we develop an exact scheme which will then be used to derive the FOFDM. We denote by n, a positive integer and let the interval [0, 1] be divided into n equal sub-intervals through the nodes xm = mh, m = 0(1)n where h = 1/n is the mesh size and the approximate solution to u(x) at the grid points xj0 s is denoted by mj0 s . Assume that all the coefficient functions in (2.4) are constants throughout the interval [0, 1]. When f = 0, the theory of finite difference methods yields the following exact finite difference scheme: pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi!     Ah h A2  4BC e Ah mm  exp   exp þ ð3:1Þ mmþ1 þ 2 cosh mm1 ¼ 0: 2C e 2C e 2C e This scheme is ‘‘exact’’ in the sense that it preserves most of the essential properties of the corresponding differential equation and therefore solves the problem exactly. (A major advantage of having an exact difference scheme for a differential equation is that questions related to the usual considerations of consistency, stability and convergence [9,10,25] do not arise.) For the reasons, to be mentioned shortly, we require the exact scheme for the reduced case. To this end, the exact scheme corresponding to the problem (with constant A and f) C e u00 ðxÞ þ Au0 ðxÞ ¼ f can easily be deduced from the scheme (3.1) and is given by Ce

mjþ1  2mj þ mj1 mjþ1  mj ¼ fj ; þ Aj 2 h w

ð3:2Þ

where

    hC e Ah w ¼ exp 1 : Ce A 2

Considering the differential equation (2.4) at the fixed grid points xj, we have Luðxj Þ  C e ðxj Þu00 ðxj Þ þ Ae ðxj Þu0 ðxj Þ þ Bðxj Þuðxj Þ ¼ f ðxj Þ:

ð3:3Þ

Inspired by Eq. (3.2), we may approximate Eq. (3.3) by the nonstandard scheme Lh m  C e;j

mjþ1  2mj þ mj1 mjþ1  mj þ Bj mj ¼ fj ; þ Aj 2 h wj

where the variable denominator function w2j is      3 hC e;j Ae;j h h w2j  w2j ðh; C e;j Þ :¼ exp  1 ¼ h2 þ O : C e;j Ae;j C e;j

ð3:4Þ

ð3:5Þ

Remark 3.1. The FOFDM for the case when the layer occurs near the right end (i.e., when Ae(x) < 0) can be obtained in an analogous way and will be given by C e;j

mjþ1  2mj þ mj1 mj  mj1 þ Bj mj ¼ fj ; þ Ae;j h w2j

ð3:6Þ

with the expression for w2j as the one given in (3.5) and Bj ¼ Bj , fj ¼ fj . From hereafter, the scheme (3.4) will be referred to as the FOFDM. Remark 3.2. One should note that the denominator function is derived for a simpler case. However, if not impossible, the extraction of the denominator function for the general constant coefficient homogeneous problem is not straightforward and the first choice is to rely on this function while retaining the essential

126

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

properties, e.g., the corresponding governing differential equations have the solutions possessing similar layer behavior, etc. Remark 3.3. The above FOFDM looks like the exponentially fitted schemes. However, unlike the ideas used in the construction of such exponentially fitted schemes, the idea used here is extendable for singularly perturbed partial differential equations. Since such an extension is beyond the scope of this paper, authors dealing them separately. The FOFDM (3.4) satisfies Lemma 3.4 (Discrete minimum principle). Let pi be any mesh function that satisfies p0 P 0, pn P 0 and Lh pi 6 0, i = 1(1)n  1, then pi P 0 "i = 0(1)n. Proof. The proof is obtained by contradiction. Let j be such that pj = mini pi and suppose that pj < 0. Clearly, j 62 {0, n}, pj+1  pj P 0 and pj  pj1 6 0. Therefore C e;j Ae;j ðpjþ1  pj Þ þ Bj pj Lh pj ¼ 2 ðpjþ1  2pj þ pj1 Þ þ h wj ¼

C e;j Ae;j ðpjþ1  pj Þ þ Bj pj P 0; ½ðpjþ1  pj Þ  ðpj  pj1 Þ þ 2 h wj

where the strict inequality holds if pj+1  pj > 0. This is a contradiction and therefore pj P 0. Since j is arbitrary, we have pi P 0 "i = 0(1)n. h Using this discrete minimum principle, we now show that the FOFDM (3.4) also satisfies the uniform stability result provided in Lemma 3.5 (Uniform stability estimate). The operator Lh is uniformly stable in the sense that if Zi is any mesh function such that Z0 = Zn = 0, then 1 max jLh Z j j 80 6 i 6 n: jZ i j 6 A 16j6n1 Proof. Introducing two mesh functions g i defined by 1 e e ¼ max jLh Z j j; where M g i ¼ M xi  Z i ; A 16j6n1 h h   e we see that g 0 ¼ 0 and gn ¼ M xn P 0. Since L xi ¼ Ae;i þ Bi xi , we have L gi 6 0 for 1 6 i 6 n  1 and there e xi  Z i P 0 for 0 6 i 6 n fore, from the discrete minimum principle, we have gi P 0 for 0 6 i 6 n. Hence M which gives the desired result. h This uniform stability estimate plays central role in proving the convergence of the FOFDM which is given in the following sub-section: 3.1. Convergence analysis of the FOFDM To begin with the convergence analysis, first and foremost, in line with Remark 2.6, we show that Theorem 3.6. Under the assumptions Ae(x) P A > 0, B(x) 6 h < 0, Ce(x) > 0 "x 2 [0, 1], the solution of the homogeneous system corresponding to the FOFDM (3.4) together with the given boundary conditions, is of bounded variation, uniformly for any mesh size and the parameters e, d, g and satisfies the estimate n X jmj  mj1 j 6 2C 1 kBkh;1 kmkh;1 þ ðj/0 j þ jc1 jÞ; ð3:7Þ j¼1

where C = A or kAekh,1. Here kÆkh,1 is the discrete l1-norm given by kxkh;1 ¼ max jxi j: 06i6n

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

127

Proof. Consider the mesh function pj = mj  mj1. The homogeneous part of (3.4) then leads to C e;j Ae;j pjþ1 þ Bj mj ¼ 0: ðpjþ1  pj Þ þ h w2j

ð3:8Þ

Replacing the coefficients Ce,j  e + (ajd + bjg)/2 and Ae,j involved in the above equation with the choices: e þ ðaj d þ bj gÞ=2 P e þ ðdkakh;1 þ gkbkh;1 Þ=2 0 < e þ ðaj d þ bj gÞ=2 < e þ ðad þ bgÞ=2 Ae;j P A > 0

if pjþ1  pj P 0;

if pjþ1  pj < 0;

if pjþ1 P 0;

0 < Ae;j 6 kAkh;1

if pjþ1 < 0:

and a simple rearrangement of the terms, followed by taking modulus and the use of triangle inequality yields jpjþ1 j 6 Kjpj j þ ð1  KÞ

hC 2 ; C

ð3:9Þ

where C2 = kBkh,1kmkh,1 and C = A or kAkh,1 depending upon the two situations considered on pj0 s . Further, w2 will represent w2j with these choices in Ae,j and K is given by 8 < ð1 þ Aw2 =hðe þ ðad þ bgÞ=2ÞÞ1 ; if pjþ1  pj < 0; pjþ1 P 0; K¼ ð3:10Þ : ð1 þ kAk w2 =hðe þ ðdkak þ gkbk Þ=2ÞÞ1 ; if p  p P 0; p < 0: jþ1 j jþ1 h;1  h;1 h;1 Changing j to j  1 successively, we obtain hC 2 jpj j 6 Kjpj1 j þ ð1  KÞ C      hC 2  hC 2 hC 2 hC 2 2 ¼ K jp j  6 K  jpj2 j þ ð1  KÞ þ ð1  KÞ þ j2 C  C C C

ð3:11Þ

.. . 6 CK j1 þ

hC 2 ; C

where C ¼ jp1 j 

This implies that   n X 1  Kn C2 jpj j 6 C þ 1K C j¼1

hC 2 : C

as K < 1 and nh ¼ 1:

ð3:12Þ

Now from (3.8) and (3.9), we have pjþ1 6 Kpj 

ð1  KÞhBj mj : C

ð3:13Þ

Summing up for j = 1, 2, . . . , n  1 and using m0 = /0, mn = c1, we obtain p1 6 Kpn þ ð1  KÞð/0  c1 Þ 

n1 ð1  KÞh X Bj m j ; C j¼1

) jp1 j 6 Kjpn j þ ð1  KÞðj/0 j þ jc1 jÞ þ

n1 ð1  KÞh X jBj jjmj j; C j¼1

ð1  KÞC 2 6 Kjpn j þ ð1  KÞðj/0 j þ jc1 jÞ þ C

as ðn  1Þh 6 1:

ð3:14Þ

128

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

Thus from (3.11), we obtain   hC 2 hC 2 ð1  KÞC 2 n1 6 K CK Cþ þ þ ð1  KÞðj/0 j þ jc1 jÞ þ C C C ) Cð1  K n Þ 6

ð1  KÞð1  hÞC 2 þ ð1  KÞðj/0 j þ jc1 jÞ; C

which gives C

ð1  K n Þ C 2 6 þ ðj/0 j þ jc1 jÞ 1K C

as 1  h < 1:

ð3:15Þ

The proof is accomplished then by using (3.15) in (3.12). h Very often during the analysis, one encounters with the exponential terms involving e divided by the power functions in e which are always the main cause of worry. For their careful consideration while proving the e-uniform convergence, we prove Lemma 3.7. For a fixed mesh and e ! 0, we have lim max

e!0 16j6n1

expðMxj =eÞ ¼0 ek

and

lim max

e!0 16j6n1

expðMð1  xj Þ=eÞ ¼0 ek

8k 2 Iþ ;

where xj = jh, h = 1/n "j = 1(1)n  1. Proof. Consider the partition ½0; 1 :¼ f0 ¼ x0 < x1 < x2 < < xn1 < xn ¼ 1g: Clearly, for the interior grid points, we have max

expðMxj =eÞ expðMx1 =eÞ expðMh=eÞ 6 ¼ ek ek ek

max

expðMð1  xj Þ=eÞ expðMð1  xn1 Þ=eÞ expðMh=eÞ 6 ¼ ek ek ek

16j6n1

and 16j6n1

(as x1 = h, 1  xn1 = 1  (n  1)h = h). An application of L’Hospital’s rule then leads to lim e!0

expðMh=eÞ pk k!  lim ¼ lim ¼0 p!1 ðMhÞk expðMhpÞ pð¼1=eÞ!1 expðMhpÞ ek

which completes the proof.

h

Now at the grid point xj, we consider a mesh function as uj  mj and observe that jLh ðuj  mj Þj 6 Mhu00 ðnj1 Þ þ Mh2 uðivÞ ðnj2 Þ

ð3:16Þ

nji

where xj < < xjþ1 , i = 1, 2. Denoting the smooth and singular components of the numerical solution by mr and ms, respectively, and using Lemma 2.5, we obtain ! j 2 expðAn2 =C e Þ h r r jL ðuj  mj Þj 6 M h þ h ð3:17Þ C 2e and jL

h

ðusj



msj Þj

! j expðAnj1 =C e Þ 2 expðAn2 =C e Þ 6M h þh : C 2e C 4e

ð3:18Þ

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

Combining the two, considering the dominant parts and using Lemma 3.5, we obtain ! 2 expðAxj =C e Þ max juj  mj j 6 M h þ max h : 16j6n1 16j6n1 C 4e

129

ð3:19Þ

Finally, using Lemma 3.7, we obtain Theorem 3.8. Assuming that a(x), b(x), a(x), f(x) and f(x) are sufficiently smooth so that u(x) 2 C4[0, 1], the FOFDM (3.4) is first order e-uniformly convergent in the sense that the numerical solution m of the problems (2.4) obtained via (3.4) (with m0 = u(0), mn = u(1)) satisfies the error estimate sup max juj  mj j 6 Mh:

0
ð3:20Þ

Remark 3.9. It should be noted from (3.19) that for a fixed (but not so small e), we do achieve second order convergence (though not e-uniform) which can be seen in the tabular results. However, when e is very small, one can see the e-uniform numerical results confirming the outcomes stated in Theorem 3.8 as above. 4. Fitted mesh finite difference method The main purpose of designing the FMFDM is to obtain a higher order method for the class of problems stated in (2.4). Whenever the high-order accuracy is/was expected, many researchers try/tried in various different ways to design the higher order methods. Out of the two main ways to obtain higher order methods, the first one is to use the methods based on the idea of acceleration of convergence. On a mesh of Shishkin type, the use of such techniques, e.g., Richardson’s extrapolation [26] or defect correction method [5,6] may lead to the order as Oðn2 ln2 nÞ. The second approach is to use the higher order methods directly. Again on a mesh of Shishkin type, the highest possible order of convergence for most of the above class of problems that was achieved via any direct methods was Oðn1 ln nÞ. However, in this paper we develop a direct method with which we have obtained the order as Oðn2 ln2 nÞ. The resulting finite difference scheme is applied on a piecewise uniform mesh of Shishkin type. 4.1. Description of the Shishkin Mesh We describe the mesh for the case when there is one boundary layer at the left end. We recall that n is a positive integer and we denote by s the width of the boundary layer. Let the interval [0, 1] be divided into two sub-intervals: ½0; 1 :¼ ½0; s [ ½s; 1: The piecewise uniform mesh (of Shishkin type) in these sub-intervals is designed as follows: Assuming that n = 2m with m P 3, the intervals (0, s) and (s, 1) are each divided into n/2 equal mesh elements which insures that there are n/2 + 1 equidistant grid points into the intervals [0, s] whereas one has n/2 equidistant grid points in (s, 1]. With C e ¼ C e ð0Þ, we define the parameter s by s ¼ minf1=2; 8C e ln ng: ð4:1Þ Let xj0 ¼ s and ½0; 1 :¼ 0 ¼ x0 < x1 < < xj0 < < xn ¼ 1; with xj = hj  hj1, where the mesh spacing is given by  j ¼ 1; . . . ; j0 ; 2sn1 ; hj ¼ 2ð1  sÞn1 ; j ¼ j0 þ 1; . . . ; n:

ð4:2Þ

If s = 1/2, i.e., 1=2 < 8C e ln n then n1 is very small relative to C e . This is very unlikely in practice and in such a case the method can be analyzed using the standard techniques. We therefore assume that s ¼ 8C e ln n:

ð4:3Þ

130

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

Remark 4.1. If a(x)  0 in (1.1), then B(x) has more influence than Ae(x) in (2.4) which leads to two boundary layers (one at each ends). In such a situation, while using FMFDM, we divide the interval [0, 1] into three subintervals: ½0; 1 :¼ ½0; s [ ½s; 1  s [ ½1  s; 1 and we consider n o pffiffiffiffiffiffi s ¼ min 1=4; 2 C e lnðn=4Þ :

ð4:4Þ

Then letting j0 ¼ n=4; xj0 ¼ s; xnj0 ¼ 1  s we obtain the following expressions for s and hj: pffiffiffiffiffiffi s ¼ 2 C e lnðn=4Þ

ð4:5Þ

and  hj ¼

4sn1 ;

j ¼ 1; . . . ; j0 ; n  j0 þ 1; . . . ; n;

2ð1  2sÞn1 ; j ¼ j0 þ 1; . . . ; n  j0 :

ð4:6Þ

Remark 4.2. Inclusion of the factor min Ae(x) in the definition of s in (4.1) will hardly affect the theoretical analysis. However, it does affect the numerical results. Instead of (4.1), one may go for the choice s ¼ minf1=2; ð8= min Ae ðxÞÞC e ln ng which looks more general but some more numerical experiments suggested us that depending on the particular problem, various such possibilities may occur. Similar observations are made with (4.4) too. Such information is specially important when one is testing the higher order schemes. 4.2. FMFDM on Shishkin Mesh Taking the Taylor series expansion for m around xj, and neglecting the terms containing third and higher order terms, we get the following approximations for mj+1 and mj1: mjþ1  mj þ hjþ1 m0j þ mj1  mj  hj m0j þ

h2jþ1 00 m ; 2! j

ð4:7Þ

h2j 00 m : 2! j

ð4:8Þ

Adding (4.7) and (4.8), we obtain mj1  2mj þ mjþ1  ðhjþ1 

hj Þm0j

h2jþ1 þ h2j 00 mj ; þ 2!

ð4:9Þ

whereas after solving them, we obtain m0j 

1 ½h2 mjþ1  ðh2j  h2jþ1 Þmj  h2jþ1 mj1 : hj hjþ1 ðhj þ hjþ1 Þ j

ð4:10Þ

Using (2.4) and (4.10) in (4.9) and simplifying, we obtain the tridiagonal system: Gm ¼ F ;

ð4:11Þ

where G is the matrix of the system and m and F are corresponding vectors. The various entries of this matrix and the components of the rhs vector are given by

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

ðsubdiagðGÞÞj ¼ r j ; j ¼ 2; 3; . . . ; n  1;

9 ðsupdiagðGÞÞj ¼ rþ j ; j ¼ 1; 2; . . . ; n  2; > > > =

ðdiagðGÞÞj ¼ rcj ; j ¼ 1; 2; . . . ; n  1; F 1 ¼ q1 f1  r F j ¼ qj fj ; j ¼ 2; 3; . . . ; n  2; 1 m0 ; m 0 ¼ g0 ;

131

F n1 ¼ qn1 fn1  rþ n1 mn ;

m n ¼ g1 ;

> > > ;

ð4:12Þ

where r j

¼

2C e;j h2jþ1 þh2j

þ

hjþ1 hj ðhjþ1 þhj Þ





2C e;j ðhjþ1 hj Þ h2jþ1 þh2j

 Ae;j ;   2C e;j ðhjþ1 hj Þ 2C e;j hj þ  Ae;j ; rj ¼ h2 þh2  hjþ1 ðhjþ1 þhj Þ h2jþ1 þh2j j jþ1   ðh2 h2j ÞAe;j 2C ðh hj Þ2 2C e;j  h e;jh ðhjþ1 þ Bj ; þ hjþ1jþ1 rcj ¼ 2 h2 þh 2 2 hj ðhjþ1 þhj Þ þh2 Þ jþ1

jþ1 j

j

jþ1

j

9 > > > > > > > = > > > > > > ; qj ¼ 1: >

ð4:13Þ

It is to be noted that when hj = hj+1 = h (i.e., uniform mesh throughout the region) then (4.13) reduces to r j ¼

C e;j Ae;j c 2C e;j ; r ¼  2 þ Bj ; qj ¼ 1:  2h j h2 h

ð4:14Þ

We call the method consisting of (4.11) and (4.13) as the fitted mesh finite difference method (FMFDM) whereas the method consisting of (4.11) and (4.14) is referred to as the standard finite difference method (SFDM). 4.3. Convergence analysis of FMFDM Theorem 4.3. Under the assumptions of Theorem 3.6, the solution of the homogeneous system corresponding to the FMFDM together with the given boundary conditions, is of bounded variation, uniformly for any mesh size and the parameters e and satisfies the estimate n X jmj  mj1 j 6 4C 1 kBkh;1 kmkh;1 þ ðj/0 j þ jc1 jÞ; ð4:15Þ j¼1

where C = A or kAekh,1. Proof. Consider the mesh function pj = mj  mj1 and denote the piecewise uniform mesh size by ~h. The proof follows then the similar lines as the proof of Theorem 4.3 with the additional knowledge of the fact that the maximum step size via the Shishkin mesh satisfies ~h 6 2n1 . h It is to be noted that the FMFDM above does not have a discrete minimum principle. Thus the theoretical analysis via discrete minimum principle and the uniform stability lemma (as is done in [24] and various other works) cannot be carried out. However, because of the special mesh in this paper, the system obtained from FMFDM is irreducibly diagonally dominant and therefore we will use another method of analysis based on the estimate (4.16) below. The local truncation error sj(u) of FMFDM is given by sj ðuÞ ¼ ðGuÞj  F j ¼ ðGðu  mÞÞj : Thus max juj  mj j 6 kG1 k max jsj ðuÞj: j

j

ð4:16Þ

The inequality (4.16) trivially holds for the case when j = 0 or n. On the other hand, when 1 6 j 6 n  1, we have sj ðuÞ ¼ T 0 uj þ T 1 u0j þ T 2 u00j þ T 3;1 u000 ðn1;j Þ þ T 3;2 u000 ðn2;j Þ;

ð4:17Þ

132

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

where n1,j 2 (xj1, xj), n2,j 2 (xj, xj+1) and c þ T 0 ¼ ðr j þ r j þ r j Þ  q j Bj ;

T 3;1 ¼

h3jþ1 þ r ; 3! j

T 3;2 ¼ 

 T 1 ¼ hjþ1 rþ j  hj rj qj Ae;j ;

T2 ¼

h2jþ1 þ h2j  r þ rj  C e;j qj ; 2! j 2!

h3j  r : 3! j

ð4:18Þ

Since the mesh is piecewise uniform with mesh spacing ~h, one obtains via (4.17) and (4.18): h2 u000 ðnj Þ; max jsj ðuÞj 6 M ~

ð4:19Þ

16j6n1

where nj 2 (xj1, xj+1). Moreover, by a result in [40], we have þ kG1 k 6 max fjrcj j  ðjr j j þ jrj jÞg

1

16j6n1

6 M~ h6M

as 0 < ~h < 1:

ð4:20Þ

Using (4.19) and (4.20) in (4.16), we obtain h2 u000 ðnj Þ; max juj  mj j 6 M ~

ð4:21Þ

16j6n1

We consider now two cases: either the grid point xj is inside the fine mesh, i.e., j 2 f1; . . . ; j0 g

ð4:22Þ

or it is inside the coarse mesh, i.e., j 2 fj0 þ 1; . . . ; n  1g:

ð4:23Þ

Further, we denote by me,r and me,s, respectively, the regular and the singular components of the numerical solution. Then using (4.13) along with (4.2) and (4.3) into (4.21) and Lemma 2.5, we see that ( in case of ð4:22Þ; n2 ln2 n½1 þ C 1 e expðAðxj0 ÞÞ; max jur;j  mr;j j 6 M ð4:24Þ 1 1 2 16j6n1 ð2n Þ ½1 þ C e expðAðxn1 ÞÞ; in case of ð4:23Þ: Similarly,

(

max jus;j  ms;j j 6 M

16j6n1

n2 ln2 n½C 3 e expðAðxj0 ÞÞ; 1 2

ð2n Þ

½C 3 e

expðAðxn1 ÞÞ;

in case of ð4:22Þ; in case of ð4:23Þ:

ð4:25Þ

Combining (4.24) and (4.25) and using Lemma 3.7, we obtain the main result: Theorem 4.4. Let a(x), b(x), a(x), f(x) and f(x) be sufficiently smooth functions so that u(x) 2 C3[0, 1]. Let mj, j = 0(1)n, be the approximate solution of (2.4), obtained using FMFDM with m0 = u(0), mn = u(1). Then, there is a constant M independent of e and the mesh size such that sup max juj  mj j 6 Mn2 ln2 n:

0
5. Test examples and numerical results In this section we present some numerical results corresponding to constant or variable coefficient problems with d = g = 0.5e. Example 5.1. Consider the problem e2 y 00 ðxÞ þ y 0 ðxÞ  2yðx  dÞ  5yðxÞ þ yðx þ gÞ ¼ 0; with yðxÞ ¼ 1; d 6 x 6 0;

yðxÞ ¼ 0; 1 6 x 6 1 þ g:

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

133

The exact solution of the problem (2.4) with the coefficient functions obtained from the above DDE is uðxÞ ¼ ðexpðm 1 x þ m2 Þ  expðm1 þ m2 xÞÞ=ðexpðm2 Þ  expðm1 ÞÞ;  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

where m1;2 ¼ Ae  A2e  4BC e 2C e , Ae = a  ad + bg, B = a + f + b, Ce = e2 + (d2a + g2b)/2, with a = 1, a = 2, f = 5, b = 1. Property: Boundary layer near left end. Example 5.2. Consider the problem e2 y 00 ðxÞ þ 0:5y 0 ðxÞ  3yðx  dÞ  2yðxÞ þ 2yðx þ gÞ ¼ 1; with yðxÞ ¼ 1; d 6 x 6 0;

yðxÞ ¼ 0; 1 6 x 6 1 þ g:

The exact solution of the problem (2.4) with the coefficient functions obtained from the above DDE is ðB  1Þðexpðm1 x þ m2 Þ  expðm1 þ m2 xÞÞ  ðexpðm2 xÞ  expðm1 xÞÞ 1 þ ; Bðexpðm2 Þ  expðm1 ÞÞ B  qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi

2C e , Ae = a  ad + bg, B = a + f + b, Ce = e2 + (d2a + g2b)/2, with ¼ Ae  A2e  4BC e

uðxÞ ¼ where m1;2

a = 0.5, a = 3, f = 2, b = 2. Property: Boundary layer near left end. Example 5.3. Consider the problem (2.4) with the coefficient functions obtained from e2 y 00 ðxÞ þ ð1 þ x2 Þy 0 ðxÞ  ð1 þ x2 Þyðx  dÞ  ð3 þ x3 ÞyðxÞ þ ð2 þ xÞyðx þ gÞ ¼ 1; with yðxÞ ¼ 1; d 6 x 6 0;

yðxÞ ¼ 0; 1 6 x 6 1 þ g:

Exact solution: Not known. Property: Boundary layer near left end. Example 5.4. Consider the problem (2.4) with the coefficient functions obtained from e2 y 00 ðxÞ  ð1 þ expðx2 ÞÞy 0 ðxÞ  x expðxÞyðx  dÞ  x2 yðxÞ  ð1  expðxÞyðx þ gÞ ¼ 1; with yðxÞ ¼ 0; d 6 x 6 0;

yðxÞ ¼ 1; 1 6 x 6 1 þ g:

Exact solution: Not known. Property: Boundary layer near right end. Example 5.5. Consider the problem e2 y 00 ðxÞ þ 3yðx  dÞ  8yðxÞ þ 2yðx þ gÞ ¼ 1; with yðxÞ ¼ 1; d 6 x 6 0;

yðxÞ ¼ 0; 1 6 x 6 1 þ g:

The exact solution of the problem (2.4) with the coefficient functions obtained from the above DDE has the expression as in Example 5.2 with Ae = ad + bg, B = a + f + b, Ce = e2 + (d2a + g2b)/2, a = 3, f = 8, b = 2. Property: Two boundary layers, one at each end. Having mj  mnj (the approx. solution obtained via FOFDMs or FMFDM) for different values of n and e, the maximum errors (denoted by Een;e in the case when the exact solution is known) at all the mesh points are evaluated using the formula

134

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

Een;e :¼ max juðxj Þ  mj j; 06j6n

whereas when the exact solution is not available, the maximum errors (denoted by Edn;e ) are evaluated using the double mesh principle [3] for FOFDMs: Edn;e :¼ max jmnj  m2n 2j j: 06j6n

Further we will tabulate the errors e n;e ; En ¼ max E 0
e n;e denotes Ee or Ed : where E n;e n;e

The numerical rates of convergence are computed using the formula [3]: e nk ;e = E e 2nk ;e Þ; rk  rk;e :¼ log2 ð E

k ¼ 1; 2; . . .

and the numerical rate of ‘‘uniform convergence’’ is computed using Rn :¼ log2 ðEn =E2n Þ: 6. Conclusions and future directions Based on some modeling rules [22] for the nonstandard finite difference methods, in this paper, we have designed a fitted operator finite difference method (FOFDM) which is first order e-uniformly convergent for the problems arising from SPDDEs. To improve the order of convergence while having just one function evaluation per-step, we have developed a fitted mesh finite difference method (FMFDM) which is almost second order e-uniformly convergent. It is to be noted that in the construction of FOFDM in Section 3, we have assumed that the coefficient functions are constant. So when we really have the constant coefficient homogeneous problems (see, e.g., Example 5.1) then we will have exact schemes. No numerical instability will occur and the resulting behavior can be seen from the results in Table 1. The FOFDMs (3.4) and (3.6) have been applied to test problems in Examples 5.2, 5.3 and 5.4, respectively. The e-uniform results can be seen in Tables 2–4. Table 5 shows the rates of convergence obtained via FOFDM 3.4. In Tables 6 and 7, we respectively present maximum errors and rates of convergence obtained via FMFDM. Finally, Table 8 contains the maximum errors for the problem in Example 5.5 whose solution has two layers. Note that the use of double mesh principle to evaluate the maximum errors via FMFDM (when exact solution is not known) is not appropriate because the grid points xnj and x2n 2j may not coincide due to the use of the factor ln n in the Shishkin mesh. It is known that the standard finite difference methods do not perform well in the layer regions. Just to confirm this further, graphs have been plotted for one example. In Fig. 1, we plot the difference of exact solution in Example 5.2 and the numerical solution obtained via FOFDM (3.4) and its standard finite difference analogue (in which the denominator function is replaced by h2) referred to as Standard FDM in the figure legends.

Table 1 Results for Example 5.1 (max. errors: using exact scheme (3.1)) e

n = 64

n = 128

n = 256

n = 512

n = 1024

n = 2048

21 24 28 29 210 211 212

0.55E14 0.14E15 0.69E17 0.27E18 0.79E22 0.30E28 0.38E42

0.21E13 0.23E14 0.28E16 0.35E17 0.54E19 0.93E22 0.95E29

0.83E13 0.31E14 0.56E16 0.28E16 0.85E21 0.27E18 0.66E22

0.30E12 0.23E13 0.33E15 0.69E17 0.28E16 0.69E17 0.93E22

0.78E12 0.64E13 0.56E15 0.28E16 0.11E15 0.28E16 0.54E19

0.62E12 0.27E12 0.22E14 0.56E16 0.22E15 0.22E15 0.54E19

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

135

Table 2 Results for Example 5.2 (max. errors: using FOFDM (3.4)) e

n = 64

n = 128

n = 256

n = 512

n = 1024

n = 2048

21 24 28 212 215 109 1010 1012

0.26E04 0.13E03 0.17E02 0.52E02 0.55E02 0.55E02 0.55E02 0.55E02

0.66E05 0.33E04 0.45E03 0.25E02 0.28E02 0.28E02 0.28E02 0.28E02

0.17E05 0.83E05 0.11E03 0.11E02 0.14E02 0.14E02 0.14E02 0.14E02

0.41E06 0.21E05 0.29E04 0.38E03 0.67E03 0.72E03 0.72E03 0.72E03

0.10E06 0.52E06 0.72E05 0.11E03 0.31E03 0.36E03 0.36E03 0.36E03

0.26E07 0.13E06 0.18E05 0.29E04 0.13E03 0.18E03 0.18E03 0.18E03

En

0.55E02

0.28E02

0.14E02

0.72E03

0.36E03

0.18E03

Table 3 Results for Example 5.3 (max. errors: using FOFDM (3.4)) n = 64

n = 128

n = 256

n = 512

n = 1024

n = 2048

2 24 28 212 215 1010 1012 1014

0.43E05 0.39E04 0.58E03 0.76E03 0.76E03 0.76E03 0.76E03 0.76E03

0.11E05 0.98E05 0.18E03 0.39E03 0.39E03 0.39E03 0.39E03 0.39E03

0.27E06 0.24E05 0.50E04 0.20E03 0.20E03 0.20E03 0.20E03 0.20E03

0.67E07 0.61E06 0.13E04 0.95E04 0.98E04 0.98E04 0.98E04 0.98E04

0.17E07 0.15E06 0.32E05 0.39E04 0.49E04 0.49E04 0.49E04 0.49E04

0.42E08 0.38E07 0.79E06 0.12E04 0.25E04 0.25E04 0.25E04 0.25E04

En

0.76E03

0.39E03

0.20E03

0.98E04

0.49E04

0.25E04

e 1

Table 4 Results for Example 5.4 (max. errors: using FOFDM (3.6)) e

n = 64

n = 128

n = 256

n = 512

n = 1024

n = 2048

21 24 28 212 106 108 1010 1012

0.26E04 0.19E03 0.13E02 0.13E02 0.13E02 0.13E02 0.13E02 0.13E02

0.66E05 0.49E04 0.56E03 0.68E03 0.68E03 0.68E03 0.68E03 0.68E03

0.16E05 0.12E04 0.18E03 0.34E03 0.34E03 0.34E03 0.34E03 0.34E03

0.41E06 0.30E05 0.49E04 0.17E03 0.17E03 0.17E03 0.17E03 0.17E03

0.10E06 0.76E06 0.13E04 0.84E04 0.85E04 0.85E04 0.85E04 0.85E04

0.26E07 0.19E06 0.32E05 0.35E04 0.43E04 0.43E04 0.43E04 0.43E04

En

0.13E02

0.68E03

0.34E03

0.17E03

0.85E04

0.43E04

Table 5 Results for Example 5.2 (rates of convergence: using FOFDM (3.4)) nk = 64 · 2k1, k = 1(1)5 r1

r2

r3

r4

r5

2 24 28 215 109 1010 1016

0.20E+01 0.20E+01 0.19E+01 0.96E+00 0.96E+00 0.96E+00 0.96E+00

0.20E+01 0.20E+01 0.20E+01 0.98E+00 0.98E+00 0.98E+00 0.98E+00

0.20E+01 0.20E+01 0.20E+01 0.99E+00 0.99E+00 0.99E+00 0.99E+00

0.20E+01 0.20E+01 0.20E+01 0.10E+01 0.99E+00 0.99E+00 0.99E+00

0.20E+01 0.20E+01 0.20E+01 0.10E+01 0.10E+01 0.10E+01 0.10E+01

Rn

0.96E+00

0.98E+00

0.99E+00

0.99E+00

0.10E+01

e 1

136

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

Fig. 2 contains the plots for the difference of same exact solution in Example 5.2 and the numerical solution obtained via FMFDM and SFDM (described in Section 4). The analysis for the FMFDM has been carried out differently. It is interesting to note (as is also pointed out in [7,8,24]) that the lack of a discrete minimum/maximum principle does not destroy the e-uniform convergence and this is due to the fact that we have used appropriate fitted mesh. An important issue concerning this article is the use of expansions (2.1) and (2.2) for the terms y(x  d) and y(x + g) in Eq. (1.1). Although we have provided the analysis in this paper for the problem (2.3), it does not imply that the solution of (2.3) is near the solution of (1.1). Intuitively, it does, when the shift arguments are of the small order of e because then, in general, the derivatives of the solutions of such SPDDEs have the bounds 0 jy ðkÞ ðxÞj 6 M=ðoðek ÞÞ; k 2 Iþ which implies that djy (x)j 6 M, d2jy00 (x)j 6 M, etc. However, detailed derivation of bounds on such derivatives is currently under progress. Further extensions of these methods to obtain e-uniform results for the nonlinear and the turning point problems related to SPDDEs are currently being investigated by the authors. Apart from these, the extensions Table 6 Results for Example 5.2 (max. errors: using FMFDM) e

n = 64

n = 128

n = 256

n = 512

n = 1024

n = 2048

21 24 28 210 212 214 215

0.86E04 0.78E03 0.11E01 0.11E01 0.11E01 0.12E01 0.14E01

0.21E04 0.20E03 0.38E02 0.38E02 0.38E02 0.38E02 0.38E02

0.53E05 0.49E04 0.12E02 0.12E02 0.12E02 0.12E02 0.12E02

0.13E05 0.12E04 0.39E03 0.39E03 0.39E03 0.39E03 0.39E03

0.33E06 0.31E05 0.12E03 0.12E03 0.12E03 0.12E03 0.12E03

0.84E07 0.77E06 0.36E04 0.36E04 0.36E04 0.36E04 0.36E04

En

0.14E01

0.38E02

0.12E02

0.39E03

0.12E03

0.36E04

Table 7 Results for Example 5.2 (rates of convergence: using FMFDM) nk = 64 · 2k1, k = 1(1)5 e

r1

r2

r3

r4

r5

21 24 26 28 210 212 214 215

0.20E+01 0.20E+01 0.20E+01 0.16E+01 0.16E+01 0.16E+01 0.17E+01 0.19E+01

0.20E+01 0.20E+01 0.20E+01 0.16E+01 0.16E+01 0.16E+01 0.16E+01 0.16E+01

0.20E+01 0.20E+01 0.20E+01 0.16E+01 0.17E+01 0.16E+01 0.16E+01 0.17E+01

0.20E+01 0.20E+01 0.20E+01 0.17E+01 0.17E+01 0.17E+01 0.17E+01 0.17E+01

0.20E+01 0.20E+01 0.20E+01 0.17E+01 0.17E+01 0.17E+01 0.17E+01 0.17E+01

Rn

0.19E+01

0.16E+01

0.17E+01

0.17E+01

0.17E+01

Table 8 Results for Example 5.5 (max. errors: using FMFDM) n = 64

n = 128

n = 256

n = 512

n = 1024

n = 2048

2 24 28 210 212 214 215

0.65E05 0.12E03 0.19E02 0.35E02 0.35E02 0.35E02 0.35E02

0.16E05 0.29E04 0.48E03 0.14E02 0.14E02 0.14E02 0.14E02

0.41E06 0.73E05 0.12E03 0.48E03 0.51E03 0.51E03 0.51E03

0.10E06 0.18E05 0.30E04 0.12E03 0.18E03 0.18E03 0.18E03

0.25E07 0.46E06 0.75E05 0.30E04 0.58E04 0.58E04 0.58E04

0.64E08 0.11E06 0.19E05 0.75E05 0.18E04 0.18E04 0.18E04

En

0.35E02

0.14E02

0.51E03

0.18E03

0.58E04

0.18E04

e 1

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

137

0.25 Diff. of Exact and Standard FDM Solutions Diff. of Exact and FOFDM Solutions

Difference of solutions

0.2

0.15

0.1

0.05

0

0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x Fig. 1. Errors of standard FDM and FOFDM solutions for Example 5.2 for n = 128, e = 103, d = g = 0.5e.

0.45 Diff. of Exact and SFDM Solutions Diff. of Exact and FMFDM Solutions

0.4

Difference of solutions

0.35 0.3 0.25 0.2 0.15 0.1 0.05 0 0

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

1

x Fig. 2. Errors of SFDM and FMFDM solutions for Example 5.2 for n = 128, e = 103, d = g = 0.5e.

of FOFDM for PDEs for SPPs and SPDDEs as well as that of FMFDM for SPDDEs is also currently under progress. Acknowledgement Authors wishes to thank Prof. M. Stynes for many valuable comments and suggestions.

138

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

References [1] N.S. Bakhvalov, Towards optimization of methods for solving boundary value problems in the presence of a boundary layer, Zh. Vychisl. Mat. Mat. Fiz. 9 (1969) 841–859 (in Russian). [2] M.W. Derstine, H.M. Gibbs, D.L. Kaplan, Bifurcation gap in a hybrid optical system, Phys. Rev. A 26 (1982) 3720–3722. [3] E.P. Doolan, J.J.H. Miller, W.H.A. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, 1980. [4] 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, New York, 2000. [5] A. Fro¨hner, H.-G. Roos, The e-uniform convergence of a defect correction method on a Shishkin mesh, Appl. Numer. Math. 37 (2001) 79–94. [6] A. Fro¨hner, T. Linss, H.-G. Roos, Defect correction on Shishkin-type meshes, Numer. Alg. 26 (2001) 281–299. [7] A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, On numerical experiments with central difference operators on special piecewise uniform meshes for problems with boundary layers, Notes on Numerical Fluid Mechanics, vol. 46, Vieweg-Verlag, 1994, pp. 167–176. [8] A.F. Hegarty, J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Use of central difference operators for the solution of singularly perturbed problems, Comm. Appl. Numer. Meth. 10 (1994) 297–302. [9] F.B. Hilderbrand, Finite Difference Equations and Simulations, Prentice-Hall, Englewood Cliffs, NJ, 1968. [10] M.K. Jain, Numerical Solution of Differential Equations, John Wiley & Sons, New York, 1984. [11] M.K. Kadalbajoo, K.C. Patidar, A survey of numerical techniques for solving singularly perturbed ordinary differential equations, Appl. Math. Comput. 130 (2–3) (2002) 457–510. [12] M.K. Kadalbajoo, K.C. Patidar, Exponentially fitted spline in compression for the numerical solution of singular perturbation problems, Comput. Math. Applicat. 46 (5–6) (2003) 751–767. [13] M.K. Kadalbajoo, K.C. Patidar, Exponentially fitted spline approximation method for solving self-adjoint singular perturbation problems, Int. J. Math. Math. Sci. 61 (2003) 3873–3891. [14] M.K. Kadalbajoo, K.K. Sharma, Numerical analysis of singularly perturbed delay differential equations with layer behavior, Appl. Math. Comput. 157 (1) (2004) 11–28. [15] M.K. Kadalbajoo, K.K. Sharma, Numerical analysis of boundary value problems for singularly perturbed differential-difference equations: small shifts of mixed type with rapid oscillations, Comm. Numer. Meth. Eng. 20 (2004) 167–182. [16] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, 1993. [17] C.G. Lange, R.M. Miura, Singular perturbation analysis of boundary-value problems for differential difference equations. V. Small shifts with layer behavior, SIAM J. Appl. Math. 54 (1994) 249–272. [18] C.G. Lange, R.M. Miura, Singular perturbation analysis of boundary-value problems for differential difference equations. VI. Small shifts with rapid oscillations, SIAM J. Appl. Math. 54 (1994) 273–283. [19] T. Linss, H.-G. Roos, R. Vulanovic´, Uniform pointwise convergence on Shishkin type meshes for quasilinear convection diffusion problems, SIAM J. Numer. Anal. 38 (2000) 897–912. [20] A. Longtin, J. Milton, Complex oscillations in the human pupil light reflex with mixed and delayed feedback, Math. Biosci. 90 (1988) 183–199. [21] J.M.-S. Lubuma, K.C. Patidar, Contributions to the theory of non-standard finite difference methods and applications to singular perturbation problems, in: R.E. Mickens (Ed.), Invited Chapter in the Book: Advances in the Applications of Nonstandard Finite Difference Schemes, World Scientific, Singapore, 2005, pp. 513–560. [22] R.E. Mickens, Nonstandard Finite Difference Models of Differential Equations, World Scientific, Singapore, 1994. [23] J.J.H. Miller, Construction of a FEM for a singularly perturbed problem in 2 dimensions, in: Numerische Behandlung von Differentialgleichungen, Band 2 (Tagung, Math. Forschungsinst., Oberwolfach, 1975), Int. Ser. Numer. Math., 31, Birkhau¨ser, Basel, 1976, pp. 165–169. [24] J.J.H. Miller, E. O’Riordan, G.I. Shishkin, Fitted Numerical Methods for Singular Perturbation Problems, World Scientific, Singapore, 1996. [25] A.R. Mitchell, D.F. Griffiths, Finite Difference Methods in Partial Differential Equations, Wiley, New York, 1980. [26] M.C. Natividad, M. Stynes, Richardson extrapolation for a convection–diffusion problem using a Shishkin mesh, Appl. Numer. Math. 45 (2) (2003) 315–329. [27] K.C. Patidar, High order fitted operator numerical method for the self-adjoint singular perturbation problems, Appl. Math. Comput. 171 (2005) 547–566. [28] K.C. Patidar, On the use of nonstandard finite difference methods, J. Diff. Eqns. Applicat. 11 (2005) 735–758. [29] K.C. Patidar, K.K. Sharma, Uniformly convergent nonstandard finite difference methods for singularly perturbed differential difference equations with delay and advance, Int. J. Numer. Meth. Eng. 66 (2006) 272–296. [30] K.C. Patidar, K.K. Sharma, e-Uniformly convergent nonstandard finite difference methods for singularly perturbed differential difference equations with small delay, Appl. Math. Comput. 175 (2006) 864–890. [31] H.-G. Roos, M. Stynes, L. Tobiska, Numerical Methods for Singularly Perturbed Differential Equations, Springer-Verlag, Berlin, 1996. [32] H.-G. Roos, Layer adapted grids for singular perturbation problems, Z. Angew. Math. Mech. 78 (1998) 291–309. [33] H.-G. Roos, T. Linss, Sufficient conditions for uniform convergence on layer adapted grids, Computing 64 (1999) 27–45.

M.K. Kadalbajoo et al. / Applied Mathematics and Computation 182 (2006) 119–139

139

[34] K.K. Sharma, Numerical analysis of boundary value problems for singularly perturbed differential difference equations with delay as well as advance, Ph.D. thesis, Indian Institute of Technology, Kanpur, India, 2003. [35] G.I. Shishkin, A difference scheme for a singularly perturbed parabolic equation with a discontinuous boundary condition, Zh. Vychisl. Mat. Mat. Fiz. 28 (1988) 1679–1692 (in Russian). [36] R.B. Stein, A theoretical analysis of neuronal variability, Biophys. J. 5 (1965) 173–194. [37] R.B. Stein, Some models of neuronal variability, Biophys. J. 7 (1967) 37–68. [38] M. Stynes, H.-G. Roos, The midpoint upwind scheme, Appl. Numer. Math. 23 (1997) 361–374. [39] M. Stynes, A jejune heuristic mesh theorem, Comput. Meth. Appl. Math. 3 (2003) 488–492. [40] J.M. Varah, A lower bound for the smallest singular value of a matrix, Linear Alg. Appl. 11 (1975) 3–5. [41] R. Vulanovic´, On a numerical solution of a type of singularly perturbed boundary value problem by using a special discretization mesh, Univ. Novom. Sadu Zb. Rad. Prirod. Mat. Fak. Ser. Mat. 13 (1983) 187–201. [42] R. Vulanovic´, A uniform method for quasilinear singular perturbation problems without turning points, Computing 41 (1989) 97– 106. [43] R. Vulanovic´, A second order numerical method for nonlinear singular perturbation problems without turning points, J. Comput. Math. Math. Phys. 31 (1991) 522–532. [44] R. Vulanovic´, A uniform numerical method for a class of quasilinear turning point problems, in: J.J.H. Miller, R. Vichnevetsky (Eds.), Proceedings of the 13th IMACS World Congress on Computation and Applied Mathematics, Dublin, IMACS, 1991, p. 493. [45] R. Vulanovic´, P.A. Farrell, Continuous and numerical analysis of a multiple boundary turning point problems, SIAM J. Numer. Anal. 30 (1993) 1400–1418. [46] R. Vulanovic´, Fourth order algorithms for a semilinear singular perturbation problem, Numer. Alg. 16 (1997) 117–128. [47] R. Vulanovic´, A priori meshes for singularly perturbed quasilinear two-point boundary value problems, IMA J. Numer. Anal. 21 (2001) 349–366. [48] R. Vulanovic´, A higher order scheme for quasilinear boundary value problems with two small parameters, Computing 67 (2001) 287– 303.