#### Abstract

In this article, singularly perturbed parabolic differential difference equations are considered. The solution of the equations exhibits a boundary layer on the right side of the spatial domain. The terms containing the advance and delay parameters are approximated using Taylor series approximation. The resulting singularly perturbed parabolic PDEs are solved using the Crank–Nicolson method in the temporal discretization and nonstandard finite difference method in the spatial discretization. The existence of a unique discrete solution is guaranteed using the discrete maximum principle. The uniform stability of the scheme is investigated using solution bound. The uniform convergence of the scheme is discussed and proved. The scheme converges uniformly with the order of convergence , where is number of subintervals in spatial discretization and is mesh length in temporal discretization. Two test numerical examples are considered to validate the theoretical findings of the scheme.

#### 1. Introduction

Differential equations play a prominent role in many disciplines for modelling real-world problems in which the evolution of the system depends on the present state of the system. Problems that depend on the past and current evolution of the system are modelled by the delay deferential equation [1].

Singularly perturbed differential equations (SPDEs) are differential equations in which their highest order derivative term is multiplied by a small perturbation parameter. Such types of differential equations model various processes in applied mathematics, for instance, in fluid mechanics, fluid dynamics, elasticity, quantum mechanics, chemical-reactor theory, reaction-diffusion processes, and other domains of fluid motion [2, 3]. Singularly perturbed differential difference equations model physical problems for which the evaluation does not only depend on the present state of the system but also on the past history. In general, when the perturbation parameter tends to zero, the smoothness of the solution of the singularly perturbed differential difference equations deteriorates and it forms a boundary layer [4]. Such types of singularly perturbed differential difference equations have an application in the modelling of neuronal variability [5, 6].

Musila and Lansky [7] give a generalization for the model in [5, 6] in the form of singularly perturbed parabolic differential difference equations (SPPDDEs):for studying time evolution of the trajectories of the membrane potential, where the first derivative term denotes the exponential decay between two consecutive jumps caused by the input processes. The membrane potential decays exponentially to the resting level with a membrane time constant . and are diffusion moments of the Wiener process characterizing the influence of dendritic synapses on the cell excitability. The excitatory input contributes to the membrane potential by an amplitude with intensity , and similarly, the inhibitory input contributes by an amplitude with intensity .

The mathematical model in (1) is in the form of SPPDDEs. There are no known analytical methods of solution for physiologically acceptable parameters [7]. So, to simulate the solution of this model, one has to develop suitable numerical methods. Numerical methods developed for regular problems become inapplicable for solving SPPDDEs when the perturbation parameter tends to zero [8]. Due to this, there is an interest in the development of numerical methods that converges uniformly (or converges independent of the perturbation parameter).

In recent years, scholars have devoted for the numerical treatment of SPPDDEs. In a series of papers [8–11], Bansal and Sharma developed different numerical schemes for solving SPPDDEs with large shift arguments. Kumar and Kadalbajoo [12] applied the implicit Euler method in temporal discretization and B-spline collocation method on the Shishkin mesh in the spatial discretization, and they obtain a first-order uniformly convergent result. Ramesh and Kadalbajoo [4] used the implicit Euler method in temporal discretization and upwind and midpoint upwind FDM on Shishkin mesh for the spatial discretization. They obtain almost first-order uniformly convergent numerical result. Rao and Chakravarthy [13] applied the implicit Euler method in temporal discretization and exponentially fitted operator FDM in the spatial discretization. They obtain a first-order uniformly convergent result. Woldaregay and Duressa [14–16] developed numerical schemes using fitted operator methods on uniform mesh. They obtained first-order uniformly convergent result. Kumar et al. [17] developed stable FDM using the implicit Euler method for temporal discretization and a class of nonstandard FDM using infinite serious approximation for the derivative terms in spatial discretization, and they apply Richardson extrapolation technique to extend the rate of convergence of the scheme to order two. Shivhare et al. [18] developed a numerical scheme using the implicit Euler method in temporal discretization and quadratic B-spline collocation method for the spatial discretization on the exponentially graded mesh. Daba and Duressa [19] used the implicit Euler method in temporal discretization and extended cubic B-spline collocation method for the spatial discretization on the Shishkin mesh.

In this article, we developed numerical scheme using the Crank–Nicolson method in the temporal discretization and nonstandard FDM in the spatial discretization. We used Rothe’s procedure for discretizing the problem (i.e., first we discretize the temporal direction and followed by spatial discretization). The proposed scheme satisfies the discrete maximum principle and uniform stability estimate. The developed scheme treats the considered problem accurately with parameter uniform convergence.

##### 1.1. Organization of the Paper

The rest of the article is organized as follows. In Section 2, the statement of the problem is given. In Section 3, the asymptotic approximation and its properties are discussed. In Section 4, the proposed numerical scheme and its uniform convergence analysis are given. In Section 5, numerical examples, results, and discussions are given. Finally, in Section 6, the conclusion of the article is included.

##### 1.2. Notation

In this article, and are denoted for the number of mesh points in spatial and temporal discretization, respectively. The symbol (in some cases indexed) denotes a positive constant independent of the perturbation parameter and . The norm denotes the supremum (or maximum) norm.

#### 2. Statement of the Problem: Aim of the Research

On the domain , for some fixed number , with the boundary , a singularly perturbed parabolic differential difference equation is given bywhere is a singular perturbation parameter and and are small delay and advance parameters satisfying , with the initial and interval conditions

The problem in (2) and (4) is the simplified form of the model in (1) by including the source function with specified initial and interval conditions. From the convection diffusion fluid flow point of view, the terms in (2) and (4) can be classified as is called the diffusion term; is called the convection term and the nonderivative terms referred to as the reaction terms.

For the existence of a unique solution, the functions , , and are assumed to be sufficiently smooth, bounded, and independent of . The coefficient functions , and are assumed to satisfyfor some constant . For small values of , the solution of (2) and (4) exhibits left or right boundary layer depending on the sign of the coefficient function of the convection term . If , a regular boundary layer occurs in the neighbourhood of , and for the case , the layer occurs in the neighbourhood of [20].

The main contribution of this article is, to develop uniformly convergent numerical scheme using nonstandard FDM for solving the problem in (2) and (4). The developed method is applicable on the uniform mesh.

#### 3. Properties of the Analytical Solution

In this section, we discuss behaviours of the analytical solution of the problem in (2) and (4). When the shift parameters and are smaller than , using Taylor series approximation for the terms with shifts is valid [21]. So, we approximate and as

Using the approximations in (6) into (2), we obtainwith initial and boundary conditionswhere , , and .

For small values of the shift parameters , (2) and (5) and (7) and (8) are asymptotically equivalent, since the difference between the two equations is . We assumed that , where and are the lower bounds for and , respectively. Furthermore, we assumed , which implies the occurrence of the boundary layer in the neighbourhood of . The other case implies the occurrence of the boundary layer in the neighbourhood of and can be treated similarly. The layer is maintained for sufficiently small shift parameters .

Lemma 1. *To avoid the presence of a singularity in the numerical solution, it is required that the given data satisfy the compatibility conditions**, and**so that the data match at the two corners**and*. Let the functions , and be continuous on *, then* (7) and (8) *have a unique solution**. In particular, when the compatibility conditions are not satisfied, a unique classical solution still exist but may not be differentiable on all of* [22]*.*

In this article, we considered the case, when boundary layer occurs in the neighbourhood of . Using the compatibility conditions, we deduce that there exists a constant independent of such that for all .

For details, interested readers can refer [22].

*Remark 1. **Note that there does not exist a constant**independent of**such that**, since the layer occurs in the neighbourhood of**.*

The problem obtained by setting in (7) and (8) is called reduced problem and given aswith the initial conditions

It is first-order hyperbolic PDEs with initial data given along sides and of the domain . For small values of , the solution of the problem in (7) and (8) is very close to the solution of (11) and (12).

Lemma 2. *The maximum principle. Let be any sufficiently smooth function defined on , which satisfies . Then, implies that .*

*Proof. *Let be such that and suppose that . It is clear that implying that . Since from the theory of extrema of a function in calculus implies, , and giving that , which is contradiction to the assumption. Therefore, . Hence, we obtain .

Lemma 3. *Uniform stability estimate. Let be the solution of the problem in (7) and (8). Then, it satisfies the following bound:where .*

*Proof. *The proof is by the construction of barrier functions as , and applying the maximum principle, we obtain the required bound. For the detail of the proof, the interested reader can refer [16].

Lemma 4. *Let be solution of (7) and (8). Then, the derivatives of the solutions with respect to and satisfy following the bound:where is lower bound of .*

*Proof. *Refer in [23].

#### 4. Formulation of the Numerical Scheme

##### 4.1. Temporal Discretization

We discretize the temporal domain into subintervals using the grid points , where . Let is denoted for the approximation of at the time level discretization. The continuous problem in (7) and (8) semidiscretized using the Crank–Nicolson method in the temporal direction resulting to a system of boundary value problemswith the following discretized boundary conditions:where, and , for and .

Lemma 5. *Semidiscrete maximum principle. Let be a sufficiently smooth function on the domain . If , and , then, .*

*Proof. *Let be such that and suppose that . From the assumption, it is clear that , which implies that . Now, using the differential operator at the point , we haveApplying the property in calculus for minimum points, we obtainUsing the estimates in (18) into (17), we obtain , which is contradiction to the assumption . Therefore, we conclude that .

Next, let us analyse the bound of the error in temporal semidiscretization. Let us denote and be the exact and approximate solution of the problem in (7) and (8) using the above discretization. Let us denote the local error for each time step by .

Lemma 6. *Local error estimate. Suppose that . The local error estimate in the temporal discretization is given by*

*Proof. *Using Taylor’s series expansion for and centered at , we obtainSubstituting (20) into (7), we obtainwhere,

Now, it can be seen that the local error is the solution ofHence, applying the maximum principle on the operator givesNext, we need to show the bound for the global error of the temporal semidiscretization. Let us denote for the global error estimate up to the time step.

Lemma 7. *Global error estimate . The global error estimate at is given by*

*Proof. *Using the local error estimate up to the time step given in the above lemma, we obtain the global error estimate at the time step aswhere is constant independent of and .

Next, we set a bound for the derivatives of solution of the problem in (15).

Lemma 8. *For each, , the solution of the boundary value problems in (15) satisfies the following bound:*

*Proof. *See in [23].

##### 4.2. Spatial Discretization

For the spatial discretization, we apply the nonstandard FDM. The construction of nonstandard FDM depends on the knowledge of the corresponding exact FDM, which is given in [24].

###### 4.2.1. Exact Finite Difference Method

For the problem of the form in (15), in order to construct exact finite difference scheme, we follow the procedures developed in [24] for BVPs. Consider the constant coefficient subequations from (15) aswhere and . Since equation (27) is constant coefficient second-order differential equations, it has two independent solutions, that is,where .

We discretize the spatial domain , using uniform mesh length such that , where is the number of mesh intervals in spatial discretization. We denote the approximate solution of at mesh point by . The target is to calculate a difference equation, which has the same general solution as (27) at the mesh point given by

Using the theory of difference equations for second-order linear difference equations in [24], we obtainand then, by simplifying and substituting the values of and , we obtainwhich is an exact difference scheme for (44). For , we use the approximation in (32). Hence, multiplying both sides by , simplifying we obtain

Rearranging, we obtain

The denominator function for the second-order derivative approximation in (34) is . For the variable coefficient problems, we modify the denominator function aswhere is a function of , and .

##### 4.3. The Discrete Scheme

Using the denominator function in (35) into the discretization of (15), we obtain the nonstandard finite difference schemes aswhere and

##### 4.4. Stability and Uniform Convergence Analysis

In this subsection, we want to show that the discrete scheme in (24) satisfies the discrete maximum principle, uniform stability estimates. and uniform convergence.

Lemma 9. *Discrete maximum principle . Let be any mesh function satisfying . Then, , which implies that .*

*Proof. *Suppose that there exists such that . From the assumption, it implies . Also, we assume that and . Now, we haveand using the assumptions made above, we obtain , for . Thus, the supposition , for is wrong. Hence, we obtain.

Next, we want to prove that the scheme in (36) satisfies the uniform stability estimate.

Lemma 10. *Uniform stability estimate . The solution of the discrete scheme in (24) satisfies the following bound:where is the lower bound of .*

*Proof. *Let and define the barrier function by . At the boundary points, we obtainOn the discretized spatial domain , we obtainUsing the maximum principle in Lemma 9, we obtain . Hence, the required bound is obtained.

Lemma 11. *Let be the denominator function defined in (35), and then, the following bound holds*

*Proof. *Let us define . Using the expression for , we obtainNext, let us set a bound for . SinceWe obtain the bounds as . Hence, from (43) and (44), we obtainLet us define the forward and backward finite differences operators in spatial direction asrespectively, and the second-order finite difference operator as

Theorem 1. *The discrete solution of the problem in (9) satisfies the truncation error bound:*

*Proof. *We consider the truncation errorUsing the estimate in Lemma 11 for the bound of , the truncation error bound becomesUsing the bound of the derivatives of the solution in Lemma 8 in (50) givesThis completes the proof.

Lemma 12. *For a fixed mesh number and for , it holdswhere .*

Theorem 2. *The discrete solution of the problem in (15) satisfies the error bound:where .*

*Proof. *By using the results of Lemma 12 in Theorem 1 and applying the bound in Lemma 10 gives the required bound.

Theorem 3. *Under the hypothesis of boundedness of discrete solution, the solution of (7) and (8) satisfies the parameter uniform bound:*

*Proof. *The parameter uniform convergence of the scheme follows from the bound in (15) and (31):

#### 5. Numerical Results and Discussion

To validate the established theoretical results, we develop an algorithm and perform experiments using the proposed numerical scheme on the problem of the form given in (2) and (4).

*Example 1. *Consider the problemwith the initial condition and the interval conditions , and .

*Example 2. *Consider the problemwith the initial condition and the interval conditions , and .

The exact solution of the considered examples is not available. So, the maximum absolute error, -uniform error and rate of convergence are calculated using the double mesh principle. Let be the computed solution of the problem using numbers of mesh points and be the computed solution on the double number of mesh points by including the midpoints and into the mesh points. The maximum absolute error is calculated as

The uniform error estimate is calculated using

The rate of convergence of the scheme is calculated usingand the uniform rate of convergence is calculated using

The solution of the problems in example 1 and example 2 exhibits a regular boundary layer on the right side of the spatial domain. In Figures 1 and 2, one can observe boundary layer formation as the perturbation parameter goes small, particularly at , one can observe in Figures 3(d) and 4(d) . In Figures 3(a)–3(d) and 4(a)–4(d), the effect of the shift parameters and the perturbation parameter is shown on the profile of the solution by considering different values for the arameters.

**(a)**

**(b)**

**(a)**

**(b)**

**(a)**

**(b)**

**(c)**

**(d)**

**(a)**

**(b)**

**(c)**

**(d)**

The numerical result in Tables 1 and 2 show the maximum absolute error and rate of convergence of the scheme. The maximum absolute errors presented in Tables 1 and 3 show the uniform convergence of the scheme. One can observe from these tables, as the value of the perturbation parameter goes small, the maximum absolute error becomes stable and constant (i.e., converges are independent of the influence of the perturbation parameters ). Tables 3 and 4 show the maximum absolute error for different values of delay and advance parameters. As it is observed in these tables, the shift parameters do not influence the uniform convergence of the scheme. Table 5 shows the maximum absolute error and rate of convergence of the scheme for temporal discretization. As one observes in this table, the temporal discretization has the second order of convergence with correct agreement to the theoretical analysis. Tables 6 and 7 present the uniform error and uniform rate of convergence of the scheme and results from published papers in literature. The proposed scheme is more accurate than the upwind scheme in the Shishkin mesh, midpoint upwind scheme in the Shishkin mesh in [4], cubic B-spline scheme in the Shishkin mesh in [12], and the exponentially fitted scheme in [13].

Figure 5 shows the log-log plot of the maximum absolute error verses number of mesh points in space of the scheme for different values of the perturbation parameter. As one observes from the figures, for all values of the perturbation parameter, the line graphs are overlapped and parallel to the reference line . This indicates that the convergence of the scheme is independent of the perturbation parameter .

**(a)**

**(b)**

#### 6. Conclusion

In this article, a parameter uniform numerical method is developed for solving a singularly perturbed parabolic differential equation involving both delay and advance on the spatial variable of reaction terms. The considered problems’ solution exhibits a boundary layer on the right side of the spatial domain. The developed scheme constitutes the Crank–Nicolson method in temporal discretization and nonstandard FDM in the spatial discretization. The uniform stability and parameter uniform convergence of the scheme are investigated and proved. The applicability of the scheme is validated by taking numerical test examples. The effect of the perturbation parameter, effect of the delay, and advance parameters on the solution are depicted using figures. The proposed scheme converges uniformly (i.e., converges independent of the perturbation parameter) with the order of convergence one in space and two in time directions. The performance of the scheme is investigated by comparing the results with some schemes in published papers. We found that the proposed scheme is more accurate, stable, and parameter uniformly convergent.

#### Data Availability

No external data are used in this article.

#### Conflicts of Interest

The authors declare there are no potential conflicts of interest.