Abstract
This paper presents a numerical method to solve singularly perturbed differential-difference equations. The solution of this problem exhibits layer or oscillatory behavior depending on the sign of the sum of the coefficients in reaction terms. A fourth-order exponentially fitted numerical scheme on uniform mesh is developed. The stability and convergence of the proposed method have been established. The effect of delay parameter (small shift) on the boundary layer(s) has also been analyzed and depicted in graphs. The applicability of the proposed scheme is validated by implementing it on four model examples. Maximum absolute errors in comparison with the other numerical experiments are tabulated to illustrate the proposed method.
1. Introduction
Differential equations with a small parameter multiplying the highest order derivative are called singularly perturbed differential equations. Mathematically, any ordinary differential equation in which the highest derivative is multiplied by a small positive parameter and containing at least one delay/advance parameter is known as a singularly perturbed differential-difference equation [1]. Such type of equations arises frequently from the mathematical modeling of various practical phenomena, for example, in the modeling of the human pupil-light reflex [2], the study of bistable devices [3], and vibrational problems in control theory [4]. When perturbation parameter is very small, most numerical methods for solving such problems may be unstable and give inaccurate results. So, it is important to develop suitable numerical methods to solve singularly perturbed delay differential equations.
Hence, in the recent times, many researchers have been trying to develop numerical methods for solving singularly perturbed delay differential equations. For example, the authors in [5] proposed a computational method of first order for singularly perturbed delay reaction-diffusion equations with layer or oscillatory behavior. The authors in [6] presented a fourth-order finite difference scheme for second-order singularly perturbed differential-difference equations with negative shift. The authors in [7] presented exponentially fitted a second-order finite difference scheme for a class of singularly perturbed delay differential equations with large delay. In [8], the numerical solution of singularly perturbed differential-difference equations with dual layer was presented. Recently, the authors in [9] presented a computational method for solving a singularly perturbed delay differential equation with twin layers or oscillatory behavior. But, still there is a lack of accuracy because the treatment of singularly perturbed problems is not trivial and the solution depends on perturbation parameter and mesh size [10–12]. Due to this, numerical treatment of singularly perturbed delay differential equations needs improvement. Therefore, it is important to develop a more accurate and convergent numerical method for solving singularly perturbed delay differential equations.
Consider the following singularly perturbed delay differential equation of the form:subject to the interval and boundary conditionswhere is the perturbation parameter, , and is a small delay parameter of Also, are bounded smooth functions and is a given constant. The layer or oscillatory behavior 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 equations (1) and (2) exhibits layer behavior, and if it exhibits oscillatory behavior. Therefore, if the solution exhibits layer behavior, there will be two boundary layers which will occur at both end points and [12].
Thus, the purpose of this study is to develop stable, convergent, and more accurate numerical method for solving singularly perturbed delay differential equations.
2. Description of the Method
By using Taylor series expansion in the neighborhood of we have
Substituting equation (3) into equation (1), we obtain an asymptotically equivalent singularly perturbed two-point boundary value problem:where under boundary conditions
Discretizing the given interval [0, 1] in to equal parts with constant mesh size we have
Using Taylor’s expansions of and up to , we get the finite difference approximations for and aswhere , for andwhere for .
Substituting equations (6) and (7) into equation (4) and simplifying, we obtainwhere is the local truncation error and , and .
By successively differentiating both sides of equation (4) and evaluating at and using into equation (8), we obtainwhere
Now introducing a fitting parameter and using central difference approximation for and in equation (9), we have
Multiplying both sides of equation (11) by and taking the limit as we obtainwhere .
From the theory of singular perturbations and [13], we have two cases for and :
Case 1. For (right-end boundary layer), we haveThus, from equation (12), we get
Case 2. For (left-end boundary layer), we haveThus, from equation (12), we getIn general, for discretization, we take a variable fitting parameter aswhere .
Simplifying equation (11), we get the tridiagonal system of the equation as follows:whereThe tridiagonal system in equation (18) can be easily solved by the Thomas algorithm with help of Matlab 2013.
3. Stability and Convergence Analysis
Case 1: Layer behavior (i.e., for ; thus, since ).
Lemma 1. If and for all , then the solution for all for equations (4) and (5).
Proof. Suppose , such that and . Since and is a point of minima, then and
Therefore, we have since (by assumption) and . But this is a contradiction. Then, it follows that for .
Theorem 1. If the solution of the problem in equations (4) and (5) satisfiesfor some constant , then the solution is stable.
Proof. We define two functions: . Then, we get andsince and for suitable choice of . Therefore, by Lemma 1, we get , for all . So, .
Hence, the stability of the solutions of the problem in equations (4) and (5) is proved for the case of the layer behavior.
Lemma 2. The finite difference operator in equation (13) satisfies the discrete minimum principle; i.e., if is any mesh function such that and for all , then for all .
Proof. Suppose there exists a positive integer k such that . Then, from equation (13), we haveFor sufficiently small and for suitable value of we obtain However, (by the assumption) and . But, this is a contradiction. Hence, .
Theorem 2. The finite difference operator in equation (13) is stable for if is any mesh function, then for some constant
Proof. We define two functions:Then, similar to Theorem 1, we get andsince and therefore, by Lemma 2, we get for all .
Thus, .
This proves the stability of the scheme for the case of layer behavior.
Case 2: oscillatory behavior (i.e., for ; thus, since ).
The continuous maximum principle and stability of the solution of equations (4) and (5) are presented as follows for the case of oscillatory behavior.
Lemma 3. If and , for all then the solution for all for equations (4) and (5).
Proof. Suppose , such that and . Since and is a point of maxima, therefore, and . Therefore, we have since (by assumption) and . But this is a contradiction. Hence, for .
Theorem 3. If the solution of the problem in equations (4) and (5) satisfiesfor some constant,
Proof. The proof is analogous to Theorem 1. Hence, the stability of the solutions of the problem in equations (4) and (5) is proved for the case of oscillatory behavior. Now, we present the stability of the discrete problem in equation (13) for the case of oscillatory behavior.
Lemma 4. The finite difference operator in equation (13) satisfies the discrete maximum principle, if is any mesh function such that and for all , then for all .
Proof. Suppose there exists a positive integer k such that Then, from equation (13), we haveFor sufficiently small and for suitable value of we obtain However, (by the assumption) and But, this is a contradiction. Hence,
Theorem 4. The finite difference operator in equation (13) is stable for if is any mesh function, then for some constant
Proof. The proof is similar to Theorem 2. This proves the scheme for the case of oscillatory behavior.
Definition 1. (uniform convergence). Let be a solution of equations (1) and (2). Consider a difference scheme for solving equations (1) and (2). If the scheme has a numerical solution that satisfies where and are independent of and of mesh size then we say the scheme uniformly convergent to with respect to the norm [14].
Lemma 5. The bound for derivative of the solution of the problem (1)–(3) when is given by
Proof. For the proof, refer [15]
For any mesh function define the following difference operators:
Theorem 5. Let and be, respectively, the analytical solution of equations (4) and (5) and numerical solutions of equation (13). Then, for sufficiently large the following parameter uniform error estimate holds:
Proof. Let us consider the local truncation error defined aswheresince In our assumption,
By considering is fixed and taking the limit for we obtain the following:From Taylor series expansion, the bound for the difference becomeswhereNow using the bounds and the assumption equation (30) reduces toHere, the target is to show the scheme convergence independent on the number of mesh points.
By using the bounds for the derivatives of the right layer solution in Lemma 5, we obtain
Lemma 6. For a fixed mesh and for it holds
Proof. Refer from [16].
By using Lemma 6 into equation (36), it results inHence, by discrete maximum principle, we obtainThus, result of equation (39) shows equation (29). Hence proved.
Remark 1. A similar analysis for convergence may be carried out for the left layer case.
4. Numerical Examples and Results
To demonstrate the applicability of the method, we implemented the method on four numerical examples, two with twin boundary layers and two with oscillatory behavior. Since those examples have no exact solution, the numerical solutions are computed using the double mesh principle. The maximum absolute errors are computed using the double mesh principle given bywhere 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 intervals [9].
Example 1. Consider the singularly perturbed delay reaction-diffusion equation with layer behaviorunder the interval and boundary conditionsThe maximum absolute errors are presented in Tables 1 and 2 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 behaviorunder the interval and boundary conditionsThe maximum absolute errors are presented in Tables 3 and 4 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 behaviorunder the interval and boundary conditionsThe maximum absolute errors are presented in Table 5 for different values of
The graph of the computed solution for and δ = 0.003 is also given in Figure 3.
Example 4. Consider the singularly perturbed delay reaction-diffusion equation with oscillatory behaviorunder the interval and boundary conditionsThe maximum absolute errors are presented in Table 6 for different values of
The graph of the computed solution for and δ = 0.003 is also given in Figure 4.
The rate of convergence (): in the same way, in equation (40), one can define by replacing by and by that is,The computational rate of convergence is also obtained by using double mesh principle defined as follows [9]:The following tables (Tables 7 and 8) show the rate of convergence of the present method for different values of the mesh size
4.1. The Effect of Delay Term on the Solution Profile
To analyze the effect of the delay term on the solution profile of the problem, the numerical solution of the problem for different values of the delay parameters has been given by the following graphs.
5. Discussion and Conclusion
Fourth-order fitted operator numerical method for solving singularly perturbed reaction-diffusion with delay 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 , and also results are presented in the tables and figures. It is observed that, from the tables, the present method approximates the solution and the stability and convergence of the method is established well. The effect of the delay on the solution of singularly perturbed delay reaction-diffusion equation is showed by plotting graphs of four model examples. Two model examples of twin layers behavior and two model examples with oscillatory layers have been considered and solved for different values of perturbation parameter delay parameter , and mesh size The numerical solutions are tabulated (Tables 1–6) in terms of maximum absolute errors and observed that the present method improves the findings in [17]. Also, it is significant that all of the maximum absolute errors decrease rapidly as N increases. The results presented in Tables 7 and 8 confirmed that computational rate of convergence as well as the theoretical estimate indicates that the method is a fourth-order convergent.
Furthermore, 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 behavior (Figure 1). When the delay parameter is of the solution maintains layer behavior although the coefficient of the delay term in the equation of and the delay increases, the thickness of the left boundary layer decreases while that of the right boundary layer increases (Figure 2). For the oscillatory behavior case, one can conclude that the solution oscillates throughout the domain for different values of delay parameter (Figures 3 and 5). In a concise manner, the present method gives more accurate solution and is uniformly convergent for solving singularly perturbed delay reaction-diffusion equations with twin layer and oscillatory behavior. Also it can see that, as mesh size decreases, the absolute errors also decrease from Figures 5–8.
Data Availability
No data were used to support the study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Authors’ Contributions
HG proposed the main idea of this paper. HG, SB, and AD prepared the manuscript and performed all the steps of the proofs in this research. All authors contributed equally and significantly in writing this paper. All authors read and approved the final manuscript.
Acknowledgments
The authors wish to express their thanks to Jimma University, College of Natural Sciences, for technical support and the authors of literatures for the provided scientific aspects and idea for this work.