In this paper, a nonasymptotic method is presented for solving singularly perturbed delay differential equations whose solution exhibits a boundary layer behavior. The second order singularly perturbed delay differential equation is replaced by an asymptotically equivalent first order neutral type delay differential equation. Then, Simpson’s integration formula and linear interpolation are employed to get three term recurrence relation which is solved easily by Discrete Invariant Imbedding Algorithm. Some numerical examples are given to validate the computational efficiency of the proposed numerical scheme for various values of the delay and perturbation parameters.
AMS Mathematics Subject Classification : 65L10.
1. Introduction
A singularly perturbed delay differential equation is an ordinary differential equation in which the highest derivative is multiplied by a small parameter and containing delay term. In recent years, Delay differential equations have attracted the attention of many researchers because of their applications in many scientific and technical fields such as in population ecology, control theory, viscoelasticity, and materials with thermal memory, El’sgolts
[4]
, and Bellman and Cooke
[3]
. Boundary value problems involving delay differential equations are also used to study signal transmission with time delays in control theory El’sgolts
[4]
, firstexit problems in neurobiology Stein
[13]
, and Tuckwell
[14]
. As a result, many researchers tried to develop and present numerical schemes for solving such problems. For example, Amiraliyev and Cimen
[1]
, presented numerical method comprising a fitted difference scheme on a uniform mesh to solve second order delay differential equations. Lange and Miura
[10
,
11]
gave an asymptotic approach for a class of boundaryvalue problems for linear secondorder differentialdifference equations. Kadalbajoo and Sharma
[7
,
8
,
9]
presented numerical approaches to solve singularly perturbed differentialdifference equations, which contains negative shift in the convention term (i.e. in the derivative term). Furthermore, Phaneendra, et al.
[12]
, presented higher order numerical method to solve the problems of the same type.
The objective of this paper is, therefore, to present a nonasymptotic method for solving singularly perturbed delay differential equations. The method is named as nonasymptotic as it does not depend on asymptotic expansion as well as on the matching of the coefficients. In this method, we first replace the second order singularly perturbed delay differential equation by an asymptotically equivalent first order neutral type delay differential equation and employ Simpson’s rule for the numerical integration of the obtained first order delay differential equation. Then, linear interpolation is used to get the three term recurrence relation which is solved easily by discrete Invariant Imbedding algorithm. The method is demonstrated by implementing several model examples by taking various values for the delay and perturbation parameters.
2. Description of the Method
Consider singularly perturbed singular boundary value problems of the form
with the interval and boundary conditions
where
ε
is small parameter, 0 <
ε
<< 1, and the delay parameter,
δ
is of
o
(
ε
) satisfying (
ε
−
δ
a
(
x
)) > 0 for all
xϵ
[0, 1];
a
(
x
),
b
(
x
),
f
(
x
) and
ϕ
(
x
) are smooth functions and
β
is a finite constant. For
δ
= 0, the problem (1)(3) becomes a boundary value problem for singularly perturbed ordinary differential equations. The layer behavior of the problem under consideration is maintained only for
δ
≠ 0 but sufficiently small (i.e
δ
is of
o
(
ε
)). When the delay parameter,
δ
exceeds the perturbation parameter,
ε
(i.e
δ
is of
O
(
ε
)), the layer behavior of the solution is no longer maintained, rather the solution exhibits an oscillatory behavior or diminished. In this paper, we consider the cases where the boundary layer is either on the left side or right side of the interval [0, 1].
2.1. Layer on the Left Side:
In general, the solution of problem (1)(3) exhibits boundary layer behavior at one end of the interval [0, 1] depending on the sign of
a
(
x
). Now, we assume that
a
(
x
) ≥
M
> 0 throughout the interval [0, 1], where
M
is positive constant. This assumption merely implies that the boundary layer will be in the neighborhood of
x
= 0. Further, it is also assumed that
b
(
x
) satisfies the condition that
b
(
x
) ≤ −
θ
< 0 , ∀
xϵ
[0, 1] where
θ
is a positive constant.
By using Taylor series expansion in the neighborhood of the point
x
, we have
By using (4) in (1) we get an asymptotically equivalent first order delay neutral type differential equation:
The transition from equation (1) to equation (5) is admitted, because of the condition that
ε
is small, 0 <
ε
<< 1; and the truncated term
ε
^{2}
/2 is very small and approaches to 0 as
ε
is small. Further details on the validity of this transition can also be found in Elsgolt’s and Norkin
[5]
.
Now we divide the interval [0, 1] into
N
equal subintervals of mesh size
h
= 1/
N
so that
x_{i}
=
ih
,
i
= 0, 1, 2, · · · ,
N
.
Integrating equation (5) with respect to
x
from
x_{i}
to
x
_{i+1}
for
i
= 1, 2, · · · ,
N
−1, we get
where
y_{i}
=
y
(
x_{i}
),
a_{i}
=
a
(
x_{i}
),
b_{i}
=
b
(
x_{i}
),
f_{i}
=
f
(
x_{i}
).
By using Simpson’s rule to evaluate the integral, we get:
Since delay term is of the small order of
ε
, to tackle the term containing delay, we can use Taylor’s approximation (Kadalbajoo and Sharma
[9]
, Lange and Miura
[11]
). Thus, by means of Taylor series expansion and then by approximating
by linear interpolation, we get:
In similar way,
Hence, by making use of equations (8)(12) in (7) we obtain:
To make equation (13) a three term recurrence relation, we can express
y_{i}
_{+1/2}
in terms of
y_{i}
_{−1}
,
y_{i}
and
y_{i}
_{+1}
using Hermite’s interpolation as follows:
In view of equations(5) and (14), we get:
By making use equations (8) (12) in (15) and finite difference approximations, we get
Finally, making use of equations (16) in (13) and rearranging as three term recurrence relation, we get:
for
i
= 1, 2, · · · ,
N
− 1, where
We solve the above tridiagonal system by using method of invariant imbedding system described in section 3.
2.2. RightEnd Boundary Layer Problem:
We now assume that
a
(
x
) ≤ −
M
< 0 throughout the interval [0, 1], where
M
is positive constant. This assumption merely implies that the boundary layer will be in the neighborhood of
x
= 1. By using Taylor series expansion in the neighborhood of the point
x
, we have:
and consequently, equation (1) is replaced by the following first order differential equation:
Now we divide the interval [0, 1] into
N
equal subintervals of mesh size
h
= 1/
N
so that
x_{i}
=
ih
,
i
= 0, 1, 2, · · · ,
N
. Integrating equation (23) with respect to
x
from
x_{i}
_{−1}
to
x_{i}
, we get:
By using Simpson’s rule to evaluate the integral, we get:
Again, by means of Taylor series expansion and then by approximating
by linear interpolation, we get:
In similar way,
Hence, by making use of equations (26)(30) in (25) we obtain:
To make equation (31) a three term recurrence relation, we can express
y_{i}
_{−1/2}
in terms of
y_{i}
_{−1}
,
y_{i}
and
y_{i}
_{+1}
using Hermite’s interpolation as follows:
In view of equations (23) and (32), we get:
By making use of equations (26) (30) in (33) and finite difference approximations, we get:
Finally, making use of equations (34) in (31) and rearranging as three term recurrence relation, we get:
for
i
= 1, 2, · · · ,
N
− 1, where
We solve the above tridiagonal system by using method of invariant imbedding described in the next section.
3. Discrete Invariant Imbedding Algorithm
We now describe the Thomas algorithm which is also called Discrete Invariant Imbedding (Angel and Bellman
[2]
) to solve the three term recurrence relation:
for
i
= 1, 2, · · · ,
N
− 1. Let us set a difference relation of the form:
where
W_{i}
=
W
(
x_{i}
) and
T_{i}
=
T
(
x_{i}
) are to be determined. From (41), we have:
substituting (42) in (40), we get:
By comparing (43) and (41), we get the recurrence relations:
To solve these recurrence relations for
i
= 1, 2, 3, · · · ,
N
−1 , we need the initial conditions for
W
_{0}
and
T
_{0}
. If we choose
W
_{0}
= 0 , then we get
T
_{0}
=
ϕ
(0) . With these initial values, we compute
W_{i}
and
T_{i}
for
i
= 1, 2, 3, · · · ,
N
− 1 from equations (44) and (45) in forward process, and then obtain
y_{i}
in the backward process from (41).
The conditions for the stability of discrete invariant imbedding algorithm are (Angel and Bellman
[2]
, Kadalbajoo and Reddy
[6]
):
In our method, one can easily show that if the assumptions
a
(
x
) > 0,
b
(
x
) < 0 and (
ε
−
δa
(
x
)) > 0 hold, then the above conditions (46) hold and thus the invariant imbedding algorithm is stable.
4. Numerical Examples
To demonstrate the applicability of the method, we have considered two numerical experiments with leftend boundary layer and two numerical experiments with rightend boundary layer. We compared the computed results with the exact solution of the problems where the exact solution of the problems is known; and we also have tested the effect of small delay parameter on computed solution of the problem for different values of
δ
of
o
(
ε
) which are presented by numerical experiments whose exact solutions are not known.
Example 4.1
. Consider an example of singularly perturbed delay differential equation with left layer:
The exact solution is given by
where
The computational results are presented in the
Tables 1
,
2
,
3
and
4
for
ε
= 0.01 and 0.001 for different values of
δ
.
Example 4.2.
Now we consider an example of variable coefficient singularly perturbed delay differential equation with left layer:
For which the exact solution is not known.This example is also considered to show the effect of the small shift on the boundary layer solution. The computational results are presented in the
Tables 5
and
6
for
ε
= 0.01 and 0.001 for different values of
δ
.
Example 4.3.
Consider a singularly perturbed delay differential equation with right layer:
The exact solution is given by
where
The computational results are presented in the
Tables 7
and
8
for
ε
= 0.01 and 0.001 for different values of
δ
.
Example 4.4.
Now we consider an example of variable coefficient singularly perturbed delay differential equation with right layer:
For which the exact solution is not known. The computational results are presented in the
Tables 9
and
10
for
ε
= 0.01 and 0.001 for different values of
δ
.
5. Discussion and Conclusion
We have presented a numerical integration method to solve singularly perturbed delay differential equations whose solutions exhibits the layer behavior. The scheme is repeated for different choices of the delay parameter,
δ
and perturbation parameter
ε
. The choice of
δ
is not unique, but can assume any number of values satisfying the condition,
δ
(
ε
) =
τε
with
τ
=
O
(1) and
τ
is not too large Lange and Miura
[10]
.
To demonstrate the efficiency of the method, some numerical examples of boundary value problems with constant and variable coefficients are considered for different values of
ε
and
δ
. Besides solving for the numerical solutions, we have considered numerical results of boundary value problems, examples 4.2 & 4.4 , whose exact solutions are not known to show the effect of delay parameter on the boundary layer solutions.
From the tabulated solutions of the numerical experiments considered, it can be observed that the present method approximates the exact solution very well (
Tables 1

4
,
7

8
) for
h
≥
ε
for which other classical finite difference methods fails to give good results; and the small shift,
δ
affects both the boundary layer solutions (left and right) in similar fashion but reversely. That is, as δ increases the thickness of the left boundary layer decreases (
Tables 5

6
) while that of the right boundary layer increases (
Tables 9

10
). This method does not depend on asymptotic expansion as well as on the matching of the coefficients. Thus, we have devised an alternative technique of solving boundary value problems for singularly perturbed delay differential equations, which is easily implemented on computer and is also practical.
Numerical results of Example 4.1 forε= 0.01,δ= 0.001, N = 100
Numerical results of Example 4.1 for ε = 0.01, δ = 0.001, N = 100
Numerical results of Example 4.1 forε= 0.01,δ= 0.008, N = 100
Numerical results of Example 4.1 for ε = 0.01, δ = 0.008, N = 100
Numerical results of Example 4.1 forε= 0.001,δ= 0.0001, N = 100
Numerical results of Example 4.1 for ε = 0.001, δ = 0.0001, N = 100
Numerical results of Example 4.1 forε= 0.001,δ= 0.0008, N = 100
Numerical results of Example 4.1 for ε = 0.001, δ = 0.0008, N = 100
Numerical results of Example 4.2 forε= 0.01,N= 100 and different values ofδ.
Numerical results of Example 4.2 for ε = 0.01, N = 100 and different values of δ.
Numerical results of Example 4.2 forε= 0.001,N= 100 and different values ofδ.
Numerical results of Example 4.2 for ε = 0.001, N = 100 and different values of δ.
Numerical results of Example 4.3 forε= 0.01,δ= 0.001, N = 100
Numerical results of Example 4.3 for ε = 0.01, δ = 0.001, N = 100
Numerical results of Example 4.3 forε= 0.01,δ= 0.008, N = 100
Numerical results of Example 4.3 for ε = 0.01, δ = 0.008, N = 100
Numerical results of Example 4.4 forε= 0.01,N= 100 and different values ofδ.
Numerical results of Example 4.4 for ε = 0.01, N = 100 and different values of δ.
Numerical results of Example 4.4 forε= 0.001,N= 100 and different values ofδ.
Numerical results of Example 4.4 for ε = 0.001, N = 100 and different values of δ.
BIO
Gemechis File received his B.Sc from Bahir Dar University, M.Sc from Addis Ababa University, and currently persuing Ph.D. from National Institute of Technology, Warangal. His research interest focuses on the numerical treatments of singular perturbation problems.
Department of Mathematics, National Institute of Technology, Warangal, 506 004 A.P, India.
email: gammeef@yahoo.com
Y. N. Reddy received his M.Sc from NIT Warangal, and Ph.D. from IITKampur. His research interest focuses on the numerical treatments of singular perturbation problems.
Department of Mathematics, National Institute of Technology, Warangal, 506 004 A.P, India.
email: ynreddy@nitw.ac.in
Amiraliyev G. M.
,
Cimen E.
(2010)
Numerical Method For A Singularly Perturbed Convection Diffusion Problem with Delay
Applied Mathematics and Computation
216
2351 
2359
DOI : 10.1016/j.amc.2010.03.080
Angel E.
,
Bellman R.
1972
Dynamic Programming and Partial differential equations
Academic Press
New York
Bellman R.
,
Cooke K. L.
1963
Differential Difference Equations
Academic Press
New York
El’sgolts L.E.
1964
Qualitative Methods in Mathematical Analysis
Trans. Math. Monogr.
American Mathematical Society, Providence, R. I.
12
Elsgolt’s L. E.
,
Norkin S. B.
1973
Introduction to the Theory and Applications of Differential Equations with Deviating Arguments
Academic Press
New York
Kadalbajoo M. K.
,
Reddy Y. N.
(1986)
A non asymptotic method for general linear singular perturbation problems
Journal of Optimization Theory and Applications
55
256 
269
Kadalbajoo M. K.
,
Sharma K. K.
(2005)
Numerical treatment of boundary value problems for secondorder singularly perturbed delay differential equations
Comput. Appl. Math.
24
151 
172
DOI : 10.1590/S010182052005000200001
Kadalbajoo M. K.
,
Sharma K. K.
(2008)
A numerical method on finite difference for boundary value problems for singularly perturbed delay differential equations
Applied Mathematics and Computation
197
692 
707
DOI : 10.1016/j.amc.2007.08.089
Kadalbajoo M. K.
,
Sharma K. K.
(2004)
Numerical analysis of singularly perturbed delay differential equations with layer behavior
Applied Mathematics and Computation
157
11 
28
DOI : 10.1016/j.amc.2003.06.012
Lange C. G.
,
Miura R. M.
(1994)
Singular perturbation analysis of boundaryvalue problems for differentialdifference equations. v. small shifts with layer behavior
SIAM J. Appl. Math.
54
249 
272
DOI : 10.1137/S0036139992228120
Lange C. G.
,
Miura R. M.
(1994)
Singular perturbation analysis of boundaryvalue problems for differentialdifference equations. vi. Small shifts with rapid oscillations
SIAM J. Appl. Math.
54
273 
283
DOI : 10.1137/S0036139992228119
Phaneendra K.
,
Reddy Y. N.
,
Soujanya GBSL.
(2011)
A seventh order numerical method for singularly perturbed differential difference equations with negative shift
Non linear Analysis: Modeling and Control
16
206 
219
Tuckwell H. C.
(1976)
On the first exit time problem for temporally homogeneous Markov processes
J. Appl. Probab.
13
39 
46
DOI : 10.2307/3212663