Exponentially Fitted Numerical Method for Singularly Perturbed Differential-Difference Equations

*is paper presents a numerical method to solve singularly perturbed differential-difference equations. *e solution of this problem exhibits layer or oscillatory behavior depending on the sign of the sum of the coefficients in reaction terms. A fourthorder exponentially fitted numerical scheme on uniform mesh is developed. *e stability and convergence of the proposed method have been established. *e effect of delay parameter (small shift) on the boundary layer(s) has also been analyzed and depicted in graphs. *e 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.


Introduction
Differential equations with a small parameter ε(0 < ε ≤ 1) 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. e authors in [6] presented a fourth-order finite difference scheme for second-order singularly perturbed differential-difference equations with negative shift. e 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 h [10][11][12]. Due to this, numerical treatment of singularly perturbed delay differential equations needs improvement. erefore, 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: y(x) � ϕ(x), for − δ ≤ x ≤ 0 and y(1) � φ, (2) where ε is the perturbation parameter, 0 < ε ≪ 1, and δ is a small delay parameter of O(ε), 0 < δ ≪ 1. Also, a(x), b(x), f(x), and ϕ(x) are bounded smooth functions and φ is a given constant. e layer or oscillatory behavior of the problem under consideration is maintained for δ ≠ 0 but sufficiently small, depending on the sign of a(x) + b(x), for all x ∈ (0, 1). If a(x) + b(x) < 0, the solution of the problem in equations (1) and (2) exhibits layer behavior, and if a(x) + b(x) < 0, it exhibits oscillatory behavior. erefore, if the solution exhibits layer behavior, there will be two boundary layers which will occur at both end points x � 0 and x � 1 [12].
us, the purpose of this study is to develop stable, convergent, and more accurate numerical method for solving singularly perturbed delay differential equations.

Description of the Method
By using Taylor series expansion in the neighborhood of x, we have Substituting equation (3) into equation (1), we obtain an asymptotically equivalent singularly perturbed two-point boundary value problem: where Discretizing the given interval [0,1] in to N equal parts with constant mesh size h, we have Using Taylor's expansions of y i+1 and y i− 1 up to O(h 5 ), we get the finite difference approximations for y i ′ and y i ″ as where where τ 2 � − (h 4 /360)y (6) Substituting equations (6) and (7) into equation (4) and simplifying, we obtain where τ i � h 4 ((p i /120)y (5) (ξ 1 ) + (ε/360)y (6) is the local truncation error and p( By successively differentiating both sides of equation (4) and evaluating at x i and using into equation (8), we obtain where Now introducing a fitting parameter σ and using central difference approximation for y i ″ and y i ′ in equation (9), we have Multiplying both sides of equation (11) by h and taking the limit as h ⟶ 0, we obtain where ρ � h/ε. From the theory of singular perturbations and [13], we have two cases for p(x) > 0 and p(x) < 0: us, from equation (12), we get Case 2. For p(x) > 0 (left-end boundary layer), we have us, from equation (12), we get In general, for discretization, we take a variable fitting parameter as where ρ i � h/ε. Simplifying equation (11), we get the tridiagonal system of the equation as follows: where

Theorem 1. If the solution of the problem in equations
for some constant K ≥ 1, then the solution is stable.
Proof. We define two functions: since q(x) < 0 and for suitable choice of K. erefore, by 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.
e finite difference operator L N in equation (13) satisfies the discrete minimum principle; i.e., if w i is any mesh function such that w 0 ≥ 0 and L N w i ≤ 0, for all x i ∈ (0, 1), Proof. Suppose there exists a positive integer k such that w k < 0 and w i � min 0≤i≤N w i . en, from equation (13), we have For sufficiently small h and for suitable value of p k , we obtain L N w k > 0. However, w k < 0 (by the assumption) and B k ⟶ p k < 0. But, this is a contradiction. Hence, Proof. We define two functions: en, similar to eorem 1, we get ψ ± (0) ≥ 0 and is proves the stability of the scheme for the case of layer behavior.
Proof. e proof is analogous to eorem 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.
e finite difference operator L N in equation (13) satisfies the discrete maximum principle, if w i is any mesh function such that w 0 ≥ 0 and L N w i ≥ 0, for all x i ∈ (0, 1), Proof. Suppose there exists a positive integer k such that w k < 0 and w i � max 0≤i≤N w i . en, from equation (13), we have For sufficiently small h and for suitable value of p k , we obtain L N w k < 0. However, w k < 0(by the assumption) and B k ⟶ p k > 0. But, this is a contradiction. Hence, w i ≥ 0 for all x i ∈ (0, 1).

Proof.
e proof is similar to eorem 2. is proves the scheme for the case of oscillatory behavior.  (1) and (2). If the scheme has a numerical solution y N that satisfies ‖y − y N ‖ ≤ Ch p , where C > 0 and p > 0 are independent of ε and of mesh size h, then we say the scheme uniformly convergent to y with respect to the norm ‖·‖ [14].

International Journal of Differential Equations
Proof. For the proof, refer [15] For any mesh function z i , define the following difference operators: □ Theorem 5. Let y(x i ) and y i be, respectively, the analytical solution of equations (4) and (5) and numerical solutions of equation (13). en, for sufficiently large N, the following parameter uniform error estimate holds: Proof. Let us consider the local truncation error defined as where By considering N is fixed and taking the limit for ε ⟶ 0, we obtain the following: From Taylor series expansion, the bound for the difference becomes Now using the bounds and the assumption ε ≤ N − 1 , equation (30) reduces to

International Journal of Differential Equations
Here, 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 (36) □ Lemma 6. For a fixed mesh and for ε ⟶ 0, it holds Proof. Refer from [16]. By using Lemma 6 into equation (36), it results in Hence, by discrete maximum principle, we obtain

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. e maximum absolute errors are computed using the double mesh principle given by where y h i is the numerical solution on the mesh x i N 1 at the nodal point x i and x i � x 0 + ih, i � 1, 2, . . . , N − 1 and y h/2 i is the numerical solution on a mesh, obtained by bisecting the original mesh with N number of intervals [9]. Example 1. Consider the singularly perturbed delay reaction-diffusion equation with layer behavior under the interval and boundary conditions e maximum absolute errors are presented in Tables 1  and 2 for different values of ε and δ. e graph of the computed solution for ε � 0.01 and different values of δ is also given in Figure 1.

Example 2. Consider the singularly perturbed delay reaction-diffusion equation with layer behavior
under the interval and boundary conditions e maximum absolute errors are presented in Tables 3  and 4 for different values of ε and δ. e graph of the computed solution for ε � 0.01 and different values of δ is also given in Figure 2.
under the interval and boundary conditions y(x) � 1, e maximum absolute errors are presented in Table 5 for different values of δ.     International Journal of Differential Equations   International Journal of Differential Equations e graph of the computed solution for ε � 0.001 and δ � 0.003 is also given in Figure 3.
under the interval and boundary conditions e maximum absolute errors are presented in Table 6 for different values of δ. e graph of the computed solution for ε � 0.001 and δ � 0.003 is also given in Figure 4.    e computational rate of convergence ρ is also obtained by using double mesh principle defined as follows [9]: e following tables (Tables 7 and 8) show the rate of convergence ρ of the present method for different values of the mesh size h.

e 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.

Discussion and Conclusion
Fourth-order fitted operator numerical method for solving singularly perturbed reaction-diffusion with delay has been presented.  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. e 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 fourthorder 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 O(1), the delay affects the boundary layer solution but maintains the layer behavior ( Figure 1). When the delay parameter is of O(ε), the solution maintains layer behavior although the coefficient of the delay term in the equation of O(1) 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 h decreases, the absolute errors also decrease from Figures 5-8.

Data Availability
No data were used to support the study.

Conflicts of Interest
e authors declare that they have no conflicts of interest.