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.

1. # 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).

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-

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

pi1

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

i1 + 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

i1

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 yi1 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 pi1 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(xi i − δ ))+ h (q yi+1 hi+1)+ (qi y )+ h ri 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

2. ### Right End boundary Layer

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 xi1 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 i1

h

(7)

xi1

xi1

xi 1

δ δ x

= 1−  yi + yi1

h h

[y(x + ε ) y(x

+ ε )][y(x ε ) + y(x

ε )] = [ p(x) y(x δ )] i

i

i1

i i1

xi1

Substituting the equations (6) and (7) in Eq. (5), we get

xi

p(x) y(x δ )dx +

xi1

xi

[q(x) y(x) + r(x)]dx

xi1

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

i1 i

i1

i h

i h

4ε h

i1

h

δ

2 i1

δ

h 2 i1

δ δ h

+ yi1

( yi yi1) = pi y(xi δ ) pi1y(xi1 δ )

Fi =

pi1 +  + pi (1 + ) + pi1 + pi1 + qi

h 2 h h h 2 2

h

xi xi

 

2ε δ δ

p(x) y(x δ )dx +

xi1

[q(x)y(x) + r(x)]dx

xi1

Gi =

h + pi h 2 Pi

h h

By Using the Trapizoidal rule to evaluate the integral

approximation, we get

Hi = ri + ri1

2 2

2ε ( y y )

h i+1 i

2ε ( y y ) = p y(x δ ) p y(x δ )

h i i1 i i i1 i1

h

2. # NUMERICAL EXAMPLES

To describe the method we consider six test examples with

(piy(xi δ ) + pi1y(xi1 δ )) 2

h h

left and right end boundary layers.

+ (qi yi + qi1yi1)+

(ri + ri1)

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

yi1 = pi y(xi δ ) pi1y(xi1 δ ) 2 piy(xi δ )

((1 - em2 )em1x + (em1 1)em2x )

h

pi1y(xi1 δ ) +

2

h

2 (qi yi ) +

h

2 (qi1yi1) +

h

(ri + ri1)

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

yi1 2

i1

i1

i1

i1

δ δ

(10)

Example 2.

h

2(ε δ )

εy(x) + 0.25 y(x δ ) y(x) = 0

under the interval

= 1 +  yi1 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 )

ε ε ε

i1

i 1

where m1

= (1

1 + 4(ε + δ ) )

2(ε + δ )

and

+

δ

p y

δ

y

h

p +

δ y

δ y

m2 = (1 +

1 + 4(ε + δ ) )

2(ε + δ ) .

i1 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 i1 h i

i 1

2

h i h

i+1

Example 4. εy(x) y(x δ ) + y(x) = 0 under the interval with

pi11 +

yi1

yi  +

h h In this paper, we have described a numerical integration

+ ri + ri1

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 yi1 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

21

22

23

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

24 1.9899e-002 1.0335e-002 6.9813e-003 5.2710e-003 4.2338e-003

0.4

25

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

21

6.5308e-004 3.2800e-004 2.1900e-004 1.6437e-004 1.3156e-004

22 1.2766e-003 6.4293e-004 4.2966e-004 3.2264e-004 2.5830e-004

23 2.3695e-003 1.1981e-003 8.0175e-004 6.0245e-004 4.8251e-004 1

δ=0.3ε

24

25

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

21 3.7887e-003 1.9299e-003 1.2948e-003 9.7419e-004 7.8082e-004

22 1.5347e-002 7.9592e-003 5.3719e-003 4.0537e-003 3.2550e-003

0.7

numerical solution

0.6

0.5

23

24

25

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

21

22

23

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

24

25

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