Delay differential equations are equations that are dependent on the previous states and have been used in various dynamical systems. For instance, in robotics, delays occur through the manipulation of information or feedback control [1]. A surface acoustic wave sensor is modeled using a delay differential equation [2]. In chemical kinetics, the reaction time and the time taken for mixing the reactants are modeled by delay differential equations [3]. Apart from these, delayed dynamical systems have also been found useful in modeling musical instruments [4], traffic dynamics [5], models of HIV infection [6], population dynamics [7], economic cycles [8], and others.

A subclass of differential equations in which the term with the highest order derivative is multiplied by a small positive parameter (ε) and involves one or more shift arguments is said to be a singularly perturbed differential equation with delay [9]. Such problems frequently arise in the modeling of various physical systems, such as the human pupil-light reflex [10], the study of bistable devices in digital electronics [11], variational problem in control theory [12], immune response modeling [13], mathematical modeling in ecology [14], models to stabilize rotating and frozen waves [15], models for the physiological process [16]. The presence of the small positive number causes boundary layers and the spatial delay gives rise to an interior layer in the solution of the problem. The layers are asymptotically narrow regions in the neighborhood of the end points of the domain, where the gradient of the solution decreases as ε approaches zero [17]. The rapidly changing behavior of the solution in the layers causes a significant error in the solution, and hence, it is almost impossible to solve the problem analytically. On the contrary, standard numerical methods do not give satisfactory solutions, unless a large number of mesh points are considered, which requires a massive computational cost [18]. To overcome such drawbacks, there is a need of developing a numerical scheme, independent of the perturbation parameter.

Various research works are available in the literature to address the aforementioned limitations. For instance, Erdogan and Cen [19] solved a singularly perturbed convection–diffusion with delay by constructing a hybrid finite difference scheme on a Shishkin-type mesh and obtained a uniformly convergent method with respect to ε. In [20], a class of turning point singularly perturbed boundary value problems is treated by constructing a numerical scheme based on the trigonometric quintic B-spline basis functions on a piece-wise uniform mesh. The method is obtained to be almost first-order convergent irrespective of ε. In [21], a time-dependent singularly perturbed differential equation with delay is solved by constructing a numerical scheme using the Euler scheme on uniform time mesh and a hybrid finite difference scheme on piece-wise uniform Shishkin mesh in the spatial direction. In [22], two-dimensional singularly perturbed semi-linear convection–diffusion problems have been treated using the nonstandard finite difference approach. The authors linearized the continuous problem and then discretized it using the nonstandard finite difference methods. Mbroh and Munyakazi [23] solved singularly perturbed one- and two-dimensional problems by constructing a scheme using the method of lines by using the fitted operator finite difference method for the space discretization and the Crank–Nicolson method for the time discretization. In [24], a singularly perturbed problem with time lag is treated by constructing a numerical method using the standard finite difference operators centered in space and implicit in time on a piece-wise uniform mesh. Sahoo and Gupta [25] solved a singularly perturbed problem involving discontinuous convective and source terms by developing a numerical scheme using a first-order accurate, simple upwind scheme on specially designed piece-wise uniform Shishkin meshes. Appadu and Tijani [26] treated a one-dimensional generalized Burgers–Huxley equation by proposing two solutions using the classical finite difference scheme and nonstandard finite difference scheme and obtained that one of the proposed solutions contains a minor error. In [27], a singularly perturbed ordinary differential equation with a large negative shift is treated by developing a numerical scheme using the fitted operator method via domain decomposition.

Singularly perturbed differential equations involving large delays in the spatial variable have been solved by few authors. In [28], a singularly perturbed problem with a large delay is solved by developing a scheme using the Crank–Nicolson method on a temporal mesh and the central difference method on nonuniform Shishkin meshes. Bansal and Sharma [29] formulated a scheme for a singularly perturbed with delay using the implicit Euler on the temporal meshes and standard central difference on nonuniform spatial meshes. Ejere et al. [30] developed a numerical scheme for a time-dependent singularly perturbed differential equation with large spatial delays using the weighted-average method in the temporal direction and the central difference method in the spatial direction on piece-wise uniform Shishkin meshes. Alam and Khan [31] proposed a new numerical algorithm for singularly perturbed differential equations involving the shift and the advance parameters. They used Crank–Nicolson in the time direction and cubic B-spline basis functions on generalized Shishkin mesh in the spatial direction, and by this, they obtained a uniformly convergent scheme of order four in time and almost of order four in space.

In this article, we proposed a robust numerical scheme to solve singularly perturbed differential equations involving spatial and temporal delays. The scheme is developed using the implicit Euler method for the temporal variable and the nonstandard finite difference method for the spatial variable on uniform meshes. To handle the temporal delay, we used the Taylor series approximation, and the spatial delay is handled by choosing special meshes, in such a way that the term with the spatial delay coincides with the mesh point x_i = 1. Error estimate and uniform convergence analysis are investigated and proved for the proposed method. Model numerical examples are also solved to support the theoretical results.

The remaining sections of the article are organized as follows: The description of the continuous problem is presented in Section 2. The time semi-discrete scheme and the fully discrete scheme are briefly discussed in Section 3. To support the validity of the proposed scheme, numerical examples, results, and discussions are provided in Section 4, and the study is concluded in Section 5.

Notations: Throughout this article, we used C as a generic positive number, independent of ε and the mesh numbers. If w is the smooth function in D̄, then we used the maximum norm as ||w||=max(s,t)∈D̄|w(s,t)|.

We considered a singularly perturbed delay differential equation given by

Where D̄=[0,2]×[0,T], 0 < ε≪1, a temporal shift τ > 0 of o(ε) and a finite time T. We assumed that the functions α(s), β(s), γ(s, t), w₀(s, t), ψ(s, t), and φ(t) are smooth enough and bounded on the considered domain. Moreover, to avoid oscillation in the solution for arbitrary positive constant λ, the coefficient functions α(s) and β(s) satisfy the conditions [29]

Setting ε = 0, the associated reduced problem is given by

Since w₀ need not satisfy the conditions w₀(0, t − τ) = ψ(0, t − τ), w₀(0, t − τ) = ψ(0, 0), w₀(2, t − τ) = φ(2, t), w0(1-,t)=w0(1+,t), and (w0)s(1-,t)=(w0)s(1+,t), the solution exhibits two boundary layers and an interior layer at s = 1 [28, 32]. The functions involved in the continuous problem (1) should satisfy Holder continuity and compatibility conditions are imposed at the corner points as

and

From the aforementioned assumptions and conditions, we observe that the solution involves layers and there exists a constant C independent of ε [33]. So for (s,t)∈D̄, we have |w(s, t) − w₀(s, 0)| = |w(s, t) − w₀(s)| ≤ Ct. Applying Taylor's series expansion, we have

Inserting Equations (5) into (1) gives

subjected to

where

Lemma 1. Suppose z(s, t) is a smooth function in D̄. If z(s, t) ≥ 0, (s,t)∈∂D and Lεz(s,t)≥0, (s,t)∈D, then Lεz(s,t)≥0, (s,t)∈D̄.

Proof. For (ŝ,t^)∈D̄, suppose that z(ŝ,t^)=minD̄z(s,t)<0. From the considered hypothesis (ŝ,t^)∉∂D. By extreme value theorem, we have zt(ŝ,t^)=0, zs(ŝ,t^)=0, and zss(ŝ,t^)>0. Then,

Case 1: For s ∈ (0, 1], Lεz(ŝ,t^)=∂z(ŝ,t^)∂t-ε∂2z(ŝ,t^)∂s2+α(ŝ,t^)<0.

Case 2: For s ∈ (1, 2), Lεz(ŝ,t^)=∂z(ŝ,t^)∂t-ε∂2z(ŝ,t^)∂s2+α(ŝ)+β(ŝ)z(ŝ-1,t^)-τβ(ŝ)∂z(ŝ-1,t^)∂t≤∂z(ŝ,t^)∂t-ε∂2z(ŝ,t^)∂s2+[α(ŝ)+β(ŝ)]z(ŝ,t^)-τβ(ŝ)∂z(ŝ,t^)∂t<0.

Thus, Lεz(ŝ,t^)<0, which contradicts the given condition. Therefore, our supposition fails, so that z(ŝ,t^)≥0, which implies that z(s, t) ≥ 0, (s,t)∈D̄.                □

Lemma 2. The solution of the continuous problem (6)-(7) is estimated as follows:

Proof. Let us define z±(s,t)=||ϑ||λ+max{|∂D|}±w(s,t). Then, we have z^±(0, t) ≥ 0 and z^±(0, t) ≥ 0. Moreover, For s ∈ (0, 1], t ∈ [0, T],

For s ∈ (1, 2), t ∈ [0, T],

Thus, Lεz±(s,t)≥0 and by Lemma 1, we have z^±(s, t) ≥ 0, (s,t)∈D̄, which yields the stability estimate.                □

Lemma 3. The derivatives of the solution w(x, t) of Equations (6), (7) are bounded as follows:

for all nonnegative integer k, l such that 0 ≤ k + 2l ≤ 4.

Proof. For k = 0 and l = 0, it is to show the bound of the solution w(s, t), which is Lemma 2. For k = 0, l ≠ 0, we show the bound of derivatives of w(s, t) with respect to t. For k ≠ 0, l = 0, it is to show the bound of derivatives of the solution w(s, t) with respect to s and for the cases when (k, l) ≠ 0, we determine the bound of its derivatives of the solution w(s, t) with respect to s and t. Let qε,1(s)=e-sλ/ε+e-(1-s)λ/ε, s ∈ (0, 1] and qε,2(s)=e-(s-1)λ/ε+e-(2-s)λ/ε, s ∈ (1, 2). For a fixed value of λ, following the approaches of Kellogg and Tsan [34], we can get that |q_{ε, ι}(s)| ≤ c, for constant c and ι = 1, 2. For the case k = 0, Equation (6) becomes

Along the sides t = 0, we have w = 0, which implies that ∂2w(s,t)∂s2=0. Then from Equation (8), we get

Without loss of generality, assuming that all the values considered in Equation (7) are zero, we have w(s − 1, 0) = ψ(s − 1, 0) = 0, s ∈ [0, 1] and w(s − 1, 0) = w₀(s − 1, 0) = 0, s ∈ (1, 2]. Then, Equation (9) becomes ∂w(s,0)∂t=γ(s,0). Since γ is a smooth function, it implies |∂w(s,t)∂t|≤C, for sufficiently chosen C on ∂D. Applying the operator Lε given in Equation (6) on ∂w(s,t)∂t, we obtain |Lεwt(s,t)|=|ϑt(s,t)|≤C on D. Thus, applying Lemma 1 gives |∂w(s,t)∂t|≤C on D̄. For l = 2, differentiating Equation (8) with respect to t gives

Along s = 0 and s = 2, we have ∂2w(s,t)∂t2=0, and along t=0, we have w = 0 and ∂2w(s,0)∂s2=0. From ∂w(s,0)∂t=γ(s,0), we have ∂3w(s,0)∂s2∂t=∂2γ(s,0)∂s2. Assuming that the initial and boundary conditions are identically zero, we have w_t(s − 1, 0) = ψ_t(s − 1, 0) = 0, s ∈ (0, 1] and w_t(s − 1, 0) = (_w0)t(s − 1, 0) = 0, s ∈ (1, 2]. Then from Equation (10), we get

From Equation (11), along t = 0, we obtain that |∂2w(s,0)∂t2|≤0 on ∂D and the operator Lε implies |Lε∂2w(s,t)∂t2|=|∂2ϑ(s,t)∂t2|≤C on D and applying Lemma 1 yields |∂2w(s,t)∂t2|≤C on D̄. The bound of derivatives of the solution w(s, t) for the cases k ≠ 0 can be determined by similar procedures as mentioned earlier. For the case k = 1, let s∈D and consider a neighborhood R=(e,e+ε), ∀s ∈ R. For some q∈R̄ and t ∈ (0, T], the mean value theorem gives

For s∈R̄, we can get

Putting Equation (12) into the aforementioned result, we obtain |∂w(s,t)∂s|≤C(1+ε-12qε,ι), (s,t)∈D̄ for sufficiently large chosen C. For k = 2, by rearranging Equation (1) we get

For a fixed time t ∈ [0, T] and for s ∈ [0, 2], since |w| ≤ C, |∂w(s,t)∂t|≤C, and γ is a smooth function, |∂2w(s,t)∂s2|≤C(1+ε-1qε,ι). The derivative of Equation (13) with respect to t gives

Since |∂w(s,t)∂t|≤C, |∂2w(s,t)∂t2|≤C, and |γ_t(s, t)| ≤ C, we get |∂3w(s,t)∂s2∂t|≤C(1+ε-1qε,ι). In a similar procedure, the bounds on the derivatives of the solution can be easily determined for the remaining values of k and l with 0 ≤ k + 2l ≤ 4.                □