1. Introduction

Delay differential equations play an important role in the mathematical modeling of various practical phenomena in the biosciences and control theory. Any system involving a feedback control will almost always involve time delays. These arise because a finite time is required to sense information and then react to it. A delay differential equation is of the retarded type if the delay argument does not occur in the highest order derivative term. If we restrict this class to a class in which the highest derivative term is multiplied by a small parameter, then it is said to be singularly perturbed delay differential equation of the retarded type. These problems depend on a small positive parameter in such a way that the solution varies rapidly in some parts of the domain and varies slowly in some other parts of the domain. So, typically there are thin transition layers where the solution varies rapidly or jumps abruptly, while away from the layers the solution behaves regularly and varies slowly. In recent years, there has been a growing interest in the numerical study of singularly perturbed delay differential equations because of their applications in many scientific and technical fields like micro scale heat transfer, hydrodynamics of liquid helium, second-sound theory, thermo elasticity, diffusion in polymers, reaction-diffusion equations, stability, control of chaotic systems, a variety of models for physiological processes, Gemechis and Reddy ^[1].

Lange and Miura ^[2] gave an asymptotic approach for a class of boundary-value problems for linear second-order singularly perturbed differential-difference equations. Ramesh and Kadalbajoo ^[3] presented the numerical approximation of singularly perturbed linear second order reaction-diffusion boundary value problems with a small shift (δ) in the reaction term. Swamy ^[4] presented the quantitative analysis of delay reaction-diffusion equations with layer or oscillatory behaviour by employing the numerical integration. Soujanya and Reddy ^[5] presented a computational technique for solving singularly perturbed delay differential equations with layer or oscillatory behaviour in which the small delay is in the reaction term. The treatment of singularly perturbed problems presents severe difficulties that have to be addressed to ensure accurate numerical solutions, Doolan et al. ^[6], Kadalbajoo and Reddy ^[7] and Roos et al. ^[8]. Kadalbajoo and Ramesh ^[9] states that, the accuracy of the problem increased by increasing the resolution of the grid which might be impractical in some cases like higher dimensions. Pratima and Sharma ^[10] states that, till date ε-uniformly convergent methods have not been sufficiently developed for a wide class of singularly perturbed delay differential equations. In this paper, we present numerical method for solving singularly perturbed delay reaction-diffusion equations with layer or oscillatory behaviour via fourth order finite difference method which is uniformly convergent and more accurate than the others.

2. Description of the Method

Consider singularly perturbed delay reaction-diffusion equation of the form:

(1)

with the interval and boundary conditions,

(2)

where is small parameter, and is also small delay parameter, ; and are bounded smooth functions in and is a given constant. The layer or oscillatory behaviour of the problem under consideration is maintained for but sufficiently small, depending on the sign of , for all . If , the solution of the problem in Eqs. (1) and (2) exhibits layer behaviour, and if , it exhibits oscillatory behaviour. Therefore, if the solution exhibits layer behaviour, there will be two boundary layers which will occur at both the end points and.

The solution should be continuous on , continuously differentiable on and also satisfies Eqs. (1) and (2).

By using Taylor series expansion in the neighborhood of the point, we have:

(3)

Substituting Eq. (3) into Eq. (1), we obtain an asymptotically equivalent singularly perturbed two point boundary value problem of the form:

(4)

under the boundary conditions,

(5)

where, , and.

The transition from Eq. (1) to Eq. (4) is admitted, because of the condition that is sufficiently small. Further details on the validity of this transition can be found in Elsgolt’s and Norkin ^[11].

Now, divide the interval into equal parts with constant mesh length. Let be the mesh points. Then, we have .

By using Taylor series expansion, we obtain:

(6)

(7)

Subtracting Eq. (6) from Eq. (7), we obtain:

(8)

Hence, the second order finite difference approximation for the first derivative of is:

(9)

where,

Similarly, adding Eqs. (6) and (7), we obtain:

(10)

Hence, the second order finite difference approximation for the second derivative of is:

(11)

where,

Substituting Eqs. (6) and (7) into Eq. (9) yields:

(12)

where,

Again, substituting Eqs. (6) and (7) into Eq. (11), we obtain:

(13)

where,

Applying to in Eq. (9), we obtain:

(14)

Substituting Eq. (14) into Eq. (12), we obtain:

(15)

where,

Applying to in Eq. (11), we obtain a fourth order finite difference scheme for Eq. (4) as:

(16)

Substituting Eq. (16) into Eq. (13), we obtain:

(17)

where,

From Eqs. (15) and (17), we have:

(18)

Evaluating Eq. (4) at and using Eq. (18), we obtain:

(19)

Simplifying Eq.(19), we get:

(20)

where,

By successively differentiating both sides of Eq. (4), evaluating at , and using into Eq. (20), we obtain:

(21)

Substituting Eqs. (9) and (11) into Eq. (21) for and and making use of and , we obtain:

(22)

where, is the local truncation error.

Multiplying both sides of Eq. (22) by , we get the three-term recurrence relation of the form:

(23)

where,

The tri-diagonal system in Eq. (23) can be easily solved by the method of Discrete Invariant Imbedding Algorithm.

Since Eq. (23) holds for , we have linear equations in the unknowns . The matrix of this set of linear equations is denoted as .

3. Stability and Convergence Analysis

Lemma 1: For all and sufficiently small , the matrix is an irreducible and diagonally dominant matrix.

Proof:

Clearly, is a tridiagonal matrix. is irreducible if its co-diagonals contain non-zero elements only. The co-diagonal contains . It is easily seen that, for sufficiently small ,

Hence, is irreducible.

Again one can observe that, and in each row of , the sum of the two off-diagonal elements less than or equal to the modulus of the diagonal element. This proves the diagonal dominant of .

Under these conditions the discrete imbedding algorithm is stable, Kadalbajoo and Reddy ^[12].

Lemma 2: Let be the analytical solution of the problem in Eqs. (4) and (5) and be the numerical solution of the discretized problem of Eq. (23). Then, for sufficiently small and is positive constant.

Proof:

Multiplying both sides of Eq. (23) by , we get:

(24)

where,

is a local truncation error, for .

Incorporating the boundary conditions in Eq. (24), we get the systems of equations of the form:

(25)

where,

are tri-diagonal matrices of order, and

are the associated vectors of Eq. (25).

Let be the solution which satisfies the Eq. (25), we have:

(26)

Let be the discretization error, then,

Subtracting Eq. (25) from Eq. (26), we get:

(27)

Let , , , ,

Let be the element of the matrix , then:

For,

For,

Thus, for sufficiently small , we have:

Hence, the matrix is irreducible, Varga ^[13].

Let be the sum of the elements of the row of the matrix , then:

Let

then:

For sufficiently small , is monotone, Varga ^[13] and Young ^[14].

Hence, exists and .

From the error Eq. (27), we have:

(28)

For sufficiently small , we have:

Let be the element of and we define,

(29)

Since , then from the theory of matrices, we have:

Hence,

(30)

(31)

Further,

(32)

Now, from Eqs. (28) (32), we get:

(33)

where,

which is independent of perturbation parameter and mesh size .

This establishes that the method is of fourth order uniformly convergent.

4. Numerical Examples

To demonstrate the applicability of the method, we implemented the method on four numerical examples, two with boundary layers and two with oscillatory behaviour. Since, those examples have no exact solution, so the numerical solutions are computed using double mesh principle. The maximum absolute errors are computed using double-mesh principle given by:

(34)

where is the numerical solution on the mesh at the nodal point and and is the numerical solution on a mesh, obtained by bisecting the original mesh with number of mesh intervals, Doolan et al. ^[6].

Example 1. Consider the singularly perturbed delay reaction-diffusion equation with layer behaviour,

under the interval and boundary conditions

The maximum absolute errors are presented in Table 1 and Table 5 for different values of and The graph of the computed solution for and different values of is also given in Figure 1.

Example 2. Consider the singularly perturbed delay reaction-diffusion equation with layer behaviour,

under the interval and boundary conditions

The maximum absolute errors are presented in Table 2 and Table 6 for different values of and The graph of the computed solution for and different values of is also given in Figure 2.

Example 3. Consider the singularly perturbed delay reaction-diffusion equation with oscillatory behaviour,

under the interval and boundary conditions

The maximum absolute errors are presented in Table 3 for different values of . The graph of the computed solution for and different values of is also given in Figure 3.

Example 4. Consider the singularly perturbed delay reaction-diffusion equation with oscillatory behaviour,

under the interval and boundary conditions

The maximum absolute errors are presented in Table 4 for different values of . The graph of the computed solution for and different values of is also given in Figure 4.

5. Numerical Results

Table 1. The maximum absolute errors of Example 1, for different values of δ with ε=0.1

Download as

PowerPoint Slide

Larger image(png format)

• Tables index

Veiw figure

View current table in a new window

View next table

Table 2. The maximum absolute errors of Example 2, for different values of δ with ε=0.1

Download as

PowerPoint Slide

Larger image(png format)

• Tables index

Veiw figure

View current table in a new window

View previous table

View next table

Table 3. The maximum absolute errors of Example 3, for different values of δ with ε=0.1

Download as

PowerPoint Slide

Larger image(png format)

• Tables index

Veiw figure

View current table in a new window

View previous table

View next table

Table 4. The maximum absolute errors of Example 4, for different values of δ with ε=0.1

Download as

PowerPoint Slide

Larger image(png format)

• Tables index

Veiw figure

View current table in a new window

View previous table

View next table

Table 5. The maximum absolute errors of Example 1, for different values of ε with δ=0.03

Download as

PowerPoint Slide

Larger image(png format)

• Tables index

Veiw figure

View current table in a new window

View previous table

View next table

Table 6. The maximum absolute errors of Example 2, for different values of ε with δ=0.03

Download as

PowerPoint Slide

Larger image(png format)

• Tables index

Veiw figure

View current table in a new window

View previous table

View next table

The Effect of Delay Parameter on the Solution Profile

The following graphs (Figure 1 – Figure 4) show the numerical solutions obtained by the present method for different values of delay parameter .

Download as

• PPT
  PowerPoint Slide
• PNG
  Larger image(png format)

Veiw figureFigures index

•  NEW
  View larger figure in new window

View next figure

Figure 1. The numerical solution of Example 1 with ε=0.1

Download as

• PPT
  PowerPoint Slide
• PNG
  Larger image(png format)

Veiw figureFigures index

•  NEW
  View larger figure in new window
• PREV
  View previous figure

View next figure

Figure 2. The numerical solution of Example 2 with ε=0.1

Download as

• PPT
  PowerPoint Slide
• PNG
  Larger image(png format)

Veiw figureFigures index

•  NEW
  View larger figure in new window
• PREV
  View previous figure

View next figure

Figure 3. The numerical solution of Example 3 with ε=0.01

Download as

• PPT
  PowerPoint Slide
• PNG
  Larger image(png format)

Veiw figureFigures index

•  NEW
  View larger figure in new window
• PREV
  View previous figure

Figure 4. The numerical solution of Example 4 with ε=0.01

The Rate of Convergence

In the same way in Eq. (34) one can define by replacing by and by , that is:

The computational rate of convergence is also obtained by using the double mesh principle defined as, Doolan et al. ^[6]:

Table 7. Rate of Convergence ρ for ε=0.1 and δ=0.05

Download as

PowerPoint Slide

Larger image(png format)

• Tables index

Veiw figure

View current table in a new window

View previous table

6. Discussion and Conclusion

Numerical solution of second order singularly perturbed delay reaction-diffusion equations with boundary layer or oscillatory behaviour via finite difference method has been presented. To demonstrate the efficiency of the method, four model examples without exact solutions have been considered for different values of the perturbation parameter and delay parameter The numerical solutions are tabulated (Table 1 - Table 6) in terms of maximum absolute errors and observed that the present method improves the findings of Soujanya and Reddy ^[5] and Swamy ^[4]. From the results, it can be observed that as the grid size decreases, the maximum absolute errors decrease, which shows the convergence to the computed solution. The stability and convergence of the method are investigated and established well. The results presented in Table 7 confirmed that computational rate of convergence as well as theoretical estimates indicate that method is a fourth order convergent.

Further, to investigate the effect of delay on the solution of the problem, numerical solutions have been presented using graphs. Accordingly, when the order of the coefficient of the delay term is of , the delay affects the boundary layer solution but maintains the layer behaviour (Figure 2). When the delay parameter is of , the solution maintains layer behaviour although the coefficient of the delay term in the equation is of and as the delay increases, the thickness of the left boundary layer decreases while that of the right boundary layer increases (Figure 1).

To demonstrate the effect on the oscillatory behavior, we consider the examples 3 and 4 when the solution of the problem exhibits oscillatory behaviour for delay parameter equal to zero and different from zero. We observe that, if the coefficient of the delay term is of o(1), the amplitude of the oscillations increases slowly as the delay increases provided the delay parameter is of (Figure 3) and if the coefficient of the delay term is of O(1), there is no oscillation in the left half of the interval while the amplitude of the oscillations increases as the delay increases in the right half of the interval provided the delay parameter is of (Figure 4).

Conflict of Interest

The authors declare that there is no conflict of interests regarding the publication of this paper.

Acknowledgments

The authors would like to thank Jimma University for the financial and material support.

References

[1]	Gemechis File and Y.N. Reddy, Numerical integration of a class of singularly perturbed delay differential equations with small shift, International Journal of Differential Equation, 572723 (2012), 1-12.
	In article		View Article
[2]	C. G. Lange and R. M. Miura, Singular perturbation analysis of boundary-value problems for differential-difference equations very small shifts with layer behavior, SIAM J. Appl. Math., 54 (1994), 249-272.
	In article		View Article
[3]	V.P. Ramesh and M.K. Kadalbajoo, Numerical Algorithm for Singularly Perturbed Delay Differential Equations with Layer and Oscillatory Behavior, Neural, Parallel, and Scientific Computations, 19 (2011), 21-34.
	In article
[4]	D. K. Swamy, Quantative analysis of delay differential equations with layer, Advance Research and Innovations in Mechanical, Material Science, Industrial Engineering and Management - ICARMMIEM- (2014), 145-150.
	In article
[5]	G. BSL. Soujanya and Y.N. Reddy, Computational method for singularly perturbed delay differential equations with layer or oscillatory behaviour, Appl. Math. Inf. Sci. 10 (2), (2016), 527-536.
	In article		View Article
[6]	E.R. Doolan, J.J.H. Miller and W.H.A. Schilders, Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, 1980.
	In article
[7]	M.K. Kadalbajoo and Y.N. Reddy, Asymptotic and numerical analysis of singular perturbation problems, Applied Mathematics and Computation, 30 (1989), 223-259.
	In article		View Article
[8]	H.G. Roos, M. Stynes and L. Tobiska, Numerical Methods for Singularly Perturbed Differential Eqns, Convection-Diffusion and Flow Problems, Springer Verlag, Berlin, 1996.
	In article		View Article
[9]	M.K. Kadalbajoo and V.P. Ramesh, Numerical methods on Shishkin mesh for singularly perturbed delay differential equations with a grid adaptation strategy, Applied Mathematics and Computation, 188 (2007), 1816-1831.
	In article		View Article
[10]	R. Pratima and K.K. Sharma, Numerical analysis of singularly perturbed delay differential turning point problem, Applied Mathematics and Computation, 218 (2011), 3483-3498.
	In article		View Article
[11]	L.E. Elsgolt’s and S.B. Norkin, Introduction to the Theory and Applications of Differential Equations with Deviating Arguments, Academic Press, New York, 1973.
	In article
[12]	M. K. Kadalbajoo and Y.N. Reddy, A non-asymptotic method for general singular perturbation problems, Journal of Optimization Theory and Applications, 55 (1986), 256 269.
	In article
[13]	R.S. Varga, Matrix Iterative Analysis. New Jersey: Prentice-Hall, Englewood Cliffs, 1962.
	In article
[14]	D.M. Young, Iterative Solution of Large Linear Systems, Academic Press, Inc. New York, 1971.
	In article