Int. J. Anal. Appl. (2022), 20:63 Numerical Simulation of Singularly Perturbed Delay Differential Equations With Large Delay Using an Exponential Spline Ramavath Omkar, K. Phaneendra∗ Department of Mathematics, University College of Science, Osmania University, India ∗Corresponding author: kollojuphaneendra@yahoo.co.in Abstract. In this study, numerical solution of a differential-difference equation with a boundary layer at one end of the domain is suggested using an exponential spline. The numerical scheme is developed using an exponential spline with a special type of mesh. A fitting parameter is inserted in the scheme to improve the accuracy and to control the oscillations in the solution due to large delay. Convergence of the method is examined. The error profiles are represented by tabulating the maximum absolute errors in the solution. Graphs are being used to show that how the fitting parameter influence the layer structure. 1. Introduction Differential-difference equations (DDEs) are ones in which a state variable’s time evolution is in- consistently dependent on a particular history, i.e., a physical system’s rate of change is dependent not only on its current condition but also on its prior history. These equations have been widely utilized in population dynamics [7], nonlinear delay differential equations relating to physiological control systems [14], red blood cell system [13], predator-prey models [15], neuronal variability problems connected to patterns of nerve action potentials formed by unit quantal inputs occurring at random are studied [22]. Bellman and Cooke [1], Doolan et al. [2] and Driver [3] just are a few of the authors who have produced papers and books in recent years explaining various methods for solving differential-difference equations with perturbation. In [4] authors solved singularly perturbed delay differential equation (SPDE) using a fitted scheme on an uniform mesh. In [5], Glizer solved linear quadratic optimal control problem with Received: Oct. 15, 2022. 2010 Mathematics Subject Classification. 65L11, 65L12. Key words and phrases. exponential spline; differential-difference equation; large delay; truncation error. https://doi.org/10.28924/2291-8639-20-2022-63 ISSN: 2291-8639 © 2022 the author(s). https://doi.org/10.28924/2291-8639-20-2022-63 2 Int. J. Anal. Appl. (2022), 20:63 delay by Hamiltonian boundary-value problem, with boundary function method. In [6], Kadalbajoo et. al. proposed numerical research utilizing finite difference techniques. The authors in [8] developed a technique for solving SPDEs with twin layers or oscillatory behaviour using a computational scheme. The authors in [9] used a nonpolynomial spline to develop a numerical solution for a DDEs having layer with a small and large delay in the differentiated terms. Using domain decomposition, the authors of [10] suggested a mixed difference technique to solve DDEs with mixed shifts. For a class of linear second-order DDE type, Lange and Miura [11, 12] developed an asymptotic method. They developed a mathematical model by random synaptic inputs in dendrites for estimating the approximate time for the activation of action potentials in nerve cells and also discussed issues with solutions that had rapid oscillations in their study. The authors of [16] utilized a quadrature method for solving the SPDE, as well as a two-point quadrature rule to obtain a tridiagonal system. The authors in [17] developed a finite difference approach for solving SPDEs with turning points and mixed shifts. A finite difference scheme on Shishkin mesh is proposed to solve singularly perturbed delay differential equations with or without a turning point in [18]. The authors of [19, 20] used tension splines and an exponentially fitted spline for the problem with convection delay-dominated diffusion equation. The authors of [23] solved a two parameter semi linear differential equation by using an exponential spline. The authors of [21] solved SPDE with a numerical integration technique. 2. Statement of the problem Consider the continuous problem of a second order SPDE εz′′ + p(u)z′(u −δ) + q(u)z(u) = r(u) (2.1) with boundary conditions z(u) = ϕ(u) −δ ≤ u ≤ 0 , z(1) = τ (2.2) where 0 < ε << 1 is a perturbation, p(u),q(u), r(u) and ϕ(u) are bounded continuous functions on (0, 1), and τ is finite constant and δ is delay parameter. It is well known that when δ = 0 the above delay differential equation is reducing to a SPDE. The solution z(u) exhibit boundary layer on left side when p(u) is positive or on right side when p(u) is negative throughout the interval [0, 1]. If δ(ε) is of order O(ε), layer profile of the solution is no longer retained and the solution oscillates. 3. An exponential spline The region [υ, Υ] is subdivided into l equal subregions of mesh size h = (υ−Υ) l so that ui = υ + ih, i = 0, 1, . . . , l are the mesh points. To manage the delay term, the mesh size is chosen as h = δ ν , where ν= mantissa of δ which is a positive integer. Let the exact solution be z(u) to Eq. (2.1) and ui be an approximation solution to the z(ui ) attained Int. J. Anal. Appl. (2022), 20:63 3 by the segment Qi (u) passing (ui,Qi ) and (ui+1,Qi+1). Each exponential spline segment Qi (u) has the form. Qi (u) = aie k(u −ui ) + bie−k(u −ui ) + ci (u −ui ) + di f or i = 0, 1, 2, . . . , l. (3.1) where ai,bi,ci and di are the constants and k is a free parameter. To find the values of the coefficients in Eq. (3.1), the interpolatory conditions Qi (ui ) = Si,Qi (ui+1) = Si+1,Q ′′ i (ui ) = Mi,Q ′′ i (ui+1) = Mi+1 are used. Using these conditions, we get the following expressions ai = h2(Mi+1−e−θMi ) 2θ2sinh(θ) , bi = h2(Mi−e−θMi+1) 2θ2sinh(θ) , ci = (Si+1−Si ) h − h(Mi+1−Mi ) θ2 and di = Si − h2Mi θ2 where θ = kh and i = 0, 1, 2, . . . , l. Using the first derivative continuity i.e., Q′i−1 (u) = Qi ′ (u) at the point (ui,Qi ), we get the relation (zi+1 − 2zi + zi−1) = h2 (αMi+1 + βMi + αMi−1) f or i = 0, 1, . . . , l − 1 (3.2) where α = ( sinh(θ)−θ θ2sinh(θ) ) and β = 2θcosh(θ)−2sinh(θ) θ2sinh(θ) . 4. Numerical Approach Eq. (2.1) can be discretized at the mesh point uj by εMj = r ( uj ) −p ( uj ) z′j−ν −q ( uj ) zj f or j = i − 1, i, i + 1 (4.1) Using the following difference approximations of z′: z′i−ν ≈ [ 1 + 2ωh2qi+1 + ωh [3 pi+1 + pi−1] 2h ] zi−ν+1 − 2ω [pi+1 + pi−1] zi−ν − [ 1 + 2 ωh2 qi−1 −ωh [pi+1 + 3 pi−1] 2h ] zi−ν−1 + ωh [ri+1 − ri−1] z′i−ν+1 ≈ 3zi−ν+1 − 4zi−ν + zi−ν−1 2h , z′i−ν−1 ≈ −zi−ν+1 + 4zi−ν − 3zi−ν−1 2h Now substituting the above approximations in Eq. (4.1) and in Eq. (3.2) we obtained.( ε h2 + αqi+1 ) zi+1 + ( −2ε h2 + βqi ) zi + ( ε h2 + αqi−1 ) zi−1+ ( −αpi−1 2h + βpi 2h ( 1 + 2ωh2qi+1 + ωh (3pi+1 + pi−1) ) + 3α 2h pi+1)zi−ν+1+ ( 2αpi−1 h − 2βωpi (pi+1 + pi−1) − 2α h pi+1)zi−ν+ ( −3αpi−1 2h − βpi 2h ( 1 + 2ωh2qi−1 −ωh (pi+1 + 3pi−1) ) + α 2h pi+1)zi−ν−1 = (αri+1 + βri + αri−1) −βωpih (ri+1 − ri−1) (4.2) 4 Int. J. Anal. Appl. (2022), 20:63 When the shift parameter is small in comparison to the perturbation parameter, the layer behavior of the solution is maintained and reliable results are achieved. The layer behavior, however, no longer holds good if δ(ε) is of order O(ε) and oscillations manifest. We are therefore constructing a numerical approach with fitting parameter based on an exponential spline method and a particular mesh type. In order to demonstrate how significant the fitting parameter is in our suggested approach, we also examine how the solution layer behaves when there are large delays. Now inserting the fitting parameter σ in the above scheme, we have (σ� + h2αqi+1)zi+1 + (−2σ� + h2βqi )zi + (σ� + h2σqi−1)zi−1+ ( −αhpi−1 2 + βhpi 2 (1 + 2ωh2qi+1 + ωh(3pi+1 + pi−1)) + 3αh 2 pi+1)zi−ν+1+ (2αhpi−1 − 2βωpih2(pi+1 + pi−1) − 2αhpi+1)zi−ν+ ( −3αhpi−1 2 − βhpi 2 (1 + 2ωh2qi−1 −ωh(pi+1 + 3pi−1)) + αh 2 pi+1)zi−ν−1 = h2(αri+1 + βri + αri−1) −βωpih3(ri+1 − ri−1) (4.3) Re write the above equation Ai+1zi+1 + Bizi + Ci−1zi−1 + Di−ν+1zi−ν+1 + Ei−νzi−ν + Fi−ν−1zi−ν−1 = Gi (4.4) where Ai+1 = σ� + αh 2qi+1, Bi = −2σ� + βh2qi, Ci−1 = σ� + αh2qi−1 Di−ν+1 = −αhpi−1 2 + βhpi 2 (1 + 2ωh2qi+1 + ωh(3pi+1 + pi−1)) + 3αh 2 pi+1 Ei−ν = 2αhpi−1 − 2βωpih2(pi+1 + pi−1) − 2αhpi+1 Fi−ν−1 = −3αhpi−1 2 − βpih 2 (1 + 2ωh2qi−1 −ωh(pi+1 + 3pi−1)) + αh 2 pi+1 Gi = h 2(αri+1 + βri + αri−1) −βωpih3(ri+1 − ri−1) p(ui ) = pi,q(ui ) = qi, r(ui ) = ri The above scheme can be written by using the boundary conditions Ai+1zi+1 + Bizi + Ci−1zi−1 = Gi −Di−ν+1zi−ν+1 −Ei−νzi−ν −Fi−ν+1zi−ν+1,∀ 1 ≤ i ≤ ν − 1 Ai+1zi+1 + Bizi + Ci−1zi−1 + Di−ν+1zi−ν+1 = Gi −Ei−νzi−ν − Fi−ν−1zi−ν−1, ∀ i = ν Ai+1zi+1 + Bizi + Ci−1zi−1 + Di−ν+1zi−ν+1 + Ei−νzi−ν = Gi −Fi−ν−1zi−ν−1, ∀ i = ν + 1 Ai+1zi+1 +Bizi +Ci−1zi−1 +Di−ν+1zi−ν+1 +Ei−νzi−ν +Fi−ν−1zi−ν−1 = Gi ∀ ν+ 2 ≤ i ≤ l−1. (4.5) By using Gauss elimination method with partial pivoting or any other method resolves the above system of l equations. Int. J. Anal. Appl. (2022), 20:63 5 4.1. Left-end layer. Assume p(u) ≥ M > 0 and q(u) ≤ −⊆ < 0 where θ and M are positive con- stants. Then the problem (2.1) - (2.2) shows layer structure at u = 0 for small values of perturbation. Now, to enhance the scheme’s effectiveness, the fitting parameter inserted in Eq. (4.4) is determined by using the following approximate solution given in [2]. limh→0z(ih) ≈ z0(ih) + (φ(0) −z0(0))exp(−p(0)iρ) + O(ε), where ρ = h � . Using this expression in Eq. (4.4) and following the procedure given in [2], we get σ = ρ(α + β/2)p(0)ep(0)νρcoth( p(0)ρ 2 ) Lemma 4.1. Suppose Ψ0 ≥ 0 and Ψl ≥ 0, then LlΨi ≤ 0 ∀i = 1, 2, 3, . . . .., l − 1 implies that Ψi ≥ 0 ∀i = 0, 1, 2, 3. . . ..l. Proof. Let j ∈ (0, 1, . . . .l) such that Ψj = max0≤i≤lΨi. Let if possible Ψj < 0. We will show that it leads to contradiction. Clearly j /∈ (0, l). Now we have Ll Ψj =   ερi Ψ ′′ i + (2α + β)qi Ψi = Gi −Di−ν+1Ψi−ν+1 −Ei−νΨi−ν −Fi−ν−1Ψi−ν−1, ∀i = 1, 2, . . . ,ν. ερi Ψ ′′ i + (2α + β)qi Ψi + Di−ν+1Ψi−ν+1 = Gi −Ei−νΨi−ν −Fi−ν−1Ψi−ν−1, ∀ i = ν + 1, . . . ., l1 − 1. ερi Ψ ′′ i + (2α + β)qi Ψi + Di−ν+1Ψi−ν+1 + Ei−νΨi−ν = Gi −Fi−ν−1Ψi−ν−1, ∀i = l1, . . . .., l − 1. (4.6) Case 1: For i = 1, 2, . . . ,ν LlΨj = (ερi Ψ ′′ i + (2α + β)qi Ψi = Gi −Di−ν+1Ψi−ν+1 −Ei−νΨi−ν −Fi−ν−1Ψi−ν−1, ∀i = 1, 2, . . . ,ν.) ≥ 0(∵ qi ) < 0, i = 1, 2, . . . .ν Case 2: For i = ν + 1, . . . ., l1 − 1 LlΨj = (ερi Ψ ′′ i + (2α + β)qi Ψi + Di−ν+1Ψi−ν+1 = Gi −Ei−νΨi−ν −Fi−ν−1Ψi−ν−1, ∀ i = ν + 1, . . . ., l1 − 1.) ≥ 0 Case 3: For i = l1, . . . .., l − 1 LlΨj = (ερi Ψ ′′ i + (2α + β)qi Ψi + Di−ν+1Ψi−ν+1 + Ei−νΨi−ν = Gi −Fi−ν−1Ψi−ν−1, ∀i = l1, . . . .., l − 1.) ≥ 0 As a result, we have LlΨj > 0 , which is contradicts the hypothesis that LlΨj ≤ 0, 1 ≤ j ≤ l − 1. As a consequence, we assume that Ψj < 0 is incorrect and that Ψj > 0 is correct. Since j was chosen an arbitrary, we have Ψj ≥ 0,∀j = 0, 1, 2, . . . , l (4.7) � 6 Int. J. Anal. Appl. (2022), 20:63 Theorem 4.1. Let p(u) ≥ M > 0 and q(u) ≤ −θ < 0 where M and θ are both positive constants. The problem (2.1) - (2.2) with boundary conditions has a uniqueness, existence and satisfying solution. ||Z||h,∞ ≤ ||r||h,∞θ−1 + C2(||ϕ||h,∞ + τ) (4.8) where C2 ≥ 1 is a positive constant. Proof. To demonstrate the uniqueness and existence, assume li=0 and l i=0 be the two solutions to the discrete problem (2.1) - (2.2). Then ti = vi −wi is a mesh function that meets the conditions t0 = 0 = tl and for 1 ≤ i ≤ l − 1 we have Llti = Llvi −Llwi. Since vi and wi satisfy Eq. (4.5), therefore Llti = 0, 1 ≤ i ≤ l−1. As a result, the mesh function ti meets the discrete minimum principle hypothesis, and when we apply it to the mesh function ti, we get ti = (vi −wi ) ≥ 0 f or 0 ≤ i ≤ l (4.9) Again, if we set ti = −(vi −wi ) , then ti is a mesh function that satisfies t0 = 0 = tl and we have Llti = 0, 1 ≤ i ≤ l − 1 along the same lines as before. As a result, the discrete minimum principle is applied to the mesh function ti gives ti = −(vi −wi ) ≥ 0 i.e., (vi −wi ) ≤ 0 f or 0 ≤ i ≤ l (4.10) From Eq. (4.9) and Eq. (4.10), we get vi − wi = 0 . This suggests that the discrete problem (2.1)-(2.2) solution is unique. The existence of linear equation is assumed by their uniqueness. Now we’ll explain how to prove the bound on < zi >li = 0. We do this by introducing two barrier functions Ψ± i Defined by Ψ± i = ||r||h,∞θ−1 + C2(||ϕ||h,∞ + |τ|) ±Zi, 0 ≤ i ≤ l where C2 ≥ 1 is a positive constant that can be chosen at random. Then we have Ψ± i = ||r||h,∞θ−1 + C2(||ϕ||h,∞ + |τ|) ±Zi, 0 ≤ i ≤ l Now we have Ψ±0 = ||r||h,∞θ −1 + C2(||ϕ||h,∞ + |τ|) ±Z0, Ψ±0 = (||r||h,∞θ −1 + C2(||ϕ||h,∞ + |τ|) ±C2ϕ0) ≥ 0 since ||ϕ||h,∞ + |τ| ≥ ϕ0 and C2 ≥ 1, Ψ±n = ||r||h,∞θ −1 + C2(||ϕ||h,∞ + |τ|) ±Zn Ψ±n = (||r||h,∞θ −1 + C2(||ϕ||h,∞ + |τ|) ±τ) ≥ 0 Case 1. For i = 1, 2, 3, . . . .,ν LlΨ± i = ερi (Ψ ± i )′′ + (2α + β)qi Ψ ± i = Gi −Di−ν+1Ψ±i−ν+1 −Ei−νΨ ± i−ν −Fi−ν−1Ψ ± i−ν−1 = −qi (||r||h,∞θ−1 + C2(||ϕ||h,∞ + |τ|) ±LlZi ) Int. J. Anal. Appl. (2022), 20:63 7 = ( −qi θ ||r||h,∞ ± ri ) −qiC2(||ϕ||h,∞ + |τ|) ≤ 0 Case 2. For i = ν + 1, . . . ., l1 − 1 LlΨ± i = ερi (Ψ ± i )′′ + (2α + β)qi Ψ ± i + Di−ν+1Ψ ± i−ν+1 = Gi −Ei−νΨ ± i−ν −Fi−ν−1Ψ ± i−ν−1 = −qi (||r||h,∞θ−1 + C2(||ϕ||h,∞ + |τ|) ±LlZi ) = ( −qi θ ||r||h,∞ ± ri ) −qiC2(||ϕ||h,∞ + |τ|) ≤ 0 Case 3. For i = l1, . . . .., l − 1 LlΨ± i = ερi (Ψ ± i )′′ + 2α + βqi Ψ ± i + Di−ν+1Ψ ± i−ν+1 + Ei−νΨ ± i−ν = Gi −Fi−ν−1Ψ ± i−ν−1 = −qi (||r||h,∞θ−1 + C2(||ϕ||h,∞ + |τ|) ±LlZi ) = ( −qi θ ||r||h,∞ ± ri ) −qiC2(||ϕ||h,∞ + |τ|) ≤ 0 Combining above cases, we have LlΨ± i ≤ 0, 1 ≤ i ≤ l. Using the discrete maximum principle Ψ± i = ||r||h,∞θ−1 + C2(||ϕ||h,∞ + |τ|) ±Zn ≥ 0, 0 ≤ i ≤ l. which proves the desired results. � 4.2. Right-end layer. Assume p(u) ≤ −M < 0 and q(u) ≤ −⊆ < 0 where θ and M are positive constants. Then the problem (2.1) - (2.2) shows layer structure at u = 1 for small values of ε. In this case, the fitting parameter inserted in Eq. (4.4) is determined by using the following approximate solution given in [2]. limh→0z(ih) ≈ z0(0) + (φ(0) −z0(1))e−p(1)( 1 ε −iρ) + O(ε), where ρ = h � . Using this expression in Eq. (4.4) and following the process given in [2], we get σ = ρ(α + β/2)p(0)e−p(0)νρcoth( p(1)ρ 2 ) (4.11) Lemma 4.2. . Suppose Ψ0 ≥ 0 and Ψl ≥ 0, then LlΨi ≤ 0 ∀ i = 1, 2, 3, . . . ..l − 1 implies that Ψi ≥ 0 ∀ i = 0, 1, 2, 3. . . ..l. Proof. Let j ∈ (0, 1, . . . .l) such that Ψj = min0≤i≤lΨi. Let if possible Ψj < 0. We will show that it leads to contradiction. Clearly j /∈ (0, l). Now we have Ll Ψj =   ερi Ψ ′′ i + (2α + β)qi Ψi = Gi −Di−ν+1Ψi−ν+1 −Ei−νΨi−ν −Fi−ν−1Ψi−ν−1, ∀i = 1, 2, . . . ,ν. ερi Ψ ′′ i + (2α + β)qi Ψi + Di−ν+1Ψi−ν+1 = Gi −Ei−νΨi−ν −Fi−ν−1Ψi−ν−1, ∀ i = ν + 1, . . . ., l1 − 1. ερi Ψ ′′ i + (2α + β)qi Ψi + Di−ν+1Ψi−ν+1 + Ei−νΨi−ν = Gi −Fi−ν−1Ψi−ν−1, ∀i = l1, . . . .., l − 1. (4.12) 8 Int. J. Anal. Appl. (2022), 20:63 Case 1: For i = 1, 2, . . . ,ν LlΨj = (ερi Ψ ′′ i + (2α + β)qi Ψi = Gi −Di−ν+1Ψi−ν+1 −Ei−νΨi−ν −Fi−ν−1Ψi−ν−1, ∀i = 1, 2, . . . ,ν.) ≥ 0 (∵ qi < 0, i = 1, 2, . . . .ν) Case 2: For i = ν + 1, . . . ., l1 − 1 LlΨj = (ερi Ψ ′′ i + (2α + β)qi Ψi + Di−ν+1Ψi−ν+1 = Gi −Ei−νΨi−ν −Fi−ν−1Ψi−ν−1, ∀ i = ν + 1, . . . ., l1 − 1.) ≥ 0 Case 3: For i = l1, . . . .., l − 1 LlΨj = (ερi Ψ ′′ i + (2α + β)qi Ψi + Di−ν+1Ψi−ν+1 + Ei−νΨi−ν = Gi −Fi−ν−1Ψi−ν−1, ∀i = l1, . . . .., l − 1.) ≥ 0 which is contradiction to LlΨj ≤ 0, 1 ≤ j ≤ l − 1 Therefore, the assumption Ψj < 0 is false and Ψj > 0 . Since j was choosing arbitrary, we have Ψj ≥ 0,∀j = 0, 1, 2, . . . l Theorem 4.2. Let p(u) ≥ M > 0 and q(u) ≥ θ > 0 where M and θ are both positive constants. The solutions to the problem (2.1) - (2.2) with boundary conditions exist, are unique and satisfy. ||Z||h,∞ ≤ ||r||h,∞θ−1 + C1(||ϕ||h,∞ + τ) (4.13) where C1 ≥ 1 is a positive constant Proof. To prove the existence and uniqueness, assume li=0 and l i=0 be the two solutions to the problem (2.1) - (2.2). Then ti = vi −wi is a mesh function that meets the conditions t0 = 0 = tl and for 1 ≤ i ≤ l − 1 we have Llti = Llvi − Llwi. Since vi and wi satisfy Eq. (4.5), therefore Llti = 0, 1 ≤ i ≤ l − 1. As a result, the mesh function ti meets the discrete minimum principle hypothesis, and when we apply it to the mesh function ti, we get ti = (vi −wi ) ≥ 0 f or 0 ≤ i ≤ l (4.14) Again, if we set ti = −(vi − wi ) , then ti is a mesh that satisfies t0 = 0 = tl and we have Llti = 0, 1 ≤ i ≤ l − 1 along the same lines as before. As a result, the discrete minimum principle is applied to the mesh function ti gives ti = −(vi −wi ) ≥ 0 i.e., (vi −wi ) ≤ 0 f or 0 ≤ i ≤ l (4.15) From Eq. (4.14) and Eq. (4.15), we get vi − wi = 0 . This suggests that the discrete problem (2.1)-(2.2) solution is unique. The existence of linear equation is assumed by their uniqueness. Now Int. J. Anal. Appl. (2022), 20:63 9 we’ll explain how to prove the bound on < zi >li = 0. We do this by introducing two barrier functions Ψ± i Defined by Ψ± i = ||r||h,∞θ−1 + C1(||ϕ||h,∞ + |τ|) ±Zi, 0 ≤ i ≤ l where C1 ≥ 1 is a positive constant that can be chosen at random. Then we have Ψ± i = ||r||h,∞θ−1 + C1(||ϕ||h,∞ + |τ|) ±Zi, 0 ≤ i ≤ l Now we have Ψ±0 = ||r||h,∞θ −1 + C1(||ϕ||h,∞ + |τ|) ±Z0, Ψ±0 = (||r||h,∞θ −1 + C1(||ϕ||h,∞ + |τ|) ±C1ϕ0) ≥ 0 since ||ϕ||h,∞ + |τ| ≥ ϕ0 and C1 ≥ 1, Ψ± l = ||r||h,∞θ−1 + C1(||ϕ||h,∞ + |τ|) ±Zl Ψ± l = (||r||h,∞θ−1 + C1(||ϕ||h,∞ + |τ|) ±τ) ≥ 0 Case 1. For i = 1, 2, 3, . . . .,ν LlΨ± i = ερi (Ψ ± i )′′ + (2α + β)qi Ψ ± i = Gi −Di−ν+1Ψ±i−ν+1 −Ei−νΨ ± i−ν −Fi−ν−1Ψ ± i−ν−1 = −qi (||r||h,∞θ−1 + C2(||ϕ||h,∞ + |τ|) ±LlZi ) = ( −qi θ ||r||h,∞ ± ri ) −qiC2(||ϕ||h,∞ + |τ|) ≤ 0 Case 2. For i = ν + 1, . . . ., l1 − 1 LlΨ± i = ερi (Ψ ± i )′′ + 2α + βqi Ψ ± i + Di−ν+1Ψ ± i−ν+1 = Gi −Ei−νΨ ± i−ν −Fi−ν−1Ψ ± i−ν−1 = −qi (||r||h,∞θ−1 + C2(||ϕ||h,∞ + |τ|) ±LlZi ) = ( −qi θ ||r||h,∞ ± ri ) −qiC2(||ϕ||h,∞ + |τ|) ≤ 0 Case 3. For i = l1, . . . .., l − 1 LlΨ± i = ερi (Ψ ± i )′′ + (2α + β)qi Ψ ± i + Di−ν+1Ψ ± i−ν+1 + Ei−νΨ ± i−ν = Gi −Fi−ν−1Ψ ± i−ν−1 = −qi (||r||h,∞θ−1 + C2(||ϕ||h,∞ + |τ|) ±LlZi ) = ( −qi θ ||r||h,∞ ± ri ) −qiC2(||ϕ||h,∞ + |τ|) ≤ 0 Using the above cases, we have LlΨi ≤ 0, 1 ≤ i ≤ l. Using the discrete minimum principle Ψ± i = (||r||h,∞θ−1 + C1(||ϕ||h,∞ + |τ|) ±Zl) ≥ 0, 0 ≤ i ≤ l which proves the desired results. As a result of the previous theorems, the problem (2.1)-(2.2) solution is uniformly bounded, regardless of the mesh size h and the perturbation ε, demonstrating that the proposed scheme is stable for all h. � 10 Int. J. Anal. Appl. (2022), 20:63 5. Truncation error Using Taylor series, the truncation error in the method Eq. (4.2) is T (h) = {(2α + β) − 1}σεz′′i h 2+ {[( −2α 3 ) + ( 2ωσε + 1 6 ) β ] piz 3 i + ( 2α− 1 12 ) σεz4i } h4+O(h6) (5.1) The truncation error is fourth-order for any arbitrary α and β with 2α + β = 1 and any value of ω. 6. Convergence analysis The matrix form of the Eq. (4.5) is Az + B + T (h) = 0 (6.1) A =   B1 A1 0 · · · · · · · · · · · · · · · · · · 0 C2 B2 A2 0 · · · · · · · · · · · · · · · 0 0 C3 B3 A3 0 · · · · · · · · · · · · 0 ... ... ... ... ... ... ... ... ... ... Dν 0 · · · Cν Bν Aν 0 · · · · · · 0 Eν+1 Dν+1 0 · · · Cν+1 Bν+1 Aν+1 0 0 Fν+2 Eν+2 Dν+2 0 · · · Cν+2 Bν+2 Aν+2 0 0 0 Fν+3 Eν+3 Dν+3 0 Cν+3 Bν+3 Aν+3 0 0 ... ... ... ... ... ... ... ... ... 0 · · · 0 Fl−2 El−2 Dl−2 0 · · · Cl−2 Bl−2 Al−2 · · · · · · 0 Fl−1 El−1 Dl−1 0 · · · Cl−1 Bl−1   and B = [κ1,κ2, ..κν,κν+1, κν+2, ...,κl−2,κl−1] where κi =   Gi −Di−ν+1zi−ν+1 −Ei−νzi−ν −Fi−ν−1zi−ν−1 −C1z0, ∀ 1 ≤ i ≤ ν − 1 Gi −Ei−νzi−ν −Fi−ν−1zi−ν−1, ∀ i = ν Gi −Fi−ν−1zi−ν−1, ∀ i = ν + 1 Gi −Al−1zl, ∀ ν + 2 ≤ i ≤ l − 1 where Ai+1 = σ� + αh 2qi+1,Bi = −2σ� + βh2qi,Ci−1 = σ� + αh2qi−1 Di−ν+1 = −αhpi−1 2 + βhpi 2 (1 + 2ωh2qi+1 + ωh(3pi+1 + pi−1)) + 3αh 2 pi+1 Ei−ν = 2αhpi−1 − 2βωpih2(pi+1 + pi−1) − 2αhpi+1 Fi−ν−1 = −3αhpi−1 2 − βpih 2 (1 + 2ωh2qi−1 −ωh(pi+1 + 3pi−1)) + αh 2 pi+1 Gi = h 2(αri+1 + βri + αri−1) −βωpih3(ri+1 − ri−1) Int. J. Anal. Appl. (2022), 20:63 11 for i = 0, 1, . . . , l − 1 and T (h) = O ( h4 ) ,z = [z1,z2, . . . ,zl−1] T ,T (h) = [T1,T2, . . . ,Tl−1] T ,O = [0, 0, . . . , 0] T are related vectors of Eq. (6.1). Let z̃ = [z̃1, z̃2, . . . , z̃l−1] T ∼= Z satisfies the equation Az̃ + B = 0 (6.2) Let ei = zi −Zi, i = 1 (1) l − 1 be the error so that E = [e1,e2, ...,el−1] T = u −U. Subtracting Eq. (6.1) from Eq. (6.2), we obtain the error equation AE = T (h) (6.3) Let the sum of the ith row elements of the matrix A be si. Then si = h 2 (αqi−1 + βqi + αqi+1) for 1 ≤ i ≤ ν − 1 si = h 2 ( −αpi−1 + βpi ( 1 + 2ωh2qi+1 + ωh (3pi+1 + pi−1) ) + 3αpi+1 ) + h2 (αqi−1 + βqi + αqi+1) for i = ν si = h 2 ( −3αpi−1 −βpi ( 1 + 2ωh2qi−1 −ωh (pi+1 + 3pi−1) ) + αpi+1 ) + h2 (αqi−1 + βqi + αqi+1) for i = ν + 1 si = h 2 (αqi−1 + βqi + αqi+1) for ν + 2 ≤ i ≤ l − 1. Let ζ1∗ = min1≤i≤l |p (ui )| , ζ∗1 = max1≤i≤l |p (ui )| ,ζ2∗ = min1≤i≤l |q (ui )| and ζ∗2 = max1≤i≤l |q (ui )|. It is verified that A is irreducible and monotone since 0 < ε � 1 and ε ∝ O (h) with sufficiently small h. Hence A−1 exists and A−1 ≥ 0. Therefore, using Eq. (6.3) we have ||E|| ≤ ||A−1|| ||T || (6.4) For small h, we have si ≥ h2 [2 (α + β/2) ζ2∗] for 1 ≤ i ≤ ν − 1 si ≥ h2 [2 (α + β/2) ζ2∗] for i = ν si ≥ h2 [2 (α + β/2) ζ2∗] for i = ν + 1 si ≥ h2 [2 (α + β/2) ζ2∗] for ν + 2 ≤ i ≤ l − 1 Let A−1 i,k be the (i,k)th element of A−1 and we define || A|| = max1≤i≤l−1 l−1∑ k=1 A−1 i,k and ||T (h)|| = max1≤i≤l−1|Ti (h)|. 12 Int. J. Anal. Appl. (2022), 20:63 Since A−1 i,k ≥ 0 and ∑l−1 k=1 A −1 i,k . sk = 1 for i = 12, 3, ......, l − 1. Hence ν−1∑ k=1 A−1 i,k ≤ 1 min1≤k≤ν−1s k < 1 h2 [2 (α + β/2) ζ2∗] , i = 1, 2, 3, . . . ,ν − 1 (6.5) A−1 i,k ≤ 1 sν < 1 h2 [2 (α + β/2) ζ2∗] , i = ν,ν + 1 (6.6) Furthermore l−1∑ k=ν+2 A−1 i,k ≤ 1 min1≤k≤ν−1sk < 1 h2 [2 (α + β/2) ζ2∗] , i = ν + 2, ν + 3, ..., l − 1. (6.7) Using Eqs. (6.5) – (6.6) and Eq. (6.7), we get ||E|| ≤ O(h2). (6.8) As a result, the proposed exponential spline strategy converges to second order. The convergence analysis for the right end boundary layer can be examined using the same procedure. 7. Numerical Examples The effectiveness of the numerical approach can be demonstrated by the following instances. MAE (maximum absolute error) in the solution is computed using the double mesh principle El = max0≤i≤l ∣∣z li −z2l2i ∣∣. The order of convergence is also investigated by r l = log2 (Ei l/Ei 2l). We exhibit graphs showing the computed solution of the issue for various δ values. Example 1: εz′′ + z′(u −δ) −z(u) = 0 with z(u) = 1, −δ ≤ u ≤ 0 and z(1) = 0. Example 2: εz′′ + exp(−0.25u) −z(u) = 0 with z(u) = 1 for −δ ≤ u ≤ 0,z (1) = 1. Example 3: εz′′ −z′(u −δ) −z(u) = 0 with z(u) = 1, −δ ≤ u ≤ 0 and z(1) = −1. Example 4: εz′′ −z′(u −δ) + z(u) = 0 with z(u) = 1 for −δ ≤ u ≤ 0,z (1) = 1. 8. Conclusion A computational approach for solving the singularly perturbed differential equation with the large delay is derived using a special type of mesh. A numerical scheme consisting of a fitting parameter is developed to minimize the error and to control the layer structure in the solution. Four examples are solved and computational results with large delay are shown in Table 1-4. In the proposed method, we also analyzed the effect of the large delay on the layer structure or oscillatory behaviour of the solutions with and without the fitting parameter in the Figures 1-8. The graphs depict the layer behaviour in the solution of the examples with and without the fitting parameter. The impact of the fitting factor Int. J. Anal. Appl. (2022), 20:63 13 in controlling the oscillations in the solution is also depicted in the graphs. We have clearly noticed that the fitting parameter controls the oscillations in the layer for the large delay values. The proposed method is simple and it works very well with small as well as large delay. 14 Int. J. Anal. Appl. (2022), 20:63 Int. J. Anal. Appl. (2022), 20:63 15 Figure 1. Layer profile of Example 1 for 𝜀 = 0.01 without fitting factor. Figure 2. Layer profile of Example 1 for 𝜀 = 0.01 with fitting factor. Figure 3. Layer profile of Example 2 for 𝜀 = 0.01 without fitting factor. 16 Int. J. Anal. Appl. (2022), 20:63 Figure 4. Layer profile of Example 2 for 𝜀 = 0.01 with fitting factor. Figure 5. Layer profile of Example 3 for 𝜀 = 0.01 without fitting factor. Figure 6. Layer profile of Example 3 for 𝜀 = 0.01 with fitting factor. Int. J. Anal. Appl. (2022), 20:63 17 Figure 7. Layer profile of Example 4 for 𝜀 = 0.01 without fitting factor. Figure 8. Layer profile of Example 4 for 𝜀 = 0.01 with fitting factor. 18 Int. J. Anal. Appl. (2022), 20:63 Conflicts of Interest: The authors declare that there are no conflicts of interest regarding the publi- cation of this paper. References [1] R. Bellman, K.L. Cooke, Differential-Difference Equations, Academic Press, New York, USA, 1963. [2] E.P. Doolan, J.J.H. Miller, W.H.A. Schilders, et al. Uniform Numerical Methods for Problems with Initial and Boundary Layers, Boole Press, Dublin, 1980. [3] R.D. Driver, Ordinary and Delay Differential Equations, Springer, Belin-Heidelberg, New York, 1977. [4] G.M. Amiraliyev, E. Cimen, Numerical Method for a Singularly Perturbed Convection–diffusion Problem With Delay, Appl. Math. Comput. 216 (2010), 2351–2359. https://doi.org/10.1016/j.amc.2010.03.080. [5] V.Y. Glizer, Asymptotic Solution of a Boundary-Value Problem for Linear Singularly-Perturbed Functional Dif- ferential Equations Arising in Optimal Control Theory, J. Optim. Theory Appl. 106 (2000), 309–335. https: //doi.org/10.1023/a:1004651430364. [6] M.K. Kadalbajoo, K.K. Sharma, A Numerical Method Based on Finite Difference for Boundary Value Problems for Singularly Perturbed Delay Differential Equations, Appl. Math. Comput. 197 (2008), 692–707. https://doi.org/ 10.1016/j.amc.2007.08.089. [7] Y. Kuang, Delay Differential Equations with Applications in Population Dynamics, Academic Press, New York, 1993. [8] D. Kumara Swamy, K. Phaneendra, A. Benerji Babu, Y.N. Reddy, Computational Method for Singularly Perturbed Delay Differential Equations With Twin Layers or Oscillatory Behaviour, Ain Shams Eng. J. 6 (2015), 391–398. https://doi.org/10.1016/j.asej.2014.10.004. [9] M. Lalu, K. Phaneendra, S.P. Emineni, Numerical Approach for Differential-Difference Equations Having Layer Behaviour With Small or Large Delay Using Non-Polynomial Spline, Soft Comput. 25 (2021), 13709-13722. https: //doi.org/10.1007/s00500-021-06032-5. [10] L. Sirisha, K. Phaneendra, Y.N. Reddy, Mixed Finite Difference Method for Singularly Perturbed Differential Difference Equations With Mixed Shifts via Domain Decomposition, Ain Shams Eng. J. 9 (2018), 647–654. https://doi.org/10.1016/j.asej.2016.03.009. [11] C.G. Lange, R.M. Miura, Singular Perturbation Analysis of Boundary Value Problems for Differential-Difference Equations. V. Small Shifts with Layer Behavior, SIAM J. Appl. Math. 54 (1994), 249–272. https://doi.org/10. 1137/s0036139992228120. [12] C.G. Lange, R.M. Miura, Singular Perturbation Analysis of Boundary-Value Problems for Differential-Difference Equations. VI. Small Shifts with Rapid Oscillations, SIAM J. Appl. Math. 54 (1994), 273–283. https://doi.org/ 10.1137/s0036139992228119. [13] A. Lasota, M. Walewska, Mathematical Models of the Red Blood Cell System, Mat. Stos. 6 (1976), 25–40. [14] M.C. Mackey, L. Glass, Oscillation and Chaos in Physiological Control Systems, Science. 197 (1977), 287–289. https://doi.org/10.1126/science.267326. [15] A. Martin, S. Ruan, Predator-Prey Models With Delay and Prey Harvesting, J. Math. Biol. 43 (2001), 247–267. https://doi.org/10.1007/s002850100095. [16] K. Phaneendra, M. Lalu, Numerical Solution of Singularly Perturbed Delay Differential Equations Using Gaussion Quadrature Method, J. Phys.: Conf. Ser. 1344 (2019), 012013. https://doi.org/10.1088/1742-6596/1344/1/ 012013. [17] P. Rai, K.K. Sharma, Numerical Analysis of Singularly Perturbed Delay Differential Turning Point Problem, Appl. Math. Comput. 218 (2011), 3483–3498. https://doi.org/10.1016/j.amc.2011.08.095. https://doi.org/10.1016/j.amc.2010.03.080 https://doi.org/10.1023/a:1004651430364 https://doi.org/10.1023/a:1004651430364 https://doi.org/10.1016/j.amc.2007.08.089 https://doi.org/10.1016/j.amc.2007.08.089 https://doi.org/10.1016/j.asej.2014.10.004 https://doi.org/10.1007/s00500-021-06032-5 https://doi.org/10.1007/s00500-021-06032-5 https://doi.org/10.1016/j.asej.2016.03.009 https://doi.org/10.1137/s0036139992228120 https://doi.org/10.1137/s0036139992228120 https://doi.org/10.1137/s0036139992228119 https://doi.org/10.1137/s0036139992228119 https://doi.org/10.1126/science.267326 https://doi.org/10.1007/s002850100095 https://doi.org/10.1088/1742-6596/1344/1/012013 https://doi.org/10.1088/1742-6596/1344/1/012013 https://doi.org/10.1016/j.amc.2011.08.095 Int. J. Anal. Appl. (2022), 20:63 19 [18] P. Rai, K.K. Sharma, Numerical Approximation for a Class of Singularly Perturbed Delay Differential Equa- tions With Boundary and Interior Layer(s), Numer Algor. 85 (2019), 305–328. https://doi.org/10.1007/ s11075-019-00815-6. [19] A.S.V. Ravikanth, P. Murali, Numerical Treatment for a Singularly Perturbed Convection Delayed Dominated Dif- fusion Equation via Tension Splines, Int. J. Pure Appl. Math. 113 (2017), 1314-3395. [20] R.K. Adivi Sri Venkata, M.M.K. Palli, A numerical approach for solving singularly perturbed convection de- lay problems via exponentially fitted spline method, Calcolo. 54 (2017), 943-961. https://doi.org/10.1007/ s10092-017-0215-6. [21] Y.N. Reddy, GBSL Soujanya, K. Phaneendra, Numerical Integration Method for Singularly Perturbed Delay Differ- ential Equations, Int. J. Appl. Sci. Eng. 10 (2012), 249-261. [22] R.B. Stein, Some Models of Neuronal Variability, Biophys. J. 7 (1967), 37-68. https://doi.org/10.1016/ s0006-3495(67)86574-3. [23] W.K. Zahra, A.M. El Mhlawy, Numerical Solution of Two-Parameter Singularly Perturbed Boundary Value Problems via Exponential Spline, J. King Saud Univ. - Sci. 25 (2013), 201-208. https://doi.org/10.1016/j.jksus.2013. 01.003. https://doi.org/10.1007/s11075-019-00815-6 https://doi.org/10.1007/s11075-019-00815-6 https://doi.org/10.1007/s10092-017-0215-6 https://doi.org/10.1007/s10092-017-0215-6 https://doi.org/10.1016/s0006-3495(67)86574-3 https://doi.org/10.1016/s0006-3495(67)86574-3 https://doi.org/10.1016/j.jksus.2013.01.003 https://doi.org/10.1016/j.jksus.2013.01.003 1. Introduction 2. Statement of the problem 3. An exponential spline 4. Numerical Approach 4.1. Left-end layer 4.2. Right-end layer 5. Truncation error 6. Convergence analysis 7. Numerical Examples 8. Conclusion References