Numerical Integration of a Class of Singularly Perturbed Delay Differential Equations with Small Shift

Abstract

We have presented a numerical integration method to solve a class of singularly perturbed delay differential equations with small shift. First, we have replaced the second-order singularly perturbed delay differential equation by an asymptotically equivalent first-order delay differential equation. Then, Simpson’s rule and linear interpolation are employed to get the three-term recurrence relation which is solved easily by discrete invariant imbedding algorithm. The method is demonstrated by implementing it on several linear and nonlinear model examples by taking various values for the delay parameter 𝛿 and the perturbation parameter 𝜀.

1. Introduction

The singularly perturbed delay differential equations with small shift arise very frequently in the modeling of various physical and biological phenomena, for example, micro scale heat transfer [1], hydrodynamics of liquid helium [2], second-sound theory [3], thermoelasticity [4], diffusion in polymers [5], reaction-diffusion equations [6], stability [7], control of chaotic systems [8], a variety of models for physiological processes or diseases [9] and so forth. Hence in the recent times, many researchers have been trying to develop numerical methods for solving these problems. Amiraliyev and Cimen [10] presented numerical method comprising a fitted difference scheme on a uniform mesh to solve second-order delay differential equations. Lange and Miura [11, 12] gave an asymptotic approach for a class of boundary-value problems for linear second-order differential-difference equations. Kadalbajoo and Sharma [13–15] presented numerical approaches to solve singularly perturbed differential-difference equations, which contains negative shift in the convention term (i.e., in the derivative term). Lange and Miura [16] considered the boundary value problem for a singularly perturbed nonlinear differential difference equation with shift and discussed the existence and uniqueness of their solutions. Furthermore, Kadalbajoo and Sharma [17] have discussed the numerical solution of the singularly perturbed nonlinear differential equations with small negative shifts.

In this paper, we have presented a numerical integration method for solving a class of singularly perturbed delay differential equations with small shift. First, the second-order singularly perturbed delay differential equation is replaced by an asymptotically equivalent first-order delay differential equation. Then we employed Simpson’s rule and linear interpolation to get three-term recurrence relation which is solved easily by discrete invariant imbedding algorithm. The method is demonstrated by implementing it on several linear and nonlinear model examples by taking various values for the delay and perturbation parameters.

2. Description of the Method

Consider a class of singularly perturbed boundary value problems of the following form: 𝐿𝑦≡𝜀𝑦^′′(𝑥)+𝑎(𝑥)𝑦^′(𝑥−𝛿)+𝑏(𝑥)𝑦(𝑥)=𝑓(𝑥),0≤𝑥≤1, (2.1) with the interval and boundary conditions𝑦(0)𝑦(1)=𝛼,−𝛿≤𝑥≤0,=𝛽, (2.2a)(2.2b)

where 𝜀 is small parameter, 0<𝜀≪1, and 𝛿 is also a small shifting parameter, 0<𝛿≪1; 𝑏(𝑥), and𝑓(𝑥) are bounded continuous functions in (0,1), and 𝛼,𝛽 are finite constants. Further, we assume that 𝑎(𝑥)≥𝑀>0 throughout the interval [0,1], where 𝑀 is positive constant. This assumption merely implies that the boundary layer will be in the neighborhood of 𝑥=0.

By using Taylor series expansion in the neighborhood of the point 𝑥, we have 𝑦𝑥−√𝜀=𝑦(𝑥)−√𝜀𝑦^′(𝑥)+𝜀2𝑦^′′(𝑥) (2.3) and consequently, (2.1) is replaced by the following first-order differential equation: 𝑦^′(𝑥)=𝑝(𝑥)𝑦^′(𝑥−𝛿)+𝑞(𝑥)𝑦𝑥−√𝜀+𝑟(𝑥)𝑦(𝑥)+𝑠(𝑥), (2.4) where 𝑝(𝑥)=−𝑎(𝑥)2√𝜀,𝑞(𝑥)=−1√𝜀(𝑥)=(𝑥),𝑟2−𝑏2√𝜀(𝑥)=𝑓(𝑥),𝑠2√𝜀. (2.5)

The transition from (2.1) to (2.4) is admitted, because of the condition that 𝜀 is small, 0<𝜀≪1. This replacement is significant from the computational point of view. Further details on the validity of this transition can be found in [18].

Now we divide the interval [0,1] into 𝑁 equal subintervals of mesh size ℎ=1/𝑁 so that  𝑥_𝑖=𝑖ℎ,  𝑖=0,1,2,…,𝑁.

Integrating (2.4) with respect to 𝑥 from 𝑥_𝑖 to 𝑥_𝑖+1 for =1,2,…,𝑁−1, we get 𝑦_𝑖+1−𝑦_𝑖=𝑝_𝑖+1𝑦𝑥_𝑖+1−𝛿−𝑝_𝑖𝑦𝑥_𝑖+−𝛿_{𝑥𝑖+1𝑥𝑖}−𝑝^′(𝑥)𝑦(𝑥−𝛿)+𝑞(𝑥)𝑦𝑥−√𝜀+𝑟(𝑥)𝑦(𝑥)+𝑠(𝑥)𝑑𝑥, (2.6) where 𝑦_𝑖=𝑦(𝑥_𝑖),  𝑝_𝑖=𝑝(𝑥_𝑖),  𝑞_𝑖=𝑞(𝑥_𝑖),  𝑟_𝑖=𝑟(𝑥_𝑖),  𝑠_𝑖=𝑠(𝑥_𝑖).

By using Simpson’s rule to evaluate the integral in (2.6), we get 𝑦_𝑖+1−𝑦_𝑖=𝑝_𝑖+1𝑦𝑥_𝑖+1−𝛿−𝑝_𝑖𝑦𝑥_𝑖−ℎ−𝛿6𝑝^′_𝑖𝑦𝑥_𝑖−𝛿+4𝑝^′_𝑖+1/2𝑦𝑥_𝑖+1/2−𝛿+𝑝^′_𝑖+1𝑦𝑥_𝑖+1+ℎ−𝛿6𝑞_𝑖𝑦𝑥_𝑖−√𝜀+4𝑞_𝑖+1/2𝑦𝑥_𝑖+1/2−√𝜀+𝑞_𝑖+1𝑦𝑥_𝑖+1−√𝜀+ℎ6𝑟_𝑖𝑦_𝑖+4𝑟_𝑖+1/2𝑦_𝑖+1/2+𝑟_𝑖+1𝑦_𝑖+1+ℎ6𝑠_𝑖+4𝑠_𝑖+1/2+𝑠_𝑖+1. (2.7) By the means of Taylor series expansion and then by approximating 𝑦^′(𝑥) by linear interpolation, we get𝑦𝑥_𝑖−𝛿=𝑦𝑥_𝑖−𝛿𝑦^′𝑥_𝑖=𝑦_𝑖𝑦−𝛿^𝑖−𝑦_𝑖−1ℎ=𝛿1−ℎ𝑦_𝑖+𝛿ℎ𝑦_𝑖−1𝑦𝑥_𝑖+1,−𝛿=𝑦𝑥_𝑖+1−𝛿𝑦^′𝑥_𝑖+1=𝑦_𝑖+1𝑦−𝛿^𝑖+1−𝑦_𝑖ℎ=𝛿1−ℎ𝑦_𝑖+1+𝛿ℎ𝑦_𝑖𝑦𝑥_𝑖,−√𝜀=𝑦𝑥_𝑖−√𝜀𝑦^′𝑥_𝑖=𝑦_𝑖−√𝜀𝑦_𝑖−𝑦_𝑖−1ℎ=√1−𝜀ℎ𝑦_𝑖+√𝜀ℎ𝑦_𝑖−1,𝑦𝑥_𝑖+1−√𝜀=𝑦𝑥_𝑖+1−√𝜀𝑦^′𝑥_𝑖+1=𝑦_𝑖+1−√𝜀𝑦_𝑖+1−𝑦_𝑖ℎ=√1−𝜀ℎ𝑦_𝑖+1+√𝜀ℎ𝑦_𝑖. (2.8a)(2.8b)(2.8c)(2.8d) In similar way, 𝑦𝑥_𝑖+1/2−𝛿=𝑦𝑥_𝑖+1/2−𝛿𝑦^′𝑥_𝑖+1/2=𝑦_𝑖+1/2𝑦−𝛿^𝑖+1−𝑦_𝑖ℎ=𝑦_𝑖+1/2−𝛿ℎ𝑦_𝑖+1+𝛿ℎ𝑦_𝑖. (2.8e)

Hence, by making use of (2.8a)–(2.8e) in (2.7) we obtain 𝑦_𝑖+1−𝑦_𝑖=−𝛿ℎ𝑝_𝑖+ℎ6𝑝^′_𝑖+√𝜀6𝑞_𝑖𝑦_𝑖−1+𝛿ℎ𝑝_𝑖+1−ℎ6𝑝^′_𝑖+1−𝛿1−ℎ𝑝_𝑖+ℎ6𝑝^′_𝑖−4𝛿6𝑝^′_𝑖+1/2+ℎ6√1−𝜀ℎ𝑞_𝑖+4√𝜀6𝑞_𝑖+1/2+√𝜀6𝑞_𝑖+1+ℎ6𝑟_𝑖𝑦_𝑖+𝛿1−ℎ𝑝_𝑖+1−ℎ6𝑝^′_𝑖+1+4𝛿6𝑝^′_𝑖+1/2+ℎ6√1−𝜀ℎ𝑞_𝑖+1−4√𝜀6𝑞_𝑖+1/2+ℎ6𝑟_𝑖+1𝑦_𝑖+1+4ℎ6−𝑝^′_𝑖+1/2+𝑞_𝑖+1/2+𝑟_𝑖+1/2𝑦_𝑖+1/2+ℎ6𝑠_𝑖+4𝑠_𝑖+1/2+𝑠_𝑖+1. (2.9) To make (2.9) a three-term recurrence relation, we can express 𝑦_𝑖+1/2 in terms of 𝑦_𝑖−1, 𝑦_𝑖 and 𝑦_𝑖+1 using Hermite’s interpolation as follows: 𝑦_𝑖+1/2=12𝑦_𝑖+𝑦_𝑖+1+ℎ8𝑦^′_𝑖−𝑦^′_𝑖+1+𝑂ℎ⁴. (2.10) In view of (2.4) and (2.10), we get 𝑦_𝑖+1/2=12𝑦_𝑖+𝑦_𝑖+1+ℎ8𝑝_𝑖𝑦^′𝑥_𝑖−𝛿+𝑞_𝑖𝑦𝑥_𝑖−√𝜀+𝑟_𝑖𝑦_𝑖+𝑠_𝑖−ℎ8𝑝_𝑖+1𝑦^′𝑥_𝑖+1−𝛿+𝑞_𝑖+1𝑦𝑥_𝑖+1−√𝜀+𝑟_𝑖+1𝑦_𝑖+1+𝑠_𝑖+1. (2.11) By making use of (2.8a)–(2.8e) in (2.11) and finite difference approximations, we get 𝑦_𝑖+1/2=𝛿𝑝8ℎ^𝑖+1−𝑝_𝑖+√𝜀8𝑞_𝑖𝑦_𝑖−1+12+𝛿𝑝8ℎ^𝑖−𝑝_𝑖+1−18𝛿1−ℎ𝑝_𝑖−𝑝_𝑖+1+ℎ8√1−𝜀ℎ𝑞_𝑖+ℎ8𝑟_𝑖−√𝜀8𝑞_𝑖+1𝑦_𝑖+12+18𝛿1−ℎ𝑝_𝑖−𝑝_𝑖+1−ℎ8√1−𝜀ℎ𝑞_𝑖+1−ℎ8𝑟_𝑖+1𝑦_𝑖+1+ℎ8𝑠_𝑖−𝑠_𝑖+1. (2.12) Finally, making use of (2.12) in (2.9) and rearranging as three-term recurrence relation, we get 𝐸_𝑖𝑦_𝑖−1−𝐹_𝑖𝑦_𝑖+𝐺_𝑖𝑦_𝑖+1=𝐻_𝑖, (2.13) for 𝑖=1,2,…,𝑁−1, where 𝐸_𝑖=𝛿ℎ𝑝_𝑖+ℎ6𝑝^′_𝑖−√𝜀6𝑞_𝑖−4ℎ6−𝑝^′_𝑖+1/2+𝑞_𝑖+1/2+𝑟_𝑖+1/2−𝛿𝑝8ℎ^𝑖−𝑝_𝑖+1+√𝜀8𝑞_𝑖,𝐹_𝑖𝛿=1+ℎ𝑝_𝑖+1−ℎ6𝑝^′_𝑖+1−𝛿1−ℎ𝑝_𝑖+ℎ6𝑝^′_𝑖−4𝛿6𝑝^′_𝑖+1/2+ℎ6√1−𝜀ℎ𝑞_𝑖+4√𝜀6𝑞_𝑖+1/2+√𝜀6𝑞_𝑖+1+ℎ6𝑟_𝑖+4ℎ6−𝑝^′_𝑖+1/2+𝑞_𝑖+1/2+𝑟_𝑖+1/2×12+𝛿𝑝8ℎ^𝑖−𝑝_𝑖+1−18𝛿1−ℎ𝑝_𝑖−𝑝_𝑖+1+ℎ8√1−𝜀ℎ𝑞_𝑖+ℎ8𝑟_𝑖−√𝜀8𝑞_𝑖+1,𝐺_𝑖𝛿=1−1−ℎ𝑝_𝑖+1−ℎ6𝑝^′_𝑖+1−4𝛿6𝑝^′_𝑖+1/2−ℎ6√1−𝜀ℎ𝑞_𝑖+1+4√𝜀6𝑞_𝑖+1/2−ℎ6𝑟_𝑖+1−4ℎ6−𝑝^′_𝑖+1/2+𝑞_𝑖+1/2+𝑟_𝑖+1/2×12+18𝛿1−ℎ𝑝_𝑖−𝑝_𝑖+1−ℎ8√1−𝜀ℎ𝑞_𝑖+1−ℎ8𝑟_𝑖+1,𝐻_𝑖=ℎ6𝑠_𝑖+4𝑠_𝑖+1/2+𝑠_𝑖+1+4ℎ6−𝑝^′_𝑖+1/2+𝑞_𝑖+1/2+𝑟_𝑖+1/2ℎ8𝑠_𝑖−𝑠_𝑖+1. (2.14) This tridiagonal system is solved by using method of discrete invariant imbedding algorithm which is described in the next section.

3. Discrete Invariant Imbedding Algorithm

We now describe the Thomas algorithm which is also called discrete invariant imbedding [19] to solve the three-term recurrence relation: 𝐸_𝑖𝑦_𝑖−1−𝐹_𝑖𝑦_𝑖+𝐺_𝑖𝑦_𝑖+1=𝐻_𝑖,for𝑖=1,2…,𝑁−1. (3.1) Let us set a difference relation of the form 𝑦_𝑖=𝑊_𝑖𝑦_𝑖+1+𝑇_𝑖for𝑖=𝑁−1,𝑁−2,…,2,1, (3.2) where 𝑊_𝑖=𝑊(𝑥_𝑖) and 𝑇_𝑖=𝑇(𝑥_𝑖) are to be determined.

From (3.2), we have 𝑦_𝑖−1=𝑊_𝑖−1𝑦_𝑖+𝑇_𝑖−1. (3.3) Substituting (3.3) in (3.1), we have 𝑦_𝑖=𝐺_𝑖𝐹_𝑖−𝐸_𝑖𝑊_𝑖−1𝑦_𝑖+1+𝐸_𝑖𝑇_𝑖−1−𝐻_𝑖𝐹_𝑖−𝐸_𝑖𝑊_𝑖−1. (3.4) By comparing (3.2) and (3.4), we get the recurrence relations 𝑊_𝑖=𝐺_𝑖𝐹_𝑖−𝐸_𝑖𝑊_𝑖−1,𝑇_𝑖=𝐸_𝑖𝑇_𝑖−1−𝐻_𝑖𝐹_𝑖−𝐸_𝑖𝑊_𝑖−1. (3.5)(3.6) To solve these recurrence relations for 𝑖=1,2,3,…,𝑁−1, we need the initial conditions for 𝑊₀ and 𝑇₀. If we choose 𝑊₀=0, then we get 𝑇₀=𝛼. With these initial values, we compute 𝑊_𝑖 and 𝑇_𝑖 for 𝑖=1,2,3,…,𝑁−1 from (3.5) and (3.6) in forward process and then obtain 𝑦_𝑖 in the backward process from (3.2).

The conditions for the discrete invariant imbedding algorithm to be stable are (see [18–21]) 𝐸_𝑖>0,𝐺_𝑖>0,𝐹_𝑖≥𝐸_𝑖+𝐺_𝑖,𝐸_𝑖≤𝐺_𝑖. (3.7) In our method, one can easily show that if the assumptions 𝑎(𝑥)>0, 𝑏(𝑥)<0 and (𝜀−𝛿𝑎(𝑥))>0 hold, then the above conditions (3.7) hold, and thus the discrete invariant imbedding algorithm is stable.

4. Numerical Experiments

To demonstrate the applicability of the method, we have implemented it on two linear and two nonlinear problems with left-end boundary layers. Computational results are compared with exact solutions wherever exact solutions are available. When exact solution is not available, we have tested the effect of small delay parameter on solution of the problem for different values of 𝛿 of 𝑜(𝜀).

4.1. Linear Problems

Example 4.1. Consider an example of singularly perturbed delay differential equation with left layer: 𝜀𝑦^′′(𝑥)+𝑦^′(𝑥−𝛿)−𝑦(𝑥)=0;𝑥∈[0,1]with𝑦(0)=1,𝑦(1)=1. (4.1) The exact solution is given by 𝑦(𝑥)=1−𝑒^𝑚2𝑒^𝑚1𝑥+𝑒^𝑚1−1𝑒^𝑚2𝑥𝑒^𝑚1−𝑒^𝑚2, (4.2) where 𝑚₁=−1−1+4(𝜀−𝛿)/2(𝜀−𝛿) and 𝑚₂=−1+1+4(𝜀−𝛿)/2(𝜀−𝛿).
The computational results are presented in Tables 1, 2, 3, and 4 for 𝜀=0.001 and 0.0001 for different values of 𝛿.

Table 1 

Numerical results of Example 4.1 for 𝜀=0.001,𝛿=0.0001, 𝑁=100.

Table 2 

Numerical results of Example 4.1 for 𝜀=0.001,𝛿=0.0008, 𝑁=100.

Table 3 

Numerical results of Example 4.1 for 𝜀=0.0001,𝛿=0.00001, 𝑁=100.

Table 4 

Numerical results of Example 4.1 for 𝜀=0.0001,𝛿=0.00008, 𝑁=100.

Example 4.2. Now we consider an example of variable coefficient singularly perturbed delay differential equation with left layer: 𝜀𝑦^′′(𝑥)+𝑒^−0.5𝑥𝑦^′(𝑥−𝛿)−𝑦(𝑥)=0with𝑦(0)=1,𝑦(1)=1. (4.3) For which the exact solution is not known. This example is considered to show the effect of the small shift on the boundary layer solution.
The computational results are presented in Tables 5 and 6 for 𝜀=0.001 and 0.0001 for different values of 𝛿.

Table 5 

Numerical results of Example 4.2 for 𝜀=0.001, 𝑁=100, and different values of 𝛿.

Table 6 

Numerical results of Example 4.2 for 𝜀=0.0001, 𝑁=100, and different values of 𝛿.

4.2. Nonlinear Problems

Nonlinear problems are linearized by the quasilinearization process. Then we have applied the present method.

Example 4.3. Consider a singularly perturbed nonlinear delay differential equation: 𝜀𝑦^′′(𝑥)+𝑦(𝑥)𝑦^′(𝑥−𝛿)−𝑦(𝑥)=0 (4.4) under the interval and boundary conditions 𝑦(𝑥)=1,−𝛿≤𝑥≤0,𝑦(1)=1. (4.5) The exact solution is not known.
The computational results are presented in Tables 7 and 8 for 𝜀=0.01 for different values of 𝛿.

Table 7 

Numerical results of Example 4.3 for 𝜀=0.001, 𝑁=100, and different values of 𝛿.

Table 8 

Numerical results of Example 4.3 for 𝜀=0.0001, 𝑁=100, and different values of 𝛿.

Example 4.4. Consider an example of singularly perturbed nonlinear delay differential equation: 𝜀𝑦^′′(𝑥)+2𝑦^′(𝑥−𝛿)+𝑒^𝑦(𝑥)=0 (4.6) under the interval and boundary conditions 𝑦(𝑥)=0,−𝛿≤𝑥≤0,𝑦(1)=0. (4.7) The exact solution is not known.
The computational results are presented in Tables 9 and 10 for 𝜀=0.01 and 0.001 for different values of  𝛿.

Table 9 

Numerical results of Example 4.4 for 𝜀=0.001, 𝑁=100, and different values of 𝛿.

Table 10 

Numerical results of Example 4.4 for 𝜀=0.0001, 𝑁=100, and different values of 𝛿.

5. Discussions and Conclusions

We have presented a numerical integration method to solve singularly perturbed delay differential equations. 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 𝜏=𝑂(1) and 𝜏 is not too large Lange and Miura [12]. To demonstrate the efficiency of the method, we have implemented it on two linear and two nonlinear model examples with the boundary layer on the left for different values of 𝜀 and 𝛿. From the computational results, it is observed that the proposed method approximates the exact solution very well (see Tables 1–4), and the small shift, 𝛿, affects the boundary layer solutions. That is, as 𝛿 increases, the size/thickness of the left boundary layer decreases (see Tables 5–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.