Integration Technique for Singularly Perturbed Delay Differential Equations
D. Kumara Swamy1, A. Benerji Babu1, Y.N. Reddy1, K. Phaneendra2*
Abstract— In this paper, we have described a numerical integration technique for solving singularly perturbed delay differential equations. The second order singularly perturbed boundary value problem is transformed into an asymptotically equivalent first order neutral differential equation. Then numerical integration and linear interpolation is used to get the tri-diagonal system. Discrete invariant imbedding algorithm is used to solve this tridiagonal system. The error analysis of the technique is discussed. To demonstrate the technique and affect of the delay argument in the layer, we have implemented the technique on several test examples.
Index Terms— Singularly perturbed delay differential equation, Boundary layer, Trapezoidal rule, Linear interpolation, Tridiagonal system, maximum absolute error.
—————————— ——————————
1 INTRODUCTION
Any ordinary differential equation, in which the highest derivative is multiplied by a small parameter and involving at least one delay term is called singularly perturbed delay dif- ferential equation. In recent years, there has been a growing interest in the numerical treatment of such differential equa- tions. The boundary value problems of delay differential equa- tions are ubiquitous in the variational problems in control the-
technique is discussed. To demonstrate the technique and affect of the delay argument in the layer, we have implement- ed the technique on several test examples.
DESCRIPTION OF THE METHOD
Consider a singularly perturbed delay differential equation
ory [3]. Lange and Miura [7, 8] gave an asymptotic approach for a class of boundary-value problems for linear second-order differential-difference equations in which the highest order
εy′ (x) + a(x) y′(x − δ ) + b(x) y(x) = f (x)
on 0 < x <1, 0 < δ << 1, with
(1)
derivative is multiplied by small parameter and shows the
y(x) = φ(x) ,
−δ ≤ x ≤ 0 , (2)
effect of very small shifts on the solution and pointed out that they drastically affect the solution and therefore cannot be neglected. Kadalbajoo and Sharma [6] presented a numerical approach to solve singularly perturbed differential-difference equation, which contains only negative shift in the differenti- ated term. In this method authors present a numerical method composed of a standard upwind finite difference scheme on a
special type of mesh shifts are either o(ε ) or O(ε ). Pratima
y(1)= γ (3)
where a(x), b(x),f(x) are smooth functions , γ is a constant and
δ is the delay. For the function y(x) be a smooth solution to the problem (1), it must satisfy: boundary value problem be continuous on [0, 1] and be continuously differentiable on (0, 1).
Layer on the left side
Rai and Sharma [10] presented a numerical method for singu-
Here, we consider the case
a(x) ≥ M > 0, ∀x ∈[0,1], M being
larly perturbed delay differential equation with turning points. Reddy et. al. [11] presented a numerical integration of a class of singularly perturbed delay differential equations with small shift, where delay is in differentiated term.
In In this paper, we have described a numerical integration technique for solving singularly perturbed delay differential equations. The second order singularly perturbed boundary value problem is transformed into an asymptotically equiva- lent first order neutral differential equation. Then numerical
positive constant. In this case the solution of the boundary
value problem exhibits boundary layer behaviour on the left
side of the interval [0, 1] i.e., at x = 0.
In this section we construct a numerical scheme for solving
the boundary value problem based on numerical integration
for the case when the solution of the problem exhibits a layer
on the left side.
By using Taylor Series expansion in the neighbourhood of
the point x , we have
integration and linear interpolation is used to get the tri- diagonal system. Discrete invariant imbedding algorithm is
y′(x − ε ) ≈ y′ − εy′
Therefore,
and
y′(x + ε ) ≈ y′ + εy′
used to solve this tri-diagonal system. The error analysis of the
1Department of Mathematics, National Institute of Technology, Warangal- 506004
2Department of Mathematics, Nizam College, Osmania University, Hydera-
εy′ = y′(x + ε ) − y′(x − ε )
2
Substituting the above equation in (1), it reduces to an as- ymptotically equivalent first order neutral differential equa-
bad-500001,email:kollojuphaneendra@yahoo.co.in
IJSER © 2013
http://www.ijser.org
tion
ε yi +1 −yi −
ε yi −yi −1
y′(x + ε ) − y′(x − ε ) = p(x) y′(x − δ ) + q(x) y(x) + r(x)
(4)
2 2
h h
where
p(x) = −2a ,
q(x) = −2b , r(x) = 2 f (x)
= p
δ
h
2
− − p′
− δ
=h q y
The transition from equation (1) to equation (4) is admitted,
because of the condition that ε is small. This replacement is
i+1 1
h
i+11
+
h 2
i+1
i+1
significant from the computational point of view. Further
δ δ h
h
δ h
δ h
+
details on the validity of this transition can be found [5, 9].
Now divide the interval [0, 1] into N equal subintervals of
+ pi+1
pi 1 −
−
h 2
pi′+1 −
h 2
pi′1 −
h
qi yi
2
+ − p δ − h p′ δ y
h r h r
mesh size h =1/N so that xi = ih, i = 0, 1, 2, …, N.
h
h
2
Integrating eq. (4) with respect to x from xi to xi +1 , we get
i i
i−1 + 2
i + 2
i+1
xi+1
∫
y′(x + ε )dx −
xi+1
∫
y′(x − ε )dx =
Rearrange the above equation into three recurrence rela- tion, we get
xi xi
p δ + δ
p′ +
2ε y
xi+1
∫
[ p(x) y′(x − δ )dx + q(x) y(x) + r(x)]dx
h
2
i i
h
i−1
xi 4ε
− +
δ − − δ − δ ′
− h ′ − δ h
h pi+1 h
pi 1 h
pi+1
2
pi 1
2
+ qi yi
h 2
2ε
+ − p
δ
−
h
2
+ p′
δ
−
h
− q y
h r h r
[ p(x) y(x − δ )] xi+1
xi+1
∫
p′(x) y(x − δ )dx
i+1 1
h h
i+1 1
h
i+1
2
i+1 = 2
i + 2
i+1
xi xi
xi+1
The above equation can be written as
+ ∫ [q(x) y(x) + r(x)]dx
xi
Ei yi−1 − Fi yi + Gi yi+1 = Hi , i = 1,2,3......n − 1
(8)
2εyi′+1 − 2εyi′ = pi+1 y(xi+1 − δ ) − pi y(xi − δ )
where
δ δ 2ε
xi+1
− ∫
xi
p′(x) y(x − δ )dx +
xi+1
∫
xi
q(x)y(x) + r(x)dx
Ei = pi h + 2
4ε
pi′ +
h
δ
δ δ
h δ h
By using the Trapezoidal rule to evaluate the integral
Fi =
h + pi+1 h − pi 1 − h − 2 pi′+1 − 2 pi′1 − h + 2 qi
approximation, we get
h
2ε
δ h
δ h
2εyi′+1 − 2εyi′ = pi+1y(xi+1 − δ ) − pi y(xi − δ ) −
(pi′+1y(xi+1 − δ ))
Gi =
h − pi+11 − h + 2 pi′+11 − h − 2 qi+1
2
− + ri
h (p′ y(x i i |
− δ ))+ h (q y i+1 |
h i+1)+ (qi y |
)+ h r i i+1 |
h | H = h r + h r |
2 | 2 | 2 | 2 | 2 |
We solve the tridiago |
i 2 i+1 2 i
(5)
Again by means of Taylor series expansion and linear
imbedding algorithm.
nal system (5) using discrete invariant
interpolation for y′(x) , we get
−
y(x
δ ) ≈ y(x
) − δ y′(x
) = y
δ yi+1
yi
i+1
δ
i+1
δ
i+1
i+1
h
(6)
Now for the right - end boundary layer, integrating Eq. (4) with respect to x from xi−1 to xi , we get
= 1 − yi+1 + yi
h h
y − y
xi
∫ y′(x + ε )dx −
xi
∫ y′(x − ε )dx =
xi
∫ [ p(x) y′(x − δ )dx + q(x) y(x) + r(x)]dx
y(xi − δ ) ≈ y(xi ) − δy′(xi ) = yi − δ
i i−1
h
(7)
xi−1
xi−1
xi −1
δ δ x
= 1− yi + yi−1
h h
[y(x + ε ) − y(x
+ ε )]− [y(x − ε ) + y(x
− ε )] = [ p(x) y(x − δ )] i
i
i−1
i i−1
xi−1
Substituting the equations (6) and (7) in Eq. (5), we get
xi
− ∫ p′(x) y(x − δ )dx +
xi−1
xi
∫ [q(x) y(x) + r(x)]dx
xi−1
y + ε yi+1 −
i
y
− y
− ε (y − y
i
)− y + ε yi+1 − y
where E
= 2ε + p
δ h
1 + + p′
(1 +
δ ) − h q
i h
ε
h
i−1 i
i−1
i h
i h
4ε h
i−1
h
δ
2 i−1
δ
h 2 i−1
δ δ h
+ yi−1 −
( yi − yi−1) = pi y(xi − δ ) − pi−1y(xi−1 − δ )
Fi =
− pi′1 + + pi (1 + ) + pi−1 + pi′−1 + qi
h 2 h h h 2 2
h
xi xi
2ε δ δ
− ∫ p′(x) y(x − δ )dx +
xi−1
∫ [q(x)y(x) + r(x)]dx
xi−1
Gi =
h + pi h − 2 Pi′
h h
By Using the Trapizoidal rule to evaluate the integral
approximation, we get
Hi = ri + ri−1
2 2
2ε ( y − y ) −
h i+1 i
2ε ( y − y ) = p y(x − δ ) − p y(x − δ )
h i i−1 i i i−1 i−1
h
NUMERICAL EXAMPLES
To describe the method we consider six test examples with
− (pi′ y(xi − δ ) + pi′−1y(xi−1 − δ )) 2
h h
left and right end boundary layers.
+ (qi yi + qi−1yi−1)+
(ri + ri−1)
Example 1.
εy′ (x) + y′(x − δ ) − y(x) = 0 ;
x ∈[0,1] under the
2
2ε 4ε 2ε
2 interval with boundary conditions y(x) = 1, −δ ≤ x ≤ 0, y(1) = 1
h The exact solution is given by
h yi+1 − h
yi + h
yi−1 = pi y(xi − δ ) − pi−1y(xi−1 − δ ) − 2 pi′ y(xi − δ )
((1 - em2 )em1x + (em1 − 1)em2x )
h
− pi′−1y(xi−1 − δ ) +
2
h
2 (qi yi ) +
h
2 (qi−1yi−1) +
h
(ri + ri−1)
2
y(x) =
(em1 − em2 )
(9)
Again by means of Taylor series expansion and then corre-
where
m1 = (−1 −
1 + 4(ε − δ ) )
2(ε − δ )
and
sponding
y′(x) by linear interpolation, we have
−
m = (−1 +
1 + 4(ε − δ ) ) .
y(x
δ ) ≈ y(x
) − δ y′(x
) = y
δ yi
yi−1 2
i−1
i−1
i−1
i−1
δ δ
(10)
Example 2.
h
2(ε − δ )
εy′ (x) + 0.25 y′(x − δ ) − y(x) = 0
under the interval
= 1 + yi−1 − yi
h h
with boundary conditions y(x) = 1, −δ ≤ x ≤ 0, y(1) = 0
y − y
Example 3.
εy′ (x) − y′(x − δ ) − y(x) = 0 ;
under the interval
y(xi − δ ) ≈ y(xi ) − δy′(xi ) = yi − δ
i+1 i
h
(11)
with boundary conditions y(x) = 1, −δ ≤ x ≤ 0, y(1) = -1 The exact solution is given by
δ δ
m m x m m x
= 1 + yi − yi+1
h h
y(x) = ((1 + e
2 )e 1
(e
1 + 1)e 2 )
m m
2 |
y − | 4 |
y + | 2 | y = p + | δ |
y − | δ | y |
h | i+1 | h | i | h | | h | i | h | i+1 |
Substituting the equations (10) and (11) in Eq.(9), we get
(e 2 − e 1 )
ε ε ε
i−1
i 1
where m1
= (1 −
1 + 4(ε + δ ) )
2(ε + δ )
and
+
δ
− p y
δ
− y −
h
p′ +
δ y
δ y
m2 = (1 +
1 + 4(ε + δ ) )
2(ε + δ ) .
i−1 1
| boundary conditions y(x) = 1, | −δ ≤ x ≤ 0, | y(1) = - 1 |
h | | δ | | δ | | h | (qi yi )+ | h (q y ) i−1 i−1 | | | |
2 | | h | | h | | 2 | | 2 | 4. DISCUSSIONS AND CONCLUSIONS | |
h i−1 h i
i 1
2
h i − h
i+1
Example 4. εy′ (x) − y′(x − δ ) + y(x) = 0 under the interval with
− pi′−11 +
yi−1 −
yi +
h h In this paper, we have described a numerical integration
+ ri + ri−1
2 2
Rearrange the above equation into three recurrence relation, we get
technique for solving singularly perturbed delay differential equations. The second order singularly perturbed boundary value problem is transformed into an asymptotically equiva- lent first order neutral differential equation. Then numerical integration and linear interpolation is used to get the tri-
Ei yi−1 − Fi yi + Gi yi+1 = Hi , i = 1,2,3......n − 1
(12)
diagonal system. Discrete invariant imbedding algorithm is used to solve this tri-diagonal system. The error analysis of the
technique is discussed. To demonstrate the technique and affect of the delay argument in the layer, we have implement- ed the technique on several test examples and we have shown the layer behaviour through graphs.
From the numerical examples presented here, we observed that as δ increases, the thickness of the left end boundary layer decreases and as δ increases, the thickness of right end boundary layer increases. As the grid size h decreases, the maximum error decreases, which shows the convergence to the computed solution.
TABLE 1. THE MAXIMUM ABSOLUTE ERRORS FOR δ = 0.03
1
0.9
0.8
delta=0.01
exact computed ******
ε / N
Example 1
100 200 300 400 500
0.7
delta=0.03 delta=0.06
2−1
2−2
2−3
1.5805e-003 7.9569e-004 5.3170e-004 3.9924e-004 3.1962e-004
4.1130e-003 2.0791e-003 1.3911e-003 1.0452e-003 8.3706e-004
9.1564e-003 4.6690e-003 3.1325e-003 2.3570e-003 1.8894e-003
0.6
0.5
delta=0.08
2−4 1.9899e-002 1.0335e-002 6.9813e-003 5.2710e-003 4.2338e-003
0.4
2−5
4.5123e-002 2.4551e-002 1.6873e-002 1.2868e-002 1.0397e-002
0 0.2 0.4 0.6 0.8 1
Example 2
Fig.1.The numerical solution of example 1 with ε = 0.1
2−1
6.5308e-004 3.2800e-004 2.1900e-004 1.6437e-004 1.3156e-004
2−2 1.2766e-003 6.4293e-004 4.2966e-004 3.2264e-004 2.5830e-004
2−3 2.3695e-003 1.1981e-003 8.0175e-004 6.0245e-004 4.8251e-004 1
δ=0.3ε
2−4
2−5
4.2895e-003 2.1844e-003 1.4653e-003 1.1024e-003 8.8353e-004
8.0121e-003 4.1358e-003 2.7884e-003 2.1028e-003 1.6880e-003
0.9
0.8
δ=0.6ε
δ=0.9ε
Example 3
2−1 3.7887e-003 1.9299e-003 1.2948e-003 9.7419e-004 7.8082e-004
2−2 1.5347e-002 7.9592e-003 5.3719e-003 4.0537e-003 3.2550e-003
0.7
numerical solution
0.6
0.5
2−3
2−4
2−5
3.4927e-002 1.8498e-002 1.2574e-002 9.5224e-003 7.6640e-003
7.5272e-002 4.3189e-002 3.0328e-002 2.3451e-002 1.9078e-002
8.0792e-002 1.2306e-001 1.4428e-001 1.5343e-001 1.5555e-001
0.4
0.3
0.2
Example 4
2−1
2−2
2−3
4.7433e-003 2.3865e-003 1.5943e-003 1.1969e-003 9.5814e-004
1.0028e-002 5.0628e-003 6.4266e-003 4.8333e-003 2.0371e-003
1.8863e-002 9.5867e-003 6.2857e-003 4.7533e-003 3.8732e-003
0.1
0
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x
2−4
2−5
3.3542e-002 1.7246e-002 1.1608e-002 8.7488e-003 7.0204e-003
5.6064e-002 2.9375e-002 1.9905e-002 1.5053e-002 1.2103e-002
Fig.2.The numerical solution of example 2 with ε = 0.01
REFERENCES
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
exact computed ******
delta=0.000 delta=0.007
delta=0.015 delta=0.025
[1] E. Angel and R. Bellman, “Dynamic Programming and Partial differen tial equations, Academic Press, New York, 1972.
[2] R. Bellman and K. L. Cooke, “Differential-Difference Equations”, Acdemaic Press, New York, USA, 1963.
[3] M.W. Derstine, F.A.H.H.M. Gibbs and D. L. Kaplan, “Bifurcation gap in a hybrid optical system”, Phys. Rev. A, vol. 26, pp. 3720–3722, 1982.
[4] R. D. Driver, “Ordinary and Delay Differential Equations”, Belin- Heidelberg, New York, Springer, 1977.
[5] L. E. El’sgol’ts and S. B. Norkin, “Introduction to the Theory and Applica- tions of Differential Equations with Deviating Arguments”, Academic Press, New York, 1973.
[6] M.K. Kadalbajoo and K.K. Sharma, “A numerical method based on finite difference for boundary value problems for singularly perturbed delay differential equations”, Applied Mathematics and Computation, vol. 197, pp. 692–707, 2008.
[7] Lange, C.G. Miura, R.M. Singular perturbation analysis of boundary- value problems for differential-difference equations. v. small shifts with layer behavior, SIAM J. Appl. Math., 54, pp. 249–272, 1994.
0 0.2 0.4 0.6 0.8 1
Fig 3.The numerical solution of example 3 with ε = 0.1
1
δ=0.3ε
[8] C.G. Lange and R.M. Miura, “Singular perturbation analysis of boundary-value problems for differential-difference equations. vi. Small shifts with rapid oscillations”, SIAM J. Appl. Math., vol. 54, pp. 273–283, 1994.
[9] R.E. O’Malley, “Singular Perturbation Methods for Ordinary Differential Equations, Springer-Verlag, New York, 1991.
[10] Pratima Rai and Kapil K. Sharma, “Numerical analysis of singularly perturbed delay differential turning point problem”, Applied Mathe- matics and Computation, vol. 218, pp. 3483- 3498, 2011.
[11] Y. N. Reddy, G. B. S. L. Soujanya and K. Phaneendra, “Numerical Integration Method for Singularly Perturbed Delay Differential Equations”, International Journal of Applied Science and Engineering,
0.8
δ=0.7ε
δ=0.9ε
vol. 10, no. 3, pp. 249-261, 2012.
0.6
0.4
numerical solution
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
x
Fig 4.The numerical solution of example 4 with ε = 0.01