convergence, finite difference method.

A differential-difference equation or delay differential equation (DDE) is an equation in which the system’s evolution at a specific time is contingent upon its state at a prior period. This differs from ordinary differential equations (ODEs), in which the derivatives rely solely on the present value of the independent variable. A delay differential equation is classified as singularly perturbed if the highest derivative is multiplied by a small parameter. Delay differential equations emerge in the mathematical modeling of diverse practical phenomena, such as microscale heat transfer[1], hydrodynamics of liquid helium[2], second-sound theory[3], thermoelasticity[4], diffusion in polymers[5], reaction-diffusion equations[6], stability[7], and various models for physiological processes or diseases[8, 9].

A delay differential equation (DDE) is classified as a retarded delay differential equation (RDDE) if the delay argument is absent from the highest order derivative term; otherwise, it is referred to as a neutral delay differential equation (NDDE). By confining it to a category where the largest derivative term is multiplied by a small parameter, we obtain singularly perturbed delay differential equations. This study focuses on second-order singularly perturbed retarded differential delay equations (RDDEs). DDEs are often transformed into ODEs with coefficients dependent on the delay by second-order accurate Taylor series expansions of the delay-involved terms. Asymptotic and numerical methods for delay differential equations (DDEs) are sometimes more complex to implement than those for ordinary differential equations (ODEs) due to the necessity of employing suitable approximations for the retarded arguments y(x− ) or y0(x− ) and their derivatives, as well as the requirement for the algorithm to address the discontinuities caused by the delay. For more details, one may suggest to go through some of the highly recognized monographs and articles: [10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23] and [24, 25, 9, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48] etc.

The literature reveals a limited number of studies addressing the treatment of a category of SPDDEs with small delays ( ) in both convection and reaction terms. The authors propose a three-term exponentially fitted upwind technique [49]. Rate of convergences are quadratic-accelerated using Richard-extrapolation. The researcher in [50] employed a cubic spline fitted method for small delay singularly perturbed differential equations. The Taylor series expansion has been employed to approximate the delay-containing terms, transforming the main problem into a singularly perturbed problem in ordinary differential equations. This is further processed to generate a tridiagonal scheme and solved by the Thomas algorithm. The theoretical study indicates that the scheme attains a second-order rate of convergence for larger values of the perturbation parameter (") and a first-order rate of convergence for smaller values. Results are expressed as maximum absolute point-wise errors. A Shishkin mesh [51] strategy has been proposed by authors of [52] for B-spline collocation techniques in SPDDEs with small delays, which are applied over the underlying domain of the problem. The proposed scheme is of almost second-order rate of convergence. A non-standard finite difference approach is described in [53] to address problems in both the left and right end boundary layers. Practical investigates demonstrate that, for fixed values of the grid points and perturbation parameter, as the delay term increases, the width of the boundary layer decreases. In the work [53], the authors devised a mid-point finite difference approach on a uniform mesh. The authors have examined the stability of the suggested numerical scheme utilizing the maximum principle and barrier functions for solution constraints. The results convergence is occurring uniformly at a linear rate of convergence. The researchers of the study [54] have developed three distinct numerical methods to address a category of SPDDEs that incorporate delayed arguments in convection and reaction terms.

Objectives: Inspired by the literature, our primary objective is to investigate a category of boundary value problems in SPDDEs with small delays in both the convection and reaction components using a straightforward and efficient exponentially fitted finite difference method. The scheme is straightforward and readily customizable by the computer.

This paper introduces an exponentially fitted, finite difference, three-term approach to address a category of Singularly Perturbed Differential-Difference Delay Equations (SPDDEs) with small delays in both the convection and response terms on a uniform mesh. The Taylor series expansion is initially utilized to approximate terms with delay arguments, followed by the introduction of a fitting factor, or artificial viscosity, to the fitted finite difference-based three-term scheme, resulting in the formation of a tridiagonal system. This system has been resolved with the renowned Thomas algorithm. The innovative aspect of the proposed strategy lies in the lack of prior research that addresses this type of problem using fitting factors on a uniform grid. This proposed method can yield highly precise results when the condition " > N is met. The results produced by the proposed method have been analyzed in terms of maximum absolute point-wise errors (MAPEs) and rates of convergence (ROCs).

The remaining portion of the paper is arranged as follows: Considered continuous problem is stated in Section:2. The methods of solution is described in the sections 3. Convergence analysis of the fitted method are analyzed in the Section: 4. Considered test example problems with the computational results are presented in the Section: 5. Conclusions are presented in the Section: 6. Paper ends with the bibliography.

We consider the following class of 2nd order singularly perturbed DDEs with a small delay in first derivative and undifferentiated terms:

Statement of the problem. To introduce the proposed exponentially fitted finite difference scheme, we examine second-order linear singularly perturbed differential equations with delay in the reaction and convection terms, which can be expressed as:

"u00 (w) + (w) u0 (w − ) + r (w) u (w − ) + ! (w) u (w) = f (w) (1) 8 w 2 (0, 1) = under the interval and the conditions:

u (w) = (w) on − w 0, andu (w) = (w) on1 w 1 + , (2) where is the small positive delay (negative shift) parameters to be order of 0 < = o(") << 1, "(0 < " << 1) is small singular perturbation parameter. The functions (w), r(w), !(w), f(w), (w)and (w) are sufficiently smooth in 8 w 2 (0, 1) = . The unique solution of equation (1) with (2) exhibits the layer behavior at the left end (neighbourhood of u = 0) of the interval if (w) − r (w) > 0 and at the right end (neighbourhood of u = 1) of the interval if (w) − r (w) < 0. If the shift parameter is smaller than epsilon, it is acceptable to use Taylor’s series expansion for the term that includes the shift argument. Hence, to obtain an accurate estimation of the term incorporating the delay parameter, the Taylor series expansion is utilized in the following manner: w0 (u − ) w0 (u) − w00 (u) + O( 2) and 0 2 00 3 (3) w (u − ) w (u) − w (u) + 2 w (u) + O( ) Substituting equation (3) into equation (1), we get an asymptotically equivalent singularly perturbed two point boundary value problem(BVP) of the following form:

µu00 (w) + (w) u0 (w) + (w) u (w) = f (w) , 8w 2 [0, 1] (4) with the BCs u (0) = (0) = 0, u (1) = (1) = 1, (5) with 0 < µ = " − (u) + 2 r(u) << 1, (w) = (w) − r (w) and (w) = r (w) + ! (w) where 0 and 1 are finite constants. Equation (1) is transitionable to equation (4) if is sufficiently small, improving computational efficiency. Further details on the validity of this transition can be found in Elsgolt’s and Norkin [55] and the articles [56, 57, 38]. Thus, the solution of equation (4) will provide a good approximation to the solution of equation (1).

For the differential equation given in (4), let us represent the differential operator L as Lu(w) = µu00(w) + (w)u0(w) + (w)u(w) [58].

Lemma 1. (Maximum Principle) Let %(w) be a differentiable function such that %(0) is greater than or equal to 0 and %(1) is greater than or equal to 0. Given L %(w) 0 for every w in the interval (0, 1), it follows that %(w) is greater than or equal to 0 for all w in the closed interval [0, 1].

Proof : See on [57]

Lemma 2. (Stability Principle) Let u(w) denote the solution to the problem described in equation (4). By combining it with equation (5), we obtain ku(w)k −1 kfk + max(| 0|, | |). (6) where k · k is the L1 norm given by ku(w)k = max |u(w)| and −1 is lower of (w) on [0, 1].

0 w 1 Proof : Assume that + and − are distinct barrier functions, and that ± (w) = −1 kfk+max(| 0|, | |)±u(w). Applying Lemma-1 to the comparison function (w) ± u(w) yields the desired estimate without delay. For more information, please refer to [57].

Lemma 3. The problem represented by equations (4) - (5) has a solution denoted as u(w) that meets

|u(w)| C 1 + exp

and (m) −m − w u (w) C 1 + µ exp , 8 m 1,

where is lower bound of (w).

Proof : See on [59] or [19].

This section describes the proposed method for the solution of the problem (2.4) with (2.5) which is equiva-lent/valid verson to the original problem (2.1) with (2.2).

To describe the method we assume that the function u(w) is smooth on (0, 1) and has the value ui at uniform mesh points wi = ih for i = 0, 1, 2, ..., N, where h = N 1 is the interval of differencing(mesh size).

Further, we let the coefficients of the equation (2.4) are evaluated at the midpoint of the each of the interval [wi−1, wi+1] and the solution of the equation (2.4) with (2.5) is given by:

− i(w µ−wi) (w − wi) u (w) = Ai + Bie − iui, wi−1 < w < wi+1 (7)

where Ai, Bi are arbitrary constants whose values are to be determined.

Using the notations u (wi−1) = ui−1, u (wi) = ui, u (wi+1) = ui+1, we have the following from the equation (3.1):

µ ih h − µih h ui−1 = Ai + Bie + iui, ui+1 = Ai + Bie − iui, ui = Ai + Bi

ih − ih

Obviously, ui−1 − 2ui + ui+1 = Bi e µ − 2 + e µ .

On solving the above relations for Ai and Bi, we get:

ih − h

u − 2u + u e µ + e µi ui − (ui+1 + ui−1)

i+1 i i−1

Bi = ih − ih , Ai = ih − ih e µ − 2 + e µ e µ − 2 + e µ

Substituting Ai and Bi in equation (3.1), we get:

ih − ih i(w−wi) i(w−wi)

e µ + e µ − 2e− µ ui e− µ − 1 (ui+1 + ui−1) w − wi u (y) = ih − ih + ih − ih − iui (8) e µ − 2 + e µ e µ − 2 + e µ i

From equation (3.2), we have:

ih ih ih

e− µ − e µ ui + e µ − 1 (ui+1 + ui−1) h u (wi−1) = ih − ih + iui (9) e µ − 2 + e µ i

ih ih ih

e µ − e− µ ui + e− µ − 1 (ui+1 + ui−1) h u (wi+1) = ih − ih − iui (10) e µ − 2 + e µ i Let = µ h , then using (3.3) and (3.4) we obtain:

(e i − e− i ) 2h u (wi−1) − u (wi+1) = (ui−1 − 2ui + ui+1) i − i + iui (11) (e − 2 + e ) i

Now using the Taylor’s series expansion procedure(u(wi+1) = ui+1 and u(wi−1) = ui−1), we have:

ui+1 = ui + hu0i + h 2! 2 u00i + h 3! 3 u000i + h 4! 4 u(4) i + h 5! 5 u(5) i + h 6! 6 u(6) i + h 7! 7 u(7) i + h 8! 8 u(8) i + O(h9) (12)

ui−1 = ui − hu0i + h 2! 2 u00i − h 3! 3 u000i + h 4! 4 u(4) i − h 5! 5 u(5) i + h 6! 6 u(6) i − h 7! 7 u(7) i + h 8! 8 u(8) i − O(h9) (13) u00 i+1 = u00 i + hu(3) i + h 2! 2 u(4) i + h 3! 3 u(5) i + h 4! 4 u(6) i + h 5! 5 u(7) i + h 6! 6 u(8) i + h 7! 7 u(9) i + h 8! 8 u(10) i + O(h11) (14)

u00 i−1 = u00 i − hu(3) i + h 2! 2 u(4) i − h 3! 3 u(5) i + h 4! 4 u(6) i − h 5! 5 u(7) i + h 6! 6 u(8) i − h 7! 7 u(9) i + h 8! 8 u(10) i − O(h11) (15)

Using the above expansion equations (3.6) to (3.9), we have

2h2 00 2h4 (4) 2h6 (6) 2h8 (8) 10 ui−1 − 2ui + ui+1 = ui + ui + ui + ui + O(h ) (16) 2! 4! 6! 8!

00 00 00 2h2 (4) 2h4 (6) 2h6 (8) 2h8 (10) 12 ui−1 − 2ui + ui+1 = ui + ui + ui + ui + O(h ) (17) 2! 4! 6! 8!

and 0 ui+1 − ui−1 2 ui = + O h (18)

h4 (6)

Substituting the expression for 12 ui from the equation (3.11) in the equation (3.10) and simplifying, we get

2 00 h2 00 00 00 2 (4) h6 (8) h4 (4) 2h8 (8) 10

ui−1 − 2ui + ui+1 = h u i + 30 u i−1 − 2u i + u i+1 − h ui − 360 ui + 12 ui + 8! ui + O(h ) Or, h2 00 00 00 h4 (4) ui−1 − 2ui + ui+1 = (u i−1 + 28u i + u i+1) + ui + R (19) 30 20

13h8 (8)

where R = − 302400 ui ( ) for 2 [wi−1, wi+1] .

Differentiating equation (2.4) twice, we obtain the following for ui :

(4) fi00 ifi0 i 2 (2 i 0 + i) 00 i ( i 0 + i) ( i 00 + 2 i 0) 0 i i 0 i 00 ui = − 2 + 2 − ui + 2 − ui + 2 − ui (20) µ µ µ µ µ µ µ µ From the equation (2.4), we get µu00 i+1 = − i+1u0 i+1 − i+1ui+1 + fi+1 (21) µu00 i = − iu0 i − iui + fi (22) µu00 i−1 = − i−1u0 i−1 − i−1ui−1 + fi−1 (23) The non-symmetric finite difference approximations for first order derivatives u0 i+1 and u0 i−1 are given by:

0 ui−1 − 4ui + 3ui+1 ui−1 − 2ui + ui+1 2 ui+1 = − + O h (24) 2h h 0 −3ui−1 + 4ui − ui+1 ui−1 − 2ui + ui+1 2 ui−1 = + + O h (25) 2h h Substituting the expressions for u00 i , u00 i+1, and u00 i−1 obtained by using the equations (3.12), (3.18) and (3.19) in the equations (3.16), (3.15) and (3.17) respectively, and the expression for ui from the equation (3.14) in the equation (3.13) and then simplifying, we get µ ui−1−2 h u 2i+ui+1 + 30 i−1 −3ui−1+4 2hui−ui+1 + ui−1−2h ui+ui+1 + 30 i−1 ui−1 + 28 30 i ui+12 − hui−1 +

28 i−1 i+1 i+1 ui−1−4ui+3ui+1 ui−1−2ui+ui+1 h2 i( i+ i) 00 0 30 ui + 30 ui+1 + 30 2h − h − 20 µ − ( i + 2 i) u −u h2 2 (2 0 i+ i) u −u h2 2 (2 0 i+ i) (26) i+12h i−1 + 20 i µ2 i − µ i+12h i−1 + 20 i µ2 i − µ −

2 0 i 2 0 2 2 2 0 + h i i 00 1 h 00 ifi h fi i ( i i) 20 µ − i ui = 30 [fi−1 + 28fi + fi+1] + 20 fi − µ + 20 µ2 − µ Now by making use of a constant fitting parameter ( ) in the above scheme (3.20), we obtain ( ) µ ui−1−2 h u 2i+ui+1 + 30 i−1 −3ui−1+4 2hui−ui+1 + ui−1−2h ui+ui+1 + 30 i−1 ui−1+

28 30 i ui+12 − hui−1 + 2830 i−1 ui + 30 i+1 ui+1 + 30 i+1 ui−1−42 uh i+3ui+1 − ui−1−2h ui+ui+1 −

0 2 0

h2 i( i+ i) 00 0 ui+1−ui−1 h2 i i (2 i+ i) ui+1−ui−1 20 µ − ( i + 2 i) 2h + 20 µ2 − µ 2h + (27)

2 2 0 0 0

h i i (2 i+ i) h2 i i 00 1 h2 00 ifi 20 µ2 − µ − 20 µ − i ui = 30 [fi−1 + 28fi + fi+1] + 20 fi − µ +

h2fi 2 i (2 0 i+ i) 20 µ2 − µ

where the values of the parameter ( ) is be determined.

Multiplying the above equation(3.21) by h and taking the limit as h ! 0, we get ( ) (0) 2 3 lim [ui+1 − 2ui + ui−1] = + ( (0)) lim [ui−1 − ui+1] (28) h!0 2 40 h!0 Now, an use of the equation (3.5) in (3.22) gives the required value of the fitting parameter ( ) as:

2 3 4 1 5 ( ) = (0) + ( (0)) (29) 2 20 tanh (0)

Finally, by making use of the equation (3.21) and the value of ( ) obtained in the equation (3.23), we get the following fitted three-term finite difference scheme:

Lh Eiui−1 − Fiui + Giui+1 = Ri, (i = 1, 2, 3, ..., N − 1) (30) with µ( ) i−1 i−1 28 i i+1 h i( 0 i+ i) 00 0 hi 2 i (2 0 i+ i) Ei = h2 − 60h + 30 − 60h − 60h + 40 µ − ( i + 2 i) − 40 µ2 − µ

2 2 2 0 + 2 h 0 i 2µ( ) 28 i h i i ( i i) h i i 00 Fi = h2 − 30 − 20 µ2 − µ + 20 µ − i

µ( ) i−1 i+1 28 i i+1 h i( 0 i+ i) 00 0 hi 2 i (2 0 i+ i) Gi = h2 + 60h + 30 + 60h + 60h − 40 µ − ( i + 2 i) + 40 µ2 − µ

0 2 2 0

1 h2 00 ifi h fi i (2 i+ i) Ri = 30 [fi−1 + 28fi + fi+1] + 20 fi − µ + 20 µ2 − µ

The scheme (3.24) with (2.5) provides a tri-diagonal algebraic system of (N − 1) equations with (N − 1) unknowns u1 to uN−1. We will solve this system by ‘Discrete Invariant Imbedding algorithm’ (also known as "Thomas Algorithm") described in [38, 60, 36]. It is known that the "Thomas algorithm” provides a numerically stable and convergent technique for solving the tridiagonal system of equations of the form(3.24)-(2.5) under the following conditions [61]:

Fi > (Ei + Gi), |Ei| |Gi|, Ei > 0, Gi > 0 (31) By considering the functions (w) = and (w) = as constant functions on [0,1] where and are finite constants, one can easily observe that these conditions (3.25) hold true for the problems whose solution either exhibits left layer or right layer behavior. In case of the problems whose solution exhibits left layer and right h i layer behaviour, all the conditions (3.25) are satisfied when − − h20 2 "2 − " > 0 where (w) = < 0.

In this section, we will discuss the convergence/error analysis of the proposed fitted method [62, 63].

Definition 1. Let i[u] Lhu(wi) − L u(wi), i = 1, 2, ..., N − 1, where u is a continuously differentiable function on I = [0, 1] with Lh as the discrete difference operator. Then the difference problem (3.24) with (2.5) is consistent with the difference problem (2.4) with (2.5) if

| i[u]| ! 0 as h ! 0.

The quantities i[u], i = 1, 2, ..., N − 1, are said to be the local truncation (or local discretization) errors.

Definition 2. The difference problem (3.24) with (2.5) has local rth− order accuracy if, for sufficiently smooth data, a positive constant C exists independent of h and " such that

max | i[u]| Chr.

1 i N

The consistency of the difference problem (3.24)-(2.5) with (2.4)-(2.5) and its locally second-order accuracy is shown by the following lemma.

Lemma 4. If u 2 C4(I), then

28h2 (3) 3 | i| max u (w) + O h ; i = 1, 2, ..., N − 1.

wi−1 w wi+1 180 proof. By definition n o n o i = ( )µ ui+1−2 h u 2i+ui−1 − u00 i + 30 i−1 −3ui−1+4 2hui−ui+1 + ui−1−2h ui+ui+1 − u0 i−1 + n o n o 28 30 i ui+12 − hui−1 − u0 i + 30 i+1 ui−1−42 uh i+3ui+1 − ui−1−2h ui+ui+1 − u0 i+1 .

n 2 4 o n 2 o

h (4) h (6) i−1 00 h (3)

) i = ( )µ 12 ui + 360 ui + ... + 30 hui − 3 ui + ... + n 2 4 o n 2 o

28 i h (3) h (5) i+1 00 h (3)

30 6 ui + 120 ui + ... + 30 −hui − 3 ui + ...

n 2 o ) | i| max h 12µ u(4)(w) + max h 30 |u00(w)| + wi−1 w wi+1 n 2 o wi−1 w wi+1 max 28 180 h u(3)(w) − max h 30 |u00(w)|

wi−1 w wi+1 wi−1 w wi+1

h2µ (4) 28h2 (3) ) | i| max u (w) + max u (w) wi−1 w wi+1 12 wi−1 w wi+1 180 Using the relation (3.23), we get

3 h(0)

h (0) coth 2µ h2 2(0)

) | i| 24 1 + 20µ2 max u(4)(w) + n 2 o wi−1 w wi+1 max 28 180 h u(3)(w)

wi−1 w wi+1

Using the expansion of coth h 2µ (0) , we get

h3 n o 28h2 ) | i| + O h5 max u(4)(w) + max u(3)(w) 24 wi−1 w wi+1 wi−1 w wi+1 180

28h2 (3) 3 ) | i| max u (w) + O(h ) wi−1 w wi+1 180

) | i| O h2 ; i = 1, 2, ..., N − 1. (32) Thus, the desired result is obtained.

Definition 3. (Uniform Convergence) Let u(w) be the exact solution and uN be the approximate solution obtained by any numerical method/scheme for the problem under consideration. If there exist a relation of the form ku − uN k Chr, where C > 0 and r > 0 are independent of µ and of the mesh h, then we say that the method/scheme is of r − th order accurate and converges uniformly to u(w) with respect to the norm k.k [64]. Theorem 1. Let u(w) be the analytical solution of the equation (2.4) with (2.5) [and hence (2.1) with (2.2)] and let uN (w) be the approximate solution obtained by the discretized scheme (3.24) with (2.5). Then for sufficiently small h, the error component satisfies the uniform error estimate: ku − uN k C h2 with C a positive constant.

Proof. Clearly the system of equations (3.24) has the following matrix-vector form [63]:

V u = D (33)

where V = (aij) , 1 i, j N − 1 is a tri-diagonal matrix of order (N − 1), with

µ( ) i−1 i+1 28 i i+1 h i( 0 i+ i) 00 0 hi 2 i (2 0 i+ i) ai,i+1 = 2 + 60h + 30 + 60h + 60h − 40 µ − ( i + 2 i) + 40 µ2 − µ

2 2 2 0 + 2 h 0 i 2µ( ) 28 i h i i ( i i) h i i 00 ai,i = − h2 − 30 − 20 µ2 − µ + 20 µ − i

µ( ) i−1 i−1 28 i i+1 h i( 0 i+ i) 00 0 hi 2 i (2 0 i+ i) ai,i−1 = h2 − 60h + 30 − 60h − 60h + 40 µ − ( i + 2 i) − 40 µ2 − µ

0 2 2 0

1 h2 00 ifi h fi i (2 i+ i) and D = (di) is a column vector with di = 30 [fi−1 + 28fi + fi+1] + 20 fi − µ + 20 µ2 − µ , for

i = 1, 2, ..., (N − 1), with the local truncation error i given by

| i| O h2 (34)

from the equation (4.1).

We also have V u − (h) = D (35) where u = (u0, u1, u2, ..., uN ) and (h) = ( 1(h), 2(h), ..., N (h)) denote the actual solution and the local truncation error respectively.

From equations (4.2) and (4.4), we have V (u − u) = (h) (36) Thus, the error equation is V E = (h) (37) where E = u − u = (e0, e1, e2, ..., eN ) .

Let Si be the sum of elements of the ith row of V, then we have N P −1 2 2 0 + h 2 i −µ( ) (28 1+ 2) ( 0+28 1+ 2) 1 ( 1 1) h 1 h1 S1 = a1,j = h2 + 30 + 60h + µ2 − µ 20 + 40 − 2 h 0 i 0 + h 1 1 00 h 1( 1 1) 00 0 20 µ − 1 − 40 µ − ( 1 + 2 1) N P −1 2 h 0 i SN−1 = aN−1,j = −µ h2( ) + (28 N−30 1+ N−2) − ( N−2+28 60h N−1+ N ) + h 20 N−1 µ N−1 − N 00 −1 −

2 2 0 + h 2 i 0 + N−1 ( N−1 N−1) h N−1 hN−1 h N−1( N−1 N−1) 00 0 µ2 − µ 20 + 40 − 40 µ − N−1 + 2 N−1 N P −1 1 2 Si = ai,j = 30 [ i−1 + 28 i + i+1] + O(h )

= bi = Bi0 ; i = 2 (1) N − 2; where Bi0 = bi = 30 1 [ i−1 + 28 i + i+1] .

Since 0 < << 1, the matrix V is irreducible and monotone for sufficiently small h. Thus V −1 exists with the non-negative elements. Hence, from equation(4.6) we obtain

E = V −1 (h) (38)

||E|| ||V −1|||| (h)|| (39) Let aki be the (ki)th elements of V −1. Since aki 0, from the operations of matrices we have.

N X −1 akiSi = 1; k = 1, 2, ..., N − 1 (40) Therefore, it follows that N X −1 1 1 1 aki = (41) i=1 min Si Bi0 |Bi0 |

0 i N−1

for some i0 between 1 and (N − 1) and Bi0 = bi.

Therefore, from equations(4.3),(4.7) and (4.9) we obtain N X −1 ej = aki i(h); j = 1(1)N − 1

which implies O h2 ej ; j = 1(1)N − 1 (42) |bi| By making use of the definitions and the equation(4.11) we obtain ||E|| = o(h2) Thus the method proposed is of second order convergence.

Here, we will demonstrate the performance of the present method by solving four test example problems; presenting the results in terms of maximum absolute pointwise errors(MAPEs) and the numerical order of convergence(NOC). We will also present the graphical solutions with varying values of the delay parameter to examine the effect of small delay on the left layer and right layer behavior of the solutions. Since the exact solution of the considered test example problems are not available, the well known double mesh principle [14] is used to calculate the MAPEs and NOC. The formula for finding MAPEs (E" N ) is given by:

E" N = max uN i − u2N i , i = 0, 1, 2, ..., N,

0 i N

where uN i and u2N i are the numerical solutions obtained on the mesh wi = w0 + ih, i = 0, 1, 2, ..., N and on the mesh obtained by halving of the original mesh size h with 2N number of mesh intervals respectively. Furthermore,

N N E" N

the NOC (R" ) is obtained using the formula: R" = log2 E2N .

Example problems exhibiting left layer behaviour. In this subsection we will solve two SPDDEs whose solutions exhibite the left layer behaviour in [0, 1].

Example 1. Firstly, we consider the following problem "u00(w) + (1 + w)u0(w − ) + exp(−2w)u(w − ) − 2 exp(−w)u(w) = 0; with u(0) = 1 and u(1) = 0 Example 2. We consider the following problem "w00(u) + (1 + u)w0(u − ) + sin(2u)w(u − ) − exp(−u)w(u) = sin(2u) + 3 exp(−u); with w(0) = 1 and w(1) = 1

We solved these two example problems for various values of , " and N, and tabulated the computational results in the Tables 1 to 2. Table 1 and 2 presents the results in terms of MAPEs for different N and " with = 0.5 ". Further, a simple look on the numerical order of convergence(NOC) (R" N ) presented in the Table 1 for example-1 and the Table 2 for example-2 clearly reveals that the method proposed is of the second order convergence rate. In order to examine the effect of small delay ( ) on the left layer behaviour of the solutions for the example-1 and 2 we plotted the graphs for " = 0.1, 0.01 in the figures 1 to 2 and 4 to 5 respectively with fixed N = 100 and the varying values of where = o("). By the simple observation of the figures 1 to 2 and 4 to 5 we easily confirm that there is no significant effect of small delay on the left layer behavior of the solution when = o(") with " small. That is, when the values of " tends to zero the effect of small delay on the left layer behavior becoming negligible for all = o("). Thus the effect of small delay can be ignored in case when = o(") with " sufficiently small. For the example-1 qnd 2, the effect of the delay is to increase and decrease the width of the left boundary layers. Plotted the Loglog graphs for the MAPEs (E" N ) for different N with various " in the Figures 3 and 6. From these Figures 3 and 6, one can easily visualize the higher accuracy of the suggested method for the increasing " more clearly.

Example problems exhibiting right layer behaviour. We performed the numerical experiments on the following two SPDDEs type whose solutions exhibite right layer behaviour in the interval [0, 1]:

Example 3. We consider the following problem −"u00(w) + (1 + w)u0(w − ) − exp(−2w)u(w − ) + exp(−w)u(w) = 0; with u(0) = 1 and u(1) = −1 Example 4. Finally, consider the following problem −"u00(w) + (1 + w)u0(w − ) − exp(−2w)u(w − ) + exp(−w)u(w) = exp(x − 1); with u(0) = 1 and u(1) = −1 We solved these two example problems for various values of , " and N, and tabulated the computational results in the Tables 3 to 4. Table 3 and 4 presents the results in terms of MAPEs for different N and " with = 0.5 ". Further, a simple look on the numerical order of convergence(NOC) (R" N ) presented in the Table 3 for example-3 and the Table 4 for example-4 clearly reveals that the method proposed is of the second order convergence rate. In order to examine the effect of small delay ( ) on the right layer behaviour of the solutions for the example-3 and 4 we plotted the graphs for " = 0.1, 0.01 in the figures 7 to 8 and 10 to 11 respectively with fixed N = 100 and the varying values of where = o("). By the simple observation of the figures 7 to 8 and 10 to 11 we example 2 for different value of example 2 for different value of

easily confirm that there is no significant effect of small delay on the right layer behavior of the solution when = o(") with " small. That is, when the values of " tends to zero the effect of small delay on the right layer behavior becoming negligible for all = o("). Thus the effect of small delay can be ignored in case when = o(") with " sufficiently small. For the example-3 and 4, the effect of the delay is to decrease and increase the width of the right boundary layers. Plotted the Loglog graphs for the MAPEs (E" N ) for different N with different " in the Figures 9 and 12. From this Figures 9 and 12, one can easily visualize the higher accuracy of the sugested method for the increasing " more clearly.

We developed a new method for finding the solution of the considered SPDDEs exhibiting left-layer or right-layer behavior and implemented it on four standard widely discussed problems. To show the applicability, accuracy, and efficiency of the method, we tabulated the results in terms of MAPEs (maximum absolute pointwise errors) and the NOC (numerical order of convergence). Both the theory and the numerical computation state that example 3 for different value of example 3 for different value of with N = 100, " = 0.1 with N = 100, " = 0.01

the present method is capable of approximating the solution very well by producing second-order convergent solutions. Furthermore, we plotted the solution graphs for various values of , while keeping N and " fixed, to examine how the delay affects the behavior of the solution in both the left and right layers. One can easily observe from figures 1 to 2 and 4 to 5, for example, problems 1 and 2, respectively, that there is no significant effect of small delay on the left layer behavior of the solution when = o(") with " sufficiently small. As the values of " decrease, the effect of a small delay on the behavior of the left layer of the solution becomes negligible. Also, the effect of the delay is to increase and decrease the width of the left boundary layers. Plotted the Loglog graphs for the MAPEs (E" N ) for different N with various " in the Figures 3 and 6. From these Figures 3 and 6, one can easily visualize the higher accuracy of the suggested method for the increasing " more clearly. In case of problems with the right layer behavior of the solutions, an examination of figures 7 to 8 and 10 to 11 for example problems 3 and 4, respectively, shows there is no significant effect of small delay on the right layer behavior of the solution when = o(") with " sufficiently small. As the values of " decrease, the effect of small delay on the right layer behavior of the solution becomes negligible. Furthermore, the effect of the delay is to decrease and increase the width of the right boundary layers. Plotted the Loglog graphs for the MAPEs (E" N ) for different N with different " in the Figures 9 and 12. From this Figures 9 and 12, one can easily visualize the higher accuracy of the sugested method for the increasing " more clearly. Finally, it is concluded that the current method seems to be one of the best options for approximating the solution of the considered class of SPDDEs.

of example 4 for different value of of example 4 for different value of

The authors would like to express their sincere gratitude to the editors and anonymous reviewers for their invaluable comments and constructive feedback, which significantly contributed to the enhancement of this paper.

A. Kumar: conceptualization and editing. B. Prasad: methodology and software. R. Ranjan: writing and editing. All authors have read and agreed to the published version of the manuscript.

The authors declare that no external funding or support was received for the research presented in this paper, including administrative, technical, or in-kind contributions. Data Availability All data supporting the reported findings in this research paper are provided within the manuscript.

The authors declare that there is no conflict of interest concerning the reported research findings. Funders played no role in the study’s design, in the collection, analysis, or interpretation of the data, in the writing of the manuscript, or in the decision to publish the results.