Accelerated Numerical Method for Singularly Perturbed Differential Difference Equations

—In this paper, accelerated finite difference method for solving singularly perturbed delay reaction-diffusion equations is presented. First, the solution domain is discretized. Then, the derivatives in the given boundary value problem are replaced by finite difference approximations and the numerical scheme that provides algebraic systems of equations is obtained, which can easily be solved by Thomas algorithm. The consistency, stability and convergence of the method have been established. To increase the accuracy of our established scheme we used Richard- son’s extrapolation techniques. To validate the applicability of the proposed method, four model examples have been considered and solved for different values of perturbation parameters and mesh sizes. The numerical results have been presented in tables and graphs to illustrate; the present method approximates the exact solution very well. Moreover, the present method gives better accuracy than the existing numerical methods mentioned in the literature.


I. INTRODUCTION
S INGULARLY perturbed ordinary differential equation with a delay is ordinary differential equations in which the highest derivative is multiplied by a small parameter and involving at least one delay term. 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 [1], the study of bi stable devices [2], and vibrational problems in control theory [3], etc. When perturbation parameter is very small, most numerical methods for solving such problems may 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, in [4] proposed computational method of first order for singularly perturbed delay reaction-diffusion equations with layer or oscillatory behavior. Authors in [5] presented fourth order finite difference scheme for second order singularly perturbed differential-difference equation with negative shift. Authors of [6] presented exponentially fitted second order finite difference scheme for a class of singularly perturbed delay differential equations with large delay. In [7], the numerical solution of singularly perturbed differential-difference equations with dual layer was presented. Recently, [4] presented computational method for solving singularly perturbed delay differential equation with twin layers or oscillatory behavior. Similarly, in [8] fourth order numerical method for singularly perturbed delay differential equations. But, still there is a lack of accuracy because of the treatment of singularly perturbed problems is not trivial and the solution depends on perturbation parameter and mesh size [9], [10], [11]. Due to this, numerical treatment of singularly perturbed delay differential equations needs improvement. Therefore, it is important to develop more accurate and convergent numerical method for solving singularly perturbed delay differential equations. Thus, the purpose of this study is to develop stable, convergent and more accurate numerical method for solving singularly perturbed delay reaction-diffusion equations.

II. DESCRIPTION OF THE METHOD
To describe the method, we first consider a linear singularly perturbed delay reaction-diffusion equation of the form: subject to the interval and boundary conditions, where ε is 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. The 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 Eqs. (1) and (2) exhibits layer behavior and if a(x) + b(x) > 0, 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 x = 0 and x = 1 [12]. By using Taylor series expansion in the neighborhood of we have: Substituting Eq. (3) into Eq. (1), we obtain an asymptotically equivalent singularly perturbed two point boundary value problem: under boundary conditions, To describe the scheme, we divide the intervalΩ = [0, 1] into N equal subintervals of uniform mesh length h. Let x i be the mesh points for i = 0, 1, 2, ..., N such that 0 = x 0 < N and x i = ih, for i = 0, 1, 2, ..., N . Let us denote a( related to a mesh point x i . Assume that the solution y(x) has continuous higher order derivatives onΩ = [0, 1] and expanding y i±1 by Taylor's series expansion about the point x i for x i−1 and x i+1 gives Subtracting Eq. (7) from Eq. (6) and adding up both Eqs. (6) and (7) gives the following two equations respectively.
For more clarity, let us re-write Eq. (10) as and τ 3 = p i τ 1 + ετ 2 . Since y(x) assumed to have continuous higher order derivatives, differentiating the continuous problem in Eq. (4) successively and considering the results of differentiation at the nodal point x i gives the following Evaluating Eq. (13) in terms of y ′ i and y ′′ i , and then substitute into the D which given in Eq. (11) written as From Eqs. (8) and (9), we have the following finite approximations for y ′ i and y ′′ i as Substituting Eq. (15) into Eq. (14) yields D written in term of the difference term y i and y i+1 . Again, putting this value of D in to Eq. (11), and then after collecting like terms, we get the following three term recurrence relation.

III. RICHARDSON EXTRAPOLATION
Extrapolation can be applied whenever it is known that an approximation technique has an error term with a predictable form, one that depends on a parameter, usually the step size h. We combine two approximations obtained from the above scheme in Eq. (16) using different values of the mesh sizes, h and h/2 to obtain a higher order approximation. Such technique is known as Richardson extrapolation. This procedure is 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. In our case, from Eq. (11) we know that truncation error for the present scheme is O(h 4 ), which implies: where y(x i ) and Y N are exact and approximate solutions respectively, C is constant independent of mesh sizes h.
Let Ω 2n be the mesh obtained by bisecting each mesh interval in Ω n and denote the approximation of the solution on Ω 2n by Y 2N . Consider Eq. (17) works for any h ̸ = 0, which yields: So that, it works for any h 2 ̸ = 0 produces: where the remainder terms, R N and R 2N are O(h 6 ).

Combination of inequalities given in Eqs. (18) and (19) clues to
15y This suggests that which is also an approximation of y(x i ).
Using this approximation to evaluate the truncation error, we obtain: Hence, Richardson extrapolation method accelerates the rate of convergence for the developed finite difference scheme from fourth order to sixth order convergent.
IV. CONSISTENCY OF THE METHOD Local truncation errors refer to the differences between the original differential equation and its finite difference approximations at the grid points. They measure how well a finite difference approximates the differential equation [13]. A finite difference scheme is called consistent if the limit of truncation error (T E ) is equal to zero as the mesh size h goes to zero. Hence, in this work, definition of consistency on the proposed method which is given in Eq. (16) Considering Eq. (22) and after incorporating the boundary conditions, y 0 = y(0) = α and y N = y(1) = β into Eq. (16), we have a system, which can be written in matrix form: where the matrices: Here, the coefficient matrix A is a tridiagonal matrix. The codiagonals of A contain E i , G i such that for sufficiently small h(i.e.h → 0), E i ̸ = 0 and G i ̸ = 0, ∀i = 1, 2, ..., N − 1.
Hence, A is irreducible [14]. Again one can observe that, ∥E i ∥ > 0, ∥G i ∥ > 0, ∥F i ∥ > 0 and the sum of the two off diagonal elements is less than or equal to the modulus of the diagonal element. i.e.∥E i + G i ∥ ≤ ∥F i ∥, ∀i = 1, 2, ..., N − 1 This proves that the diagonal dominance of A. Hence, A is diagonally dominant. Under these conditions, the Thomas Algorithm is stable for sufficiently small h, as shown in [15]. Moreover, as proved by Smith [16] the Eigenvalues of a tridiagonal matrix (N − 1) × (N − 1) of matrix A are: .., N −1 (24) Also, from trigonometric identity, we have 1 − cos sπ N = 2sin 2 sπ 2N , Hence the Eigenvalues of matrix A can be re-written as: A finite difference method for the boundary value problems is stable if A is invertible and where C and r are two constants that are independent of h [16].
Since the matrix A is symmetric, also its inverse A −1 is symmetric and the Eigenvalues of A −1 is given by 1 λs , we have Thus the developed scheme in Eq. (16) is stable. A consistent and stable finite difference method is convergent by Lax's equivalence theorem [16]. Hence, as we have shown above, the proposed method is satisfying the criteria for both consistency and stability which are equivalent to the convergence of the method.
Example 3: Consider the singularly perturbed delay reaction-diffusion equation with oscillatory behavior, εy ′′ (x) + 0.25y(x − δ) + y(x) = 1 under the interval and boundary conditions  1.57e-02 1.84e-03 1.20e-04 7.79e-06 4.89e-07 2 10 3.38e-02 6.20e-03 4.93e-04 3.15e-05 1.99e-06 The maximum absolute error are presented in Table V for different values of δ. The graph of the computed solution for ε = 0.001 and different values of δ is also given in the Figure  5. Example 4: Consider the singularly perturbed delay reaction-diffusion equation with oscillatory behavior,   under the interval and boundary conditions The maximum absolute error are presented in Table 6 for different values of δ. The graph of the computed solution for ε = 0.001 and different values of δ is also given in the Figure  6.

A. 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 have been given by the following graphs.       (Tables (I) to (VI)) in terms of maximum absolute errors and observed that the present method improves the findings of [8]. Also, it is significant that all of the maximum absolute errors decrease rapidly as N increases.
Further, 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  4 and 5). In a concise manner, the present method gives more accurate solution for solving singularly perturbed delay reaction-diffusion equations with twin layer and oscillatory behavior. Also it can see that as mesh size h decrease the absolute errors also decrease from (Figures 3, 4, 7, 8).