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-

      image

      ε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

      image

      image

      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

      image

      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

        image

        image

        image

        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

        image

        image

        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

        image

        image

        h

        2

        i i

        image

        h

        i1

        xi 4ε

        − +

        δ δ δ

        h δ h

        h pi+1 h

        pi 1 h

        pi+1

        2

        pi 1

        2

         + qi yi

        h 2


        2ε

        image

        + − p


        δ

          h

        image

        image

        2

          + p

         


        δ

          h

        image

        image

          − q y


        image

        image

        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

      image

      image

      Ei = pi h + 2

      4ε

      image

      pi′ +

      h

      δ


      δ δ


      h δ h

      By using the Trapezoidal rule to evaluate the integral

      Fi =

      image

      image

      image

      image

      image

      image

      image

      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

      image

      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

         


      image

      image

      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

      image

      h


      (6)


      Now for the right - end boundary layer, integrating Eq. (4) with respect to x from xi1 to xi , we get

      image

      image

      = 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 δ

      image

      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

      image

      i h


      ε

      image

      h

      i1 i

      image

      i1

      i h

      image

      i h

      4ε h

      image

      image

      i1

      h

      δ

      2 i1

      δ

      image

      image

      h 2 i1

      δ δ h

      + yi1

      image

      ( 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

      image

      image

      Hi = ri + ri1

      2 2

      image

      2ε ( y y )

      h i+1 i


      image

      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

image

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


h h

left and right end boundary layers.

image

+ (qi yi + qi1yi1)+

image

(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

image

image

h yi+1 h

image

yi + h

image

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


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

h

image

pi1y(xi1 δ ) +

2

h

image

2 (qi yi ) +

h

image

2 (qi1yi1) +

image

h

(ri + ri1)

2

y(x) =


image

(em1 em2 )

(9)

Again by means of Taylor series expansion and then corre-


where

m1 = (1

image

1 + 4(ε δ ) )


2(ε δ )


and

sponding

y(x) by linear interpolation, we have


m = (1 +


image

1 + 4(ε δ ) ) .

y(x

1 + 1)e 2 )

image

  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

image

(e 2 e 1 )


ε ε ε


i1


i 1


where m1

image

= (1

1 + 4(ε + δ ) )


2(ε + δ )


and


+

δ

image

p y


δ

image

y


h

image

p +


image

δ y


image

δ 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 )

i1 i1

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

image

image

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

image

ε / N


image

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

image

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

image

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


image

-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


image

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