Higher-Order Numerical Method for Singularly Perturbed Delay Reaction-Diffusion Problems

In this paper, a higher-order numerical method is presented for solving the singularly perturbed delay differential equations. Such kind of equations have a delay parameter on reaction term and exhibits twin boundary layers or oscillatory behavior. Recently, different numerical methods have been developed to solve the singularly perturbed delay reaction-diffusion problems. However, the obtained accuracy and its rate of convergence are satisfactory. Thus, to solve the considered problem with more satisfactory accuracy and a higher rate of convergence, the higher-order numerical method is presented. First, the given singularly perturbed delay differential equation is transformed to asymptotically equivalent singularly perturbed twopoint boundary value convection-diffusion differential equation by using Taylor series approximations. Then, the constructed singularly perturbed boundary value differential equation is replaced by three-term recurrence relation finite difference approximations. The Richardson extrapolation technique is applied to accelerate the fourth-order convergent of the developed method to the sixth-order convergent. The consistency and stability of the formulated method have been investigated very well to guarantee the convergence of the method. The rate of convergence for both the theoretical and numerical have been proven and are observed to be in accord with each other. To demonstrate the efficiency of the method, different model examples have been considered and simulation of numerical results have been presented by using MATLAB software. Numerical experimentation has been done and the results are presented for different values of the parameters. Further, The obtained numerical results described that the finding of the present method is more accurate than the findings of some methods discussed in the literature.


Introduction
The particularity of noticing the relation between causes and effects arises when the cause is small and the effect is large. In the philosophy of perturbation for mathematics and physical systems, the study of this relation got a significant amount of attention in past and recent years. A singularly perturbed differential equation is a differential equation in which the highest order derivative is multiplied by a small parameter that is recognized as a perturbation parameter. The solution of singularly perturbed differential equations varies rapidly in the regions called boundary layers. The study of the solution of these equations is of great significance due to the formation of sharp boundary layers when the perturbation parameter approaches zero. Singularly perturbed delay differential equation is an equation in which the evolution of the system at a convinced time depends on the rate at an earlier time.
The delay in the process rises due to the requirement of a definite time to sense the instruction and react to it. The delay differential equation can be classified as retarded delay differential equation and neutral differential equation. The applications of delay differential equations arise in the modeling of neural variables, variational problems in control theory, description of human pupil reflex, physical and biological phenomena like optically bi-stable devices, numerical modeling in biosciences, and HIV infection [1][2][3][4][5].
In singular perturbation problems, when the perturbation parameter equal to zero and if the order is reduced by one then the problem is called convection-diffusion type, whereas if the order is reduced by two it is called reaction-diffusion type. Hence, second-order singularly perturbed delay differential equations are may be convection-diffusion or reaction-diffusion types. The former type of problem exhibits left or right boundary layer only depending on the sign of coefficient in diffusion and convection terms, while reaction-diffusion delay differential equations have dual boundary layers (both at the left and right side of the domain) or oscillatory behavior of solution depending on the sign of the sum of coefficients in reaction terms. Researchers have been developed and analyzed different numerical methods for solving singularly perturbed delay differential equations. For instance, numerical methods proposed by various authors for solving singularly perturbed delay differential equations of convection-diffusion type are; Parameter-robust numerical method based on defect-correction technique [6], Fitted mesh numerical method [7], Numerical integration method [8], using B-Spline collocation method [9], Fitted Method [10], Terminal boundary-value technique [11], Numerical integration using exponential integrating factor [12], Fitted fourth-order scheme [13] and A non-asymptotic method for general singular perturbation problems [20] so on. Further, more recently, various authors have developed numerical methods for solving singularly perturbed delay reactiondiffusion problems like; Computational method [14][15][16][17], a fourth-order numerical method [18], trigonometric B-spline [4], and Computational method for singularly perturbed delay differential equations of the reaction-diffusion type with negative shift [22].
All these works concern second-order singularly perturbed delay differential equations in which the developed methods are analyzed in different approaches and produce good accurate numerical solutions corresponding with a diverse rate of convergence to demonstrate the efficiency of the methods. However, the obtained approximate solution and the corresponding order of convergence are not more satisfactory which indicates that yet to solve singularly perturbed delay reaction-diffusion problems demands to develop of other numerical methods to produce a more accurate numerical solution. Therefore, the main objective of this paper is to present a higher-order numerical method to solve the singularly perturbed delay reaction-diffusion problems.

Formulation of the Method
Consider the singularly perturbed delay reaction-diffusion problem: subject to the interval and boundary conditions, where ε is perturbation parameter satisfies 0  (1) and (2) exhibits two boundary layer behavior which will occur at both endpoints 0 x = and 1 x = , while ( ) ( ) 0 a x b x + > , it exhibits oscillatory behavior, [4,5,16].
Taylor's series is familiar to handle the delay term as Substituting equation (3) into equation (1), gives: under the boundary conditions To develop the finite difference method for the problem in equation (4) Assume that ( ) y x has continuous higher-order derivatives on [0,1], using Taylor's series expansion, we have: Subtracting equation (7) from equation (6) or adding the two equations, we yield Substituting equation (8) into the discrete form of equation (4) gives: where 0 Once and twice successively differentiating both sides of equation (4), concerning the independent variable and considering at the nodal point i x , yields (4) and i i y y ′′′ respectively. Then substituting these values into equation (9), gives the three-term recurrence relation of the form:

Richardson Extrapolation
The purpose of this section is to improve the accuracy and the order of convergence by convergence acceleration technique which involves a combination of two computed approximate solutions. The linear combination turns out to be a better approximation, [24,26].
Ever since from equation (8) or equation (9), we know that the truncation error of the formulated method is where ( ) and So, it works for any mesh size 0 2 h ≠ leads to 4 2 2 2 ( ) , 16 where the remainders, N R and 2N R are Assume that C is a constant in both equations (12) and (13), we multiply equation (13) by 16, and then combining the two inequalities in equations (12) and (13) gives 16 16 .
This can be re-written as and further by notation is also an approximation to the exact solution ( ) i y x . Hence, using the approximate solution in equation (14), and the truncation terms in equations (12) and (13), we further obtain the truncation error: This indicates that the Richardson extrapolation method accelerates the order of convergence from fourth-order to sixth-order.

Consistency of the Method
Local truncation error refers to the difference between the differential equation and its finite difference approximations. A finite difference scheme is consistent if the limit of truncation error (TE) is equal to zero as the mesh size h goes to zero. Further, local truncation errors measure how well a finite difference discretization approximates the differential equation [19,21,[23][24][25][26].
Truncation error ( TE ) from equations (8) and (9) then, it is possible to show that (5) This implies that the developed method is accurate in order four. Hence, from the definition of consistency, both the truncation errors in equations (15) and (16) Thus, the proposed method is consistent, (See [24]).

Stability of the Method
Consider the developed scheme in equation (10) which is given by: If we multiply both sides of equation (10) where the matrices: This proves the diagonal dominant of M . Under these conditions, any tri-diagonal solver is stable for sufficiently small h, as shown in the book "Numerical solution of differential equations, Introduction to finite difference and finite element methods", [19]. As proved by Smith [21], the eigenvalues of a tri-diagonal ( ) ( ) A finite difference method for solving differential equations is stable if M is invertible and where C and 0 h are two constants that are independent of h , [19,21].
Since matrix M is symmetric also its inverse matrix Thus, the developed scheme in equation (10) is stable. A consistent and stable finite difference method is convergent by Lax's equivalence theorem. Hence, as we have shown above the proposed method is satisfying the criteria for both consistency and stability which are equivalents to the convergence of the method.

Numerical Examples and Results
To confirm the applicability of the formulated method, we have been realized the method on four examples, the first two with twin boundary layers, while the remaining with oscillatory behavior. Since these examples have no exact solution; consequently the numerical solutions are computed using the double mesh principle. The maximum absolute errors are calculated using the double-mesh principle before and after applying Richardson extrapolation is given by The corresponding order of convergence is determined by Example 1: Consider the singularly perturbed delay differential equation ( ) 2 ( ) ( ) 1, (0,1), The maximum absolute errors are presented in Tables 1  and 3 for different values of parameters and numerical computation in graphs as Figures 1 and 2.
Example 2: Consider the differential equation   Comparisons of maximum absolute errors are presented in Tables 2 and 3   (1) 0.
The maximum absolute errors are presented in Tables 4 and 7 with the corresponding rate of convergence given in Table 6 and graphical representation in Figures 3 and 4.
Example 4: Consider the delay differential equation ( ) ( ) 2 ( ) 1, (0,1), The maximum absolute errors are presented in Tables 5 and 7 with the corresponding rate of convergence given in Table 6 and graphical representation in Figures 3 and 4.

Discussion and Conclusion
Numerical solution of the second-order singularly perturbed delay reaction-diffusion equations that exhibits boundary layer or oscillatory behavior via higher-order numerical method have been presented. To achieve the higher-order method, the fourth-order finite difference method is accelerated to the sixth-order one using the Richardson extrapolation technique. The efficiency of the method is validated by numerical examples and results for different parameters. The obtained numerical results have been compared with the results obtained by the methods in "Solution of second-order singular perturbed delay differential equation using Trigonometric B-Spline" and "Exponentially Fitted Numerical Method for Singularly Perturbed Differential-Difference Equations", [4,18], (See Tables 1, 2, 3, 7). Numerical confirmation to the contribution of applying the Richardson extrapolation technique is presented in Tables 4, 5, 6. Besides, the obtained maximum absolute errors decrease rapidly as the number of mesh points N increases which indicates the convergence of the formulated method. Also, the consistency and stability of the method were investigated to guarantee convergence analysis.
Furthermore, Figures 1, 2, 3 demonstrate the effects of the parameters on the numerical solutions. Figure 4, to verify that the consequence before and after applying Richardson extrapolation and effect of decreasing number of mesh sizes of the domain on the numerical solution in case of different delay parameters and the number of mesh points. Accordingly, after applying Richardson extrapolation accuracy of the solution improved with an accelerated rate of convergence, and as the number of mesh point's N increases, the accuracy of the numerical solution increases.
Overall, a higher-order numerical method for solving the singularly perturbed delay differential equation is presented. This method is stable, consistent, and produces a more accurate solution than some existing methods for the differential equation under consideration [4,18]. The interested researcher will be formulating the eighth or higher-order convergent numerical methods to obtain a more accurate solution.