A NON-ASYMPTOTIC METHOD FOR SINGULARLY PERTURBED DELAY DIFFERENTIAL EQUATIONS

Abstract

In this paper, a non-asymptotic method is presented for solving singularly perturbed delay differential equations whose solution exhibits a boundary layer behavior. The second order singularly perturbed delay differential equation is replaced by an asymptotically equivalent first order neutral type delay differential equation. Then, Simpson's integration formula and linear interpolation are employed to get three term recurrence relation which is solved easily by Discrete Invariant Imbedding Algorithm. Some numerical examples are given to validate the computational efficiency of the proposed numerical scheme for various values of the delay and perturbation parameters.

1. Introduction

A singularly perturbed delay differential equation is an ordinary differential equation in which the highest derivative is multiplied by a small parameter and containing delay term. In recent years, Delay differential equations have attracted the attention of many researchers because of their applications in many scientific and technical fields such as in population ecology, control theory, viscoelasticity, and materials with thermal memory, El’sgolts [4], and Bellman and Cooke [3]. Boundary value problems involving delay differential equations are also used to study signal transmission with time delays in control theory El’sgolts [4], first-exit problems in neurobiology Stein [13], and Tuckwell [14]. As a result, many researchers tried to develop and present numerical schemes for solving such problems. For example, Amiraliyev and Cimen [1], presented numerical method comprising a fitted difference scheme on a uniform mesh to solve second order delay differential equations. Lange and Miura [10,11] gave an asymptotic approach for a class of boundary-value problems for linear second-order differential-difference equations. Kadalbajoo and Sharma [7,8,9] presented numerical approaches to solve singularly perturbed differential-difference equations, which contains negative shift in the convention term (i.e. in the derivative term). Furthermore, Phaneendra, et al.[12], presented higher order numerical method to solve the problems of the same type.

The objective of this paper is, therefore, to present a non-asymptotic method for solving singularly perturbed delay differential equations. The method is named as non-asymptotic as it does not depend on asymptotic expansion as well as on the matching of the coefficients. In this method, we first replace the second order singularly perturbed delay differential equation by an asymptotically equivalent first order neutral type delay differential equation and employ Simpson’s rule for the numerical integration of the obtained first order delay differential equation. Then, linear interpolation is used to get the three term recurrence relation which is solved easily by discrete Invariant Imbedding algorithm. The method is demonstrated by implementing several model examples by taking various values for the delay and perturbation parameters.

 

2. Description of the Method

Consider singularly perturbed singular boundary value problems of the form

with the interval and boundary conditions

where ε is small parameter, 0 < ε << 1, and the delay parameter, δ is of o(ε) satisfying (ε − δa(x)) > 0 for all xϵ[0, 1]; a(x), b(x), f(x) and ϕ(x) are smooth functions and β is a finite constant. For δ = 0, the problem (1)-(3) becomes a boundary value problem for singularly perturbed ordinary differential equations. The layer behavior of the problem under consideration is maintained only for δ ≠ 0 but sufficiently small (i.e δ is of o(ε)). When the delay parameter, δ exceeds the perturbation parameter, ε (i.e δ is of O(ε)), the layer behavior of the solution is no longer maintained, rather the solution exhibits an oscillatory behavior or diminished. In this paper, we consider the cases where the boundary layer is either on the left side or right side of the interval [0, 1].

2.1. Layer on the Left Side: In general, the solution of problem (1)-(3) exhibits boundary layer behavior at one end of the interval [0, 1] depending on the sign of a(x). Now, we assume that a(x) ≥ M > 0 throughout the interval [0, 1], where M is positive constant. This assumption merely implies that the boundary layer will be in the neighborhood of x = 0. Further, it is also assumed that b(x) satisfies the condition that b(x) ≤ −θ < 0 , ∀xϵ[0, 1] where θ is a positive constant.

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

By using (4) in (1) we get an asymptotically equivalent first order delay neutral type differential equation:

The transition from equation (1) to equation (5) is admitted, because of the condition that ε is small, 0 < ε << 1; and the truncated term ε2/2 is very small and approaches to 0 as ε is small. Further details on the validity of this transition can also be found in Elsgolt’s and Norkin [5].

Now we divide the interval [0, 1] into N equal subintervals of mesh size h = 1/N so that xi = ih, i = 0, 1, 2, · · · ,N.

Integrating equation (5) with respect to x from xi to xi+1 for i = 1, 2, · · · ,N −1, we get

where yi = y(xi), ai = a(xi), bi = b(xi), fi = f(xi).

By using Simpson’s rule to evaluate the integral, we get:

Since delay term is of the small order of ε, to tackle the term containing delay, we can use Taylor’s approximation (Kadalbajoo and Sharma [9], Lange and Miura [11]). Thus, by means of Taylor series expansion and then by approximating by linear interpolation, we get:

In similar way,

Hence, by making use of equations (8)-(12) in (7) we obtain:

To make equation (13) a three term recurrence relation, we can express yi+1/2 in terms of yi−1, yi and yi+1 using Hermite’s interpolation as follows:

In view of equations(5) and (14), we get:

By making use equations (8)- (12) in (15) and finite difference approximations, we get

Finally, making use of equations (16) in (13) and rearranging as three term recurrence relation, we get:

for i = 1, 2, · · · ,N − 1, where

We solve the above tridiagonal system by using method of invariant imbedding system described in section 3.

2.2. Right-End Boundary Layer Problem: We now assume that a(x) ≤ −M < 0 throughout the interval [0, 1], where M is positive constant. This assumption merely implies that the boundary layer will be in the neighborhood of x = 1. By using Taylor series expansion in the neighborhood of the point x, we have:

and consequently, equation (1) is replaced by the following first order differential equation:

Now we divide the interval [0, 1] into N equal subintervals of mesh size h = 1/N so that xi = ih, i = 0, 1, 2, · · · ,N. Integrating equation (23) with respect to x from xi−1 to xi, we get:

By using Simpson’s rule to evaluate the integral, we get:

Again, by means of Taylor series expansion and then by approximating by linear interpolation, we get:

In similar way,

Hence, by making use of equations (26)-(30) in (25) we obtain:

To make equation (31) a three term recurrence relation, we can express yi−1/2 in terms of yi−1, yi and yi+1 using Hermite’s interpolation as follows:

In view of equations (23) and (32), we get:

By making use of equations (26)- (30) in (33) and finite difference approximations, we get:

Finally, making use of equations (34) in (31) and rearranging as three term recurrence relation, we get:

for i = 1, 2, · · · ,N − 1, where

We solve the above tridiagonal system by using method of invariant imbedding described in the next section.

 

3. Discrete Invariant Imbedding Algorithm

We now describe the Thomas algorithm which is also called Discrete Invariant Imbedding (Angel and Bellman [2]) to solve the three term recurrence relation:

for i = 1, 2, · · · ,N − 1. Let us set a difference relation of the form:

where Wi = W(xi) and Ti = T(xi) are to be determined. From (41), we have:

substituting (42) in (40), we get:

By comparing (43) and (41), we get the recurrence relations:

To solve these recurrence relations for i = 1, 2, 3, · · · ,N −1 , we need the initial conditions for W0 and T0 . If we choose W0 = 0 , then we get T0 = ϕ(0) . With these initial values, we compute Wi and Ti for i = 1, 2, 3, · · · ,N − 1 from equations (44) and (45) in forward process, and then obtain yi in the backward process from (41).

The conditions for the stability of discrete invariant imbedding algorithm are (Angel and Bellman [2], Kadalbajoo and Reddy [6]):

In our method, one can easily show that if the assumptions a(x) > 0, b(x) < 0 and (ε − δa(x)) > 0 hold, then the above conditions (46) hold and thus the invariant imbedding algorithm is stable.

 

4. Numerical Examples

To demonstrate the applicability of the method, we have considered two numerical experiments with left-end boundary layer and two numerical experiments with right-end boundary layer. We compared the computed results with the exact solution of the problems where the exact solution of the problems is known; and we also have tested the effect of small delay parameter on computed solution of the problem for different values of δ of o(ε) which are presented by numerical experiments whose exact solutions are not known.

Example 4.1. Consider an example of singularly perturbed delay differential equation with left layer:

The exact solution is given by

where

The computational results are presented in the Tables 1, 2, 3 and 4 for ε = 0.01 and 0.001 for different values of δ.

Example 4.2. Now we consider an example of variable coefficient singularly perturbed delay differential equation with left layer:

For which the exact solution is not known.This example is also considered to show the effect of the small shift on the boundary layer solution. The computational results are presented in the Tables 5 and 6 for ε = 0.01 and 0.001 for different values of δ.

Example 4.3. Consider a singularly perturbed delay differential equation with right layer:

The exact solution is given by

where

The computational results are presented in the Tables 7 and 8 for ε = 0.01 and 0.001 for different values of δ.

Example 4.4. Now we consider an example of variable coefficient singularly perturbed delay differential equation with right layer:

For which the exact solution is not known. The computational results are presented in the Tables 9 and 10 for ε = 0.01 and 0.001 for different values of δ.

 

5. Discussion and Conclusion

We have presented a numerical integration method to solve singularly perturbed delay differential equations whose solutions exhibits the layer behavior. 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 τ = O(1) and τ is not too large Lange and Miura[10].

To demonstrate the efficiency of the method, some numerical examples of boundary value problems with constant and variable coefficients are considered for different values of ε and δ. Besides solving for the numerical solutions, we have considered numerical results of boundary value problems, examples 4.2 & 4.4 , whose exact solutions are not known to show the effect of delay parameter on the boundary layer solutions.

From the tabulated solutions of the numerical experiments considered, it can be observed that the present method approximates the exact solution very well (Tables 1-4, 7-8) for h ≥ ε for which other classical finite difference methods fails to give good results; and the small shift, δ affects both the boundary layer solutions (left and right) in similar fashion but reversely. That is, as δ increases the thickness of the left boundary layer decreases (Tables 5-6) while that of the right boundary layer increases (Tables 9-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.

TABLE 1.Numerical results of Example 4.1 for ε = 0.01, δ = 0.001, N = 100

TABLE 2.Numerical results of Example 4.1 for ε = 0.01, δ = 0.008, N = 100

TABLE 3.Numerical results of Example 4.1 for ε = 0.001, δ = 0.0001, N = 100

TABLE 4.Numerical results of Example 4.1 for ε = 0.001, δ = 0.0008, N = 100

TABLE 5.Numerical results of Example 4.2 for ε = 0.01, N = 100 and different values of δ.

TABLE 6.Numerical results of Example 4.2 for ε = 0.001, N = 100 and different values of δ.

TABLE 7.Numerical results of Example 4.3 for ε = 0.01, δ = 0.001, N = 100

TABLE 8.Numerical results of Example 4.3 for ε = 0.01, δ = 0.008, N = 100

TABLE 9.Numerical results of Example 4.4 for ε = 0.01, N = 100 and different values of δ.

TABLE 10.Numerical results of Example 4.4 for ε = 0.001, N = 100 and different values of δ.