# A Numerical Scheme for Singularly Perturbed Delay Differential Equations of Convection-Diffusion Type on an Adaptive Grid.

1 Introduction

A singularly perturbed delay differential equation is a differential equation in which the highest order derivative is multiplied by a small parameter and involving at least one delay term. Singular perturbation problems are generally the first approximation of the considered physical model. Hence in such cases, more realistic model should include some of the past and the future states of the system and hence, a real system should be modelled by differential equations with delay or advance. Such type of equation arises frequently in the mathematical modelling of various practical phenomena, for example, in the modelling of the human pupil-light reflex [26], model of HIV infection [10,33], the study of bistable devices in digital electronics [11], variational problem in control theory [17, 18], first exist time problems in modelling of activation of neuronal variability [37], immune response [30], evolutionary biology [32], dynamics of networks of two identical amplifiers [9], mathematical ecology [19], population dynamics [20] and in a variety of models for physiological process [27,28,29]. The theory and numerical treatment of delay differential equations can be found in [7,12].

The numerical solution of singularly perturbed delay differential equations with large delays can be found in Lange and Miura [22], Amiraliyev and Erdogan [2], Amiraliyev and Cimen [1], Amiraliyeva et al. [3], Erdogan and Amiraliyev [13], Subburayan and Ramanujam [34], Bansal et. al. [5,6].

It is well-known that numerical methods for singularly perturbed boundary value problems have to be very carefully created, due to boundary or interior layers in the solution. Finite difference methods are always a convenient choice for solving boundary value problems because of their simplicity. Whenever central finite difference schemes are applied to solve the singularly perturbed differential equations numerically on uniform meshes, oscillations are observed in the numerical solution and their magnitude increases in layers regions. The presence of oscillations in the approximated solution shows that central finite difference operators are unstable. To eliminate these oscillations while retaining the order of accuracy, one needs a fine mesh at the layer regions. This may be done either via uniformly fine meshing or via adaptive mesh strategy. The former strategy increases significantly in computational cost as the perturbation parameter e decreases. The use of adaptive mesh refinement techniques is nowa-days a standard component in numerical computation. In fact, Bakhvalov [4] was the first to use layer adapted meshes in the context of reaction-diffusion problems. In the late 1970s and early 1980s, the special meshes for convection-diffusion problems have been investigated by Gartland [16], Liseikin [25], Vulanovic [35, 36] and others in order to achieve uniform convergence. In early 1990s, special piecewise-uniform meshes have been proposed by Shishkin [31]. Because of their simple structure, they have attracted much attention and are widely used for numerical approximation of SPPs. The limitation with Shishkin meshes is the requirement of a priori knowledge of the location of the layer regions. However, the performance of Shishkin meshes is inferior to that of Bakhvalov meshes, which has prompted efforts to improve them while retaining some of their simplicity. A detailed survey about layer adapted meshes for convection diffusion problems can be found in [24]. The aim of this paper is to improve the performance of the adaptive mesh such that there will be more mesh points in the layer region. In this paper, we generate the layer adaptive meshes to suppress the oscillations, via an entropy production operator which is positive automatically whenever the solution is unphysical, which effectively corresponds to insufficient resolution. We restrict ourselves to the class of delay differential equations which correspond to second order boundary value problems with a non-zero coefficient of the convection term.

The outline of this paper is as follows: In Section 2, we state singularly perturbed delay differential equation of convection-diffusion type. In Section 3, we introduce an adaptive mesh and use classical finite difference scheme to solve singularly perturbed delay differential equation of convection-diffusion type. In Section 4, we provide the error analysis for the proposed method. In Section 5, three numerical examples have been solved to demonstrate the applicability and efficiency of the proposed method and confirms the theoretical estimates. The paper ends with Section 6 with a brief conclusion.

2 Statement of the problem

In this paper, we shall study the following singularly perturbed delay differential equation of convection-diffusion type:

Lu(x) = -[epsilon]u"(x) + a(x)u'(x) + b(x)u(x - 1) = f(x), 0 < x < 2 (2.1)

with the interval and boundary conditions,

u(x) = [phi](x), x [member of] [-1, 0], u(2) = [beta], (2.2)

where 0 < [epsilon] [much less than] 1 and a(x) [less than or equal to] [alpha] > 0, [[beta].sub.0] < b(x) < 0, 2[alpha] + 5[[beta].sub.0] [greater than or equal to] [gamma] > 0, a(x),b(x), f (x) are given sufficiently smooth functions on [bar.[OMEGA]], [OMEGA] = (0, 2), [bar.[OMEGA]] = [0, 2], [[OMEGA].sup.-] = (0,1), [[OMEGA].sup.+] = (1, 2), [phi](x) is a smooth function on [-1, 0] and [beta] is a given constant which is independent of [epsilon].

To generate the adaptive mesh, we followed the steps of Kumar and Srinivasan [21]. We now define an auxiliary equation i.e., the entropy production equation by multiplying with an appropriate test function. From the case of scaler hyperbolic conservation laws in Leveque [23], we know that for scalar conservation laws, [u.sup.2] is always an appropriate entropy variable and, therefore, 2u(x) is a suitable multiplying test function. On multiplying equation (2.1) (componentwise), we obtain

Lu * 2u = f * 2u, (3.1)

after simplifying, we can write the above equation (3.1) as

-[epsilon]([xi]" - 2[(u').sup.2]) + a[xi]' + 2buu(x - 1) - 2uf = 0, where [xi] = [u.sup.2].

Rewriting the above equation, we get

-[epsilon][xi]" + a[xi]' - 2uf + 2buu(x - 1) = -2[epsilon][(u').sup.2]. (3.2)

We label the linear operator on the left hand side (LHS) in the above equation (3.2), as the entropy production operator with analogy to similar operators in the hyperbolic conservation laws in Leveque [23]. The continuous operator should obviously be negative for all x [member of] [0, 2] (as the right hand side (RHS) is always negative for all x [member of] [0, 2]).

3.1 Discretization of entropy production operator

We know that if we solve the discrete problem of equation (2.1)-(2.2) computationally using central difference method, we get non-physical oscillations inside and near the boundary layer region. Similarly, if we calculate the discrete problem of the left hand side equation in (3.2) using the same central difference operator by taking [[xi].sub.i] = [u.sup.2.sub.i], where u is the central difference computed solution for equation (2.1), we observe that LHS is negative wherever the solution is smooth enough and positive where we have boundary layers (or oscillation in the computed solution of equation (2.1)). After investigation, we have found that if we write the RHS term - 2[epsilon][(u').sup.2] of equation (3.2) in the difference operator at the ith mesh point, we get

-2[epsilon]([u.sub.i] - [u.sub.i-1]/[x.sub.i] - [x.sub.i-1]) ([u.sub.i+1] - [u.sub.i-1]/[x.sub.i+1] - [x.sub.i]) (3.3)

by applying a combination of forward and backward difference operator for u which is of opposite sign wherever the oscillations occur and hence results in positive value of the difference operator as in equation (3.2) for the oscillatory numerical solution [u.sub.i]. We exploit this property of the entropy production operator to precisely determine regions of insufficient resolution.

3.2 Mesh selection strategy

The following algorithm is proposed to generate an adaptive mesh for solving second order singularly perturbed delay differential equations:

Step 1: Choose an initial grid with few uniform mesh points.

Step 2: Calculate the entropy on the mesh points using (3.3).

Step 3: Identify mesh points where entropy is positive and then choose the mesh point with the maximum entropy.

Step 4: Add one mesh point on the left and the right side of the location found in Step 3 and then calculate the entropy using (3.3) at each point on newly generated non-uniform mesh.

Step 5: If the entropy is positive for at least one mesh point, go to the Step

3. Otherwise stop the iterative process. The resulting mesh is our adapted mesh which is represented by [N.sub.g].

3.3 Finite difference scheme

Initially, we discretize the interval [0,2] into few equal parts with mesh spacing h. In this case, we discretize in such a way that the delay x = 1 is a mesh point. As mentioned in the previous section, keep on adding the points on both sides of the mesh point where the entropy is maximum and positive in each iteration, we get a stage in which the entropy is negative at each mesh point through out the interval. We assume the final non-uniform (or adaptive mesh [N.sub.g]) mesh as [[bar.[OMEGA]].sup.N] = {0 = [x.sub.0] < [x.sub.1] < .... < [x.sub.2N] = 2}. In each iteration, we find the index m such that xm = m.hj = 1. We discretize the delay differential equation (2.1)-(2.2), using central finite difference scheme on a non-uniform mesh as follows:

[mathematical expression not reproducible], (3.4)

where,

[mathematical expression not reproducible],

and the second order centered difference operator [[delta].sup.2] is defined by

[mathematical expression not reproducible].

Since, [h.sub.i] = [x.sub.i+1] - [x.sub.i], [h.sub.i-1] = [x.sub.i] - [x.sub.i-1] and [h.sub.i] + [h.sub.i-1] = [x.sub.i+1] - [x.sub.i-1], the equation (3.4), can be written in the form

[mathematical expression not reproducible].

Multiplying both sides by (hi + hi-1)hihi-1, we get

[mathematical expression not reproducible]. (3.5)

The equation (3.5) can be rewritten in the form

[mathematical expression not reproducible], (3.6)

with the boundary conditions,

[u.sub.0] = [phi](0), [U.sub.2N] = [beta], (3.7)

where

[mathematical expression not reproducible],

and

[mathematical expression not reproducible].

We have solved the system of equations (3.6) with the boundary conditions (3.7) by Gauss elimination method with partial pivoting.

4 Error analysis

4.1 Truncation error

By Taylor series expansion we have,

[mathematical expression not reproducible]

The truncation error at each nodal point is given by

[mathematical expression not reproducible].

After simplification, we get the truncation error as

[[tau].sub.i] = ([h.sub.i] - [h.sub.i-1]) (-[epsilon]/3u'"([[zeta].sub.i) + [a.sub.i]/2 u"([[eta].sub.i])).

As [h.sub.i] [right arrow] 0, the truncation error tends to zero, which shows that scheme is consistent. The order of the truncation error is given by O([h.sub.i] - [h.sub.i-1]).

4.2 Discretization error

Our aim is to express, if possible, the discretization error [e.sub.i] = [u.sub.i] - [U.sub.i] at the ith mesh point in terms of hi. The difference equations approximating the problem are given by

[mathematical expression not reproducible], (4.1)

[mathematical expression not reproducible]. (4.2)

Subtracting equation (4.1) from equation (4.2), we get

[mathematical expression not reproducible].

Multiplying both sides by ([h.sub.i] + [h.sub.i-1])[h.sub.i][h.sub.i-1], we get

[mathematical expression not reproducible], (4.3)

The error equation (4.3) can be represented as

[mathematical expression not reproducible], (4.4)

where

[mathematical expression not reproducible].

It can be easily verified that

[absolute value of [F.sub.i] - [E.sub.i] - [G.sub.i] - [H.sub.i]] = [absolute value of [b.sub.i]]([h.sub.i] + [h.sub.i- 1])[h.sub.i][h.sub.i-1] > 0.

Now from equation (4.4)

[mathematical expression not reproducible], (4.5)

Where

[g.sub.i] = - [epsilon]/u'"([[zeta].sub.i) + [a.sub.i]/2 u"([[eta].sub.i]).

Let [lambda] = [parallel]e[parallel][sub.[infinity]] and select the index i for which [absolute value of [e.sub.i]] = [parallel]e[parallel][sub.[infinity]] = [lambda]. Now from equation (4.5), we have

[mathematical expression not reproducible],

where

[mathematical expression not reproducible],

which is bounded and independent of [h.sub.i]. Thus [parallel]e[parallel][sub.[infinity]] = O([h.sub.i] - [h.sub.i-1]) as [h.sub.i] [right arrow] 0. This shows the convergence of the proposed numerical scheme.

5 Numerical examples and results

In this section three illustrative examples will be studied. To demonstrate the applicability and efficiency of the proposed method, we have implemented the present method on three examples with right-end boundary layer in the interval [0, 2]. These problems were widely discussed in the literature.

Since the exact solutions of the problems are not known, the maximum absolute errors for the examples are calculated using the following double mesh principle

[E.sup.N.sub.[epsilon]] = max [absolute value of [U.sup.N.sub.i] - [U.sup.2N.sub.2i]].

To use the double mesh principle, we incorporated the mesh points at the middle of each mesh and applied the central finite difference scheme on variable mesh. The numerical rate of convergence is estimated by the formula

[R.sup.N] = log [absolute value of [E.sup.N.sub.[epsilon]]/ [E.sup.2N.sub.[epsilon]]]/log 2.

Example 1. We consider singularly perturbed delay differential equation of convection-diffusion type as [8]

-[epsilon]u"(x) + 3u'(x) - u(x - 1)=0, u(x) = 1, -1 [less than or equal to] x [less than or equal to] 0, u(2) = 2.

The exact solution of this problem is given by

[mathematical expression not reproducible],

where

[mathematical expression not reproducible].

The maximum absolute errors, number of adaptive mesh generated points and the rate of convergence are presented in Table 1 for different values of perturbation parameter [epsilon].

Example 2. We consider singularly perturbed delay differential equation of convection-diffusion type as [8]

[mathematical expression not reproducible].

The maximum absolute errors, number of adaptive mesh generated points and the rate of convergence are presented in Table 2 for different values of perturbation parameter [epsilon].

Example 3. We consider singularly perturbed delay differential equation of convection-diffusion type [8]

[mathematical expression not reproducible].

In the previous two examples, it is assumed that f (x) is continuous on [0, 2]. Motivated by the works of [14,15], in this example, we suppose that f (x) has a simple discontinuity at x = 1, that is f(1-) [not equal to] f (1+), where f (1-) and f (1+) are left and right limits respectively. This boundary value problem exhibits a strong boundary layer at x = 2 and an interior week layer at x = 1.

The maximum absolute errors, number of adaptive mesh generated points and the rate of convergence are presented in Table 3 for different values of perturbation parameter e. It is observed from the Table 3 that after [epsilon] = [2.sup.-9] the rate of convergence has been reduced to almost 1, which is due to the effect of weak interior layer at x = 1. From the Figure 3, it is clear that some mesh points are added around x = 1 to resolve the weak boundary layer and some mesh points are added around x = 2 to resolve the strong boundary layer.

6 Conclusions

In this paper, we have presented an adaptive mesh method to deal with oscillation produced by central finite difference method, when applied to convection dominated singularly perturbed delay differential equations. The method is analysed for convergence. The proposed method is an [epsilon]-uniform convergent. An extensive amount of computational work has been carried out to demonstrate the proposed method. The maximum absolute error is tabulated in the form of Tables 1-3 for the considered examples in support of the predicted theory. The graphs of the solution of the considered examples for certain values of perturbation parameter are plotted in Figures 1-3. On the basis of the numerical results of a variety of examples, it is concluded that the present method offers significant advantage for the linear singularly perturbed delay differential equations.

https://doi.org/10.3846/mma.2018.041

Acknowledgements

The authors are grateful to the referees for their valuable suggestions and comments.

References

[1] G.M. Amiraliyev and E. Cimen. Numerical method for a singularly perturbed convection-diffusion problem with delay. Appl. Math. Comp., 216:2351-2359, 2010. https://doi.org/10.1016/j.amc.2010.03.080.

[2] G.M. Amiraliyev and F. Erdogan. Uniform numerical method for singularly perturbed delay differential equations. Comput. Math. Applicat., 53:1251-1259., 2007. https://doi.org/10.1016/jxamwa.2006.07.009.

[3] I.G. Amiraliyeva, F. Erdogan and G.M. Amiraliyev. A uniform numerical method for dealing with a singularly perturbed delay initial value problem. Appl. Math. L., 23:1221-1225, 2010. https://doi.org/10.1016/j.aml.2010.06.002.

[4] N.S. Bakhvalov. On the optimization of the methods for solving boundary value problems in the presence of boundary layers. Zh. Vychisl. Mater. Fiz., 9:841-859, 1969.

[5] K. Bansal, P. Rai and K.K. Sharma. Numerical treatment for the class of time dependent singularly perturbed parabolic problems with general shift arguments. Diff. Equ. Dyn. Sys., 25(2):327-346, 2017. https://doi.org/10.1007/s12591-015-0265-7.

[6] K. Bansal and K.K. Sharma. Parameter uniform numerical scheme for time dependent singularly perturbed convection-diffusion-reaction problems with general shift arguments. Num. Algor., 75(1):113-145, 2017.

[7] A. Bellen and M. Zennaro. Numerical methods for delay differential Equations. Oxford Science Publications, New York, 2003. https://doi.org/10.1093/acprof:oso/9780198506546.001.0001.

[8] P.P. Chakravarthy, S.D. Kumar and R.N. Rao. Numerical solution of second order singularly perturbed delay differential equations via cubic spline in tension. Int. J. Appl. Comput. Math, 3(3):1703-1717, 2016. https://doi.org/10.1007/s40819-016-0204-5.

[9] Y. Chen and J. Wu. The asymptotic shapes of periodic solutions of a singular delay differential system. J. Diff. Equ., 169:614-632, 2001. https://doi.org/10.1006/jdeq.2000.3910.

[10] R.V. Culshaw and S. Ruan. A delay differential equation model of hiv infection of CD4+ T-cells. Math. Biosci, 165(1):27-39, 2000. https://doi.org/10.1016/S0025-5564(00)00006-7.

[11] M.W. Derstine, H.M. Gibbs, F.A. Hopf and D.L. Kaplan. Bifurcation gap in a hybrid optical system. Phys. Rev. A, 26:3720-3722, 1982. https://doi.org/10.1103/PhysRevA.26.3720.

[12] R.D. Driver. Ordinary and delay differnential equations. Springer, New York, 1977. https://doi.org/10.1007/978-1-4684-9467-9.

[13] F. Erdogan and G.M. Amiraliyev. Fitted finite difference method for singularly perturbed delay differential equations. Numer. Algor., 59:131-145, 2012. https://doi.org/10.1007/s11075-011-9480-7.

[14] P.A. Farrell, A.F. Hegarty, J.J.H. Miller, E. O'Riordan and G.I. Shishkin. Singularly perturbed convection-diffusion problems with boundary and weak interior layers. J. Comput. Appl. Math., 166:133-151, 2004. https://doi.org/10.1016/j.cam.2003.09.033.

[15] P.A. Farrell, J.J.H. Miller, E. O'Riordan and G.I. Shishkin. Singularly perturbed differential equations with discontinuous source terms. In J.J.H. Miller G.I. Shishkin and L. Vulkov(Eds.), Analytical and Numerical Methods for Convection-Dominated and Singularly Perturbed Problems, pp. 23-32, New York, USA, 2000. Nova Science Publishers. Inc.

[16] E.C. Gartland. Graded-mesh difference schemes for singularly perturbed two-point boundary value problems. Math. comput., 51(184):631-657, 1988. https://doi.org/10.1090/S0025-5718-1988-0935072-1.

[17] V.Y. Glizer. Asymptotic solution of a boundary-value problem for linear singularly-perturbed functional differential equations arising in optimal control theory. J. Opt. Th. and Applicat, 106:49-85, 2000.

[18] V.Y. Glizer. Block wise estimate of the fundamental matrix of linear singularly perturbed differential systems with small delay and its application to uniform asymptotic solution. J. Math. Anal. Applicat., 278:409-433, 2003. https://doi.org/10.1016/S0022-247X(02)00715-1.

[19] M. Kot. Elements of Mathematical Ecology. Cambridge University Press, 2001. https://doi.org/10.1017/CBO9780511608520.

[20] Y. Kuang. Delay Differential Equations with Applications in Population Dynamics. Academic Press, New York, 1993.

[21] V. Kumar and B. Srinivasan. An adaptive mesh strategy for singularly perturbed convection diffusion problems. Appl. Math. Modell., 39:2081-2091, 2015. https://doi.org/10.1016/j.apm.2014.10.019.

[22] C.G. Lange and R.M. Miura. Singular perturbation analysis of boundary-value problems for differential-difference equations. SIAM J. Appl. Math., 42:502-531, 1982. https://doi.org/10.1137/0142036.

[23] R.J. Leveque. Finite volume methods for hyperbolic problems. Cambridge University Press, 2002. https://doi.org/10.1017/CBO9780511791253.

[24] T. Linb. Layer-adapted meshes for convection diffusion problems. Comput. Methods Appl. Mech. Eng., 192:1061-1105, 2003. https://doi.org/10.1016/S0045-7825(02)00630-8.

[25] V.D. Liseikin. The use of special transformations in the numerical solution of boundary value problems. Comput. Math. Math. Phys., 30:43-53, 1990. https://doi.org/10.1016/0041-5553(90)90006-E.

[26] A. Longtin and G.J. Milton. Complex oscillations in the human pupil light reflex with mixed and delayed feedback. Math. Biosci, 90:183-199, 1988. https://doi.org/10.1016/0025-5564(88)90064-8.

[27] M.C. Mackey and L. Glass. Oscillation and chaos in physiological control systems. Science, 197:287-289, 1977. https://doi.org/10.1126/science.267326.

[28] J.M. Mahaffy, J. Belair and M.C. Mackey. Hematopoietic model with moving boundary condition and state dependent delay: application in erythropoiesis. J. Theoret. Biol., 190:135-146, 1998. https://doi.org/10.1006/jtbi.1997.0537.

[29] J. Mallet-Paret and R.D. Nussbaum. A differential-delay equations arising in optics and physiology. SIAM J. Math. Anal., 20:249-292, 1989. https://doi.org/10.1137/0520019.

[30] H. Mayer, K.S. Zaenker and U. an der Heiden. A basic mathematical model of immune response. Chaos, 5:155-161, 1995. https://doi.org/10.1063/1.166098.

[31] J.J.H. Miller, E.O. Riordan and I.G. Shishkin. Fitted numerical methods for singular perturbation problems. Word Scientific, Singapore, 1996. https://doi.org/10.1142/2933.

[32] J.D. Murray. Mathematical Biology. An Introduction, Springer, 2002.

[33] P.W. Nelson and A.S. Perelson. Mathematical analysis of delay differential equation models of HIV-1 infection. Math. Biosci., 179:73-94, 2002. https://doi.org/10.1016/S0025-5564(02)00099-8.

[34] V. Subburayan and N. Ramanujam. An initial value technique for singularly perturbed convection - diffusion problems with a negative shift. J. Opt. Th. Applicat., 158:234-250, 2013. https://doi.org/10.1007/s10957-012-0200-9.

[35] R. Vulanovic. On a numerical solution of a type of singularly perturbed boundary value problem by using a special discretization mesh. Zb. Rad.,Prir.-Mat. Fak., Univ. Novom Sadu,Ser. Mat., 13:187-201, 1983.

[36] R. Vulanovic. Mesh construction for discretization of singularly perturbed Boundary value problems. Ph.d thesis, University of Novi sad, 1986.

[37] W.J. Wilbur and J. Rinzel. An analysis of Stein's model for stochastic neuronal excitation. Biol. Cyber., 45:107-114, 1982. https://doi.org/10.1007/BF00335237.

P. Pramod Chakravarthy (a), Trun Gupta (a) and R. Nageshwar Rao (b)

(a) Department of Mathematics, Visvesvaraya National Institute of Technology Nagpur, 440010, India.

(b) School of Advanced Sciences, VIT University Vellore, Tamil Nadu 632014, India.

E-mail(corresp.): pramodpodila@yahoo.co.in

E-mail: trun936@gmail.com

E-mail: nraojragi@yahoo.co.in

Received January 10, 2018; revised September 18, 2018; accepted September 19, 2018

Caption: Figure 1. Numerical solution for Example 1 with [epsilon] = [2.sup.-8], N = 23 (initially)

Caption: Figure 2. Numerical solution for Example 2 with [epsilon] = [2.sup.-8], N = 23 (initially)

Caption: Figure 3. Numerical solution for Example 3 with [epsilon] = [2.sup.-8], N = 23 (initially)
```Table 1. Maximum absolute error for Example 1 for different values
of e with N = 23 (initially)

[epsilon]      Max. error      Generated              Rate of
mesh ([N.sub.g])   Convergence ([R.sup.N])

[2.sup.-5]       0.0089            29                  2.1924
[2.sup.-6]       0.0089            31                  2.1953
[2.sup.-7]       0.0089            33                  2.2041
[2.sup.-8]       0.0090            35                  2.2195
[2.sup.-9]       0.0091            37                  2.2377
[2.sup.-10]      0.0092            39                  2.2496
[2.sup.-11]      0.0092            41                  2.2516
[2.sup.-12]      0.0092            43                  2.2486
[2.sup.-13]      0.0092            47                  2.2454
[2.sup.-14]      0.0092            49                  2.2435
[2.sup.-15]      0.0092            51                  2.2425
[2.sup.-16]      0.0092            53                  2.2421
[2.sup.-17]      0.0092            55                  2.2419
[2.sup.-18]      0.0092            57                  2.2419
[2.sup.-19]      0.0092            59                  2.2418
[2.sup.-20]      0.0092            63                  2.2418

Table 2. Maximum absolute error for Example 2 for different values
of e with N = 23 (initially)

[epsilon]      Max. error    Generated       Rate of
mesh        Convergence
([N.sub.g])     ([R.sup.N)]

[2.sup.-5]       0.0257          39           2.2519
[2.sup.-6]       0.0256          41           2.2508
[2.sup.-7]       0.0254          43           2.2431
[2.sup.-8]       0.0252          45           2.2358
[2.sup.-9]       0.0250          47           2.2453
[2.sup.-10]      0.0250          47           2.2695
[2.sup.-11]      0.0250          49           2.2917
[2.sup.-12]      0.0250          49           2.3044
[2.sup.-13]      0.0250          51           2.3100
[2.sup.-14]      0.0250          53           2.3121
[2.sup.-15]      0.0250          53           2.3129
[2.sup.-16]      0.0250          55           2.3132
[2.sup.-17]      0.0250          57           2.3133
[2.sup.-18]      0.0250          59           2.3134
[2.sup.-19]      0.0250          61           2.3134
[2.sup.-20]      0.0250          63           2.3134

Table 3. Maximum absolute error for Example 3 for different values
of e with N = 23 (initially)

[epsilon]     Max. error      Generated              Rate of
mesh ([N.sub.g])   Convergence ([R.sup.N])

[2.sup.-5]      0.0720            43                  2.1185
[2.sup.-6]      0.0719            47                  2.1180
[2.sup.-7]      0.0718            49                  2.1180
[2.sup.-8]      0.0718            53                  2.1177
[2.sup.-9]      0.0717            53                  2.1177
[2.sup.-10]     0.0414            39                  1.1016
[2.sup.-11]     0.0392            41                  0.9597
[2.sup.-12]     0.0379            43                  0.9034
[2.sup.-13]     0.0376            45                  0.8961
[2.sup.-14]     0.0380            47                  0.9167
[2.sup.-15]     0.0382            49                  0.9281
[2.sup.-16]     0.0383            51                  0.9341
[2.sup.-17]     0.0384            53                  0.9372
[2.sup.-18]     0.0384            55                  0.9387
[2.sup.-19]     0.0384            57                  0.9395
[2.sup.-20]     0.0384            59                  0.9399
```