# Parameter-uniform fitted mesh method for singularly perturbed delay differential equations with layer behavior.

Abstract. Boundary value problems for singularly perturbed differential difference equations containing delay with layer behavior are considered. There are a number of realistic models in the literature where one encounters BVPs for singularly perturbed differential difference equations with small delay, such as in variational problems in control theory and first exit time problems in modeling of activation of neurons. In some recent papers, the terms negative shift for 'delay' and positive shift for 'advance' are used. In this paper, a numerical method based on the fitted mesh approach to approximate the solution of these types of boundary value problems is presented. In this method the piecewise-uniform meshes are constructed and fitted to the boundary layer regions to adapt singular behavior of the operator in these narrow regions. Both the cases, layer on the left side boundary and layer on the right side boundary, are discussed. It is shown that the method composed of an upwind difference operator on the piecewise uniform mesh is parameter-uniform by establishing a robust error estimate. The effect of small delay on the boundary layer solution is shown by plotting the graphs of the solution for different delay values for several numerical examples. Numerical results in terms of maximum absolute error are tabulated to demonstrate the efficiency of the method.Key words. Fitted mesh, finite difference, singular perturbation, differential-difference equation, delay, boundary layer, action potential

AMS subject classifications. 34K28, 34K26, 34K10

1. Introduction and examples. Boundary value problems for singularly perturbed differential difference equations are ubiquitous in the mathematical modeling of various practical phenomena in biology and physics, such as in variational problems in control theory and first exit time problems in the modeling of the determination of expected time for the generation of action potentials in nerve cells by random synaptic inputs in dendrites.

In 1965, Stein [12] gave a realistic model for the stochastic activity of neurons. Stein's model contains assumptions concerning random synaptic inputs that excitatory impulses occur randomly and arrive according to the Poisson process [pi](fe, t), each event of which leads to an instantaneous increase in depolarization V(t) by [a.sub.e], whereas inhibitory impulses arrive at event times in the second Poisson process [pi](fi, t), which is independent of [pi](fe, t) and causes V(t) to decrease by [a.sub.e]. The neuron fires an impulse when V(t) reaches or exceeds a threshold value r. After each neuronal firing, there is a refractory period of duration [t.sub.0], during which the impulses have no effect and the membrane depolarization V (t) is reset to zero, etc. In 1967, Stein [13] generalized his model to handle a distribution of postsynaptic potential amplitudes. Various theoretical models of nerve membrane potential in the presence of random synaptic activations are given and many are discussed in Holden [4]. One avenue of attack has been through computer simulation, e.g., Segundo et al. (1968) [11], which provides a useful first step in the light of the analytic difficulties encountered in solving any realistic model.

The determination of the expected time for the generation of action potentials in nerve cells by random synaptic inputs in the dendrites can be modeled as a first-exit time problem. In Stein's model, the distribution representing inputs is taken as a Poisson process with exponential decay. If in addition, there are inputs that can be modeled as a Wiener process with variance parameter [sigma] and drift parameter [mu], then the problem for expected first-exit time y, given initial membrane potential x [member of] ([x.sub.1], [x.sub.2]), can be formulated as a general boundary value problem for the linear second order differential difference equation [8]

[[sigma].sup.2]/2 y"(x) + ([mu] - x)y'(x) + [[lambda].sub.E]y(x + [a.sub.E]) + [[lambda].sub.I]y (x - [a.sub.I]) - ([[lambda].sub.E] + [[lambda].sub.I])y(x) = -1,

where the values x = [x.sub.1] and x = [x.sub.2] correspond to the inhibitory reversal potential and to the threshold value of membrane potential for action potential generation, respectively. [sigma] and [mu] are the variance and drift parameters, respectively, y is the expected first-exit time and the first order derivative term -xy'(x) corresponds to exponential decay between synaptic inputs. The undifferentiated terms correspond to excitatory and inhibitory synaptic inputs, modeled as a Poisson process with mean rates [[lambda].sub.E] and [[lambda].sub.I], respectively, and produce jumps in the membrane potential of amounts [a.sub.E] and [a.sub.I], respectively, which are small quantities and could depend on the voltage. The boundary condition is

y(x) [equivalent to] 0, x [not member of] ([x.sub.1], [x.sub.2]).

This biological problem motivates the study of boundary value problems for singularly perturbed differential difference equations with delay as well as advance, which was initiated by Lange and Miura [8, 9], where they introduced the new terminology "negative shift" for "delay" and "positive shift" for "advance".

In this paper, we consider the boundary value problems for a certain simpler class of singularly perturbed differential difference equations where there is only a small delay in the convection term and there is no term containing advance. The objective of the present paper is to continue the numerical study for solving boundary value problems for a class of singularly perturbed differential difference equations which contain small delays (negative shift) in the convection term, provided the delays are of small order of the singular perturbation parameter (0 < [member of] [much less than] 1), and to provide a more sophisticated and robust numerical scheme to approximate the solution of such type of boundary value problems which works for a wide class of problems where the existing numerical schemes fails. Since the delay is of the small order of the singular perturbation parameter, to tackle the term containing delay, we use Taylor's approximation as pointed out by Cunningham [1, p. 222], Lange et al. [9, p. 275] and Tian in his thesis work [14]. Then we construct three numerical schemes based on finite difference methods, namely, i) the standard upwind finite difference scheme which is discussed in detail in [6], ii) the fitted operator finite difference scheme which is discussed in detail in [5], and iii) the fitted mesh finite difference scheme which is considered in this paper.

The standard upwind finite difference scheme is not robust with respect to the singular perturbation parameter, while the fitted operator and fitted mesh finite difference schemes are robust with respect to the singular perturbation parameter. Now it is justified to switch over from a standard upwind finite difference scheme to the schemes adopting fitted operator or fitted mesh approach to achieve parameter uniform convergence. It remains to justify the need of construction of the fitted mesh method over the fitted operator method or, more briefly, the superiority of the fitted mesh method to the fitted operator method. There are some problems for which no parameter uniform numerical scheme can be constructed on a uniform mesh using the fitted operator approach, while for the same problems, a parameter uniform numerical scheme is constructed based on the fitted mesh approach [10, 3]. Moreover, the numerical schemes constructed using the fitted mesh approach are usually simpler than the numerical schemes based on fitted operators for singularly perturbed differential equations. These schemes are also easier to generalize to problems in more than one dimension and to nonlinear problems.

When the delay is zero, the above boundary value problem reduces to a singularly perturbed ordinary differential equation of the convection-diffusion type and the solution of the differential equation so obtained exhibits layer behavior. The effect of small delay on the layer behavior of the solution is analyzed by plotting the graphs of the solutions of the considered examples for different values of the delay.

2. Description of the problem. In this section, we state the boundary value problems for a class of singularly perturbed differential difference equations of the convection-diffusion type with small delay

[epsilon]y''(x) + a(x)y'(x - [delta]) + b(x)y(x) = f(x), (2.1)

on [OMEGA] = (0, 1), under the interval and boundary conditions

y(x) = [phi](x) on - [delta] [less than or equal to] x [less than or equal to] 0, y(1) = [gamma], (2.2)

where [epsilon] is a small parameter, 0 < [epsilon] [much less than] 1 and [delta] is of o([epsilon]) satisfying ([epsilon] - [delta]a(x)) > 0 for all x [member of] [0, 1]; a(x), b(x), f(x) and [phi](x) are smooth functions and [gamma] is a constant. For a function y(x) to be a smooth solution to the problem (2.1), (2.2), it must satisfy (2.1), (2.2), be continuous on [bar.[OMEGA]] = [0, 1] and be continuously differentiable on [OMEGA] = (0, 1). It is also assumed that b(x) satisfies the condition

b(x) [less than or equal to] -[theta] < 0 [for all]x [member of] [bar.[OMEGA]],

where [theta] is a positive constant. The boundary value problems for the above class of singularly perturbed differential difference equations contain delay only in the first-order derivative term. For [delta] = 0, the problem (2.1), (2.2) is converted into a boundary value problem for the singularly perturbed ordinary differential equation. The reduced problem corresponding to the singularly perturbed differential equation obtained by setting [epsilon] = 0 in the problem (2.1), (2.2) for [delta] = 0 is the problem a(x)y'(x) + b(x)y(x) = f(x). Due to the loss of order of the differential equation by one, the solution of the reduced problem cannot be made to satisfy both arbitrary preassigned boundary conditions simultaneously at the boundary points {0, 1} of the domain of consideration [??]. Thus in general, the solution y(x) exhibits boundary layer behavior at one of the end points of the interval [0, 1] depending on the sign of a(x), i.e., the boundary layer will be on the left side or on the right side of the interval [0, 1] according as a(x) > 0 or a(x) < 0 throughout the interval [0, 1], respectively. The layer is maintained for [delta] [not equal to] 0 but sufficiently small. In this paper, we consider both the cases where either the boundary layer is on the left side or on the right side of the interval [0, 1]. If a(x) changes sign throughout the domain of consideration then the solution can exhibit more complicated turning point behavior. The problems whose solutions exhibit turning point behavior are not discussed here.

Throughout the paper, C, M and [theta] denote generic positive constants that are independent of [epsilon] and in the case of discrete problems, also independent of the mesh parameter N. Where the value of C may change from result to result but remains constant in each. [parallel]x[parallel] denotes the global maximum norm over the appropriate domain of the independent variable, i.e.,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

3. Construction of numerical scheme. Upon using the Taylor approximation to the term containing the delay, the boundary value problem (2.1), (2.2) reduces to

([epsilon] - [delta]a(x))u"(x) + a(x)u'(x) + b(x)u(x) = f(x), (3.1),

u(0) = [[phi].sub.0], [[phi].sub.0] = [phi](0), u(1) = [gamma]. (3.2)

Since y [member of] [C.sup.2][0, 1] and the delay argument is sufficiently small, the solution u(x) of the problem (3.1), (3.2) provides a good approximation to the solution y(x) of the problem (2.1), (2.2). We denote by [L.sub.[epsilon]] the differential operator for the above problem (3.1), (3.2) which is defined for any function [PSI] [member of] [C.sup.2]([bar.[OMEGA]]) as

[L.sub.[epsilon]][PSI](x) = ([epsilon] - [delta]a(x))[PSI]"(x) + a(x)[PSI]'(x) + b(x)[PSI](x).

3.1. Layer on the left side. Here, we assume that a(x) [greater than or equal to] M > 0 and ([epsilon] - [delta]a(x)) > 0 throughout the interval [0, 1], where M is a positive constant. In this case, the solution of the problem exhibits layer behavior on the left side of the interval [0, 1].

3.1.1. Analytical results.

LEMMA 3.1. (Minimum Principle [2]) Suppose [PSI] is as mooth function satisfying [PSI](0) [greater than or equal to] 0, [PSI](1) [greater than or equal to] 0. Then [L.sub.[epsilon]][PSI](x) [less than or equal to] 0 for all x [member of] [OMEGA] implies [PSI](x) [greater than of equal to] 0 for all x [member of] [bar.[OMEGA]].

Proof. Let [x.sup.*] [member of] [bar.[OMEGA]] be such that [PSI]([x.sup.*]) = [min.sub.x[member of][bar.[OMEGA]]][PSI](x) and assume that [PSI]([x.sup.*]) < 0. Clearly [x.sup.*] [not member of] {0, 1}, therefore [PSI]'([x.sup.*]) = 0 and [PSI]''([x.sup.*]) [greater than or equal to] 0.

Now consider

[L.sub.[epsilon][PSI]([x.sup.*]) = ([epsilon] - [gamma]a([x.sup.*]))[PSI]''([x.sup.*]) + a([x.sup.*])[PSI]'([x.sup.*]) + b([x.sup.*])[PSI]([x.sup.*]) > 0,

which is a contradiction. It follows that [PSI]([x.sup.*]) [greater than or equal to] 0 and thus [PSI](x) [greater than or equal to] 0 [for all]x [member of] [bar.[OMEGA]].

LEMMA 3.2. Let u(x) be the solution of the problem (3.1), (3.2), then we have

[parallel]u[parallel] [less than or equal to] [[theta].sup.1] [parallel]f[parallel] + max([absolute value of [[phi].sub.0]|, [absolute value of [gamma]]).

Proof. Let us construct the two barrier functions [[PSI].sup.[+ or -]] defined by

[[PSI].sup.[+ or -]](x) = [[theta].sup.-1] [parallel]f[parallel] + max([absolute value of [[phi].sub.0]|, [absolute value of [gamma]]) [+ or -] u(x).

Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

since b(x) [less than or equal to] -[theta] < 0, so we have b(x)[[theta].sup.-1] [greater than or equal to] -1. Using this inequality in the above inequality, we get

[L.sub.[epsilon]][[PSI].sup.[+ or -]](x) [less than or equal to] (-[parallel]f[parallel] [+ or -] f(x)) max([absolute value of [[phi].sub.0]], [absolute value of [gamma]]) [less than or equal to] 0 [for all]x [member of] [OMEGA], since [parallel]f[parallel] [greater than or equal to] f(x).

Therefore by the minimum principle, we obtain [[PSI].sup.[+ or -]](x) [greater than or equal to] 0 for all x [member of] [bar.[OMEGA]], which gives the required estimate.

THEOREM 3.3. Let u be the solution of the problem (3.1), (3.2). Then for k = 1, 2, 3

[parallel][u.sup.(k)][parallel] [greater than or equal to] C[([epsilon] - [delta]M).sup.-k].

Proof. Let x [member of] [OMEGA] and construct a neighborhood [N.sub.x] = (c, c + ([epsilon] - [delta][parallel]a[parallel])), where c is a positive constant chosen so that x [member of] [N.sub.x] and [N.sub.x] [subset] [OMEGA]. Then by the Mean Value Theorem, there exists a point z [member of] [N.sub.x] such that

u'(z) = u(c + ([epsilon] - [delta]M)) - u(c)/([epsilon] - [delta]M)

and so

[absolute value of ([epsilon] - [delta]M)u'(z)] [less than or equal to] 2[parallel]u[parallel. (3.3)

Integrating the differential equation (3.1) from z to x we get

([epsilon] - [delta]M)u'(x) - ([epsilon] - [delta]M)u'(z) = [[integral].sup.x.sub.z] (f - a(t)u'(t) - b(t)u(t))dt,

taking modulus on both sides and using the fact that the maximum norm of a function is never smaller than the value of the function over the domain of consideration, we get

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.4)

We have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Once again taking modulus on both sides and using the fact that the maximum norm of a function is always greater than the value of the function over the domain of consideration, we get

[[integral].sup.x.sub.z [absolute value of a(t)u'(t)]dt [less than or equal to] (2[parallel]a[parallel] + [parallel]a'[parallel])[parallel]u[parallel] (3.5)

Using inequalities (3.3), (3.5), 0 [less than or equal to] [absolute value of x - z] [less than or equal to] 1 and Lemma 3.2 for the bound on u in inequality (3.4), we get

[absolute value of u'(x)] [less than or equal to] C[([epsilon] - [delta]M).sup.-1],

which gives [parallel]u'[parallel] [less than or equal to] C[([epsilon] - [delta]M).sup.-1], where C = [parallel]f[parallel] + (2 + 2)[parallel]a[parallel] + [parallel]a'[parallel] + [parallel]b[parallel])([[theta]. sup.-1][parallel]f[parallel] + max([absolute value of [[theta].sub.0], [absolute value of [gamma]])) is independent of [epsilon]. Similarly one can easily find out the bounds on u" and u"' by using the differential equation (3.1) and the bounds on u and u'.

These bounds for derivatives of u were first obtained by O'Riordan et al. [10], using techniques based on Kellogg et al. [7]. However in order to prove that the numerical method is [epsilon]-uniform, one needs more precise information about the behavior of the exact solution of the problem (3.1), (3.2). This is obtained by decomposing the solution y into a smooth component v and a singular component w as follows

u = v + w

where v can be written in the form v(x) = [v.sub.0](x) + ([epsilon] - [delta]M)[v.sub.1](x) + [([epsilon] - [delta]M).sup.2][v.sub.2](x) and [v.sub.0](x), [v.sub.1](x) and [v.sub.2](x) are defined respectively to be the solutions of the problems

a(x)[v'.sub.0](x) + b(x)[v.sub.0](x) = f(x), x [member of] [OMEGA], [v.sub.0](1) = u(1) (3.6)

a(x)[v'.sub.1](x) + b(x)[v.sub.l](x) = -([epsilon] - [delta]a(x))[v".sub.0](x)/([epsilon] - [delta]M), x [member of] [OMEGA], [v.sub.1](l) = 0 (3.7)

[L.sub.[epsilon]][v.sub.2](x) = -([epsilon] - [delta]a(x))[v".sub.1](x)/([epsilon] - [delta]M), x [member of] [OMEGA], [v.sub.2](0) = 0, [v.sub.2](1) = 0. (3.8)

The smooth component v(x) is the solution of

[L.sub.[epsilon]]v(x) = f(x), x [member of] [OMEGA], v(0) = [v.sub.0](0) + ([epsilon] - [delta]M)[v.sub.l](0), v(1) = u(1)

and consequently the singular component w(x) is the solution of the homogeneous problem

[L.sub.[epsilon]]w(x) = 0, x [member of] [OMEGA], w(0) = u(0) - v(0), w(1) = 0. (3.9)

LEMMA 3.4. The bounds on [v.sub.0] and its derivatives for 0 [less than or equal to] k [less than or equal to] 3 satisfy

[parallel][v.sup.(k).sub.0][parallel] [less than or equal to] C.

Proof. The problem (3.6) can be written as

[v'.sub.0](t) + p(t)[v.sub.0](t) = f(t)/a(t), [v.sub.0](1) = [gamma], (3.10)

where p(t) = b(t)/a(t).

The problem (3.10) is a first order linear differential equation in [v.sub.0]. Therefore to obtain the solution of the problem, we multiply (3.10) by exp([integral] p(t)dt), and simplification yields

(exp ([integral] p(t)dt) [v.sub.0])' = f(t) exp ([integral] p(t)dt) /a(t). (3.11)

Now integrating (3.11) from x to 1, for some x [member of] (0, 1), we obtain

[([v.sub.0](t) exp ([integral] p(t)dt)).sup.l.sub.t=x] = [integral.sup.1.sub.x] [f(t) exp ([integral] p(t)dt) /a(t)]dt,

which on simplification reduces to

[v.sub.0](x) =-[gamma]A/B(x) - [[integral].sup.1.sub.x] [f(t) exp ([integral] p(t)dt) /a(t)]dt/B(x), (3.12)

where A = exp [([integral] p(t)dt)|.sub.t=1] and B(x) = exp [(integral p(t)dt)|.sub.t=x]. Now since a(t) and b(t) are sufficiently smooth, i.e., bounded for all t [member of] [0, 1] and we have a(t) [greater than or equal to] M > 0 [for all]t [member of] [0, 1], i.e., a(t) [not equal to] 0 for all t [member of] [0, 1], therefore p(t) = b(t)/ a(t) < 0 is bounded for all t [member of] [0, 1] and also f(t) is bounded for all t [member of] [0, 1]. Thus combining all the facts, we conclude that A, B and the second term on the right side of equation (3.12) are bounded which implies the required result that [v.sub.0] is bounded. Again from equation (3.6), we have

[v'.sub.0](x) = f(x)/a(x) - p(x)[v.sub.0](x),

and the boundedness of [v.sub.0] implies that [v'>sub.0] is bounded. Using the boundedness of [v.sub.0] and [v'.sub.0] and differentiating the differential equation (3.6) successively, we obtain the bounds on [v".sub.0] and [v"'.sub.0]'. Thus for 0 [less than or equal to] k [less than or equal to] 3, we have

[parallel][v.sup.(k).sub.0][parallel] [less than or equal to] C.

THEOREM 3.5. Let u be the solution of the problem (3.1), (3.2) and let u = v + w. For 0 [less than or equal to] k [less than or equal to] 3 and for sufficiently small [epsilon]; v, w and their derivatives satisfy the following bounds

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.13)

Proof. Since [v.sub.1] is the solution of the first order linear differential equation (3.7) and [v".sub.0](x), ([epsilon] - [delta]a(x))/([epsilon] - [delta]M), a(x) and b(x) are bounded on [0, 1] therefore the right side of equation (3.7) is bounded. Thus using similar steps as we have used in the proof of Lemma 3.4, we obtain

[parallel][v.sub.l][parallel] [less than or equal to] C,

where C is a constant. Now using this bound on [v.sub.1] and the differential equation (3.7), we get [parallel][v'.sub.1][parallel] [less than or equal to] C. After differentiating the (3.7) successively and using the bound on [v.sub.1] and [v'.sub.1], one can easily obtain the bounds on [v".sub.1]" and [v"'.sub.1]".

The quantities [v".sub.1] and ([epsilon] - [delta]a(x))/([epsilon] - [delta]M) are bounded by a constant independent of [epsilon], so the right side of equation (3.8) is bounded independently of [epsilon]. Thus [v.sub.2] is the solution of a boundary value problem similar to the problem (3.1), (3.2). Hence by Theorem 3.3, we have for 0 [less than or equal to] k [less than or equal to] 3

[parallel][v.sup.(k).sub.2][parallel] [less than of equal to] C[([epsilon] - [delta]M).sup.-k],

which gives the required estimate for the regular component v. Now to obtain the required bounds on the singular component w and its derivatives we consider the two barrier functions [[PSI].sup.[+ or -] defined by

[[PSI].sup.[+ or -](x) = ([absolute value of u(0)] + [absolute value of [v.sub.0]]) exp(-xM/([epsilon] - xM)) [+ or -] w(x).

Then we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Therefore by the minimum principle, we have

[[PSI].[+ or -](x) = ([absolute value of u(0)] + [absolute value of [v.sub.0]]) exp(-xM/([epsilon] - [delta]M)) [+ or -] w(x) [greater than or equal to] 0, x [member of] [OMEGA],

which on simplification gives

[absolute value of w(x)] [less than or equal to] C exp(-xM[([epsilon] - [delta]M).sup.-1]), x [member of] [bar.[OMEGA]],

where C = ([absolute value of u(0)] + [absolute value of [v.sub.0](0)]).

Now to find out the bounds on the derivatives of the singular component w of the solution y, we use the same technique as we did in the proof of Theorem 3.3. For x [member of] [OMEGA], construct a neighborhood [N.sub.x] = (x, x + ([epsilon] - [delta]M)). Therefore by the Mean Value Theorem, there exists a point z [epsilon] [N.sub.x] such that

w'(z) = w(x + ([epsilon] - {delta]M)) - w(x)/([epsilon] -[delta]M),

which implies that

[absolute value of ([epsilon] - [delta]M)w'(z)] [less than or equal to] 2[[parallel]w[parallel]. (3.14)

Now we have

[[integral].sup.x.sub.z] w"(t)dt = w'(x) - w'(z),

i.e.,

w'(x) = w'(z) + [[integral].sup.x.sub.z] w"(t)dt.

Using equation (3.9) for w"(t) in the above equation, we obtain

w'(x) = w'(z) + [[integral].sup.x.sub.z] -[([epsilon] - [delta]a(t)).sup.-1][a(t)w'(t) + b(t)w(t)]dt.

Taking modulus on both sides, we obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.15)

We have, by integration by parts,

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Using the fact that the maximum norm of a function is always greater than the value of the function over the domain of consideration and the inequality 0 < [absolute value of x - z] [less than or equal to] 1 followed by a simplification yields

[absolute value of [[integral].sup.x.sub.z] a(t)w'(t)dt] [less than or equal to] (2[parallel]a[parallel] + [parallel]a'[parallel])[parallel]w[parallel]. (3.16)

Using inequalities (3.14) and (3.16) in the inequality (3.15), we obtain

[absolute value of w'(x)] [less than of equal to] C[([epsilon] - [delta]M).sup.-1][parallel]w[parallel]. (3.17)

For x [member of] [N.sub.x], we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Using this value of w in inequality (3.17), we obtain

[absolute value of w'(x)] [less than or equal to] C[([epsilon] - [delta]M).sup.-1] exp(-xM/([epsilon] - [delta]M)),

which gives the required result. The estimate for w" can be easily obtained from the differential equation and the bounds on w and w'.

3.1.2. Standard finite difference method. To discretize the boundary value problem (3.1), (3.2), we place a uniform mesh of size h = 1/N on the interval [0, 1]. Denote [x.sub.i] = ih, [u.sub.i] = u([x.sub.i]), [a.sub.i] = a([x.sub.i]), [b.sub.i] = b([x.sub.i]) and fi = f([x.sub.i]), i = 0, 1, ... N. In the problem (3.1), (3.2), we approximate u"(x) and u'(x) by central difference and forward difference, respectively.

[L.sup.N.sub.l, l] = ([epsilon] - [delta]a(xi))D+[D.sup.-][u.sub.i] + a([x.sub.i])[D.sup.+][u.sub.i] + b([x.sub.i])[u.sub.i] = f([x.sub.i]),

u(0) = [[phi].sub.0], u(1) = [gamma],

where [D.sup.+][D.sup.-1][u.sub.i] = ([u.sub.i-1] - 2[u.sub.i+1] + [u.sub.i+1])/[h.sup.2], [D.sup.+][u.sub.i] = ([u.sub.i+1] - [u.sub.i])h and [D.sup.-][u.sub.i] = ([u.sub.i] - [u.sub.i-1])/h, which on simplification gives a three point difference scheme

[L.sup.N.1, l] = [E.sub.i][u.sub.i-1] - [F.sub.i][u.sub.i] + [G.sub.i][u.sub.i+1] = [H.sub.i] (3.18)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The difference equations (3.18) form a tridiagonal system of N - 1 equations with N + 1 unknowns [u.sub.0], [u.sub.1], ... [u.sub.N]. The N - 1 equations together with the given two boundary conditions are sufficient to solve the system. The stability and convergence analysis of the scheme is discussed in [6].

3.1.3. Fitted mesh finite difference method. In this section, we discretize the boundary value problem (3.1), (3.2) using the fitted mesh finite difference method composed of a standard upwind finite difference operator on a fitted piecewise uniform mesh, condensing in the boundary layer region at x = 0. The fitted piecewise-uniform mesh [[bar.[OMEGA]].sup.N.sub.l] on the interval [0, 1] is constructed by partitioning the interval into two subintervals [0, [tau]], [[tau], 1], where the transition parameter is chosen such that [tau] = min[0.5, C(epsilon] - [delta]M) In N] with C = 1/[theta]. It is assumed that N = [2.sup.k], where k [greater than or equal to] 2 is an integer, which guarantees that there is at least one point in the boundary layer region. On each of these subintervals, a uniform mesh is placed. A fitted finite difference method for the problem (3.1), (3.2) on the piecewise uniform mesh [[bar.[OMEGA]].sup.N.sub.l] is defined by

[L.sup.N.sub.2, l][u.sub.i] = f([x.sub.i]), [x.sub.i] [member of] [[OMEGA].sup.N.sub.l], (3.19)

[u.sub.0 = [[theta].sub.0],

[u.sub.N] = [gamma], (3.20)

where the discrete operator [L.sup.N.sub.2, l] is defined as

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

with

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Discrete Minimum Principle. Assume that [[PSI].sub.0] [greater than of equal to] 0 and [[PSI].sub.N] [greater than or equal to] 0. The [L.sup.N.sub.2, l][[PSI].sub.i] [less than or equal to] 0, 1 [less than or equal to] i [less than or equal to] N - 1 implies that [[PSI].sub.i] [greater than or equal to] 0 for all i, 0 [less than or equal to] i less than or equal to] N.

Proof. Let k be such that [[PSI].sub.k] = [min.sub.0[less than or equal to]I[less than or equal to]N] [[PSI].sub.i] and suppose [[PSI].sub.k] < 0. Then [[PSI].sub.k] - [[PSI].sub.k-1] [less than or equal to] 0, [[PSI].sub.k+1] - [[PSI].sub.k] [greater than or equa to] 0 and we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which contradicts the hypothesis that [L.sup.N.sub.2, l][[PSI].sub.i] [less than or equal to] 0, 1 [less than or equal to] i [less than or equal to] N - 1. Therefore [[PSI].sub.k] [greater than or equal to] 0 and we have chosen k fixed but arbitrary, so [[PSI].sub.i] [greater than or equal to] 0 for all i, 0 [less than or equal to] i [less than or equal to] N.

LEMMA 3.6. Let [Z.sub.i] be any mesh function such that [Z.sub.0] = [Z.sub.N] = 0. Then for all i, 0 [less than or equal to] I [less than or equal to] N

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Proof. Let us introduce two mesh functions [[PSI].sup.[+ or -].sub.i] defined by

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

and for 1 [less than or equal to] I [less than or equal to] N - 1

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Therefore by the discrete minimum principle, we have [[PSI].sup.[+ or -].sub.i] [greater than or equal to] 0 [for all]i, 0 [less than or equal to] I [less than or equal to] N

which gives the required estimate.

THEOREM 3.7. Let u(x) be the solution of the continuous problem (3.1), (3.2) and [U.sup.N] = [<[u.sub.i]>.sup.N.sub.i=0] be the solution of the corresponding discrete problem (3.19), (3.20). The fitted mesh finite difference method with standard upwind finite difference operator on the piecewise uniform mesh [[bar.[OMEGA]].sup.N.sub.l], condensing at the boundary layer at x = 0, is [epsilon]-uniform. Moreover y and [U.sub.N] = [<[u.sub.i]>.sup.N.sub.i=0] satisfy the following error estimate

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

where C is a constant independent of [epsilon].

Proof. As in the case of the continuous problem, the solution [U.sup.N] = [<[u.sup.i]>.sup.N.sub.i] = 0] of the discrete problem (3.19), (3.20) can be decomposed into regular and singular components. Thus

[U.sup.N] = [V.sub.N] + [W.sub.N],

where [V.sup.N] is the solution of the inhomogeneous problem

[L.sup.N.sub.2. l][V.sup.N]([x.sub.i]) = f([x.sub.i]) for all [x.sub.i] [member of] [[OMEGA].sup.N.sub.1], [V.sub.N](0) = v(0), [V.sup.N](1) = v(1)

and [W.sup.N] is the solution of the homogeneous problem

[L.sup.N.sub.2, l][W.sup.N](xi) = 0 for all [x.sub.i] [member of] [[OMEGA].sup.N.sub.1], [W.sup.N](0) = w(0), [W.sup.N](1) = w(1).

Then the error can be written in the form

[U.sup.N] - u = ([V.sup.N] - v) + ([W.sup.N] - w).

Thus the errors in the regular and singular components of the solution can be estimated separately. To estimate the error for the regular component, from differential and difference equations, we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (3.21)

Let [x.sub.i] [member of] [[OMEGA].sup.N.sub.l]. Then for any [psi] [member of] [C.sup.2] ([[bar.[OMEGA]].sub.l]),

[absolute value of [D.sup.+] - d/dx)[psi]([x.sub.i])] [less than or equal to] ([x.sub.i+1] [x.sub.i])[parallel][[psi].sup.(2)][parallel]/2

and for any [psi] [member of] [C.sup.3] ([[bar.[OMEGA].sub.l]),

[absolute value of ([D.sup.[+ or -]] - [d.sup.2]/[dx.sup.2])[psi]([x.sub.i])] [less than or equal to] ([x.sub.i+1] - [x.sub.i])[parallel][[psi].sup.(3)][parallel]/3.

For a proof of these results, see [10, Lemma 1, p. 21]. Using these results in (3.21), we obtain

[absolute value of [L.sub.2]([V.sup.N] - v)([x.sub.i])] [less than or equal to] ([x.sub.i+1] - [x.sub.i-1]) (([epsilon] - [delta]a([x.sub.i]))/3 [parallel][v.sup.(3)[parallel] + a([x.sub.i])/2 [parallel][v.sup.(2)][parallel].

Since [x.sub.i+1] - [x.sub.i-1] [less than or equal to] 2[N.sup.-1] and using Theorem 3.5 for bounds on [v.sup.(2)] and [v.sup.(3)], we obtain

[absolute value of [L.sup.N.sub.2, l]([V.sup.N] - v)([x.sub.i])] [less than or equal to] C[N.sup.-1], [x.sub.i] [member of] [[OMEGA].sup.N.sub.l] (3.22)

Now an application of Lemma 3.6 for the mesh function ([V.sup.N] - v) ([x.sub.i]) gives

[absolute value of [L.sup.N.sub.2, l]([V.sup.N] - v]([x.xub.i])] [less than or equal to] C[N.sub.-1], [x.sub.i] [member of [[OMEGA].sup.N.sub.l] (3.23)

Using inequality (3.22) in inequality (3.23), we obtain

[absolute value of ([V.sup.N] - v)([x.sub.i])] [less than or equal to] C[N.sup.-1]. (3.24)

Arguments for the estimation of the singular component of the error [L.sup.N.sub.2, l]([W.sup.N] - w) depends on the value of the transition parameter [tau], whether [tau[ = 1/2 or [tau] = C([epsilon] - [delta]M) ln N, where C = 1/[theta].

Case i) C([epsilon] - [delta]M) ln N [greater than or equal to] 1/2, i.e., when the mesh is uniform

In this case we go through the same arguments as we did in the case of the estimation of the regular part of the error which gives

[absolute value of [L.sup.N.sub.2, l]([W.sup.N] - w)] [less than or equal to] C([x.sub.i+1] - [xi.sub.i]) (([epsilon] - [delta]a([x.sub.i]))[parallel][w.sup.(3)] + a([x.sub.i])[parallel][w.sup.(2)][parallel])

Using Theorem 3.5 for bounds on [w.sup.(2)] and [w.sup.(3)] and the fact that ([x.sub.i+1] - [x.sub.i-1]) = 2[N.sup.-1] for the uniform mesh, we obtain

[absolute value of [L.sup.N.sub.2, l]([W.sup.N] - w)([x.sub.i])] [greater than or equal to] C[([epsilon] - [delta]M).sup.-1] [N.sup.-1]. (3.25)

In this case, we have [([epsilon] - [delta][parallel]a[parallel]).sup.-1] [less than or equal to] 2C ln N. Using this inequality in the above inequality (3.25), we obtain

[absolute value of [L.sup.N.sub.2, l]([W.sup.N] - w)([x.sub.i])] [less than or equal to] C[N.sup.-1] [(ln N).sup.2]. (3.26)

Now an application of Lemma 3.6 for the mesh function ([W.sup.N] - w) ([x.sub.i]) gives

[absolute value of ([W.sup.N] - w)([x.sub.i])] [less than or equal to] [absolute value of [L.sup.N.sub.2, l] ([W.sup.N] - w)([x.sub.i])] [for all][x.sub.i] [member of] [[OMEGA].sup.N.sub.l] (3.27)

Using (3.26) in (3.27), we obtain

[absolute value of ([W.sup.N] - w)([x.sub.i])] C[N.sup.-1][(ln N).sup.2] [for all][x.sub.i] [member of] [[OMEGA].sup.N.sub.l].

Case ii) C([epsilon] - [delta]M) ln N < 1/2, i.e., when the mesh is piecewise uniform with mesh spacing 2[tau]/N in the subinterval [0, [tau]] and 2(1 - [tau])/N in the subinterval [[tau], 1]

We give separate proofs for the bounds on the singular component of the error in the coarse and fine mesh subintervals. First we obtain the bound on the singular component in the outer region, i.e., in the subinterval [[tau], 1]. Using the triangular inequality, we have

[absolute value of (W - w)([X.sub.i])] [less than or equal to] [absolute value of W([x.sub.i])] + w([x.sub.i])|.

From inequality (3.13), we have

[absolute value of w([x.sub.i])] [less than or equal to] C exp(-M[x.sub.i]/([epsilon] - [delta]M)) (3.28)

for all [x.sub.i] [member of] [[tau], 1]. Exp(-M[x.sub.i]/([epsilon] - [delta]M)) is a decreasing function and [x.sub.i] [greater than or equal to] [tau]. Using these facts in above inequality (3.28) we have

[absolute value of w([x.sub.i])] [less than or equal to] C exp(-M[tau]/([epsilon] - [delta]M))

In this case we have [[tau].sub.i] = C([epsilon] - [delta]M) ln N. Using this value of [tau] in the above inequality, we get

[absolute value of w([x.sub.i])] [less than or equal to] C[N.sup.-1]. (3.29)

Now to obtain the bound on [W.sup.N], we construct a mesh function [[??].sup.N] defined as the solution of the following problem

([epsilon] - [delta]a([x.sub.i]))[D.sup.[+ or -][[??].sup.N]([x.sub.i]) + M [D.sup.+][[??].sup.N]([x.sub.i]) + b([x.sub.i])[[??].sup.N]([x.sub.i]) = 0,

1 [less than or equal to] i [less than or equal to] N - 1, under the same boundary conditions as for W. Then an application of [10, Lemma 5, p. 53] yields

[absolute value of [W.sup.N]([x.sub.i])] [less than or equal to] [absolute value of [[??].sup.N]([x.sub.i]], 0 =[less than or equal to] i [less than or equal to] N. (3.30)

Applying [10, Lemma 3, p. 51 ] for [[??].sup.N.sub.[epsilon] gives

[absolute value of [[??].sup.N]([x.sub.i])] [less than or equal to] CN-1, N/2 [less than or equal to] i [less than or equal to N.

Using this estimate for [[??].sup.N.sub.[epsilon] in inequality (3.30), we obtain

[absolute value of [W.sup.N]([x.sub.i])] [less than or equal to] C[N.sup.-1], N/2 [less than or equal to] i [less than or equal to] N. (3.31)

Thus from inequalities (3.29) and (3.31), we obtain the bound on the singular component of error in the outer region [[tau], 1]

[absolute value of [W.sup.N] - w([x.sub.i])] [less than or equal to] C[N.sup.-1], N/2 [less than or equal to] i [less than or equal to] N. (3.32)

Now it remains to prove the result for [x.sub.i] [member of] [0, [tau]], i.e., in the boundary layer region. For i = 0, there is nothing to prove. For [x.sub.i] [member of] (0, [tau]) the proof follows on the same lines as for the case i) except that we use the discrete minimum principle on [0, [tau]] and the already established bounds W([X.sub.N/2]) [less than or equal to] C[N.sup.-1] and w([X.sub.N/2]) [less than or equal to] C[N.sup.-1]. Thus by using similar arguments as we have used in the estimation of the regular component of the error, we get

[absolute value of [L.sup.N.sub.2, l](W.sup.N] - w)] [less than or equal to] 2[tau][N.sup.-1] [([epsilon] - [delta][parallel]a[parallel]).sup.-2], [absolute value of [W.sup.N](0) - w(0)] = 0,

and

[absolute value of [W.sup.N]([x.sub.N/2] - w([x.sub.N/2]] [less than or equal to] [absolute value of [W.sup.N]([x.sub.N/2]] + [absolute value of w([x.sub.N/2]] [less than or equal to] C[N.sup.-1].

Now let us introduce comparison functions [[psi].sup.[+ or -].sub.i] defined by

[[psi].sup.[+ or -].sub.i] = ([tau] - [x.sub.i])[C.sub.1][([epsilon] - [delta][parallela[parallel]). sup.-2][tau][N.sup.-1] + [C.sub.2[N.sup.-1] [+ or -] ([W.sup.N] - w)([x.sub.i],

where [C.sub.1] and [C.sub.2] are arbitrary constants. Then we have

[[psi].sup.[+ or -].sub.0] = [C.sub.1][theta][N.sup.-1] [([epsilon] - [delta[parallel]a[parallel]).sup.-2] + [C.sub.2][N.sup.-1] [greater than or equal to] 0, [[psi].sup.[+ or -].sub.N/2] = [C.sub.2][N.sup.-1] [+ or -] ([W.sup.N] - w)([x.sub.N/2]).

We choose [C.sub.2] so that the first term dominate the second term on the right of the above equation which gives [[psi].sup.[+ or -].sub.N/2] [greater than or equal to] 0 and consider [MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

Now we choose the constant [C.sub.1] so that ([a.sub.i][C.sub.1] [not equal to] 2) [greater than or equal to] 0, thus all the terms on the right side in the above inequality are negative which gives

[L.sup.N.sub.2, l][[psi].sup.[+ or -].sub.i] [less than or equal to] 0, 1 [less than or equal to] I [less than or equal to] N/2 - 1.

Then by the discrete minimum principle, we have

[[psi].sup.[+ or -].sub.i] [greater than or equal to] 0, 0 [less than or equal to] I [less than or equal to] N/2,

which on simplification gives

[absolute value of ([W.sup.N] - w)([x.sub.i])] [less than or equal to] [C.sub.1][([epsilon] - [delta]M). sup.-2][[tau].sup.2][N.sup.-1] + [C.sub.2][N.sup.-1],

Since [tau] = C([epsilon] - [gamma]M) ln N, where C = 1/[theta], we obtain

[absolute value of ([W.sup.N] - w)([x.sub.i])] [less than or equal to] C[N.sup.-1] [(ln N).sup.2]. (3.33)

Now combining inequalities (3.32) and (3.33) to obtain the bound on the singular component of error throughout the interval [0, 1], we obtain

[absolute value of ([W.sup.N] - w)([x.sub.i])] [less than or equal to] C[N.sup.-1] [(ln N).sup.2], 0 [less than or equal to] i [less than or equal to] N. (3.34)

Now after combining the two inequalities, inequality (3.24) to bound the regular error component and inequality (3.34) to bound the singular error component, we obtain the required error estimate.

3.2. Layer on the right side. Now we assume that a(x) [less than or equal to] -M < 0 throughout the interval [0, 1], where M is a constant. This assumption implies that the boundary layer will be in the neighborhood of 1, i.e., on the right side of the interval [0, 1].

3.2.1. Analytical results. As we have established the estimates in section 3.1 for the solution of the continuous problem and its derivatives in the case when the solution of the problem exhibits boundary layer behavior on the left side of the interval [0, 1], one can easily obtain similar estimates in this case on the same lines. The key difference is that in this case we approximate the first derivative by the backward finite difference operator in place of the forward finite difference operator as we did in case of left side boundary layer.

3.2.2. Standard finite difference method. We use the standard backward finite upwind difference scheme to the discretize the continuous problem (3.1), (3.2) and for i = 1, 2, ... N to obtain

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

which on simplification gives

[E.sub.i][u.sub.i-1] - [F.sub.i][u.sub.i] + [G.sub.i][u.sub.i+l] = [H.sub.i], (3.35)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

The difference equations (3.35) form a tridiagonal system of N - 1 equations with N + 1 unknowns [u.sub.0], [u.sub.1], ..., [u.sub.N]. The N - 1 equations together with the given two boundary conditions are sufficient to solve the system. The stability and convergence analysis of the scheme is discussed in [6].

3.2.3. Fitted mesh finite difference method. In this case, we discretize the boundary value problem (3.1), (3.2) using the fitted mesh finite difference method composed of a standard backward upwind finite difference operator on a fitted piecewise uniform mesh, condensing at the boundary x = 1. The fitted piecewise-uniform mesh [[bar.[OMEGA]].sup.N.sub.[tau] on the interval [0, 1] is constructed by partitioning the interval into two subintervals [0, (1 - [tau])] and [(1 - [tau]), 1], where the transition parameter is chosen such that [tau] [equivalent to] min[0.5, C([epsilon] - [delta]M) ln N] with C = 1/[theta] and it is assumed that N = [2.sup.k], where k [greater than or equal to] 2 is an integer, which guarantees that there is at least one point in the boundary layer region. On each of these subintervals, a uniform mesh is placed. A fitted finite difference method for the problem (3.1), (3.2) on the piecewise uniform mesh [[bar.[OMEGA]].sup.N.[tau]], is defined by

[L.sup.N.sub.2, [tau]][u.sub.i] = f([x.sub.i]), x [member of] [[OMEGA].sup.N.sub.[tau], [u.sub.0] = [[phi].sub.0], [u.sub.N] = [gamma],

where the discrete operator [L.sup.N.sub.2, [tau]] is defined as

[L.sup.N.sub.2, [tau]][u.sub.i] = ([epsilon] - [delta]a([x.sub.i]))[D.sup.[+ or -][u.sub.i] + a([x.sub.i]) [D.sup.-][u.sup.i] + b([x.sub.i])[u.sub.i].

Also one can easily show that in this case, the solution of the discretized problem converges [epsilon]-uniformly to the solution of the continuous problem. One can obtain the error estimate in this case on the same lines as we have done in section 3.1.3 for the case of left side boundary layer.

4. Computational results. Some numerical examples are considered and solved using the methods proposed in this paper. The exact solution of the boundary value problem given by (3.1), (3.2) for constant coefficients (i.e., a(x) = a and b(x) = b are constant), with f (x) = 0 is

u(x) = [([gamma]-[[phi].sub.0] exp([m.sub.2])) exp([m.sub.1]x)-([gamma]-[[phi].sub.0] exp([m.sub.i])) exp([m.sub.2]x)]/(exp([m.sub.i])-exp([m.sub.2])),

where

[m.sub.1] = [-a + [square root of [a.sup.2] - 4([epsilon] - a[delta]b]/2([epsilon] - a[delta]),

[m.sub.2] = [-a + [square root of [a.sup.2] - 4([epsilon] - a[delta]b]/2([epsilon] - a[delta].

Example 4.1. [epsilon]y"(x) + y'(x - [delta]) - y(x) = 0, under the interval and boundary conditionsy(x) = 1, -[delta] [less than or equal to] x [less than or equal to] 0, y(1) = 1.

Example 4.2. [epsilon]y"(x) + exp(-x)y'(x - [delta]) - xy(x) = 0, under the interval and boundary conditions y(x) = 1, -[delta] [less than or equal to] x [less than or equal to] 0, y(1)=1.

Example 4.3. [epsilon]y"(x) - y'(x - [delta]) - y(x) = 0, under the interval and boundary conditions y(x) = 1, -[delta] [less than or equal to] x [less than or equal to] 0, y(1) = -1.

Example 4.4. [epsilon]y"(x) - (1 + x)y'(x - [delta]) - exp(-x)y(x) = 0, under the interval and boundary conditions y(x) = 1, -[delta] [less than or equal to] x [less than or equal to] 0, y(1) = -1.

Example 4.5. [epsilon]y"(x) - (1 + x)y'(x - [delta]) - exp(-x)y(x) = 1, under the interval and boundary conditions y(x) = 1, -[delta] [less than or equal to] x [less than or equal to] 0, y(1) = -1.

[FIGURE 4.1 OMITTED]

[FIGURE 4.2 OMITTED]

5. Discussion. A parameter uniform numerical difference scheme based on finite difference on piecewise uniform mesh is presented to solve the boundary value problems for singularly perturbed differential difference equations of the convection-diffusion type with small delay. The theoretical analysis is presented to show that the proposed method is parameter-uniform, i.e., the method converges independently of the singular perturbation parameter [epsilon].

In support of the predicted theory in the paper and to see how the delay affects the boundary layer solution of the problem, a set of numerical experiments is carried out in section 4. The proposed difference scheme is compared with the classical finite difference scheme via tabulating the maximum absolute error for the considered examples in the Tables 4.1-4.4, where [E.sup.N] = [max.sub.0<[epsilon]<1] [E.sup.N.sub.[epsilon], with [E.sup.N.sub.[epsilon] = [max.sub.0[less than or equal to]I[less than or equal to]] [absolute value of [e.sub.i]] and [e.sub.i] = ([u.sup.N.sub.i] - u([x.sub.i])). Tables 4.1 and 4.3 show that the standard upwind difference scheme on the uniform mesh works nicely till h < [epsilon] but does not behave uniformly with respect to the singular perturbation parameter [epsilon], i.e., the method is not parameter uniform; and as the condition h < [epsilon] is violated, the scheme fails in the sense that the absolute maximum error increases as the mesh parameter h decreases. Tables 4.2 and 4.4 show that the proposed fitted mesh method is parameter uniform and works nicely independent of the mesh parameter h and the singular perturbation parameter [epsilon].

[FIGURE 4.3 OMITTED]

To demonstrate the effect of delay on the boundary layer solution of the problem, the graphs of the solution for different values of [delta] are plotted in the form of Figures 4.1-4.8. From Figures 4.1-4.8, we observe that as the delay increases, the thickness of the boundary layer decreases in the case when the solution exhibits boundary layer behavior on the left side, while it increases in the case when the solution exhibits boundary layer behavior on the right side of the interval (0,1).

[FIGURE 4.4 OMITTED]

REFERENCES

[1] W. J. CUNNINGHAM, Introduction to Nonlinear Analysis, McGraw-Hill Book Company, Inc, New York, 1958.

[2] P. A. FARRELL, A. F. HEGARTY, J. J. H. MILLER, E. O'RIORDAN, AND G. I. SHISHKIN,Robust Computational Techniques for Boundary Layers, CRC Press LLC, Florida, 2000.

[3] P. A. FARRELL, E. O'RIORDAN, AND J. J. H. MILLER, Parameter-uniform fitted mesh method for quasilinear differential equations with boundary layers, Comput. Methods Appl. Math., 1 (2001), pp. 154-172.

[4] A. V. HOLDEN, Models of the Stochastic Activity of Neurons, Belin-Heidelberg, New York, Spinger-Verlag, 1976.

[5] M. K. KADALBAJOO AND K. K. SHARMA, An [epsilon]-uniform fated operator method for solving boundary value problems for singularly perturbed delay differential equations: layer behavior, Int. J. Comput. Math., 80 (2003), pp. 1261-1276.

[FIGURE 4.5 OMITTED]

[FIGURE 4.6 OMITTED]

[6]--Numerical analysis of singularly perturbed delay differential equations with layer behavior, Appl. Math. Comput., 157 (2004), pp. 11-28.

[7] R. B. KELLOGG AND A. TSAN, Analysis of some difference approximations for a singular perturbation problem without turning points, Math. Comp., 32 (1978), pp. 1025-1039.

[8] C. G. LANGE AND 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), pp.249-272.

[9]--Singular perturbation analysis of boundary-value problems for differential-difference equations. vi. small shifts with rapid oscillations, SIAM J. Appl. Math., 54 (1994), pp. 273-283.

[FIGURE 4.7 OMITTED]

[10] J. J. H. MILLER, E. O'RIORDAN, AND G. SHISHKIN, Fitted Numerical Methods for Singularly Perturbed Problems: Error Estimates in the Maximum Norm for Linear Problems in One and Two Dimension, World Scientific Publication, Singapore, 1996.

[11] J. P. SEGUNDO, D. H. PERKEL, H. WYMAN, H. HEGSTAD, AND G. P. MOORE, Input-output relations in computer-simulated nerve cell: influence predependence of excitatory pre-synaptic terminals, Kybernetik, 4 (1968), pp. 157-171.

[12] R. B. STEIN, A theoretical analysis of neuronal variability, Biophys. J, 5 (1965), pp. 173-194.

[13]--Some models of neuronal variability, Biophys. J, 7 (1967), pp. 37-68.

[14] H. TIAN, Numerical treatment of singularly perturbed delay differential equations, PhD thesis, University of Manchester, 2000. Submitted.

* Received June 7, 2004. Accepted for publication January 29, 2006. Recommended by V. Olshevsky.

M.K. KADALBAJOO ([dagger]) AND K.K. SHARMA ([double dagger])

([dagger]) Department of Mathematics, IIT Kanpur, India, (kadal@iitk.ac.in). ([double dagger]) Institut de Mathematiques-MAB, INRIA Futurs-Bordeaux, Universite Bordeaux 1, 351, Cours de la Liberation, 33405 Talence Cedex, France (kapil.sharma@math.u-bordeaux1.fr).

TABLE 4.1 The maximum absolute error for Example 4.1, [delta] = 0.5 [epsilon] (Standard Finite Difference Method) [epsilon] Number of Mesh Points (N) 64 128 [2.sup.-1] 0.00613354 0.00310678 [2.sup.-2] 0.01302363 0.00666864 [2.sup.-3] 0.02558967 0.01339603 [2.sup.-4] 0.04767658 0.02595431 [2.sup.-5] 0.08292925 0.04802051 [2.sup.-6] 0.12457740 0.08321306 [2.sup.-7] 0.11509818 0.12486142 [2.sup.-8] 0.07139804 0.11497041 [2.sup.-9] 0.03927297 0.07071252 [2.sup.-10] 0.02162100 0.03823164 [2.sup.-11] 0.01238366 0.02039224 [2.sup.-12] 0.00765658 0.01105879 [2.sup.-13] 0.00526520 0.00628303 [2.sup.-14] 0.00406245 0.00386715 [2.sup.-15] 0.00345930 0.00265211 [E.sup.N] 0.12457740 0.12486142 [epsilon] Number of Mesh Points (N) 256 512 [2.sup.-1] 0.00156363 0.00078441 [2.sup.-2] 0.00337691 0.00169922 [2.sup.-3] 0.00686154 0.00347366 [2.sup.-4] 0.01359511 0.00696570 [2.sup.-5] 0.02615483 0.01370389 [2.sup.-6] 0.04820459 0.02626121 [2.sup.-7] 0.08336200 0.04830016 [2.sup.-8] 0.12500829 0.08343839 [2.sup.-9] 0.11490790 0.12508302 [2.sup.-10] 0.07036844 0.11487702 [2.sup.-11] 0.03770840 0.07019610 [2.sup.-12] 0.01977472 0.03744615 [2.sup.-13] 0.01039297 0.01946520 [2.sup.-14] 0.00559274 0.01005922 [2.sup.-15] 0.00316455 0.00524673 [E.sup.N] 0.12500829 0.12508302 TABLE 4.2 The maximum absolute error for Example 4.1, [delta] = 0.5[epsilon] (Fitted Mesh Finite Difference Method) [epsilon Number of Mesh Points (N) 64 128 [2.sup.-1] 0.00613354 0.00310678 [2.sup.-2] 0.01302363 0.00666864 [2.sup.-3] 0.02813705 0.01589305 [2.sup.-4] 0.03203328 0.02064657 [2.sup.-5] 0.03149705 0.01954023 [2.sup.-6] 0.03088492 0.01850837 [2.sup.-7] 0.03050696 0.01791699 [2.sup.-8] 0.03030218 0.01760203 [2.sup.-9] 0.03019609 0.01744008 [2.sup.-10] 0.03014214 0.01735804 [2.sup.-11] 0.03011495 0.01731675 [2.sup.-12] 0.03010130 0.01729604 [2.sup.-13] 0.03009446 0.01728567 [2.sup.-14] 0.03009104 0.01728048 [2.sup.-15] 0.03008933 0.01727789 [E.sup.N] 0.03203328 0.02064657 [epsilon Number of Mesh Points (N) 256 512 [2.sup.-1] 0.00156363 0.00078441 [2.sup.-2] 0.00337691 0.00169922 [2.sup.-3] 0.00709501 0.00347366 [2.sup.-4] 0.01388512 0.01008758 [2.sup.-5] 0.01249807 0.00851910 [2.sup.-6] 0.01121125 0.00705171 [2.sup.-7] 0.01047647 0.00622051 [2.sup.-8] 0.01009093 0.00579052 [2.sup.-9] 0.00989678 0.00557297 [2.sup.-10] 0.00979846 0.00546403 [2.sup.-11] 0.00974899 0.00540920 [2.sup.-12] 0.00972418 0.00538170 [2.sup.-13] 0.00971175 0.00536793 [2.sup.-14] 0.00970554 0.00536104 [2.sup.-15] 0.00970243 0.00535759 [E.sup.N] 0.01388512 0.01008758 TABLE 4.3 The maximum absolute error for Example 4.3, [delta] = 0.5[epsilon] (Standard Finite Difference Method) Number of Mesh Points (N) [epsilon 64 128 [2.sup.-1] 0.00251716 0.00126635 [2.sup.-2] 0.00727328 0.00367700 [2.sup.-3] 0.01719480 0.00875821 [2.sup.-4] 0.03566613 0.01848879 [2.sup.-5] 0.07013438 0.03726925 [2.sup.-6] 0.12705085 0.07188758 [2.sup.-7] 0.22136722 0.12938377 [2.sup.-8] 0.27214642 0.22347704 [2.sup.-9] 0.20393719 0.27504772 [2.sup.-10] 0.11285694 0.20663541 [E.sup.N] 0.27214642 0.27504772 Number of Mesh Points (N) [epsilon 256 512 [2.sup.-1] 0.00063514 0.00031806 [2.sup.-2] 0.00184859 0.00092683 [2.sup.-3] 0.00442065 0.00222089 [2.sup.-4] 0.00940797 0.00474620 [2.sup.-5] 0.01924935 0.00979803 [2.sup.-6] 0.03820138 0.01973081 [2.sup.-7] 0.07284344 0.03870772 [2.sup.-8] 0.13060579 0.07334352 [2.sup.-9] 0.22456093 0.13123155 [2.sup.-10] 0.27652225 0.22511038 [E.sup.N] 0.27652225 0.22511038 TABLE 4.4 The maximum absolute error for Example 4.3, [delta] = 0.5[epsilon] (Fitted Mesh Finite Difference Method) Number of Mesh Points (N) [epsilon] 64 128 [2.sup.-1] 0.00251716 0.00126635 [2.sup.-2] 0.00727328 0.00367700 [2.sup.-3] 0.01719480 0.00875821 [2.sup.-4] 0.03566613 0.01848879 [2.sup.-5] 0.03272907 0.01887111 [2.sup.-6] 0.03373696 0.01970300 [2.sup.-7] 0.03498209 0.02118177 [2.sup.-8] 0.03571756 0.02211210 [2.sup.-9] 0.03610636 0.02263908 [2.sup.-10] 0.03630522 0.02291006 [E.sup.N] 0.03630522 0.02291006 Number of Mesh Points (N) [epsilon] 256 512 [2.sup.-1] 0.00063514 0.00031806 [2.sup.-2] 0.00184859 0.00092683 [2.sup.-3] 0.00442065 0.00222089 [2.sup.-4] 0.00940797 0.00474620 [2.sup.-5] 0.01014654 0.00704474 [2.sup.-6] 0.01027200 0.00813917 [2.sup.-7] 0.01175482 0.00588529 [2.sup.-8] 0.01275544 0.00686774 [2.sup.-9] 0.01331548 0.00743558 [2.sup.-10] 0.01360524 0.00773978 [E.sup.N] 0.01360524 0.00813917

Printer friendly Cite/link Email Feedback | |

Author: | Kadalbajoo, M.K.; Sharma, K.K. |
---|---|

Publication: | Electronic Transactions on Numerical Analysis |

Date: | Apr 1, 2006 |

Words: | 10036 |

Previous Article: | On fast factorization pivoting methods for sparse symmetric indefinite systems. |

Next Article: | Numerical study of normal pressure distribution in entrance flow between parallel plates, I. Finite difference calculations. |