#### 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.

#### 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 . 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 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 Section 2, 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 used 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

The interval and boundary conditions are given aswhere and the functions , , and are sufficiently smooth on such that

Also, is a smooth function on and is a given constant, which is independent of . Problems (1)-(2) can be rewritten as

With on , , , and , where and represent the left and right side limit of at , respectively. Due to the influence of the perturbation parameter, the solution of the boundary value problems (1)-(2) exhibits strong boundary layers at and , and due to the large negative shift, strong interior layers occur on left and right sides of [24]. Moreover, such type of problem has a unique solution [25].

#### 3. Analytical Results

Lemma 1 (continuous maximum principle). *For the smooth functions and satisfying (3), let be in such that , , and on . Then, on .*

*Proof. *Let be given such that , and for the contrary, suppose . By the given conditions, . It follows from calculus that and . Consequently, we havewhich is a contradiction to the hypothesis. Therefore, it follows that and , for all .

Lemma 2 (stability result). *Let the conditions in (3) hold true and be any function in . Then, for all , we have .*

*Proof. *Let us define two barrier functions as . By these functions, we haveAnd by (4a4b), for all , we haveApplying Lemma 1, it follows that , for all , and hence, , for all .

Lemma 3 (bBound for the derivatives of the solution). *Let the conditions in (3) hold true and be the solution of (1)-(2). Then, we have*

*Proof. *The case for is the result of Lemma 2. To handle the case for , let and construct an associated neighborhood , such that and . By the mean value theorem, for some , we haveThen, for , we haveIn a similar procedure, for , we construct neighborhood . Then, by the mean value theorem as above, we obtain .

The case, for , follows rearranging terms in (1) asFrom this, we have . Differentiating both sides of (11) once and twice, we obtainFor 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 ) 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) aswhere , , , and . Then, the corresponding zero-order asymptotic expansions of (13a)-(13b) are given bywhere and are the asymptotic solutions of the reduced problem on and respectively, which do not satisfy the conditions in (2). From the solution of the reduced problem, the value of (say) can be obtained. In (14a)-(14b), is the left boundary layer function, is the left side of interior layer function, is the right side of interior layer function, and is the right boundary layer function.

The functions , , , and satisfy, respectively, the following transformed homogeneous differential equations:

Solving each subequations in (15a15b15c15d) at corresponding terminal points and by (14a)-(14b), we obtain the zeros order solution of (1)-(2) asso that applying the interval and boundary conditions, we can get , , , and as

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

We divide the interval into equal parts with uniform mesh length . Suppose be the mesh points. Then, we have . Let us choose the terminal points as , , , and . The left boundary layer is in the interval , the right boundary layer is in the interval , the left interior layer is in the interval , the right interior layer is in the interval , the left outer region is in the interval , and the right outer region is in the interval .

Now, at , (1) can be written as . Approximating the differential operator by the central finite difference as , where , , and introducing a fitting factor to obtain an -independently convergent solution, we obtain

From this, we havewhere , , , , and for .

*Case 1. *Left boundary layer.

On the interval , we introduce a fitting factor in (19a) asTo determine the fitting factor on the left boundary layer, we use the left boundary layer asymptotic solution with the left outer solution asWe assume that the solution converges uniformly to the solution of (1)-(2). From (20), we have , where . Taking the limit as , it becomesFrom (21), we have and . Inserting these in (22) and simplifying giveswhich is the fitting factor in the interval . Having this fitting factor in (20), we obtain

*Case 2. *Left outer region.

On the interval , we have the left outer region. Setting , (19a) becomeswhich is the left outer region scheme.

*Case 3. *Left-side interior layer.

On the interval , the interior layer will be on the left-hand side of . Introducing a fitting factor in (19a), we haveTo determine the fitting factor on the left interior layer, we use the corresponding interior layer asymptotic solution with the left outer region solution asAssuming that the solution converges uniformly to the solution of (1)-(2), from (26), we havewhere . From (27), we have and . Substituting these in (28) and simplification gives a fitting factor in left side of interior layer asWith this fitting factor, we have

*Case 4. *Right-side interior layer.

On the interval , the interior layer will be on the right side of . If we introduce a fitting factor in (19b), then we obtainTo determine , we use the right-side interior layer asymptotic solution with the right outer solution asAgain, assuming that the solution converges uniformly to the solution of (1)-(2) from (31), we havewhere . Using (32), we have and . Putting these in (33) and simplifying gives the fitting factor in the right side of interior layer asHaving this fitting factor, (31) becomes

*Case 5. *Right outer region.

On the interval , we have the right outer region. Setting , (19b) is reduced towhich is the right outer region scheme.

*Case 6. *Right boundary layer.

On the interval , introducing a fitting factor in (19b), we haveTo determine , use the right boundary layer asymptotic solution with left outer solution asAssuming that the solution converges uniformly to the solution of (1)-(2), we havewhere . Using (38), we have and , and by putting these into (39) and simplifying, we get the fitting factor:With this fitting factor, we have the right boundary layer scheme asThe 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:where and with as in (23) or (29) and subject to the boundary conditions and . Now, we setwhere and are determined as follows. From (43), we have , and substituting in (42), we obtain

From (43) and (44), we can obtain

Using the initial condition , we can solve (45a) and (45b). If we choose , it follows that , and using these and , we can compute and , for all . Suppose that, in computing from (45a), a small error has been occurred. Then, we have and . However, we are computing . So, we have

From the assumption in (3), we have , and since and , it follows that , for . Then, from (45a), we have and as and successively; if we continue in this manner, we get , for . Thus, under the assumption that there is small initial error, we have , which implies that . 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 aswhere is as in (23) or (29), , and . Using the boundary conditions and in (47), we can get a system of equations in the matrix form aswhere is a tri-diagonal matrix of order and is a diagonal matrix of order .

And the associated vectors of (48) are , and . Now, let us consider as the approximation of satisfying the system of equations.

Suppose that the truncation error be , , such that . Taking the difference between (48) and (49) yields

Let , for arbitrary constant , and let be the entry of matrix . Then, since for a sufficiently small , we have and , for . These indicate that is an irreducible matrix.

In matrix , if is the sum of entries on the row, then we have , and . Let us define and , which imply that . Thus, for a sufficiently small value of , is monotone so that it is nonsingular and . From (50), we have

In each row of , we can see that , for , and , for . Suppose that be the entry of , and let us define and . Since , and , it follows that , and . From (51), using the row sum, we have

Thus, for sufficiently chosen large value of , 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 , and if the exact solution of a problem is not known, we apply the double mesh principle [31] as and maximum absolute error is , where is the exact solution and and are numerical solutions computed at common mesh points. The point-wise rate of convergence is calculated as and the parameter uniform rate of convergence is .

*Example 1 (see [22]). *Consider a constant coefficient boundary value problem .

*Example 2 (see [21]). *Consider a variable coefficient boundary value problem .

*Example 3 (see [20]). *Consider a constant coefficient and zero source function boundary value problem .

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 Tables 1–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.

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(a)**

**(b)**

**(c)**

**(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 . 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.

#### Data Availability

No data were used to support the findings of the study.

#### Conflicts of Interest

The authors declare there are no potential conflicts of interest.