Accelerated fitted operator finite difference method for singularly perturbed delay differential equations with non-local boundary condition

In this paper, accelerated fitted finite difference method for solving singularly perturbed delay differential equation with non-local boundary condition is considered. To treat the non-local boundary condition, Simpson’s rule is applied. The stability and parameter uniform convergence for the proposed method are proved. To validate the applicability of the scheme, two model problems are considered for numerical experimentation and solved for different values of the perturbation parameter ε and mesh size h. The numerical results are tabulated in terms of maximum absolute errors and rate of convergence, and it is observed that the present method is more accurate and ε-uniformly convergent for h ≥ ε where the classical numerical methods fails to give good result, and it also improves the results of the methods existing in the literature.


Introduction
A differential equation is said to be singularly perturbed delay differential equation, if it includes at least one delay term, involving unknown functions occurring with different arguments, and also, the highest derivative term is multiplied by a small parameter. Such type of delay, differential equations play a very important role in the mathematical models of science and engineering, such as, the human pupil light reflex with mixed delay type [1], variational problems in control theory with small state problem [2], models of HIV infection [3], and signal transition [4]. Any system involving a feedback control almost involves time delay. The delay occurs because a finite time is required to sense the information and then react to it. Finding the solution of singularly perturbed delay differential equations, whose application mentioned above, is a challenging problem. In response to these, in recent years, there has been a growing interest in numerical methods on singularly perturbed delay differential equations. The authors of [5][6][7] have developed various numerical schemes on uniform meshes for singularly perturbed first and second order differential equations with integral boundary conditions. The authors of [8][9][10] have proved that the problem of differential equations with integral boundary conditions is well posed. The authors of [11] proposed a first order uniform convergent fitted finite difference scheme for singularly perturbed boundary value problem for a linear second order delay differential equation with large delay in reaction term.
The standard numerical methods used for solving singularly perturbed differential equation are sometime ill posed and fail to give analytical solution when the perturbation parameter ε is small. Therefore, it is necessary to develop suitable numerical methods which are uniformly convergent to solve this type of differential equations. In [12][13][14][15], finite difference and finite element methods are proposed to solve this kind of equations with large and small shifts.
As far as the researchers' knowledge is concerned, numerical solution of singularly perturbed boundary value problem containing integral boundary condition via accelerated exponential fitted operator method is first being considered. The basic essence of accelerated fitted operator finite difference method is fitting an operator into a finite difference scheme, determining the value for the operator, formulating the intended scheme, and then applying Richardson extrapolation that can be explained as whenever the leading term in the error for an approximation formula is known; we can combine two approximations obtained from the formula using different values of the mesh sizes h and 0.5h to obtain a higher order approximation, and the technique is known as Richardson extrapolation. This procedure is a convergence acceleration technique which consists of a linear combination of two computed approximations of a solution (applied on two nested meshes). The linear combination turns out to be a better approximation. Therefore, the main objective of this study is to develop ε-uniformly convergent and more accurate numerical method for solving singularly perturbed delay differential equations with non-local boundary condition. Hence, in the present paper, motivated by the works of [16], we developed a fitted operator finite difference scheme on uniform mesh for the numerical solution for second order singularly perturbed convection-diffusion equations with negative shift and non-local boundary condition.
The present paper is organized as follows. Statement of the problem is given in the "Statement of the problem" section. In the "Properties of continuous solution" section, properties of continuous solution are presented. The "Formulation of the numerical scheme" section describes formulation of the numerical scheme. Convergence analysis for approximate solution is given in the "Convergence analysis" section. Numerical results are given in the "Numerical examples and results" section. Discussion and conclusion is given in the "Discussion and conclusion" section.
Suppose the theorem does not hold true, then μ > 0.
it is a contradiction: Hence, the proof of the theorem.
Proof: This theorem can be proved by using Lemma 1, and the barrier functions are the test function as in Lemma 1. Lemma 3: Let y(x) be the solution for (1)-(3). Then we have the following bounds: Proof: For the proof, refer to [16]. Lemma 4: The bound for derivative of the solution y(x) of the problems (1)-(3) when x ∈ Ω 1 = (0, 1) is given by: :::; N−1: Proof: For the proof, refer to [17].

Formulation of the numerical scheme
For small values of ε, the boundary value problem, (1)-(3) exhibit strong boundary layer at x = 2 and interior layer at x = 1 (see [16]). The linear ordinary differential Eq. (1) cannot, in general, be solved analytically because of the dependence of a(x), b(x), and c(x) on the spatial coordinate x. We divide the interval [0, 2] into 2N equal parts with constant mesh length h. Let 0 = x 0 , x 2 , ..., x N = 1, x N + 1 , x N + 2 , ..., x 2N = 2 be the mesh points. Then, we have x i = ih, i = 0, 1, 2, ...2N. If we consider, the interval x ∈ (0, 1) and the coefficients of (1) are evaluated at the midpoint of each interval; then, we will obtain the differential equation: Now, the domain [0, 1] is discretized into N equal number of subintervals, each of length h. Let 0 = x 0 < x 1 < ... < x N = 1 be the points such that x i = ih, i = 0, 1, 2, ..., N. For the discretization, we apply a exponentially fitted operator finite difference method (FOFDM).
From (9), we have To find the numerical solution of (10), we use the theory applied in asymptotic method for solving singularly perturbed BVPs. In the considered case, the boundary layer is in the right side of the domain, i.e., near x = 1. From the theory of singular perturbations given by O'Malley [18] and using Taylor's series expansion for a(x) about x = 1 and restriction to their first terms, we get the asymptotic solution as follows: where y 0 (x) is the solution of the reduced problem (obtained by setting ε = 0) of (10) which is given by: Considering h is small enough, the discretized form for (11) becomes which is simplified to To handle the effect of the perturbation parameter, artificial viscosity (an exponentially fitting factor σ(ρ)) is multiplied on the term containing the perturbation parameter as follows: with boundary conditions y 0 (0) = ϕ(0) and y(N) = θ, where y(N) is evaluated by Runge-Kutta method from the reduced solution of (12).
Next, we consider the difference approximation of (9) on a uniform grid Ω For any mesh function z i , define the following difference operators: by applying the central finite difference scheme on (14) takes the form: with the boundary conditions y 0 (0) = ϕ(0) and y(N) = θ.
Debela and Duressa Journal of the Egyptian Mathematical Society (2020) 28:16 multiplying (18) by h and considering h is small and truncating the term Now, by using Taylor's series for y i − 1 and y i + 1 up to first term and substituting the results in (19) into (16) and simplifying, the exponential fitting factor is obtained as follows: Assume that Ω 2N denotes partition of [0, 2] into 2N subintervals such that 0 = x 0 < Case 1: Consider (4) on the domain Ω 1 = (0, 1) which is given by: Hence, the required finite difference scheme becomes for i = 0, 1, 2, ..., N.
The numerical scheme in (22) can be written in three-term recurrence relation as follows: where  (1, 2), for right layer in the domain Ω 2 using exponentially fitted finite difference method, which is given by: Similarly, this equation can be written as follows: where y j ¼ y x i −1 ð Þ; j ¼ 1; 2; :::; N Substituting (25) into (3) gives: (2), this equation can be re-written as follows: Therefore, on the whole domain Ω ¼ ½0; 2, the basic schemes to solve (1)
It is a contradiction. Hence the proof of the theorem. Lemma 6: Let ψ(x i ) be any mesh function then for 0 ≤ i ≤ 2N, Proof: For the proof, refer to [16].
The following theorem shows the parameter uniform convergence of the scheme developed.
Theorem 1: Let y(x i ) and y i be respectively the exact solution of (1)-(3) and numerical solutions of (17). Then, for sufficiently large N, the following parameter uniform error estimate holds: Proof: Let us consider the local truncation error defined as follows: where εσðρÞ ¼ að1Þ N −1 2 cothðað1Þ N −1 2ε Þ since ρ ¼ N −1 ε . In our assumption, ε ≤ h = N −1 . By considering N is fixed and taking the limit for ε → 0, we obtain the following: Now, using the bounds and the assumption ε ≤ N −1 , (33) reduces to: 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 solution in Lemma 4, we obtain: Lemma 7: For a fixed mesh and for ε → 0, it holds: Proof: Refer to [19]. By using Lemma 7 into (35), results to Hence, by discrete maximum principle, we obtain: Thus, result of (38) shows (32). Hence, the proof. Remark: A similar analysis for convergence may be carried out for the finite difference scheme (24).

Richardson Extrapolation
This technique is acceleration technique which involves combination of two computed approximations of a solution. The combination goes out to be an improved approximation. From the local truncation term, we have: where y(x i ) and y i are exact and approximate solutions respectively, and C is constant free from mesh size h.
Let Ω 4N be the mesh found by dividing each mesh interval in Ω 2N and symbolize the calculation of the solution on Ω 4N by y i . Consider (39) works for any h ≠ 0, which implies: So that it works for any h 2 ≠0 yields: where the remainders R 2N and R 4N are O(h 2 ). Combination of inequalities in (40) and (41) leads to yðx i Þ−ð2y i −y i Þ ≈ Oðh 2 Þ which proposes that is also a rough calculation of y(x i ). By means of this approximation to estimate the truncation error, we obtain: where C is free of mesh size h. Thus, using Richardson extrapolation first order convergent method is accelerated into second order convergent as provided in (43). Thus, we can say that the proposed method is second order convergent.

Numerical examples and results
In this section, two examples are considered to illustrate the applicability of the numerical method discussed above. The exact solutions of these test problems are not known. Therefore, double mesh principle is used to estimate the errors and compute the numerical rate of convergence to the computed solution. The double mesh formula to determine maximum absolute error is defined as follows: where Y N i and Y 2N 2i are the i th components of the numerical solutions for N and 2N, respectively. We compute the uniform error and the rate of convergence using the formula: The numerical results are presented for the values of the perturbation parameter ε ∈ {10 −4 , 10 −8 , ..., 10 −20 }.

Discussion and conclusion
This study introduces accelerated fitted operator numerical method for solving singularly perturbed delay differential equations with integral boundary condition. The behavior of the continuous solution of the problem is studied and shown that it satisfies the continuous stability estimate, and the derivatives of the solution are also bounded. The numerical scheme is developed on uniform mesh using fitted operator finite difference method in the given differential equation. The integral boundary condition is treated by using Simpson's rule. The stability of the developed numerical method is established, and its uniform convergence is proved. To validate the applicability of the method, two model problems are considered for numerical experimentation for different values of the perturbation parameter and mesh points. The numerical results are tabulated in terms of maximum absolute errors, numerical rate of convergence, and uniform errors (see Tables 1, 2 layer at x = 2, and interior layer at x = 1. Figures 3 and 4 show that as the mesh size decrease or as the number of mesh point increase, the absolute error decreases. The loglog scale plot in Figs. 5 and 6 depicted the ε-uniformly convergence of the method for h ≥ ε where the classical numerical methods fail to converge. The method is shown to be ε-uniformly convergent with order of convergence O(h 2 ). The performance of the proposed scheme is investigated by comparing with prior study (Tables 2 and 4). The proposed method is stable, more accurate, and convergent independent of the values of the perturbation parameter and the mesh size. The authors suggested that one can extend the work or solve the problem under consideration by applying higher order fitted operator numerical methods or Bakhavlov-type fitted mesh numerical method to obtain more accurate numerical results.