Fitted Numerical Scheme for Second-Order Singularly Perturbed Differential-Difference Equations with Mixed Shifts

This paper presents the study of singularly perturbed differential-difference equations of delay and advance parameters. The proposed numerical scheme is a fitted fourth-order finite difference approximation for the singularly perturbed differential equations at the nodal points and obtained a tridiagonal scheme. This is significant because the proposed method is applicable for the perturbation parameter which is less than the mesh size
 
 ,
 
 where most numerical methods fail to give good results. Moreover, the work can also help to introduce the technique of establishing and making analysis for the stability and convergence of the proposed numerical method, which is the crucial part of the numerical analysis. Maximum absolute errors range from 
 
 
 
 10
 
 
 −
 03
 
 
 
 up to 
 
 
 
 10
 
 
 −
 10
 
 
 
 , and computational rate of convergence for different values of perturbation parameter, delay and advance parameters, and mesh sizes are tabulated for the considered numerical examples. Concisely, the present method is stable and convergent and gives more accurate results than some existing numerical methods reported in the literature.


Introduction
Singularly perturbed differential-difference equations (SPDDEs) occur frequently in the mathematical modeling of various physical and biological phenomena, for example, control theory, viscous elasticity, and population dynamics [1]. Recently, many researchers have started developing different numerical methods for solving differential equations. Reference [2] investigated that the nonlinear thermal radiation and dissipation with the Darcy-Forchheimer equation in the porous medium analysis by using the fifth-order Runge-Kutta method, [3] also discussed the cross-fluid flow containing gyrotactic microorganisms and nanoparticles on a horizontal and three-dimensional cylinder by using the Runge-Kutta Fehlberg fifth-order technique, [4] studied the three-dimensional convective heat transfer of magnetohydrodynamics nanofluid flow through a rotating cone by using the fifth-order Runge-Kutta method. Reference [5] discussed that the heat transfer hybrid nanofluid contains 1-butanol as the base fluid and MoS 2 -Fe 3 O 4 hybrid nanoparticles by using the finite element method. In these papers, the influence of various parameters on velocity profile and temperature has been investigated. Reference [6] presented a finite difference numerical method for solving singularly perturbed delay differential equations, [7] introduced a Galerkin method for solving this problem, [8] also introduced a fitted tension spline method for solving such problem, [9] presented a fitted second-order numerical method for singularly perturbed problems, [10][11][12][13][14] also developed some numerical methods of different orders for solving singularly perturbed delay differential equations, and so on. However, the issue of convergence and accuracy still needs attention and improvement. In this paper, we present a stable, convergent, and more accurate exponentially fitted fourth-order numerical scheme for solving SPDDEs and investigate the influence of delay and advance parameters on the solution profile.
Using the uniform mesh technique over the domain, we have By using the Taylor series expansion, we obtain Subtracting Equation (8) from Equation (7) gives the approximation δ 1 c y i , for the first derivative of y i as where Similarly, adding Equations (7) and (8) provides the approximation δ 2 c y i , for the second derivative of y i as where T 2 = −ðh 2 /12Þy ð4Þ i . Substituting Equations (7) and (8) into Equation (9) yields where T 3 = ðh 4 /120Þy i . Again, substituting Equations (7) and (8) into Equation (10) gives where T 4 = ðh 4 /360Þy Using Equation (13) into Equation (11), we obtain where T 5 = ð13h 4 /360Þy i . Applying δ 2 c to y i ″ in Equation (10), we get a four-order finite difference scheme for Equation (5): Substituting Equation (15) into Equation (12), we get where T 6 = ð7h 4 /720Þy i . From Equations (14) and (16), we get After evaluating Equation (5) at nodal point x i and using Equation (17), we obtain Rearranging Equation (5) and successively differentiating, evaluating at x i , and substituting into Equation (19), we get Introducing the fitting factor ðσÞ into Equation (20), we have Using the central difference approximation for δ 2 c y i and δ 1 c y i , multiplying both sides of Equation (21) by h, and evaluating the limit as h ⟶ 0, we obtain where ρ = ðh/εÞ. From the theory of singular perturbations and O'Malley [17], we have two cases for pðxÞ > 0 andpðxÞ < 0.
Thus, from Equation (22), we get Case 2. For pðxÞ > 0 (left-end boundary layer), we have Thus, from Equation (22), we get In general, for discretization, we take a variable fitting parameter as Now using Equations (9) and (10) into Equation (21) for δ 1 c y i and δ 2 c y i and making use of δ c

Abstract and Applied Analysis
where Simplifying Equation (24), we get a tridiagonal system: where ; , ; , We have used a discrete invariant imbedding algorithm to solve Equation (25).

Stability Analysis
The continuous minimum principle, continuous maximum principle, and stability of the solution of Equations (5) and (6) are presented in [13]. We present the stability of the scheme in Equation (25) for both cases.
Proof. Suppose ∃k ∈ ℤ + such that w k < 0 and w k = min 0≤i≤N w i .
Then, from Equation (25), we have where

Abstract and Applied Analysis
For sufficiently small h (i.e., h ⟶ 0) and for suitable values of p k , we obtain L N w k > 0. Since, w k < 0 (by assumption), ε, σ k > 0 and F k ⟶ q k < 0:.
However, this is a contradiction. Hence, w i ≥ 0 for all x i ∈ ð0, 1Þ. jLw i jg, for some constant λ ≥ 1.

Lemma 3 (discrete maximum principle).
If w i is any mesh function such that w 0 ≥ 0 and L N w i ≥ 0, then w i ≥ 0 for all x i ∈ ð0, 1Þ.
Proof. Suppose ∃k ∈ ℤ + such that w k < 0 and w k = max 0≤i≤N w i .
Then, from Equation (25), we have For sufficiently small h and for suitable values of p k , we obtain L N w k < 0. Since, w k < 0 (by assumption), ε, σ k > 0 and F k ⟶ q k > 0:.
However, this is a contradiction. Hence, w i ≥ 0 for all x i ∈ ð0, 1Þ.     Proof. The proof is analogous to Theorem 2.
This proves the stability of the scheme for the case of oscillatory behaviour. ☐

Convergence Analysis
Writing the scheme in Equation (25) in matrix form, we obtain where M = ðm i j Þ, i, j = 1, 2, ⋯, N − 1, is a tridiagonal matrix of order N − 1.

Multiplying both sides of Equation (25) by ð−h 2 Þ, we have
and V = ðv i Þ is a column vector, where with a local truncation error: where K = ðp i /3Þy i . Equation (34) can also be written in error form as where Y = ð y 1 , y 2 , ⋯, y N−1 Þ t stands for the exact solution and TðhÞ = ðT 1 ðhÞ, T 2 ðhÞ, ⋯, T N−1 ðhÞÞ t is a local truncation error.
From Equations (34) and (38), we obtain which implies where E = Y − Y = ðe 1 , e 2 , ⋯, e N−1 Þ t . Let S i be the sum of elements of the i th row of a matrix M, then we get  Abstract and Applied Analysis Table 5: Numerical results of Example 7 for ε = 0:01, N = 100.
Thus, Equation (40), gives Let m k,i ≥ 0 be the ðk, iÞ th element of M −1 . From the theory of matrices, we have Therefore, where which implies where i 0 is some number between i and N. Therefore, kEk = Oðh 4 Þ: Hence, the present method is of fourth-order convergence.

Numerical Examples and Results
To show the applicability of the method, three model examples have been considered. The exact solution of Equations (1) and (2) with constant coefficients is where ; ; For the variable coefficients, the maximum absolute errors are computed using the double mesh principle [13].
The following graphs (Figure 3) show the pointwise absolute errors for different values of mesh size h.

Discussion and Conclusion
This study introduces an exponentially fitted fourth-order numerical method for solving singularly perturbed differential-difference equations. The results observed from the tables demonstrate that the present method approximates the solution very well and depicts the betterment over some existing numerical methods reported in the literature. The stability and convergence of the scheme have been established. The solutions presented in Table 1 confirm that the numerical rate of convergence as well as theoretical error estimates indicates that the present method is of fourth-order convergence.
To demonstrate the effect of delay and advance parameters on the left and right boundary layers of the solution, the graphs for different values of delay parameter δ and advance parameter η are plotted in Figures 1 and 2. Accordingly, based on the sign of pðxÞ, one can see that, from Figure 1, as δ increases, the width of the right boundary layer decreases for a fixed value of η, but as η increases, the width of the right boundary layer increases for a fixed value of δ while the width of the left boundary layer decreases when δ or η increases ( Figure 2). Furthermore, as h decreases, the absolute error also decreases (see Tables 2-6 and Figure 3). In general, the present method is stable, convergent, and more accurate.

Data Availability
No data were used to support this study.

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