Robust numerical method for singularly perturbed differential equations having both large and small delay

Purpose – The purpose of this study is to develop stable, convergent and accurate numerical method for solving singularly perturbed differential equations having both small and large delay. Design/methodology/approach – This study introduces a fitted nonpolynomial spline method for singularly perturbed differential equations having both small and large delay. The numerical scheme is developed on uniform mesh using fitted operator in the given differential equation. Findings – The stability of the developed numerical method is established and its uniform convergence is proved. To validate the applicability of the method, one model problem is considered for numerical experimentation for different values of the perturbation parameter and mesh points. Originality/value – In this paper, the authors consider a new governing problem having both small delay on convection term and large delay. As far as the researchers ’ knowledge is considered numerical solution of singularly perturbed boundary value problem containing both small delay and large delay is first being considered.


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 second-order differential equations having small delay on convection term.The authors of [8][9][10][11][12] have presented second-order differential equations having large delay.
In this paper, we consider a new governing problem having both small delay on convection term and large delay.Additionally, in recent years the correlative physical phenomena in-depth of the problem under consideration have been done by the authors [13][14][15][16][17].As far as the researchers' knowledge is considered numerical solution of singularly perturbed boundary value problem containing both small delay and large delay is first being considered.Thus, the purpose of this study is to develop stable, convergent and accurate numerical method for solving singularly perturbed differential-difference equations having both small and large delay.
Throughout our analysis C is generic positive constant that are independent of the parameter ε and number of mesh points 2N.We assume that Ω ¼ ½0; 2, Ω ¼ ð0; 2Þ, L 1 and L 2 are the linear operator associated to the domain Ω 1 and Ω 2 , respectively.
The boundary value problem 1-2 exhibits strong boundary layer at x ¼ 2 and interior layer at x ¼ 1. Expand y 0 ðx − δÞ about x using the Taylor's expansion and discard higher order terms.Then, the above problem can be approximated by where c ε;δ ðxÞ ¼ −ε − δdðxÞ and pðxÞ ¼ aðxÞ þ dðxÞ; As we observed from Eqns (3) and ( 4), the values of yðx − 1Þ is known for the domain Ω 1 and unknown for the domain Ω 2 due to the large delay at x ¼ 1.So, it is impossible to treat the problem throughout the domain ðΩÞ.Thus, we have to treat the problem at Ω 1 and Ω 2 separately.Eqns (3)-( 4) are equivalent to where LyðxÞ ¼
For each segment ½x i ; x iþ1 ; i ¼ 1; 2; . . .; N − 1 the non-polynomial cubic spline S Δ ðxÞ has the following form where a i ; b i ; c i and d i are unknown coefficients, and w ≠ 0 arbitrary parameter which will be used to increase the accuracy of the method.
To determine the unknown coefficients in (10) in terms of y i ; y iþ1 ; M i and M iþ1 first we define

Robust numerical method
The coefficients in (10) are determined as where θ ¼ wh.
Using the continuity condition of the first derivative at Reducing indices of Eqn ( 12) by one and substituting into Eqn (13), we obtain where If h → 0, then θ ¼ hk → 0. Thus, using L'Hopitals rule we have lim Using Taylorâ™s series expansions of y i−1 ; y iþ1 ; y 0 i−1 and y 0 iþ1 simplifying, we have where T 1 ¼ − h 2 6 y 000 ðξÞ and T 2 ¼ h 2 12 y 000 ðξÞ, for ξ Substituting Eqn (17) into Eqn ( 14) and rearranging, we get where, 12 y 000 ðξÞ is the local truncation error.From the theory of singular perturbations described in [18] and the Taylorâ™s series expansion of yðxÞ about the point '0' in the asymptotic solution of the problem in Eq. ( 9), we have Introducing a fitting factor σðρÞ in to Eq. ( 18), we get σðρÞc ε;δ ðxÞ Multiplying Eqn ( 19) by h and taking a limit as h → 0 we get σ ρ Thus, we consider two cases of the boundary layers.

Stability and convergence analysis 5.1 Truncation error
Let expand the terms y i±1 and M i±1 from Eqn ( 14), using Taylor's series as The local truncation error T i ðhÞ obtained from Eqn ( 14) as Substituting the series of y i±1 and M i±1 from Eqn (31)-( 32) and collecting like terms gives But from the values of α ¼ 1 6 and β ¼ 1 3 , Eqn (33) becomes where i .This establishes that the developed method is second order accurate or its order of convergence is Oðh 2 Þ.

Convergence analysis
Local truncation errors refer to the differences between the original differential equation and its finite difference approximation at a mesh points.Finite difference scheme is called consistent if the limit of truncation error ðT i ðhÞÞ is equal to zero as the mesh size h goes to zero.Hence, the proposed method in Eqn (26) with local truncation error in Eqn (34) satisfies the definition of consistency as Thus, the proposed scheme is consistent.

Stability analysis
Consider the developed scheme in Eqn (26), where the coefficients E i , F i and G i are as in Eqn (26).If we multiply both sides of Eqn (26) by h 2 and consider the values of E i , F i and G i for sufficiently small h, we get Considering Eqn (37) into Eqn (26) the one which is multiplied by h 2 the developed scheme can be written in a matrix form where the matrices Here, the coefficient matrix A is a tri-diagonal matrix with size ðN − 1Þ3ðN − 1Þ.Matrix A is irreducible if its codiagonals contain nonzero elements only.The codiagonals contain E i and G i .It is clearly seen that, for sufficiently small h both E i ≠ 0 and G i ≠ 0 for i ¼ 1; 2; . . .; N − 1.Hence, A is irreducible.
Again we can see that all jE i j, jFj i , jG i j > 0 for i ¼ 1; 2; . . .; N − 1 and in each row of A, the modulus of diagonal element is greater than or equal to the sum of the modulus of the two codiagonal elements (i.e.jF i j ≥ jE i j þ jG i j).This implies that A is diagonally dominant.Under this condition, the Thomas algorithm is stable for sufficiently small h.
As discussed in [19] the eigenvalues of a tri-diagonal matrix A are given by Hence, the eigenvalues of matrix A in Eqn (38) are But from trigonometric identity, we have 1 − cos sπ N ¼ 2sin 2 sπ 2N .Thus, the eigenvalues of A

Robust numerical method
A finite difference method for the boundary value problems is stable if A is nonsingular and A −1 ≤ C, for 0 < h < h 0 .where,C and h 0 are two constants that are independent of h.
Since A is real and symmetric it follows that A −1 is also real and symmetric so that, its eigenvalues are real and given by 1 λ s .Hence, as [19] the stability condition of the method will be satisfied when; , where, C is independent of h .Thus the developed scheme in Eqn (26) is stable.A consistent and stable finite difference method is convergent by [20].Hence as we have shown above, the proposed method is satisfying both the criteria of consistency and stability which are equivalent to convergence of the method.

Numerical examples and results
In this section, one example is given to illustrate the numerical method discussed above.The exact solutions of the test problem is not known.Therefore, we use the double mesh principle to estimate the error and compute the experiment rate of convergence to the computed solution.For this we put where Y N i and Y 2N 2i are the i th components of the numerical solutions on meshes of N and 2N, respectively.We compute the uniform error and the rate of convergence as The numerical results are presented for the values of the perturbation parameter ε ∈ f10 −4 ; 10 −8 ; . . .; 10 −20 g. −εy 00 ðxÞ þ 10y 0 ðxÞ À yðx À 1Þ þ y 0 ðx À εÞ ¼ x x ∈ ð0; 1Þ ∪ ð1; 2Þ; Subject to the boundary conditions yðxÞ ¼ 1; x ∈ ½−1; 0; yð2Þ ¼ 2:

Discussion and conclusion
This study introduces a fitted nonpolynomial spline method for singularly perturbed differential equations having both small and large delay.The numerical scheme is developed on uniform mesh using fitted operator in the given differential equation.The stability of the developed numerical method is established and its uniform convergence is proved.To validate the applicability of the method, one model problem is 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 Table 1).Further, behavior of the numerical solution (Figure 1), point-wise absolute error (Figure 2) and the ε -uniform convergence of the method is shown by the log-log plot (Figure 3).The method is shown to be ε -uniformly convergent with order of convergence Oðh 2 Þ.The proposed method gives accurate, stable and ε -uniform numerical result.

Example 6 . 1 .
Consider the model singularly perturbed boundary value problem: Figure 3. ε-uniform convergence with fitted operator in log-log scale for Example 6.1