# An Exponentially Fitted Numerical Scheme via Domain Decomposition for Solving Singularly Perturbed Differential Equations with Large Negative Shift

Authors: Ababi Hailu Ejere; Gemechis File Duressa; Mesfin Mekuria Woldaregay; Tekle Gemechu Dinka
Journal: Journal of Mathematics
 (2022)

---

## Abstract

In this study, we focus on the formulation and analysis of an exponentially fitted numerical scheme by decomposing the domain into subdomains to solve singularly perturbed differential equations with large negative shift. The solution of problem exhibits twin boundary layers due to the presence of the perturbation parameter and strong interior layer due to the large negative shift. The original domain is divided into six subdomains, such as two boundary layer regions, two interior (interfacing) layer regions, and two regular regions. Constructing an exponentially fitted numerical scheme on each boundary and interior layer subdomains and combining with the solutions on the regular subdomains, we obtain a second orderε-uniformly convergent numerical scheme. To demonstrate the theoretical results, numerical examples are provided and analyzed.


---

## Body

## 1. Introduction

In science and engineering, many phenomena can be modeled and described by observing the relation between causes and effects. When the cause is small and its effect is large, the relation has considerable physical system[1]. A mathematical equation associated with differential equations involving small parameter causing large effect in the problem is said to be singularly perturbed differential equation, whereas the simplified differential equation (the one that does not include the small parameter) is called the unperturbed model.Depending on the influence of the small parameter, perturbation problems can be categorized as regular and singular perturbation problems. A regular perturbation problem is the one in which the perturbed problem for small and nonzero values of the perturbation parameterε is qualitatively the same as the unperturbed problem for ε=0. A singular perturbation problem is a problem, which is qualitatively different from the unperturbed one. In this case, we can obtain an asymptotic, but possibly divergent expansion of the solution, which depends singularly on the perturbation parameter [2, 3]. A subclass of differential equation in which the second-order derivative term is multiplied by the perturbation parameter and involves at least one delay parameter is said to be singularly perturbed delay differential equation. Such types of equations are used in the mathematical modeling of various physical phenomena, for instance, for the modeling of human pupil-light reflex [4], in studying bi-stable devices [5], in neuronal variability [6], and in variational problems in control theory [7].The solution of a singularly perturbed problem varies rapidly in some region and varies slowly in other parts of the problem’s domain. The region where the solution varies rapidly is called the inner region and the region where the solution varies slowly is known as the outer region. In physical system, we observe several phenomena characterized by such rapid variations of quantities, for instance, occurrence in shock waves in gas motions, in boundary layer flow along the surface of a body, and in edge effects in the deformation of elastic plates [1]. Due to the rapid variation, one encounters difficulties to obtain satisfactory solution to the problem.For the response of the above difficulties, many research articles are available in literature, which mostly cover differential equations with small shift or without shifting parameter or convection-diffusion type problems. For instance, Bansal and Sharma [8] proposed a parameter uniform numerical method based on the nonstandard finite difference method to capture the significant properties of singularly perturbed parabolic partial differential equations with general shift arguments and obtained parameter-uniformly convergent result. In [9], singularly perturbed nondelayed reaction-diffusion problem was treated applying Galerkin finite element technique on a piece-wise uniform Shishkin mesh using linear basis functions. Arora and Kaur [10] solved singularly perturbed differential equation with small shift applying the collocation method with the modified B-spline basis functions.In [11], the authors treated singularly perturbed nonlinear differential-difference equations with negative shift. To simplify the difficulties due to the presence of nonlinearity, they obtained linear differential equation from the nonlinear one using the quasi-linearization process. To tackle the shift term, Taylor’s series expansion is used, and the fitted mesh method is applied to resolve the problem due to perturbation parameter. In [12], singularly perturbed differential equation involving both positive and negative shift arguments was treated using the fitted operator method as well as the fitted mesh method and obtained ε-uniform numerical results. In [13], a singularly perturbed differential-difference equation with small negative and positive shifts was solved by using modified Numerov’s method. In [14], a singularly perturbed reaction-diffusion problem was treated by developing a fourth-order exponentially fitted numerical scheme on a uniform mesh. Kadalbajoo and Sharma [15] treated singularly perturbed differential equation with small shifts of mixed typed using numerical approach on a uniform mesh. In [16], the authors constructed a fitted operator finite difference method using the nonstandard finite difference method to solve singularly perturbed differential equation involving both small negative and positive shifts. Gupta et al. [17] solved time-dependent singularly perturbed differential-difference convection-diffusion equations by proposing a numerical scheme using implicit Euler method in time direction and hybrid finite difference scheme on piece-wise uniform spatial mesh.However, there are only few research works on singularly perturbed differential equations with large delay. In [18], the authors solved such problem by suggesting an initial value technique and obtained almost second-order convergence with respect to ε. The same type of problem was also studied in [19], where hp finite element method was applied and obtained as a convergent result. In [20], the problem is treated by constructing a numerical method on a piece-wise uniform Shishkin mesh using classical finite difference methods. By this method, it was indicated that the boundary layers and interior layers were resolved and first-order convergent result was obtained. In [21], an exponentially fitted numerical method is constructed applying the Numerov finite difference method. The method resolves the boundary and interior layers and converges uniformly with respect to ε. By suggesting iterative techniques for the boundary value problem, a convergent numerical result was also reported by Selvi and Ramanujam [22]. In [23], a numerical method was constructed by defining a fitting comparison problem and replacing ε by a fitting factor, and by such technique, a convergent result was reported.Our aim in this study is to develop and analyze an exponentially fitted numerical scheme to solve singularly perturbed differential equations with large negative shift by decomposing the domain into subdomains. Since the solution of the considered problem involves two boundary layers and interior layer, we decompose the original domain into two boundary layer subdomains, two interior (interfacing) layer subdomains, and two regular subdomains. Then, on each boundary and interior layer subdomains, we formulate an exponentially fitted numerical scheme and on each regular subdomains, and we solve the reduced problem by settingε to be zero. Combining and analyzing these results give us a second-order ε-uniformly convergent method.The study is organized as follows. In Section2, the considered model problem is presented. In Section 3, properties of the analytic results are described briefly. Description of the numerical methods and derivation of the schemes are discussed in detail in Section 4. Demonstration of the proposed method by numerical examples is presented in Section 5, and our study is concluded in Section 6.


### 1.1. Notations

Throughout this study, we usedC as a generic constant, which is independent of the perturbation parameter and the mesh elements. The norm ⋅Ω is used to represent a continuous maximum norm.


## 1.1. Notations

Throughout this study, we usedC as a generic constant, which is independent of the perturbation parameter and the mesh elements. The norm ⋅Ω is used to represent a continuous maximum norm.


## 2. Statement of the Problem

Consider a second-order singularly perturbed delay differential equations of the form(1)Lu≔−εd2udx2+rxux+sxux−1=wx,x∈0,2.The interval and boundary conditions are given as(2)ux=ψx,x∈−1,0,u2=γ,where u∈Ω=C00,2∩C10,2∩C40,1∪1,2 and the functions rx, sx, and wx are sufficiently smooth on 0,2 such that(3)rx+sx≥2α>0,forrx>0andsx<0.Also,ψx is a smooth function on −1,0 and γ is a given constant, which is independent of ε. Problems (1)-(2) can be rewritten as(4a)L1u≔−εd2udx2+rxux=wx−sxψx−1,x∈0,1,(4b)L2u≔−εd2udx2+rxux+sxux−1=wx,x∈1,2.Withux=ψx on −1,0, u1−=u1+, u′1−=u′1+, and u2=γ, where u− and u+ represent the left and right side limit of u at x=1, respectively. Due to the influence of the perturbation parameter, the solution ux of the boundary value problems (1)-(2) exhibits strong boundary layers at x=0 and x=2, and due to the large negative shift, strong interior layers occur on left and right sides of x=1 [24]. Moreover, such type of problem has a unique solution [25].


## 3. Analytical Results

Lemma 1 (continuous maximum principle).
For the smooth functionsrx and sx satisfying (3), let ϕ be in Ω such that ϕ0≥0, ϕ2≥0, and Lϕx≥0 on 0,2. Then, ϕx≥0 on 0,2.Proof.
Lety∈Ω be given such that ϕy=minx∈0,2ϕx, and for the contrary, suppose ϕy<0. By the given conditions, y∉0,2. It follows from calculus that dϕ/dy=0 and d2ϕ/dy2≥0. Consequently, we have(5)L1ϕy=−εd2ϕdy2+ryϕy<0,y∈0,1,L2ϕy=−εd2ϕdy2+ryϕy+syϕy−1≤−εd2ϕdy2+ry+syϕy<0,y∈1,2,which is a contradiction to the hypothesis. Therefore, it follows that ϕy≥0 and ϕx≥0, for all x∈0,2.Lemma 2 (stability result).
Let the conditions in (3) hold true and ux be any function in Ω. Then, for all x∈0,2, we have ux≤1/αLu+maxu0,u2.Proof.
Let us define two barrier functions asϕ±x=1/αLu+maxu0,u2±ux. By these functions, we have(6)ϕ±0=1/αLu+maxu0,u2±u0≥1/αLu≥0,ϕ±2=1/αLu+maxu0,u2±u2≥1/αLu≥0.
And by (4a4b), for all x∈0,2, we have(7)L1ϕ±x=−εd2ϕ±dx2+rxϕ±x=∓εd2udx2+rxαL1u+rxmaxu0,u2±rxux=±wx−sxψx−1+rxαwx−sxψx−1+rxmaxu0,u2≥rxmaxu0,u2≥0,L2ϕ±x=−εd2ϕ±dx2+rxϕ±x+sxϕ±x−1=∓εd2udx2+rxαL2ux+sxαL2ux−1+rx+sxmaxu0,u2±rxux±sxux−1=±wx+rxαL2ux+sxαL2ux−1+rx+sxmaxu0,u2≥rx+sxmaxu0,u2≥2αmaxu0,u2≥0.
Applying Lemma1, it follows that ϕ±x≥0, for all x∈0,2, and hence, ux≤1/αLu+maxu0,u2, for all x∈0,2.Lemma 3 (bBound for the derivatives of the solution).
Let the conditions in (3) hold true and ux be the solution of (1)-(2). Then, we have(8)ukx≤Cε−k/2u+w,ifk=0,1,ukx≤Cε−k/2u+w+wk−2εk−2/2,ifk=2,3,4.Proof.
The case fork=0 is the result of Lemma 2. To handle the case for k=1, let x∈0,1 and construct an associated neighborhood Nx=ι,ι+ε, such that x∈Nx and Nx⊂0,1. By the mean value theorem, for some a∈Nx¯, we have(9)u′a=uι+ε−uιι+ε−ι=ε−1/2uι+ε−uι≤2ε−1/2u.
Then, forx∈Nx, we have(10)u′x=u′a+u′x−u′a=u′a+∫axu″tdt=u′a+1ε∫ax−wt+rtut+stψt−1dt,t∈a,x≤u′a+1εw+ru+sψε≤Cε−1/2u+w.
In a similar procedure, forx∈1,2, we construct neighborhood Nx=ι,ι+ε⊂1,2. Then, by the mean value theorem as above, we obtain u′x≤Cε−1/2u+w.
The case, fork=2, follows rearranging terms in (1) as(11)u″x=ε−1rxux+sxux−1−wx.
From this, we haveu″x≤Cε−1u+w. Differentiating both sides of (11) once and twice, we obtain(12)u′″x=ε−1r′xux+rxu′x+s′xux−1+sxu′x−1−w′x≤ε−1C+Cε−1/2u+w+w′≤Cε−3/2u+w+w′ε1/2,u4x=ε−1r′xux′+rxu′x′+s′xux−1′+sxu′x−1′−w″x≤Cε−2u+w+w″ε.
For more sharp bound on the derivative of the solution via Shishkin decomposition, we refer the approaches in [20, 26].


## 4. Numerical Method

### 4.1. Description of the Numerical Method

Since the solution of the problem in (1) exhibits two boundary layers and interior layer, we divide the domain into six subdomains, such as two boundary layer subdomains, two interior layer (left and right side of x=1) subdomains, and two outer region subdomains. The boundary and interior layer problems can be transformed to regular problems by appropriate transformations using stretching of variables. Consider the asymptotic expansion solution of the problem in (1)-(2) as(13a)ux,ε=∑i=0∞ui1x+vi1τ1+vi2τ2εi,x∈0,1,(13b)ux,ε=∑i=0∞ui2x+vi3τ3+vi4τ4εi,x∈1,2,where τ1=x/ε, τ2=1−x/ε, τ3=x−1/ε, and τ4=2−x/ε. Then, the corresponding zero-order asymptotic expansions of (13a)-(13b) are given by(14a)ux=u01x+v01τ1+v02τ2,x∈0,1,(14b)ux=u02x+v03τ3+v04τ4,x∈1,2,where u01x=wx−sxψx−1/rx and u02x=wx−sxu0x−1/rx are the asymptotic solutions of the reduced problem on 0,1 and 1,2 respectively, which do not satisfy the conditions in (2). From the solution of the reduced problem, the value of ux=1=μ (say) can be obtained. In (14a)-(14b), v01 is the left boundary layer function, v02 is the left side of x=1 interior layer function, v03 is the right side of x=1 interior layer function, and v04 is the right boundary layer function.The functionsv01, v02, v03, and v04 satisfy, respectively, the following transformed homogeneous differential equations:(15a)−d2v01τ1dτ12+r0v01τ1=0;τ1∈0⋅∞,(15b)−d2v02τ2dτ22+r1v02τ2=0;τ2∈0⋅∞,(15c)−d2v03τ3dτ32+r1v03τ3=0;τ3∈0⋅∞,(15d)−d2v04τ4dτ42+r2v04τ4=0;τ4∈0⋅∞.Solving each subequations in (15a15b15c15d) at corresponding terminal points and by (14a)-(14b), we obtain the zeros order solution of (1)-(2) as(16a)ux=u01x+λ1e−r0/εx+λ2e−r1/ε1−x,x∈0,1,(16b)ux=u02x+λ3e−r1/εx−1+λ4e−r2/ε2−x,x∈1,2,so that applying the interval and boundary conditions, we can get λ1, λ2, λ3, and λ4 as(17a)λ1=ψ0−u010−μ−u011e−r1/ε1−e−r0/ε+r1/ε,(17b)λ2=μ−u011−ψ0−u010e−r0/ε1−e−r0/ε+r1/ε,(17c)λ3=μ−u021−γ−u022e−r2/ε1−e−r1/ε+r2/ε,(17d)λ4=γ−u022−μ−u021e−r1/ε1−e−r1/ε+r2/ε.


### 4.2. Derivation and Properties of the Numerical Scheme

We divide the interval0,2 into N equal parts with uniform mesh length h. Suppose 0=x0,x1,x2,…,xN/2=1,xN/2+1,xN/2+2,…,xN=2 be the mesh points. Then, we have xi=ih,i=0,1,2,…,N. Let us choose the terminal points as xn1=ε, xn2=1−ε, xn3=1+ε, and xn4=2−ε. The left boundary layer is in the interval 0,ε, the right boundary layer is in the interval 2−ε,2, the left interior layer is in the interval 1−ε,1, the right interior layer is in the interval 1,1+ε, the left outer region is in the interval ε,1−ε, and the right outer region is in the interval 1+ε,2−ε.Now, atx=xi, (1) can be written as −εd2u/dxi2+rxiuxi+sxiuxi−1=wxi. Approximating the differential operator by the central finite difference as d2u/dxi2=ui−1−2ui+ui+1/h2−h2/12ui4+R, where R=−h4/360ui6t, t∈xi−h,xi+h, and introducing a fitting factor σ to obtain an ε-independently convergent solution, we obtain(18)Luxi≔−εσh2uxi−1+2εσh2+rxiuxi−εσh2uxi+1+εσh212u4xi−εσR+sxiuxi−1=wxi,∀i=1,2,…,N−1.From this, we have(19a)L1ui≔Aiui−1+Biui+Aiui+1=Di+Ti,∀i=1,2,…,N2−1,(19b)L2ui≔Aiui−1+Biui+Aiui+1=Ei+Ti,∀i=N2+1,N2+1,…,N−1,where Ai=−εσ/h2, Bi=2εσ/h2+ri, Di=wi−siψi−N/2, Ei=wi−siui−N/2, and Ti=−εσh2/12u4t+εσOh4 for t∈xi−h,xi+h.Case 1.
Left boundary layer.
On the interval0,ε, we introduce a fitting factor σ1 in (19a) as(20)−εσ1ui−1−2ui+ui+1h2+riui=Di+Tifori=1,2,…,N/8.
To determine the fitting factorσ1 on the left boundary layer, we use the left boundary layer asymptotic solution with the left outer solution as(21)ui=u01xi+λ1e−r0/εxi.
We assume that the solution converges uniformly to the solution of (1)-(2). From (20), we have −σ1/ρ2ui−1−2ui+ui+1=Di+Ti−riui, where ρ=h/ε. Taking the limit as h⟶0, it becomes(22)−limh⟶0σ1ρ2ui−1−2ui+ui+1=limh⟶0Di+Ti−riui.
From (21), we have ui−1=u01xi−1+λ1e−r0/εxier0ρ and ui+1=u01xi+1+λ1e−r0/εxie−r0ρ. Inserting these in (22) and simplifying gives(23)σ1=r0ρ/2sinhr0ρ/22,which is the fitting factor in the interval 0,ε. Having this fitting factor in (20), we obtain(24)−εσ1h2ui−1+2εσ1h2+riui−εσ1h2ui+1=Di+Ti,fori=1,2,…N8.Case 2.
Left outer region.
On the intervalε,1−ε, we have the left outer region. Setting ε=0, (19a) becomes(25)riui=Di,fori=N8+1,N8+2,…,3N8,which is the left outer region scheme.Case 3.
Left-side interior layer.
On the interval1−ε,1, the interior layer will be on the left-hand side of x=1. Introducing a fitting factor σ2 in (19a), we have(26)−εσ2ui−1−2ui+ui+1h2+riui=Di+Ti,fori=3N8+1,3N8+2,…,N2.
To determine the fitting factorσ2 on the left interior layer, we use the corresponding interior layer asymptotic solution with the left outer region solution as(27)ui=u01xi+λ2e−r1/ε1−xi.
Assuming that the solution converges uniformly to the solution of (1)-(2), from (26), we have(28)−limh⟶0σ2ρ2ui−1−2ui+ui+1=limh⟶0Di+Ti−riui,where ρ=h/ε. From (27), we have ui−1=u01xi−1+λ2e−r1/ε1−xie−r1ρ and ui+1=u01xi+1+λ2e−r1/ε1−xier1ρ. Substituting these in (28) and simplification gives a fitting factor in left side of x=1 interior layer as(29)σ2=r0ρ/2sinhr1ρ/22.
With this fitting factor, we have(30)−εσ2h2ui−1+2εσ2h2+riui−εσ2h2ui+1=Di+Ti,fori=3N8+1,3N8+2,…N2.Case 4.
Right-side interior layer.
On the interval1,1+ε, the interior layer will be on the right side of x=1. If we introduce a fitting factor σ3 in (19b), then we obtain(31)−εσ3ui−1−2ui+ui+1h2+riuxi=Ei+Ti,fori=N2+1,N2+2,…,5N8.
To determineσ3, we use the right-side interior layer asymptotic solution with the right outer solution as(32)ui=u02xi+λ3e−r1/εxi−1.
Again, assuming that the solution converges uniformly to the solution of (1)-(2) from (31), we have(33)−limh⟶0σ3ρ2ui−1−2ui+ui+1=limh⟶0Ei+Ti−riui,where ρ=h/ε. Using (32), we have ui−1=u02xi−1+λ3e−r1/εxi−1er1ρ and ui+1=u02xi+1+λ3e−r1/εxi−1e−r1ρ. Putting these in (33) and simplifying gives the fitting factor σ3 in the right side of x=1 interior layer as(34)σ3=r0ρ/2sinhr1ρ/22.
Having this fitting factor, (31) becomes(35)−εσ3h2ui−1+2εσ3h2+riui−εσ3h2ui+1=Ei+Ti,fori=N2+1,N2+2,…,5N8.Case 5.
Right outer region.
On the interval1+ε,2−ε, we have the right outer region. Setting ε=0, (19b) is reduced to(36)riui=Ei,fori=5N8+1,5N8+2,…,7N8,which is the right outer region scheme.Case 6.
Right boundary layer.
On the interval2−ε,2, introducing a fitting factor σ4 in (19b), we have(37)−εσ4ui−1−2ui+ui+1h2+riui=Ei+Ti,fori=7N8+1,7N8+2,…,N.
To determineσ4, use the right boundary layer asymptotic solution with left outer solution as(38)ui=u02xi+λ4e−r2/ε2−xi.
Assuming that the solution converges uniformly to the solution of (1)-(2), we have(39)−limh⟶0σ4ρ2ui−1−2ui+ui+1=limh⟶0Ei+Ti−riui,where ρ=h/ε. Using (38), we have ui−1=u02xi−1+λ4e−r2/ε2−xier2ρ and ui+1=u02xi+1+λ4e−r2/ε2−xie−r2ρ, and by putting these into (39) and simplifying, we get the fitting factor:(40)σ4=r0ρ/2sinhr2ρ/22.
With this fitting factor, we have the right boundary layer scheme as(41)−εσ4h2ui−1+2εσ4h2+riui−εσ4h2ui+1=Ei+Ti,fori=7N8+1,7N8+2,…,N.
The solution of the considered singularly perturbed differential equation with large negative shift can be obtained by solving the three term recurrence relations (24), (30), (35), and (42) together with the outer layer schemes (25) and (26).


### 4.3. Discrete Stability Analysis

To establish the computational stability of the proposed numerical scheme, we follow the approaches of [27, 28]. From (19a), consider the recurrence relation:(42)Aiui−1+Biui+Aiui+1=Di,i=1,2,…,N2,where Ai=−εσ/h2 and Bi=2εσ/h2+ri with σ as in (23) or (29) and subject to the boundary conditions u0=ψ0 and uN/2=uxN/2. Now, we set(43)ui=Piui+1+Qi,fori=N−1,N−2,…,2,1,where Pi=Pxi and Qi=Qxi are determined as follows. From (43), we have ui−1=Pi−1ui+Qi−1, and substituting in (42), we obtain(44)ui=−AiAiPi−1+Biui+1+Di−AiQi−1AiPi−1+Bi.From (43) and (44), we can obtain(45a)Pi=−AiAiPi−1+Bi,(45b)Qi=Di−AiQi−1AiPi−1+Bi.Using the initial conditionu0=ψ0=P0u1+Q0, we can solve (45a) and (45b). If we choose P0=0, it follows that Q0=u0, and using these and uN/2=uxN/2, we can compute Pi and Qi, for all i=1,2,…,N/2. Suppose that, in computing Pi from (45a), a small error ei has been occurred. Then, we have Pi∗=Pi+ei and Pi−1∗=Pi−1+ei−1. However, we are computing Pi∗=−Ai/AiPi−1∗+Bi. So, we have(46)ei=Pi∗−Pi=−AiAiPi−1∗+Bi+AiAiPi−1+Bi=AiPiAiPi−1+Aiei−1+Biei−1.From the assumption in (3), we have rx>0, and since Ai=−εσ/h2 and Bi=2εσ/h2+ri, it follows that Bi>2Ai, for i=1,2,…,N−1. Then, from (45a), we have P1=−A1/A1P0+B1<1 and P2=−A2/A2P1+B2<−A2/A2+B2<1 as P1<1 and successively; if we continue in this manner, we get Pi<1, for i=0,1,2,…,N−1. Thus, under the assumption that there is small initial error, we have ei=AiPi/AiPi−1+Aiei−1+Biei−1, which implies that ei<ei−1. This indicates that (45a) is stable. In a similar procedure, we can show that (45b) is also stable, and following the same procedure for the numerical scheme from (19b), we can obtain a similar result so that the stability estimate of the proposed scheme is established.


### 4.4. Convergence Analysis

To provide the convergence analysis of the proposed numerical scheme, we follow the approaches in [29, 30]. From (19a), we have a system of equations as(47)−εσh2ui−1+2εσh2+riui−εσh2ui+1+fi+Ti=0,∀i=1,2,…,N2,where σ is as in (23) or (29), fi=−wi+siψi−N/2, and Ti=εσh2/12ui4t,t∈xi−h,xi+h. Using the boundary conditions u0=ψ0 and u1=uxN/2 in (47), we can get a system of equations in the matrix form as(48)G+Hu+M+Th=0,where G=εσ/h22−10…0−12−1…0⋮⋮⋮…⋮⋮⋮⋮…−100…−12 is a tri-diagonal matrix of order N/2−1 and H=r100…00r20…0⋮⋮⋮…⋮⋮⋮⋮…000…0r1−N/2 is a diagonal matrix of order N/2−1.And the associated vectors of (48) are M=f1−εσψ0/h2,f2,…,fN/2−2,fN/2−1−εσ/h2μt, and Th=T1,T2,…,TN/2−1t. Now, let us consider U as the approximation of u satisfying the system of equations.(49)G+HU+M=0.Suppose that theith truncation error be ei=Ui−ui, i=1,2,…,N/2−1, such that E=e1,e2,…,eN/2t=U−u. Taking the difference between (48) and (49) yields(50)G+HE=Th.Letri≤K,∀i=1,2,…,N/2−1, for arbitrary constant K, and let Pi,i be the i,ith entry of matrix P. Then, since ri>0 for a sufficiently small h, we have Pi,i≠0 and 2εσ/h2+Pi,i≠0, for i=1,2,…,N/2−1. These indicate that G+H is an irreducible matrix.In matrixG+H, if Si is the sum of entries on the ith row, then we have Si=εσ/h2+ri,i=1,N/2−1, and Si=ri,i=2,3,…,N/2−2. Let us define K1=min1≤i≤N/2−1ri and K2=max1≤i≤N/2−1ri, which imply that 0<K1≤K≤K2. Thus, for a sufficiently small value of h, G+H is monotone so that it is nonsingular and G+H−1≥0. From (50), we have(51)E≤G+H−1Th.In each row ofG+H, we can see that Si>h2εσ/h2+ri, for i=1,N/2−1, and Si>h2ri, for i=2,3,…,N/2−2. Suppose that G+Hi,j−1 be the i,jth entry of G+H−1, and let us define G+H−1=max1≤i≤N/2−1∑j=1N/2−1G+Hi,j−1 and Th=max1≤i≤N/2−1Ti. Since ∑j=1N/2−1G+Hi.j−1Sj=1,i=1,2,…,N/2−1, and G+H−1≥0, it follows that G+Hi,j−1≤1/Sj<1/h2εσ/h2+ri,j=1,N/2−1, and ∑j=2N/2−2G+Hi,j−1≤1/min2≤j≤N/2−2≤1/h2ri,i=1,2,…,N/2−2. From (51), using the ith row sum, we have(52)E≤G+H−1Th=max1≤i≤N/2−1∑j=1N/2−1G+Hi,j−1×max1≤i≤N/2−1Ti≤1εσ/h2+ri+1ri+1εσ/h2+ri×εσui4t12h2≤Ch2.Thus, for sufficiently chosen large value ofC, this result shows the second-order convergence of the proposed scheme on (0, 1], and by a similar analysis, the convergence of the numerical scheme can be provided on (1, 2).


## 4.1. Description of the Numerical Method

Since the solution of the problem in (1) exhibits two boundary layers and interior layer, we divide the domain into six subdomains, such as two boundary layer subdomains, two interior layer (left and right side of x=1) subdomains, and two outer region subdomains. The boundary and interior layer problems can be transformed to regular problems by appropriate transformations using stretching of variables. Consider the asymptotic expansion solution of the problem in (1)-(2) as(13a)ux,ε=∑i=0∞ui1x+vi1τ1+vi2τ2εi,x∈0,1,(13b)ux,ε=∑i=0∞ui2x+vi3τ3+vi4τ4εi,x∈1,2,where τ1=x/ε, τ2=1−x/ε, τ3=x−1/ε, and τ4=2−x/ε. Then, the corresponding zero-order asymptotic expansions of (13a)-(13b) are given by(14a)ux=u01x+v01τ1+v02τ2,x∈0,1,(14b)ux=u02x+v03τ3+v04τ4,x∈1,2,where u01x=wx−sxψx−1/rx and u02x=wx−sxu0x−1/rx are the asymptotic solutions of the reduced problem on 0,1 and 1,2 respectively, which do not satisfy the conditions in (2). From the solution of the reduced problem, the value of ux=1=μ (say) can be obtained. In (14a)-(14b), v01 is the left boundary layer function, v02 is the left side of x=1 interior layer function, v03 is the right side of x=1 interior layer function, and v04 is the right boundary layer function.The functionsv01, v02, v03, and v04 satisfy, respectively, the following transformed homogeneous differential equations:(15a)−d2v01τ1dτ12+r0v01τ1=0;τ1∈0⋅∞,(15b)−d2v02τ2dτ22+r1v02τ2=0;τ2∈0⋅∞,(15c)−d2v03τ3dτ32+r1v03τ3=0;τ3∈0⋅∞,(15d)−d2v04τ4dτ42+r2v04τ4=0;τ4∈0⋅∞.Solving each subequations in (15a15b15c15d) at corresponding terminal points and by (14a)-(14b), we obtain the zeros order solution of (1)-(2) as(16a)ux=u01x+λ1e−r0/εx+λ2e−r1/ε1−x,x∈0,1,(16b)ux=u02x+λ3e−r1/εx−1+λ4e−r2/ε2−x,x∈1,2,so that applying the interval and boundary conditions, we can get λ1, λ2, λ3, and λ4 as(17a)λ1=ψ0−u010−μ−u011e−r1/ε1−e−r0/ε+r1/ε,(17b)λ2=μ−u011−ψ0−u010e−r0/ε1−e−r0/ε+r1/ε,(17c)λ3=μ−u021−γ−u022e−r2/ε1−e−r1/ε+r2/ε,(17d)λ4=γ−u022−μ−u021e−r1/ε1−e−r1/ε+r2/ε.


## 4.2. Derivation and Properties of the Numerical Scheme

We divide the interval0,2 into N equal parts with uniform mesh length h. Suppose 0=x0,x1,x2,…,xN/2=1,xN/2+1,xN/2+2,…,xN=2 be the mesh points. Then, we have xi=ih,i=0,1,2,…,N. Let us choose the terminal points as xn1=ε, xn2=1−ε, xn3=1+ε, and xn4=2−ε. The left boundary layer is in the interval 0,ε, the right boundary layer is in the interval 2−ε,2, the left interior layer is in the interval 1−ε,1, the right interior layer is in the interval 1,1+ε, the left outer region is in the interval ε,1−ε, and the right outer region is in the interval 1+ε,2−ε.Now, atx=xi, (1) can be written as −εd2u/dxi2+rxiuxi+sxiuxi−1=wxi. Approximating the differential operator by the central finite difference as d2u/dxi2=ui−1−2ui+ui+1/h2−h2/12ui4+R, where R=−h4/360ui6t, t∈xi−h,xi+h, and introducing a fitting factor σ to obtain an ε-independently convergent solution, we obtain(18)Luxi≔−εσh2uxi−1+2εσh2+rxiuxi−εσh2uxi+1+εσh212u4xi−εσR+sxiuxi−1=wxi,∀i=1,2,…,N−1.From this, we have(19a)L1ui≔Aiui−1+Biui+Aiui+1=Di+Ti,∀i=1,2,…,N2−1,(19b)L2ui≔Aiui−1+Biui+Aiui+1=Ei+Ti,∀i=N2+1,N2+1,…,N−1,where Ai=−εσ/h2, Bi=2εσ/h2+ri, Di=wi−siψi−N/2, Ei=wi−siui−N/2, and Ti=−εσh2/12u4t+εσOh4 for t∈xi−h,xi+h.Case 1.
Left boundary layer.
On the interval0,ε, we introduce a fitting factor σ1 in (19a) as(20)−εσ1ui−1−2ui+ui+1h2+riui=Di+Tifori=1,2,…,N/8.
To determine the fitting factorσ1 on the left boundary layer, we use the left boundary layer asymptotic solution with the left outer solution as(21)ui=u01xi+λ1e−r0/εxi.
We assume that the solution converges uniformly to the solution of (1)-(2). From (20), we have −σ1/ρ2ui−1−2ui+ui+1=Di+Ti−riui, where ρ=h/ε. Taking the limit as h⟶0, it becomes(22)−limh⟶0σ1ρ2ui−1−2ui+ui+1=limh⟶0Di+Ti−riui.
From (21), we have ui−1=u01xi−1+λ1e−r0/εxier0ρ and ui+1=u01xi+1+λ1e−r0/εxie−r0ρ. Inserting these in (22) and simplifying gives(23)σ1=r0ρ/2sinhr0ρ/22,which is the fitting factor in the interval 0,ε. Having this fitting factor in (20), we obtain(24)−εσ1h2ui−1+2εσ1h2+riui−εσ1h2ui+1=Di+Ti,fori=1,2,…N8.Case 2.
Left outer region.
On the intervalε,1−ε, we have the left outer region. Setting ε=0, (19a) becomes(25)riui=Di,fori=N8+1,N8+2,…,3N8,which is the left outer region scheme.Case 3.
Left-side interior layer.
On the interval1−ε,1, the interior layer will be on the left-hand side of x=1. Introducing a fitting factor σ2 in (19a), we have(26)−εσ2ui−1−2ui+ui+1h2+riui=Di+Ti,fori=3N8+1,3N8+2,…,N2.
To determine the fitting factorσ2 on the left interior layer, we use the corresponding interior layer asymptotic solution with the left outer region solution as(27)ui=u01xi+λ2e−r1/ε1−xi.
Assuming that the solution converges uniformly to the solution of (1)-(2), from (26), we have(28)−limh⟶0σ2ρ2ui−1−2ui+ui+1=limh⟶0Di+Ti−riui,where ρ=h/ε. From (27), we have ui−1=u01xi−1+λ2e−r1/ε1−xie−r1ρ and ui+1=u01xi+1+λ2e−r1/ε1−xier1ρ. Substituting these in (28) and simplification gives a fitting factor in left side of x=1 interior layer as(29)σ2=r0ρ/2sinhr1ρ/22.
With this fitting factor, we have(30)−εσ2h2ui−1+2εσ2h2+riui−εσ2h2ui+1=Di+Ti,fori=3N8+1,3N8+2,…N2.Case 4.
Right-side interior layer.
On the interval1,1+ε, the interior layer will be on the right side of x=1. If we introduce a fitting factor σ3 in (19b), then we obtain(31)−εσ3ui−1−2ui+ui+1h2+riuxi=Ei+Ti,fori=N2+1,N2+2,…,5N8.
To determineσ3, we use the right-side interior layer asymptotic solution with the right outer solution as(32)ui=u02xi+λ3e−r1/εxi−1.
Again, assuming that the solution converges uniformly to the solution of (1)-(2) from (31), we have(33)−limh⟶0σ3ρ2ui−1−2ui+ui+1=limh⟶0Ei+Ti−riui,where ρ=h/ε. Using (32), we have ui−1=u02xi−1+λ3e−r1/εxi−1er1ρ and ui+1=u02xi+1+λ3e−r1/εxi−1e−r1ρ. Putting these in (33) and simplifying gives the fitting factor σ3 in the right side of x=1 interior layer as(34)σ3=r0ρ/2sinhr1ρ/22.
Having this fitting factor, (31) becomes(35)−εσ3h2ui−1+2εσ3h2+riui−εσ3h2ui+1=Ei+Ti,fori=N2+1,N2+2,…,5N8.Case 5.
Right outer region.
On the interval1+ε,2−ε, we have the right outer region. Setting ε=0, (19b) is reduced to(36)riui=Ei,fori=5N8+1,5N8+2,…,7N8,which is the right outer region scheme.Case 6.
Right boundary layer.
On the interval2−ε,2, introducing a fitting factor σ4 in (19b), we have(37)−εσ4ui−1−2ui+ui+1h2+riui=Ei+Ti,fori=7N8+1,7N8+2,…,N.
To determineσ4, use the right boundary layer asymptotic solution with left outer solution as(38)ui=u02xi+λ4e−r2/ε2−xi.
Assuming that the solution converges uniformly to the solution of (1)-(2), we have(39)−limh⟶0σ4ρ2ui−1−2ui+ui+1=limh⟶0Ei+Ti−riui,where ρ=h/ε. Using (38), we have ui−1=u02xi−1+λ4e−r2/ε2−xier2ρ and ui+1=u02xi+1+λ4e−r2/ε2−xie−r2ρ, and by putting these into (39) and simplifying, we get the fitting factor:(40)σ4=r0ρ/2sinhr2ρ/22.
With this fitting factor, we have the right boundary layer scheme as(41)−εσ4h2ui−1+2εσ4h2+riui−εσ4h2ui+1=Ei+Ti,fori=7N8+1,7N8+2,…,N.
The solution of the considered singularly perturbed differential equation with large negative shift can be obtained by solving the three term recurrence relations (24), (30), (35), and (42) together with the outer layer schemes (25) and (26).


## 4.3. Discrete Stability Analysis

To establish the computational stability of the proposed numerical scheme, we follow the approaches of [27, 28]. From (19a), consider the recurrence relation:(42)Aiui−1+Biui+Aiui+1=Di,i=1,2,…,N2,where Ai=−εσ/h2 and Bi=2εσ/h2+ri with σ as in (23) or (29) and subject to the boundary conditions u0=ψ0 and uN/2=uxN/2. Now, we set(43)ui=Piui+1+Qi,fori=N−1,N−2,…,2,1,where Pi=Pxi and Qi=Qxi are determined as follows. From (43), we have ui−1=Pi−1ui+Qi−1, and substituting in (42), we obtain(44)ui=−AiAiPi−1+Biui+1+Di−AiQi−1AiPi−1+Bi.From (43) and (44), we can obtain(45a)Pi=−AiAiPi−1+Bi,(45b)Qi=Di−AiQi−1AiPi−1+Bi.Using the initial conditionu0=ψ0=P0u1+Q0, we can solve (45a) and (45b). If we choose P0=0, it follows that Q0=u0, and using these and uN/2=uxN/2, we can compute Pi and Qi, for all i=1,2,…,N/2. Suppose that, in computing Pi from (45a), a small error ei has been occurred. Then, we have Pi∗=Pi+ei and Pi−1∗=Pi−1+ei−1. However, we are computing Pi∗=−Ai/AiPi−1∗+Bi. So, we have(46)ei=Pi∗−Pi=−AiAiPi−1∗+Bi+AiAiPi−1+Bi=AiPiAiPi−1+Aiei−1+Biei−1.From the assumption in (3), we have rx>0, and since Ai=−εσ/h2 and Bi=2εσ/h2+ri, it follows that Bi>2Ai, for i=1,2,…,N−1. Then, from (45a), we have P1=−A1/A1P0+B1<1 and P2=−A2/A2P1+B2<−A2/A2+B2<1 as P1<1 and successively; if we continue in this manner, we get Pi<1, for i=0,1,2,…,N−1. Thus, under the assumption that there is small initial error, we have ei=AiPi/AiPi−1+Aiei−1+Biei−1, which implies that ei<ei−1. This indicates that (45a) is stable. In a similar procedure, we can show that (45b) is also stable, and following the same procedure for the numerical scheme from (19b), we can obtain a similar result so that the stability estimate of the proposed scheme is established.


## 4.4. Convergence Analysis

To provide the convergence analysis of the proposed numerical scheme, we follow the approaches in [29, 30]. From (19a), we have a system of equations as(47)−εσh2ui−1+2εσh2+riui−εσh2ui+1+fi+Ti=0,∀i=1,2,…,N2,where σ is as in (23) or (29), fi=−wi+siψi−N/2, and Ti=εσh2/12ui4t,t∈xi−h,xi+h. Using the boundary conditions u0=ψ0 and u1=uxN/2 in (47), we can get a system of equations in the matrix form as(48)G+Hu+M+Th=0,where G=εσ/h22−10…0−12−1…0⋮⋮⋮…⋮⋮⋮⋮…−100…−12 is a tri-diagonal matrix of order N/2−1 and H=r100…00r20…0⋮⋮⋮…⋮⋮⋮⋮…000…0r1−N/2 is a diagonal matrix of order N/2−1.And the associated vectors of (48) are M=f1−εσψ0/h2,f2,…,fN/2−2,fN/2−1−εσ/h2μt, and Th=T1,T2,…,TN/2−1t. Now, let us consider U as the approximation of u satisfying the system of equations.(49)G+HU+M=0.Suppose that theith truncation error be ei=Ui−ui, i=1,2,…,N/2−1, such that E=e1,e2,…,eN/2t=U−u. Taking the difference between (48) and (49) yields(50)G+HE=Th.Letri≤K,∀i=1,2,…,N/2−1, for arbitrary constant K, and let Pi,i be the i,ith entry of matrix P. Then, since ri>0 for a sufficiently small h, we have Pi,i≠0 and 2εσ/h2+Pi,i≠0, for i=1,2,…,N/2−1. These indicate that G+H is an irreducible matrix.In matrixG+H, if Si is the sum of entries on the ith row, then we have Si=εσ/h2+ri,i=1,N/2−1, and Si=ri,i=2,3,…,N/2−2. Let us define K1=min1≤i≤N/2−1ri and K2=max1≤i≤N/2−1ri, which imply that 0<K1≤K≤K2. Thus, for a sufficiently small value of h, G+H is monotone so that it is nonsingular and G+H−1≥0. From (50), we have(51)E≤G+H−1Th.In each row ofG+H, we can see that Si>h2εσ/h2+ri, for i=1,N/2−1, and Si>h2ri, for i=2,3,…,N/2−2. Suppose that G+Hi,j−1 be the i,jth entry of G+H−1, and let us define G+H−1=max1≤i≤N/2−1∑j=1N/2−1G+Hi,j−1 and Th=max1≤i≤N/2−1Ti. Since ∑j=1N/2−1G+Hi.j−1Sj=1,i=1,2,…,N/2−1, and G+H−1≥0, it follows that G+Hi,j−1≤1/Sj<1/h2εσ/h2+ri,j=1,N/2−1, and ∑j=2N/2−2G+Hi,j−1≤1/min2≤j≤N/2−2≤1/h2ri,i=1,2,…,N/2−2. From (51), using the ith row sum, we have(52)E≤G+H−1Th=max1≤i≤N/2−1∑j=1N/2−1G+Hi,j−1×max1≤i≤N/2−1Ti≤1εσ/h2+ri+1ri+1εσ/h2+ri×εσui4t12h2≤Ch2.Thus, for sufficiently chosen large value ofC, this result shows the second-order convergence of the proposed scheme on (0, 1], and by a similar analysis, the convergence of the numerical scheme can be provided on (1, 2).


## 5. Numerical Examples and Discussion

To show the validity and applicability of the proposed numerical scheme, we solve examples of singularly perturbed delay differential equations of the type in (1)-(2). For a problem whose exact solution is known or can be determined, we compute the maximum point-wise error as eεN=max0≤i≤Nuxi−UNxi, and if the exact solution of a problem is not known, we apply the double mesh principle [31] as eεN=max0≤i≤NUNxi−U2Nxi and maximum absolute error is eN=maxeεN, where uxi is the exact solution and UNxi and U2Nxi are numerical solutions computed at common mesh points. The point-wise rate of convergence is calculated as rεN=logeεN−logeε2N/log2 and the parameter uniform rate of convergence is rN=logeN−loge2N/log2.Example 1 (see [22]).
Consider a constant coefficient boundary value problem−εu″x+5ux−ux−1=1,ux=1,for−1≤x≤0,andu2=2.Example 2 (see [21]).
Consider a variable coefficient boundary value problem−εu″x+x+5ux−ux−1=1,ux=1,for−1≤x≤0,andu2=2.Example 3 (see [20]).
Consider a constant coefficient and zero source function boundary value problem−εu″x+2ux−ux−1=0,ux=1for−1≤x≤0,andu2=1.
Since the exact solutions are not given for the considered examples, we used the double mesh principle to show the maximum absolute error by the proposed method. In Tables1–3, the maximum point-wise errors and convergence rates of Examples 1–3 are given, respectively. From these tables, we observe that the maximum point-wise convergence is of second order and increases when mesh elements increase, which confirm the theoretical results. In Table 4, results obtained by the proposed method are compared with some published results in the literature. From this table, we can see that the numerical results are improved by the present method. The solutions of the three examples are given in Figures 1–3, respectively. In each figure, we observe that minimizing the perturbation parameter decreases the width of the boundary layers and the interior layer, which is the desired result. In Figure 4, we observe the Log-Log plot of the maximum absolute error verses the number of mesh points. In this figure, we observe that the rate of convergence of the scheme is two, which is in right agreement with the theoretical finding.Table 1
eεN, eN, rεN, and rN of Example 1, for different values of ε and N.
↓ε|N⟶292102112122132−001.8437e − 044.6688e − 051.1747e − 052.9463e − 067.3777e − 071.98151.99081.99531.99772−021.3535e − 043.4390e − 058.6676e − 062.1757e − 065.4704e − 071.97661.98831.99421.99182−044.3100e − 051.1071e − 052.8057e − 067.0623e − 071.8821e − 071.96091.98041.99021.90782−061.1424e − 052.9831e − 067.6189e − 071.9250e − 074.8824e − 081.93721.96921.98471.97922−084.8890e − 051.3402e − 053.5042e − 068.9565e − 072.2639e − 071.86711.93531.96811.98412−101.6110e − 044.8891e − 051.3402e − 053.5042e − 068.9566e − 071.72031.86711.93531.96812−124.2044e − 041.0110e − 044.8891e − 051.3402e − 053.5042e − 061.38391.72031.86711.93532−146.0346e − 044.2044e − 041.6110e − 044.8891e − 051.3402e − 050.52141.38391.72031.86712−165.9842e − 043.9522e − 041.7203e − 045.0582e − 051.3516e − 050.59851.20001.76601.90402−185.9840e − 044.0946e − 041.7380e − 045.0619e − 051.3520e − 050.54741.23631.77971.90462−205.9841e − 044.0943e − 041.7408e − 045.0620e − 051.3530e − 050.54751.23391.78201.90352−225.9841e − 044.0941e − 041.7420e − 045.0632e − 051.3584e − 050.54761.23281.78261.89812−245.9841e − 044.0940e − 041.7426e − 045.0644e − 051.3636e − 050.54761.23231.78281.89302−265.9841e − 044.0940e − 041.7428e − 045.0646e − 051.3638e − 050.54761.23211.78291.89282−285.9841e − 044.0940e − 041.7428e − 045.0648e − 051.3640e − 050.54761.23211.78281.89232−305.9841e − 044.0940e − 041.7428e − 045.0648e − 051.3640e − 050.54761.23211.78281.8923eN6.0346e − 044.2044e − 041.7428e − 045.0648e − 051.3640e − 05rN0.52141.27051.78281.8927Table 2
eεN, eN, rεN, and rN of Example 2, for different values of ε and N.
↓ε|N⟶292102112122132−001.8665e − 044.7283e − 051.1900e − 052.9948e − 067.4744e − 071.98091.99041.99531.99762−021.2319e − 043.1343e − 057.9051e − 061.9850e − 064.9735e − 071.97471.98731.99361.99682−042.7505e − 057.0856e − 061.7983e − 064.5299e − 071.1368e − 071.95671.97831.98911.99452−061.0979e − 052.9195e − 067.5245e − 071.9099e − 074.8108e − 081.91101.95611.97811.98912−084.1803e − 051.1706e − 053.1648e − 067.9473e − 072.0141e − 071.83641.88711.99361.98032−101.3501e − 044.2281e − 051.2033e − 053.0998e − 067.9521e − 071.67501.81301.95671.96282−123.3699e − 041.3657e − 044.3426e − 051.1798e − 053.1031e − 061.30311.65301.88001.92682−144.3315e − 041.3802e − 045.7230e − 052.1108e − 056.9812e − 061.65001.27001.43901.59622−164.3226e − 041.4013e − 045.6930e − 051.9876e − 056.8120e − 061.62511.29951.51821.54492−184.3218e − 041.4109e − 045.7004e − 052.0812e − 056.8120e − 061.61151.30751.45361.61132−204.3221e − 041.4113e − 045.6816e − 052.1200e − 056.8091e − 061.61471.31271.42221.63852−224.3221e − 041.4112e − 045.6819e − 052.1108e − 056.8112e − 061.61481.31251.42861.63182−244.3221e − 041.4112e − 045.6820e − 052.1112e − 056.8106e − 061.61481.31251.42831.63222−264.3221e − 041.4112e − 045.6820e − 052.1110e − 056.8109e − 061.61481.31251.42851.63202−284.3221e − 041.4112e − 045.6820e − 052.1110e − 056.8109e − 061.61481.31251.42851.63202−304.3221e − 041.4112e − 045.6820e − 052.1110e − 056.8109e − 061.61481.31251.42851.6320eN4.3315e − 041.4113e − 045.7230e − 052.1200e − 056.9812e − 05rN1.61781.30221.43271.6025Table 3
eεN, eN, rεN, and rN of Example 3, for different values of ε and N.
↓ε|N⟶292102112122132−008.5338e − 052.1595e − 055.4317e − 061.3621e − 063.4114e − 071.99851.99121.99561.99742−027.4882e − 051.8975e − 054.7761e − 061.1981e − 063.0003e − 071.98051.99021.99511.99762−045.0149e − 051.2755e − 053.2163e − 068.0754e − 072.0261e − 071.97521.98761.99381.99482−063.7180e − 059.4821e − 062.3941e − 066.0149e − 071.5058e − 071.97121.98571.99291.99802−081.1403e − 042.9523e − 057.5077e − 061.8928e − 064.7517e − 071.94951.97541.98781.99402−104.2339e − 041.1395e − 042.9504e − 057.5028e − 061.8915e − 061.89361.94941.97541.98792−121.4397e − 034.2339e − 041.1395e − 042.9504e − 057.5028e − 061.76571.89361.94941.97542−143.9298e − 031.4397e − 034.2339e − 041.1395e − 042.9504e − 051.44871.76571.89361.94942−165.9127e − 033.9298e − 031.4397e − 034.2339e − 041.1395e.040.58941.44861.76571.89362−185.6718e − 033.4217e − 031.2628e − 034.0845e − 041.2004e.040.72911.43811.62841.76662−205.6293e − 033.4020e − 031.2708e − 034.1240e − 041.1820e.040.72661.42061.62361.80282−225.9272e − 033.4112e − 031.2684e − 034.1826e − 041.1660e.040.72211.42731.60051.84282−245.6268e − 033.4120e − 031.2713e − 034.1724e − 041.1689e.040.72171.42431.60741.83572−265.6268e − 033.4118 − 031.2720 − 034.1720e − 041.1676e.040.72181.42341.60831.83722−285.6268 − 033.4118 − 031.2716 − 034.1722 − 041.1670.040.72181.42391.60781.83802−305.6268 − 033.4118 − 031.2717 − 034.1722 − 041.1668.040.72181.42381.60791.8383eN5.9127e − 033.9298e − 031.4397e − 034.2339e − 041.2004e − 04rN0.58941.44871.76571.8185Table 4
Comparison of the proposed scheme with some other published works in literature.
N⟶210211212213Example1forε=2−15Proposed schemeeN4.2044e − 041.6110e − 044.8891e − 051.3402e − 05rN1.38391.72031.8671Results in [22]eN4.7371e − 041.6440e − 045.5672e − 05rN1.52681.5622Example 2forε=2−15Proposed schemeeN2.2629e − 047.9546e − 052.2635e − 056.0860e − 06rN1.50831.81321.8950Results in [21]eN9.7137e − 044.0988e − 041.8682e − 048.8062e − 05rN1.04301.02171.01091.0056Example 3forε=2−12N⟶272829210Proposed schemeeN3.9298e − 031.4397e − 034.2339e − 041.1395e − 04rN1.44861.76571.89361.8936Results in [20]eN0.5770e − 020.3410e − 020.2010e − 020.1140e − 03rN0.61000.73200.79200.8170Figure 1
Numerical solutions of Example1 for N=210 and different values of ε. (a) ε = 2−4, (b) ε = 2−9, and (c) ε = 2−14.
(a)(b)(c)Figure 2
Numerical solutions of Example2 for N=210 and different values of ε. (a) ε = 2−4, (b) ε = 2−8, and (c) ε = 2−13.
(a)(b)(c)Figure 3
Numerical solutions of Example3 for N=210 and different values of ε. (a) ε = 2−5, (b) ε = 2−10, and (c) ε = 2−16.
(a)(b)(c)Figure 4
Log-Log scale plots of N vs. Max. pointwise errors for (a) Example1, (b) Example 2, and (c) Example 3.
(a)(b)(c)


## 6. Conclusion

In this study, a second-order singularly perturbed differential equation with large negative shift is considered. The differential equation is transformed to difference equations via decomposing the domain into its boundary, interior, and outer regions. In each boundary and interior layer subdomains, the problem is discretized and exponentially fitted numerical schemes are formulated. And in each outer layer regions, reduced problems are obtained by settingε=0. Combining these schemes, we obtained a second-order ε-uniformly convergent numerical scheme. To validate the theoretical results, three model examples are considered and solved. The numerical results of the examples confirm the theoretical analysis of the proposed numerical scheme, and hence, the proposed numerical scheme is convergent, independent of the perturbation parameter.



---

*Source: 7974134-2022-06-16.xml*