# Fitted Parameter Exponential Spline Method for Singularly Perturbed Delay Differential Equations with a Large Delay

Authors: E. Siva Prasad; R. Omkar; Kolloju Phaneendra
Journal: Computational and Mathematical Methods
 (2022)

---

## Abstract

In this paper, we present a new computational method based on an exponential spline for solving a class of delay differential equations with a negative shift in the differentiated term. When the shift parameter isOε, the proposed method works well and also controls the oscillations in the solution’s layer region. To accomplish this, we included a parameter in the proposed numerical scheme that is based on a special type of mesh, and the parameter is evaluated using the theory of singular perturbation. Maximum absolute errors and convergences of numerical solutions are tabulated to demonstrate the efficiency of the proposed computational method and to support the convergence analysis of the presented scheme.


---

## Body

## 1. Introduction

In this article, we examine a section of differential equations with a delay in the first-order derivative term and a small parameter multiplied by the second-order derivative. These equations are widely used in the simulation of a wide range of physical applications, such as construction of logical circuits by means of a network built up by 1D excitable media [1], hybrid optical systems with delayed feedback [2], population dynamics [3], nonlinear delay differential equations relating to physiological control systems [4], red blood cell system [5], predator-prey models [6], a Hamiltonian boundary value problem associated with optimal control theory [7], and neuronal variability problems connected to patterns of nerve action potentials formed by unit quantal inputs occurring at random which are studied [8].Bellman et al. [9], Doolan et al. [10], Driver [11], Norkin [12], Kokotovic et al. [13], Miller et al. [14], Smith [15], and O’Malley [16] are just a few of the books in the collection that can be used for further investigations of precise characteristics of the aforementioned class of models and perturbation problems. Lange and Miura [17, 18] devised an asymptotic method for solving equations with layer behaviour. Researchers presented a mathematical model in [17] for predicting the time it takes nerve cells with arbitrary synaptic inputs in dendrites to produce potential action. The researchers in [18] demonstrated the issues by using solutions that exhibit rapid oscillations. Kadalbajoo and his colleagues began widespread numerical work using finite differentiation techniques as well as the B-spline technique with fitted mesh [19–22]. The researchers in [23] suggested a numerical approach via nonpolynomial spline for the solution of differential-difference equations with small and large delay.In [24], researchers proposed an unconditionally stable collocation method for solving singularly perturbed one-dimensional time dependent problems. The method employs the implicit finite difference method on a uniform mesh in the temporal direction and the B-spline collocation method on a nonuniform Shishkin mesh in the spatial direction. Kadalbajoo and Gupta [25] discussed a B-spline collocation approach on a piecewise uniform Shishkin mesh for the numerical solution of a singularly perturbed problem with a boundary layer on the right side of the domain. Furthermore, it is demonstrated that the approach is uniformly convergent in the maximum norm of the second order. A brief survey of various computational techniques on different type of singularly perturbed problems is given in [26]. Gupta and Kadalbajoo [27] proposed a parameter uniform method for solving a singularly perturbed one-dimensional parabolic equation with multiple boundary turning points on a rectangular domain that combines the Euler method for time variables on a uniform mesh and the B-spline collocation method for space variables on a Shishkin mesh. The authors of [28] developed a numerical scheme for time-dependent singularly perturbed problems based on Mickens nonstandard finite difference method. The method is shown to be unconditionally stable and parameter uniformly convergent for both large and small shift arguments. Das [29] presented an a posteriori-based convergence analysis on an adaptive mesh for a nonlinear singularly perturbed system of delay differential equations. The layer-adaptive solution on the a posteriori-generated mesh has been shown to converge uniformly to the exact solution with optimal order accuracy.The authors of [30] considered a moving mesh-adaptive algorithm that adapts meshes to boundary layers for the numerical solution of singularly perturbed convection-diffusion-reaction problems with two small parameters. A special positive monitor function is used to generate the mesh and demonstrate that the proposed method is parameter independent and convergent for a class of model problems. Das and Natesan [31] used cubic spline approximation on a piecewise uniform Shishkin mesh to solve a system of reaction diffusion differential equations for Robin type boundary-value problems. This method is a second-order uniform convergence method that uses a central difference scheme for the outer region of the boundary layers and a cubic spline approximation for the inner region. The authors of [32] presented two adaptive methods for a system of reaction diffusion problems using the r-refinement strategy, which moves the mesh points toward the boundary layers. On the adaptively generated equidistributed mesh generated by a curvature-based monitor function, first-order uniform solutions are obtained using forward and backward schemes. To improve the discrete solution to second-order uniform accuracy, a cubic spline-based scheme is used. The authors of [33] investigated a numerical scheme for solving singularly perturbed parabolic partial differential equations using the Euler method and the upwind scheme on an adaptive mesh. They also examined the convergence analysis of the proposed scheme. The authors of [34] investigated the approximation of solutions to a class of fractional order Volterra-Fredholm integro-differential equations. The detailed analysis of the existence and uniqueness of the solution as well as the solution’s bounds was addressed. A perturbation approach based on homotopy analysis is used to solve the same problem. The authors of [35] considered the initial and boundary value problems of a class of fractional order Volterra integro-differential equations of the first kind. Using Leibniz’s rule, the same problem is reduced to a Volterra integro-equation of the second kind, and an operative-based method is used to approximate their solutions. Das et al. [36] investigate the homotopy perturbation method for approximating solutions of Volterra-Fredholm integro-differential equations using semianalytical approaches.The following is the paper’s outline: the problem is described in Section2 along with the conditions for the layer’s behaviour, and the exponential spline function is defined in Section 3. Section 4 discusses the numerical scheme with a large delay. Section 5 consists of truncation error and convergence analysis. Numerical experiments are addressed in Section 6. Finally, in Section 7, the article’s discussions and conclusions are addressed.


## 2. Statement of the Problem

Consider a singularly perturbed delay differential equation with a delay term, that is, a negative shift in the first-order derivative(1)Lε,δνs≡εν″s+psν′s−δ+qsνs=rs,0≤s≤1,0<ε<<1,in accordance with the conditions.
(2)νs=ϕs,−δ≤s≤0andν1=φ,where ps,qs,rs,andϕs are expected to be smooth, bounded functions in [0, 1], δ is the delay parameter, and φ is a finite constant. When δ=0, Equation (1) is reduced to a singularly perturbed equation with boundary layer behaviour and turning points that depend on the sign of ps. Throughout the interval [0, 1], the solution νs has a boundary layer on the left end side when ps is positive and on the right end side if ps is negative. When the delay parameter δε is of Oε, the behaviour of the layer can change and even be ruined, or the solution exhibits oscillatory behaviour. To monitor these oscillations in the solution profile, we devise a numerical method that combined a special mesh introduced in [21] with the exponential spline method’s fitting parameter.


### 2.1. Properties of Continuous Problem

LetLε,δ be the differential operator for the problem Equation (1) which is defined for any smooth function Ωs∈C2 as Lε,δΩs=εΩ′′s+psΩ′s−δ+qsΩs.Lemma 1 (Continuous maximum principle).
LetΩs be the smooth function satisfying Ω0≥0 and ΩN≥0. Then, Lε,δΩs≥0,∀i=1,2,⋯,N−1 implies that Ωs≥0,∀i=0,1,2,⋯,N.Proof.
Lets∗ be such that Ωs∗=min0≤s≤1Ωs and assume that Ωs∗<0. Clearly, s∗∉0,1; therefore, Ω′s∗=0 andΩ′′s∗≥0. Now, in order to prove Lε,δΩs<0, we consider different cases.
Case 1:0<s∗≤δ,Lε,δΩs∗=εΩ′′s∗+psϕ′s∗−δ+qsΩs∗<0.
Case 2:δ<s∗≤1,Lε,δΩs∗=εΩ′′s∗+psΩ′s∗−δ+qsΩs∗<0
Combining the above two cases, we getLε,δΩs<0 which contradicts the hypothesis that Lε,δΩs≥0,∀s∈0,1. Hence, our assumption that Ωs∗<0 is wrong, and thus, Ωs≥0,∀s∈01.Lemma 2.
Under the assumption thatps≥α∗>0, qs≥τ>0 where α∗, τ are positive constants, the solution of Equation (1) with boundary conditions Equation (2) exists and satisfies
(3)ν≤τ−1r+C1maxϕ0,φ.Proof.
Let us construct the two-barrier functionπ± defined by
(4)π±s=τ−1r+C1maxϕ0,φ±νs.
Then, we have(5)π±0=τ−1r+C1maxϕ0,φ±ν0=τ−1r+C1maxϕ0,φ±ϕ0≥0,π±1=τ−1r+C1maxϕ0,φ±ν1=τ−1r+C1maxϕ0,φ±φ≥0.Case 1.
For0<s≤δ,
(6)Lε,δπ±s=επ±s′′+psπ±s−δ′+qsπ±s=qsτ−1r+C1maxϕ0,φ±Lε,δνs=qsτ−1r+C1maxϕ0,φ±rs−psϕs,since qs≥τ>0, i.e., qsτ−1≥1.
(7)Lε,δπ±s≥r±rs+C1qsmaxϕ0,φ−psϕs.
Sincer>rs, and we choose the constant C1 so that the sum of the first two terms dominates the third term in the above inequality (7), we then obtain
(8)Lε,δπ±s≥0.Case 2.
Forδ<s≤1,
(9)Lε,δπ±s=επ±s′′+psπ±s−δ′+qsπ±s=qsτ−1r+C1maxϕ0,φ±Lε,δνs=qsτ−1r+C1maxϕ0,φ±rs,since qs≥τ>0, i.e., qsτ−1≥1.
(10)Lε,δπ±s≥r±rs+C1qsmaxϕ0,φ≥0.
From the inequality (8) and (10), we get Lε,δπ±s≥0,0≤s≤1.
Therefore, using Lemma1, we obtain π±s≥0, for all s∈0,1, which gives required estimate.
Lemma1 implies that the solution is unique, and since the problem under consideration is linear, the existence of the solution is implied by its uniqueness. Additionally, Lemma 2 gives the boundedness of the solution [22].Lemma 3.
The bounds on the derivatives of the solutionνs of Equation (1) with respect to x is provided by
(11)diνdsi≤C1+ε−ie−α∗1−s/ε,∀s∈01,0≤i≤4.Proof.
One can refer to [28].


## 2.1. Properties of Continuous Problem

LetLε,δ be the differential operator for the problem Equation (1) which is defined for any smooth function Ωs∈C2 as Lε,δΩs=εΩ′′s+psΩ′s−δ+qsΩs.Lemma 1 (Continuous maximum principle).
LetΩs be the smooth function satisfying Ω0≥0 and ΩN≥0. Then, Lε,δΩs≥0,∀i=1,2,⋯,N−1 implies that Ωs≥0,∀i=0,1,2,⋯,N.Proof.
Lets∗ be such that Ωs∗=min0≤s≤1Ωs and assume that Ωs∗<0. Clearly, s∗∉0,1; therefore, Ω′s∗=0 andΩ′′s∗≥0. Now, in order to prove Lε,δΩs<0, we consider different cases.
Case 1:0<s∗≤δ,Lε,δΩs∗=εΩ′′s∗+psϕ′s∗−δ+qsΩs∗<0.
Case 2:δ<s∗≤1,Lε,δΩs∗=εΩ′′s∗+psΩ′s∗−δ+qsΩs∗<0
Combining the above two cases, we getLε,δΩs<0 which contradicts the hypothesis that Lε,δΩs≥0,∀s∈0,1. Hence, our assumption that Ωs∗<0 is wrong, and thus, Ωs≥0,∀s∈01.Lemma 2.
Under the assumption thatps≥α∗>0, qs≥τ>0 where α∗, τ are positive constants, the solution of Equation (1) with boundary conditions Equation (2) exists and satisfies
(3)ν≤τ−1r+C1maxϕ0,φ.Proof.
Let us construct the two-barrier functionπ± defined by
(4)π±s=τ−1r+C1maxϕ0,φ±νs.
Then, we have(5)π±0=τ−1r+C1maxϕ0,φ±ν0=τ−1r+C1maxϕ0,φ±ϕ0≥0,π±1=τ−1r+C1maxϕ0,φ±ν1=τ−1r+C1maxϕ0,φ±φ≥0.Case 1.
For0<s≤δ,
(6)Lε,δπ±s=επ±s′′+psπ±s−δ′+qsπ±s=qsτ−1r+C1maxϕ0,φ±Lε,δνs=qsτ−1r+C1maxϕ0,φ±rs−psϕs,since qs≥τ>0, i.e., qsτ−1≥1.
(7)Lε,δπ±s≥r±rs+C1qsmaxϕ0,φ−psϕs.
Sincer>rs, and we choose the constant C1 so that the sum of the first two terms dominates the third term in the above inequality (7), we then obtain
(8)Lε,δπ±s≥0.Case 2.
Forδ<s≤1,
(9)Lε,δπ±s=επ±s′′+psπ±s−δ′+qsπ±s=qsτ−1r+C1maxϕ0,φ±Lε,δνs=qsτ−1r+C1maxϕ0,φ±rs,since qs≥τ>0, i.e., qsτ−1≥1.
(10)Lε,δπ±s≥r±rs+C1qsmaxϕ0,φ≥0.
From the inequality (8) and (10), we get Lε,δπ±s≥0,0≤s≤1.
Therefore, using Lemma1, we obtain π±s≥0, for all s∈0,1, which gives required estimate.
Lemma1 implies that the solution is unique, and since the problem under consideration is linear, the existence of the solution is implied by its uniqueness. Additionally, Lemma 2 gives the boundedness of the solution [22].Lemma 3.
The bounds on the derivatives of the solutionνs of Equation (1) with respect to x is provided by
(11)diνdsi≤C1+ε−ie−α∗1−s/ε,∀s∈01,0≤i≤4.Proof.
One can refer to [28].


## 3. Exponential Spline

We divide [0, 1] intoN equal subdomains of grid size h=1/N, so that si=ihi=0,1,2,⋯,N with 0=s0,1=sN. Let νs be the precise solution and νi be an estimate to νsi by the exponential spline Ψis passing through the points si,νi and si+1,νi+1. The exponential spline function Ψis has the following form for each ith segment, satisfying the condition of first derivative continuity at the joint nodes si,νi.
(12)Ψis=Ai+Bis−si+Kieτs−si+Lie−τs−si,i=0,1,2,⋯,N−1,where Ai,Bi,Ki, and Li are constants and τ is a parameter. Here, Ψis is of the class C20,1 interpolates νis at mesh points si,i=0,1,2,⋯,N depends on a parameter τ and reduces to the cubic spline Ψis in [0,1] as τ⟶0.DefineΨisi=νi,Ψisi+1=νi+1,Ψi′′si=Mi,Ψi″si+1=Mi+1, to calculate the coefficient values in Equation (12) in terms of νi,νi+1,Mi, and Mi+1. Following the exercise in algebra, we get
(13)Ai=νi−Miτ2,Bi=νi+1−νih−Mi+1−MiτΘ,Ki=Mi+1−e−ΘMi2τ2sinhΘ,Li=eΘMi−Mi+12τ2sinhΘ,whereΘ=τh, for i=0,1,2,⋯,N−1.The continuity condition of the first derivative, i.e.,Ψi−1′si=Ψi′si at si,νi, allows us to obtain the following relations for i=1,2,⋯,N−1.(14)αMi−1+2βMi+αMi+1=νi−1−2νi+νi+1h2,where =1/Θ2−1/ΘsinhΘ,β=cothΘ/Θ−1/Θ2,Mj=ν′′sj,j=i,i±1andΘ=τh.The local truncation errorTi of the scheme in Equation (14) is Ti=h21−2α−2βνi2+h41−12α/12νi4+h61−30α/360νi6+Oh8 for i=1,2,⋯,N−1.As a result, we obtain different orders forTi for different values of α and β in Equation (14): (i) fourth order for any arbitrary α and β with α+β=1/2; (ii) six orders for α=1/12,β=5/12.


## 4. Numerical Approach with Fitting Parameter for Large Delay

The solution’s layer behaviour is preserved, and precise results are obtained when compared to the perturbation parameter; the shift parameter is smaller [22, 23]. However, if δε is of order Oε, the solution’s layer behaviour is no longer well-preserved, and oscillations becomes visible. As a result, we are developing a numerical scheme with fitting parameters that is based on an exponential spline method and a special type of mesh developed in [21]. In addition, we examine the solution layer behaviour graphically for large delays in order to demonstrate the significance of the fitting parameter in our proposed scheme.


### 4.1. Numerical Scheme for Left-End Boundary Layer

Assume thatps≥M¯>0 and ε≥0 in [0, 1], where M¯ is a positive constant. As a result of this theory, Equations (1) and (2) show layer behaviour at s=0 for small values of ε. Now, we consider a mesh developed in [21], i.e., chooses the mesh as h=δ/m, where m=ωℓ, ℓis the mantissa of δ, and ω is a positive integer. The discretization of the boundary value problem Equations (1) and (2) with a special mesh result in
(15)ενi″=rsi−psiνi−m′−qsiνi,with(16)νi=ϕifori=−m,−m+1,⋯,0andνN=φ.Using the difference scheme in Equation (14) and the difference estimates of the first-order derivatives listed below,
(17)νi−1′=−νi+1+4νi−3νi−12h+Oh2,νi+1′=νi+1−νi−12h+Oh2,νi+1′=3νi+1−4νi+νi−12h+Oh2,we get
(18)εh2νi+1−2νi+νi−1=−αpi+12h+βpih+3αpi−12hνi−m−1+2αpi+1h−2αpi−1hνi−m+−3αpi+12h−βpih+αpi−12hνi−m+1−αqi−1νi−1−2βqiνi−αqi+1νi+1+αri+1+2βri+αri−1.Now, to enhance the scheme’s effectiveness, we incorporate a fitting parameter, which is determined using the theory of singular perturbations.(19)σρεh2νi+1−2νi+νi−1=−αpi+12h+βpih+3αpi−12hνi−m−1+2αpi+1h−2αpi−1hνi−m+−3αpi+12h−βpih+αpi−12hνi−m+1−αqi−1νi−1−2βqiνi−αqi+1νi+1+αri+1+2βri+αri−1fori=1,2,⋯,N−1.By multiplying Equation (19) by h and assuming limit as h⟶0 [4], we get
(20)limh⟶0σρνi+1h−2νih+νi−1h=α+βp0limh⟶0νi−m−1h−νi−m+1h.Based on singular perturbation theory pertaining to the layer at the interval’s left end [16], we have
(21)limνh⟶0ih≈ν0ih+ϕ0−ν00exp−p0iρ+Oε,whereρ=hϵ.Now, by plugging Equation (21) into Equation (20) and exercising, we have
(22)σρ=α+βp0e−p0i−1−mρ−e−p0i−m+1ρe−p0i+1ρ−2e−p0iρ+e−p0i−1ρ,(23)σρ=α+βp0ep0mρep0ρ−e−p0ρep0ρ/2−e−p0ρ/22=α+βp0ep0mρep0ρ/2+e−p0ρ/2ep0ρ/2−e−p0ρ/2,(24)σ=ρα+βp0ep0mρcothp0ρ2,which is the parameter that has been fitted for the layer on the left-end of the domain.


### 4.2. Numerical Scheme for Right-End Boundary Layer

Assume thatps≤M¯<0 and ε≥0 in [0, 1], where M¯ is a positive constant. As a result of this theory, Equations (1) and (2) shows layer behaviour at s=1 for small values of ε. Based on singular perturbation theory for the layer at the right end of the interval [16], we have
(25)limh⟶0νih=ν00+ϕ0−ν01e−P11/ε−iρ+Oε,whereρ=hϵ.In the exponential spline difference scheme Equation (18), we insert a fitting parameter σρ and obtain the fitting parameter
(26)σ=ρα+βp0e−p0mρcothp1ρ2,employing the similar procedure as in the case of the left boundary layer. Now, by simplifying the difference scheme in Equation (19), we get
(27)Piνi−1+Qiνi+Riνi+1+Uiνi−m+1+Viνi−m+Wiνi−m−1=Gi,i=1,2,⋯,N−1,where Pi=σε+αh2qi−1,Qi=−2σε+2βh2qi,Ri=σε+αh2qi+1,Ui=−αhpi−1/2+2βhpi/2+3αhpi+1/2,Vi=2αhpi−1−2αhpi+1,Wi=−3αhpi−1/2−2βhpi/2+αhpi+1/2,Gi=h2αri−1+2βri+αri+1.Following Equation (16), the scheme of Equation (27) can be expressed as
(28)Piνi−1+Qiνi+Riνi+1=Gi−Uiνi−m+1−Viνi−m−Wiνi−m−1,for1≤i≤m−1,Piνi−1+Qiνi+Riνi+1+Uiνi−m+1=Gi−Viνi−m−Wiνi−m−1,fori=m,Piνi−1+Qiνi+Riνi+1+Uiνi−m+1+Viνi−m=Gi−Wiνi−m−1,fori=m+1,Piνi−1+Qiνi+Riνi+1+Uiνi−m+1+Viνi−m+Wiνi−m−1=Gi,form+2≤i≤N−1.The Gauss elimination scheme with partial pivoting is used to solve the aforementioned system of equations.


## 4.1. Numerical Scheme for Left-End Boundary Layer

Assume thatps≥M¯>0 and ε≥0 in [0, 1], where M¯ is a positive constant. As a result of this theory, Equations (1) and (2) show layer behaviour at s=0 for small values of ε. Now, we consider a mesh developed in [21], i.e., chooses the mesh as h=δ/m, where m=ωℓ, ℓis the mantissa of δ, and ω is a positive integer. The discretization of the boundary value problem Equations (1) and (2) with a special mesh result in
(15)ενi″=rsi−psiνi−m′−qsiνi,with(16)νi=ϕifori=−m,−m+1,⋯,0andνN=φ.Using the difference scheme in Equation (14) and the difference estimates of the first-order derivatives listed below,
(17)νi−1′=−νi+1+4νi−3νi−12h+Oh2,νi+1′=νi+1−νi−12h+Oh2,νi+1′=3νi+1−4νi+νi−12h+Oh2,we get
(18)εh2νi+1−2νi+νi−1=−αpi+12h+βpih+3αpi−12hνi−m−1+2αpi+1h−2αpi−1hνi−m+−3αpi+12h−βpih+αpi−12hνi−m+1−αqi−1νi−1−2βqiνi−αqi+1νi+1+αri+1+2βri+αri−1.Now, to enhance the scheme’s effectiveness, we incorporate a fitting parameter, which is determined using the theory of singular perturbations.(19)σρεh2νi+1−2νi+νi−1=−αpi+12h+βpih+3αpi−12hνi−m−1+2αpi+1h−2αpi−1hνi−m+−3αpi+12h−βpih+αpi−12hνi−m+1−αqi−1νi−1−2βqiνi−αqi+1νi+1+αri+1+2βri+αri−1fori=1,2,⋯,N−1.By multiplying Equation (19) by h and assuming limit as h⟶0 [4], we get
(20)limh⟶0σρνi+1h−2νih+νi−1h=α+βp0limh⟶0νi−m−1h−νi−m+1h.Based on singular perturbation theory pertaining to the layer at the interval’s left end [16], we have
(21)limνh⟶0ih≈ν0ih+ϕ0−ν00exp−p0iρ+Oε,whereρ=hϵ.Now, by plugging Equation (21) into Equation (20) and exercising, we have
(22)σρ=α+βp0e−p0i−1−mρ−e−p0i−m+1ρe−p0i+1ρ−2e−p0iρ+e−p0i−1ρ,(23)σρ=α+βp0ep0mρep0ρ−e−p0ρep0ρ/2−e−p0ρ/22=α+βp0ep0mρep0ρ/2+e−p0ρ/2ep0ρ/2−e−p0ρ/2,(24)σ=ρα+βp0ep0mρcothp0ρ2,which is the parameter that has been fitted for the layer on the left-end of the domain.


## 4.2. Numerical Scheme for Right-End Boundary Layer

Assume thatps≤M¯<0 and ε≥0 in [0, 1], where M¯ is a positive constant. As a result of this theory, Equations (1) and (2) shows layer behaviour at s=1 for small values of ε. Based on singular perturbation theory for the layer at the right end of the interval [16], we have
(25)limh⟶0νih=ν00+ϕ0−ν01e−P11/ε−iρ+Oε,whereρ=hϵ.In the exponential spline difference scheme Equation (18), we insert a fitting parameter σρ and obtain the fitting parameter
(26)σ=ρα+βp0e−p0mρcothp1ρ2,employing the similar procedure as in the case of the left boundary layer. Now, by simplifying the difference scheme in Equation (19), we get
(27)Piνi−1+Qiνi+Riνi+1+Uiνi−m+1+Viνi−m+Wiνi−m−1=Gi,i=1,2,⋯,N−1,where Pi=σε+αh2qi−1,Qi=−2σε+2βh2qi,Ri=σε+αh2qi+1,Ui=−αhpi−1/2+2βhpi/2+3αhpi+1/2,Vi=2αhpi−1−2αhpi+1,Wi=−3αhpi−1/2−2βhpi/2+αhpi+1/2,Gi=h2αri−1+2βri+αri+1.Following Equation (16), the scheme of Equation (27) can be expressed as
(28)Piνi−1+Qiνi+Riνi+1=Gi−Uiνi−m+1−Viνi−m−Wiνi−m−1,for1≤i≤m−1,Piνi−1+Qiνi+Riνi+1+Uiνi−m+1=Gi−Viνi−m−Wiνi−m−1,fori=m,Piνi−1+Qiνi+Riνi+1+Uiνi−m+1+Viνi−m=Gi−Wiνi−m−1,fori=m+1,Piνi−1+Qiνi+Riνi+1+Uiνi−m+1+Viνi−m+Wiνi−m−1=Gi,form+2≤i≤N−1.The Gauss elimination scheme with partial pivoting is used to solve the aforementioned system of equations.


## 5. Truncation Error

The local truncation error for the scheme Equation (27) is obtained using Taylor’s series expansion, as shown below.
(29)Tih=σε1−2α+βh2νi2ξi+σε121−12ανi4ξi+13−2α+βpiξiνi3ξih4+Oh6,where si≤ξi≤si+1. Clearly, Tih=Oh4 for any random choice of α and β whose sum is equal to 1/2 and α=1/12.


### 5.1. Convergence Analysis

The matrix equation for the left-end boundary layer problem, including the specified boundary conditions in Equation (24), is as follows:
(30)Aν+B+Tih=0,(31)A=Q1R10⋯⋯⋯⋯⋯⋯⋯0P2Q2R20⋯⋯⋯⋯⋯⋯00P3Q3R30⋯⋯⋯⋯⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮Um0⋯PmQmRm0⋯⋯⋯0Vm+1Um+10⋯Pm+1Qm+1Rm+10⋯⋯0Wm+2Vm+2Um+20⋯Pm+2Qm+2Rm+20⋯00Wm+3Vm+3Um+30⋯Pm+3Qm+3Rm+300⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯⋯0WN−2VN−2UN−20⋯PN−2QN−2RN−2⋯⋯⋯0WN−1VN−1UN−10⋯PN−1QN−1,and B=r1,r2,r,⋯rm,rm+1,rm+2,⋯,rN−2,rN−1, where ri=Gi−Uiνi−m+1−Viνi−m−Wiνi−m−1−P1ν0for1≤i≤m−1Gi−Viνi−m−Wiνi−m−1fori=mGi−Wiνi−m−1Gi−RN−1νNfori=m+1form+2≤i≤N−1.
(32)Pi=σε+αh2qi−1,Qi=−2σε+2βh2qi,Ri=σε+αh2qi+1,Ui=−αhpi−12+2βhpi2+3αhpi+12,Vi=2αhpi−1−2αhpi+1,Wi=−3αhpi−12−2βhpi2+αhpi+12,Gi=h2αri−1+2βri+αri−1,fori=1,2,⋯,N−1 and Tih=Oh4, ν=ν1,ν2,⋯,νN−1T,Tih=T1,T2,⋯,TN−1T,andO=0,0,⋯,0T are related vectors of Equation (19).Letν~=ν~1,ν~2,⋯,ν~N−1T≅ν satisfies the equation
(33)Aν~+B=0.Letei=νi−νi,i=11N−1 be the error so that E=e1,e2,⋯,eN−1T=u−U.Subtracting Equation (30) from Equation (33), we obtain the error equation
(34)AE=Tih.LetS̄i be the sum of the ith row elements of the matrix A. Then, we have
(35)S̄i=h2αqi−1+2βqi+αqi+1,for1≤i≤m−1,S̄i=h2−αpi−1+2βpi+3αpi+1+h2αqi−1+2βqi+αqi+1,fori=m,S̄i=h23αpi−1+2βpi−αpi+1+h2αqi−1+2βqi+αqi+1,fori=m+1,S̄i=h2αqi−1+2βqi+αqi+1,form+2≤i≤N−1.Letζ1∗=min1≤i≤Nptiandζ1∗=max1≤i≤Npti, ζ2∗=min1≤i≤Nqti, and ζ2∗=max1≤i≤Nqti.Since0<ε≪1 and with sufficiently small h, it is verified that A is irreducible and monotone [37, 38]. Therefore, A−1existsandA−1≥0.So, we can deduce from Equation (34) that
(36)E≤A−1T.For smallh, we have S̄i≥h22α+βζ2∗,for1≤i≤m−1.
(37)S̄i≥h22α+βζ2∗,fori=m,S̄i≥h22α+βζ2∗,fori=m+1,S̄i≥h22α+βζ2∗,form+2≤i≤N−1.LetAi,k−1bethei,kthelementofA−1, and we define A−1=max1≤i≤N−1∑k=1N−1Ai,k−1andTh=max1≤i≤N−1Tih, sinceAi,k−1≥0and∑k=1N−1Ai,k−1.S̄k=1fori=12,3,⋯,N−1.Hence,(38)∑k=1m−1Ai,k−1≤1min1≤k≤m−1S̄k<1h22α+βζ2∗,i=1,2,3,⋯,m−1.Ai,k−1≤1S̄m<1h22α+βζ2∗,k=m,m+1.Furthermore,(39)∑k=m+2N−1Ai,k−1≤1min1≤k≤m−1S̄k<1h22α+βζ2∗,i=m+2,m+3,⋯,N−1.Using Equations (36), (38), and (39), we get
(40)E≤Oh2.As a result, it demonstrates that the proposed exponential spline scheme is second-order convergent. The same procedure can also be used to examine the convergence analysis for the right end boundary layer.


## 5.1. Convergence Analysis

The matrix equation for the left-end boundary layer problem, including the specified boundary conditions in Equation (24), is as follows:
(30)Aν+B+Tih=0,(31)A=Q1R10⋯⋯⋯⋯⋯⋯⋯0P2Q2R20⋯⋯⋯⋯⋯⋯00P3Q3R30⋯⋯⋯⋯⋯0⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮Um0⋯PmQmRm0⋯⋯⋯0Vm+1Um+10⋯Pm+1Qm+1Rm+10⋯⋯0Wm+2Vm+2Um+20⋯Pm+2Qm+2Rm+20⋯00Wm+3Vm+3Um+30⋯Pm+3Qm+3Rm+300⋮⋮⋮⋮⋮⋮⋮⋮⋮⋮0⋯⋯0WN−2VN−2UN−20⋯PN−2QN−2RN−2⋯⋯⋯0WN−1VN−1UN−10⋯PN−1QN−1,and B=r1,r2,r,⋯rm,rm+1,rm+2,⋯,rN−2,rN−1, where ri=Gi−Uiνi−m+1−Viνi−m−Wiνi−m−1−P1ν0for1≤i≤m−1Gi−Viνi−m−Wiνi−m−1fori=mGi−Wiνi−m−1Gi−RN−1νNfori=m+1form+2≤i≤N−1.
(32)Pi=σε+αh2qi−1,Qi=−2σε+2βh2qi,Ri=σε+αh2qi+1,Ui=−αhpi−12+2βhpi2+3αhpi+12,Vi=2αhpi−1−2αhpi+1,Wi=−3αhpi−12−2βhpi2+αhpi+12,Gi=h2αri−1+2βri+αri−1,fori=1,2,⋯,N−1 and Tih=Oh4, ν=ν1,ν2,⋯,νN−1T,Tih=T1,T2,⋯,TN−1T,andO=0,0,⋯,0T are related vectors of Equation (19).Letν~=ν~1,ν~2,⋯,ν~N−1T≅ν satisfies the equation
(33)Aν~+B=0.Letei=νi−νi,i=11N−1 be the error so that E=e1,e2,⋯,eN−1T=u−U.Subtracting Equation (30) from Equation (33), we obtain the error equation
(34)AE=Tih.LetS̄i be the sum of the ith row elements of the matrix A. Then, we have
(35)S̄i=h2αqi−1+2βqi+αqi+1,for1≤i≤m−1,S̄i=h2−αpi−1+2βpi+3αpi+1+h2αqi−1+2βqi+αqi+1,fori=m,S̄i=h23αpi−1+2βpi−αpi+1+h2αqi−1+2βqi+αqi+1,fori=m+1,S̄i=h2αqi−1+2βqi+αqi+1,form+2≤i≤N−1.Letζ1∗=min1≤i≤Nptiandζ1∗=max1≤i≤Npti, ζ2∗=min1≤i≤Nqti, and ζ2∗=max1≤i≤Nqti.Since0<ε≪1 and with sufficiently small h, it is verified that A is irreducible and monotone [37, 38]. Therefore, A−1existsandA−1≥0.So, we can deduce from Equation (34) that
(36)E≤A−1T.For smallh, we have S̄i≥h22α+βζ2∗,for1≤i≤m−1.
(37)S̄i≥h22α+βζ2∗,fori=m,S̄i≥h22α+βζ2∗,fori=m+1,S̄i≥h22α+βζ2∗,form+2≤i≤N−1.LetAi,k−1bethei,kthelementofA−1, and we define A−1=max1≤i≤N−1∑k=1N−1Ai,k−1andTh=max1≤i≤N−1Tih, sinceAi,k−1≥0and∑k=1N−1Ai,k−1.S̄k=1fori=12,3,⋯,N−1.Hence,(38)∑k=1m−1Ai,k−1≤1min1≤k≤m−1S̄k<1h22α+βζ2∗,i=1,2,3,⋯,m−1.Ai,k−1≤1S̄m<1h22α+βζ2∗,k=m,m+1.Furthermore,(39)∑k=m+2N−1Ai,k−1≤1min1≤k≤m−1S̄k<1h22α+βζ2∗,i=m+2,m+3,⋯,N−1.Using Equations (36), (38), and (39), we get
(40)E≤Oh2.As a result, it demonstrates that the proposed exponential spline scheme is second-order convergent. The same procedure can also be used to examine the convergence analysis for the right end boundary layer.


## 6. Numerical Experiments

Three problems are selected to validate the competence of the proposed scheme. The maximum absolute errors are obtained using the principle of double meshEiN=max0≤i≤NνiN−ν2i2N and the order of convergence rN=log2EiN/Ei2N for the examples, which are tabulated in Tables 1–3. Furthermore, the calculated arithmetic solutions to the considered experiments are depicted by graphs, with and without fitting parameters for various values of ε and δ.Example 1.
εν′′+ν′s−δ−νs=0withνs=1,−δ≤s≤0andν1=0.Example 2.
εν′′+0.25ν′s−δ−νs=0withνs=1,−δ≤s≤0andν1=−1.Example 3.
εν′′−ν′s−δ−νs=0withνs=1,−δ≤s≤0andν1=−1.Table 1
The maximum absolute errorsEiN of Example 1 for δ=0.03.
ε↓N=100N=200N=400N=8002−13.3803e-0068.4512e-0072.1128e-0075.2821e-0081.99992.00002.0000—2−22.1392e-0055.3490e-0061.3373e-0063.3433e-0071.99971.99992.0000—2−31.3632e-0043.4100e-0058.5264e-0062.1317e-0061.99921.99981.9999—2−47.8174e-0041.9588e-0044.8997e-0051.2251e-0051.99671.99921.9998—2−53.6000e-0039.1639e-0042.2978e-0045.7476e-0051.97401.99501.9992—2−61.2200e-0023.1000e-0037.9274e-0041.9858e-0041.97651.96731.9971—2−72.6800e-0027.5000e-0031.9000e-0034.8675e-0041.83731.98091.9647—Table 2
The maximum absolute errorsEiN of Example 2 for δ=0.03.
ε↓N=100N=200N=400N=8002−11.0966e-0062.7415e-0076.8538e-0081.7134e-0082.00002.00002.0000—2−23.5638e-0068.9099e-0072.2275e-0075.5686e-0081.99992.00002.00002−39.9549e-0062.4886e-0066.2215e-0071.5554e-0072.00012.00002.0000—2−42.2985e-0055.7462e-0061.4365e-0063.5913e-0072.00002.00012.0000—2−51.3263e-0043.3181e-0058.2966e-0062.0743e-0061.99901.99981.9999—2−67.2235e-0041.8102e-0044.5281e-0051.1322e-0051.97651.99651.9992—2−73.3000e-0038.2019e-0042.0550e-0045.1402e-0052.00841.99681.9992—Table 3
The maximum absolute errorsEiN of Example 3 for δ=0.03.
ε↓N=100N=200N=400N=8002−12.3934e-0065.9833e-0071.4958e-0073.7396e-0082.00002.00002.0000—2−24.1084e-0061.0271e-0062.5682e-0076.4205e-0082.00001.99972.0000—2−36.6395e-0051.6604e-0054.1514e-0061.0379e-0061.99951.99991.9999—2−44.2021e-0041.0515e-0042.6293e-0056.5735e-0061.99871.99971.9999—2−52.1000e-0035.1586e-0041.2911e-0043.2288e-0052.02531.99841.9995—2−69.0000e-0032.3000e-0035.7426e-0041.4375e-0041.96832.00191.9981—2−72.9000e-0028.0000e-0032.1000e-0035.1971e-0041.83731.85802.0146—


## 7. Discussions and Conclusion

The delay differential equation of order two is considered, which has a large delay in the convention term and layers at the left and right ends of the interval. An exponential spline finite difference scheme is constructed using the continuity of its first-order derivative condition at the joint nodes. We considered the case whenδε is of order Oε, layer’s behaviour can change and even be destroyed, or the solution can exhibit oscillatory behaviour. To address these drawbacks in solutions, we experimented with a new scheme based on the special mesh introduced in [21], which involves adding a parameter to the exponential spline method. Tables 1–3 illustrate the maximum errors as well as the order of convergence for Examples 1–3. Table 4 compares the maximum absolute errors in Example 1 to other methods for different shift parameter values. Figures 1–6 depict diagrams of the solutions to the experiments for various values of δ, and we compared the graphs with the fitting parameter for various values of δ=Oε to the graphs without the fitting parameter. When the δ value is larger than the ε, the solution includes oscillations, as shown in Figures 2 and 4 (without fitting parameter), whereas in Figures 1 and 3 (with parameter), the oscillations are controlled in solutions, and behaviour of the layer is well-preserved. However, when dealing with a right-end layer, as shown in Figures 5 and 6, layer behaviour of the solution is preserved even for δ=Oε. Finally, we found that exponential spline fitting parameter method with special mesh had a significant improvement in terms of regulating the oscillations in the solutions of delay differential equations.Table 4
The maximum absolute errorsEiNin Example 1 for ε=0.1.
N↓δ=0.03δ=0.05δ=0.08Our method1002.4313e-0043.0524e-0043.3301e-0042006.0837e-0057.6379e-0058.3324e-0053002.7043e-0053.3952e-0053.7039e-0054001.5213e-0051.9099e-0052.0836e-0055009.7364e-0061.2224e-0051.3335e-005Results in [21]1001.7830e-0022.5306e-0023.5989e-0022009.5140e-0031.3580e-0021.9250e-0023006.4860e-0039.2760e-0031.3132e-0024004.9190e-0037.0420e-0039.9650e-0035003.9620e-0035.6740e-0038.0280e-003Figure 1
Solution profile for Example1 with ε=0.01 and diverse values of δ with fitting parameter.Figure 2
Solution profile for Example1 with ε=0.01 and diverse values of δ without fitting parameter.Figure 3
Solution profile for Example2 with ε=0.01 and diverse values of δ with fitting parameter.Figure 4
Solution profile for Example2 with ε=0.01 and diverse values of δ without fitting parameter.Figure 5
Solution profile for Example3 with ε=0.01 and diverse values of δ with fitting parameter.Figure 6
Solution profile for Example3 with ε=0.01 and diverse values of δ without fitting parameter.



---

*Source: 9291834-2022-08-16.xml*