Skip to main content

ORIGINAL RESEARCH article

Front. Appl. Math. Stat., 24 March 2023
Sec. Numerical Analysis and Scientific Computation
This article is part of the Research Topic Singularly Perturbed Problems: Asymptotic Analysis and Numerical Solutions View all 6 articles

A robust numerical scheme for singularly perturbed differential equations with spatio-temporal delays

\nAbabi Hailu Ejere
Ababi Hailu Ejere1*Gemechis File Duressa&#x;Gemechis File Duressa2Mesfin Mekuria Woldaregay&#x;Mesfin Mekuria Woldaregay1Tekle Gemechu DinkaTekle Gemechu Dinka1
  • 1Department of Applied Mathematics, Adama Science and Technology University, Adama, Ethiopia
  • 2Department of Mathematics, Jimma University, Jimma, Ethiopia

In this article, we proposed and analyzed a numerical scheme for singularly perturbed differential equations with both spatial and temporal delays. The presence of the perturbation parameter exhibits strong boundary layers, and the large negative shift gives rise to a strong interior layer in the solution. The abruptly changing behaviors of the solution in the layers make it difficult to solve the problem analytically. Standard numerical methods do not give satisfactory results, unless a large mesh number is considered, which needs a massive computational cost. We treated such problem by proposing a numerical scheme using the implicit Euler method in the temporal variable and the nonstandard finite difference method in the spatial variable on uniform meshes. The stability and uniform convergence of the proposed scheme have been investigated and proved. To demonstrate the theoretical results, numerical experiments are carried out. From the theoretical and numerical results, we observed that the method is uniformly convergent of order one in time and of order two in space.

1. Introduction

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 xi = 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)|.

2. Description of the continuous problem

We considered a singularly perturbed delay differential equation given by

{Lεw(s,t)=(tε2s2+α(s))w(s,t)         =γ(s,t)β(s)w(s1,tτ), (s,t)D,w(s,t)=w0(s,t), (s,t)D0={(s,t):s[0,2],t[τ,0]},w(s,t)=ψ(s,t), (s,t)DL={(s,t):s[1,0], t[0,T]},w(2,t)=φ(t), (2,t)DR={(2,t), t[0,T]},    (1)

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), w0(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]

α(s)+β(s)2λ>0 and β(s)<0,s[0,2].    (2)

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

L0w(s,t)={(t+α(s))w0(s,t)=γ(s,t)β(s)ψ0(s1,tτ),   (s,t)DL,(t+α(s))w0(s,t)=γ(s,t)β(s)w0(s1,tτ),   (s,t)DR.

Since w0 need not satisfy the conditions w0(0, t − τ) = ψ(0, t − τ), w0(0, t − τ) = ψ(0, 0), w0(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

{w0(0,0)=ψ(0,0),w0(2,0)=φ(0)    (3)

and

{ψ(0,0)tε2w0(0,0)s2+α(0)w0(0,0)+β(0)w0(1,τ)=γ(0,0),φ(2,0)tε2w0(2,0)s2+α(2)w0(2,0)+β(2)w0(1,τ)=γ(2,0).    (4)

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) − w0(s, 0)| = |w(s, t) − w0(s)| ≤ Ct. Applying Taylor's series expansion, we have

w(s-1,t-τ)=w(s-1,t)-τw(s-1,t)t+O(τ2)    (5)

Inserting Equations (5) into (1) gives

Lεw(s,t)=ϑ(s,t)(s,t)D̄    (6)

subjected to

{w(s,t)=w0(s,t), (s,t)D¯w(s,t)=ψ(s,t), (s,t)DL={(s,t):s[1,0], t[0,T]}w(2,t)=φ(t), (2,t)DR={(2,t), t[0,T]},    (7)

where

Lεw(s,t)=w(s,t)tε2w(s,t)s2+α(s)w(s,t)ϑ(s,t)={γ(s,t)β(s)ψ(s1,t)+τβ(s)ψ(s1, t)t, s(0,1]γ(s,t)β(s)w(s1,t)+τβ(s)w(s1, t)t, s(1,2).

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^)tz(ŝ,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:

|w(s,t)|||ϑ||λ+max{|D|}

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],

Lεz±(s,t)=(t-εs2+α(s))z±(s,t)                  =α(s)(||ϑ||λ+max{|D|})±ϑ(x,t)0.

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

Lεz±(s,t)=(t-εs2+α(s))z±(s,t)                  +β(s)z±(s-1,t)-τβ(s)z±(s-1,t)t              [α(s)+β(s)]max{|D|}0.

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:

|k+lw(s,t)sktl|  {C(1+εk/2[esλ/ε+e(1s)λ/ε]), (s,t)DL, C(1+εk/2[e(s1)λ/ε+e(2s)λ/ε]),   (s,t)DR

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

w(s,t)t=ε2w(s,t)s2α(s)w(s,t)+{γ(s,t)β(s)ψ(s1,t)+τβ(s)ψ(s1, t)t, s(0,1]γ(s,t)β(s)w(s1,t)+τβ(s)w(s1, t)t, s(1,2).    (8)

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

w(s,0)t=ε2w(s,0)s2α(s)w(s,0)+{γ(s,0)β(s)ψ(s1,0)+τβ(s)ψ(s1, 0)t, s(0,1]γ(s,0)β(s)w(s1,0)+τβ(s)w(s1, 0)t, s(1,2).    (9)

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) = w0(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

2w(s,t)t2=ε3w(s,t)s2tα(s)w(s,t)t+{γ(s,t)tβ(s)ψ(s1,t)t+τβ(s)2ψ(s1, t)t2, s(0,1]γ(s,t)tβ(s)w(s1,t)t+τβ(s)2w(s1, t)t2, s(1,2).    (10)

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)s2t=2γ(s,0)s2. Assuming that the initial and boundary conditions are identically zero, we have wt(s − 1, 0) = ψt(s − 1, 0) = 0, s ∈ (0, 1] and wt(s − 1, 0) = (w0)t(s − 1, 0) = 0, s ∈ (1, 2]. Then from Equation (10), we get

2w(s,0)t2=ε2γ(s,0)s2α(s)γ(s,0)t+{γ(s,0)t, s(0,1]γ(s,0)t+τβ(s)2w(s1, 0)t2, s(1,2).    (11)

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 sD and consider a neighborhood R=(e,e+ε), ∀sR. For some qR̄ and t ∈ (0, T], the mean value theorem gives

|w(q,t)s|=ε-12|w(e+ε,t)-w(e,t)|2ε-12||w||.    (12)

For sR̄, we can get

w(s,t)s=w(q,t)s+w(s,t)s-w(q,t)s            =w(q,t)s+qsw(x,t)xdx=  w(q,t)s               +ε-1qs(w(x,t)x+α(x)w(x,t)               +β(x)w(x-1,t)-γ(x,t))dx.         |w(s,t)s||w(q,t)s|+Cε-1qs(||w||+||γ||)dx       =|w(q,t)s|+Cε-1(||w||+||γ||)ε12.

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

2w(s,t)s2=ε-1[w(s,t)t+α(s)w(s,t)+β(s)w(s-1,t)-γ(s,t)].    (13)

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

3w(s,t)s2t=ε-1(2w(s,t)t2+α(s)w(s,t)t+β(s)w(s-1,t)t-γ(s,t)t).

Since |w(s,t)t|C, |2w(s,t)t2|C, and |γt(s, t)| ≤ C, we get |3w(s,t)s2t|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.

3. Numerical scheme

3.1. Semi-discrete scheme in time direction

Let m be a uniform mesh number on [0, T] with step size Δt. Then, the uniform temporal discretization is given as D̄tm={tj=jΔt,Δt=T/m,j=0(1)m}. Using the implicit Euler method, we obtain a semi-discrete scheme as

LεmWj+1(s)=ϑj+1(s),    (14)

where

LεmWj+1(s)={εd2Wj+1(s)ds2+(1Δt+α(s))Wj+1(s), 0<s1,εd2Wj+1(s)ds2+(1Δt+α(s))Wj+1(s)+(1τΔt)β(s)Wj+1(s1), 1<s<2

and

ϑj+1(s)={γj+1(s)+Wj(s)Δt(1τΔt)β(s)ψj+1(s1)τβ(s)Δtψj(s1), s(0,1],γj+1(s)+Wj(s)Δtτβ(s)ΔtWj(s1), s(1,2)

with the initial and boundary conditions

{W0(s)=w0(s), s[0,2],Wj+1(s)=ψj+1(s), s(0,1],Wj+1(2)=φj+1(2), s(1,2).    (15)

Lemma 4. For a smooth function zj+1(s), suppose that zj+1(s) ≥ 0, sD and Lεmzj+1(s)0, sD. Then, zj+1(s) ≥ 0, sD̄.

Proof. Let s* ∈ [0, 2] and suppose that zj+1(s*)=minD̄zj+1(s)<0. By the considered hypothesis, s*D and by the extreme value theorem, we have dzj+1(s*)ds=0 and d2zj+1(s*)ds20.

Case 1: For s* ∈ (0, 1], Lεmzj+1(s*)=-εd2zj+1(s*)ds2+(1Δt+α(s*))zj+1(s*)<0.

Case 2: For s* ∈ (1, 2), Lεmzj+1(s*)=(-εd2ds2+1Δt+α(s*))zj+1(s*)+(1-τΔt)β(s*)zj+1(s*-1)-εd2zj+1(s*)ds2+(1Δt+α(s*)+(1-τΔt)β(s*))zj+1(s*)<0. The two cases contradict the given condition. Therefore, our assumption is wrong, which implies that zj+1(s) ≥ 0, sD̄. As a result for the operator Lεm, we have

||(Lεm)-1||(1+λΔt)-1,    (16)

which is used in estimating the truncation error.

Lemma 5. (Semi-discrete stability estimate) The solution Wj+1(s) of Equations (14), (15) can be estimated as |Wj+1(s)|||LεmW||1+λΔt+max{|Dm|}, s ∈ [0, 2].

Proof. Define two barrier functions as Z±j+1(s)=||LεmW||1+λΔt+max{|Dm|}±Wj+1(s), from which we can obtain that Z±j+1(0)0 and Z±j+1(2)0.

Case 1: For 0 < s ≤ 1, we have

LεmZ±j+1(s)=-εd2Z±j+1ds2+(1Δt+α(s))Z±j+1(s)                 =±ϑj+1(s)+(1Δt+α(s))[||LεmW||1+λΔt+max{|Dm|}]                 (1Δt+α(s))max{|Dm|}0.

Case 2: For 1 < s < 2, we have

LεmZ±j+1(s)=-εd2Z±j+1ds2+(1Δt+α(s))Z±j+1(s)                       +(1-τΔt)β(s)Z±j+1(s-1)                   -εd2Z±j+1ds2+(1Δt+α(s))Z±j+1(s)                       +(1-τΔt)β(s)Z±j+1(s)                   =±ϑj+1(s)+(1Δt+α(s)+[1-τΔt]β(s))                       [||LεmW||1+λΔt+max{|Dm|}]                   (1Δt+α(s)+[1-τΔt]β(s))max{|Dm|}0.

Thus, applying Lemma 4 yields Z±j+1(s)0, s ∈ [0, 2], which implies the required estimation.

The local error committed in the semi-discrete scheme is the difference between the exact solution w(s, tj+1) and the approximate solution Wj+1(s) of Equation (14). That is, ej+1(s)=w(s,tj+1)-Wj+1(s) and the global error Ej+1 is the contribution of the local error up to the (j + 1)th time level. The bound of error for the semi-discrete scheme is estimated as follows.

Lemma 6. Suppose that |kw(s,t)tk|C, (s,t)D̄, k = 0, 1, 2. The local error is estimated as ||ej+1||C(Δt)2, and with this condition, the global error is estimated as ||Ej+1|| ≤ Ct), j + 1 = 1, 2, …, m.

Proof. From the Taylor series expansion, we have wj+1(s)=wj(s)+Δtdw(s,tj)dt+O((Δt)2) and from this we obtain

wj+1(s)-wj(s)Δt=wt(s,tj)t+O((Δt)2).    (17)

Using Equations (17) into (6) gives

Lεmwj+1(s)+O((Δt)2)=ϑj+1(s)    (18)

From Equations (14), (18), we can obtain the boundary value problem of the form

Lεmej+1(s)=O((Δt)2),ej+1(0)=0=ej+1(2)    (19)

Using Equations (16) into (19) yields ||ej+1|| ≤ Ct)2. Using the estimation of ej+1, we have

||Ej+1||=||ξ=1jeξ||,j(Δt)T                 =||e1+e2+...+ej||                 ||e1||+||e2||+...+||ej+1||                 CT(Δt)C(Δt),j+1=1,2,...,m.

Thus, the semi-discrete scheme is convergent of order one in time.

Lemma 7. Let the solution of Equation (14) be Wj+1(s). Then, its derivatives can be bounded as follows:

|dkWj+1dsk|{C[1+εk/2(exp(λ/εs)+exp(λ/ε(1s)))], s[0,1],C[1+εk/2(exp(λ/ε(s1))+exp(λ/ε(2s)))], s(1,2],

where k = 0, 1, 2, 3, 4.

Proof. The proof can be calculated by applying the procedures in the proof of Lemma 3 for the spatial domain. Furthermore, we refer to Clavero and Gracia [35].

3.2. Fully-discrete scheme

Let us sub-divide the domain [0, 2] into n uniform meshes of size h, such that DsN={0=s0,s1,...,sn/2=1,sn/2+1,...,sn=2,si=s0+ih,i=0(1)n,h=2/n}. According to the procedures in [36], we consider a constant coefficient sub-equation from Equation (14) as

-εd2Wj+1(s)ds2+λWj+1(s)=0,    (20)

Where 1Δt+α(s)λ>0. The characteristic of Equation (20) has two distinct roots r1=eλεs and r2=e-λεs. Denoting Wij+1 as the approximation of Wj+1(s) at the grid point si, we have Wij+1=c1eλεsi+c2e-λεsi and

|Wi-1j+1eλεsi-1e-λεsi-1Wij+1eλεsie-λεsiWi+1j+1eλεsi+1e-λεsi+1|=0.

After certain manipulation, we get

Wi-1j+1-2cosh(λεh)Wij+1+Wi+1j+1=0.    (21)

The finite difference approximation of Equation (20) is given as

-εWi-1j+1-2Wij+1+Wi+1j+1σi2+λWij+1=0,    (22)

where σi is a denominator function. From Equatiosn (21), (22), the denominator function for the variable coefficients is given as

σi=2ωisinh(ωih2), where ωi=1+ΔtαiεΔt    (23)

Using the denominator function (Equation 23), we can obtain a fully discrete scheme as

Lεn,mWij+1=ϑij+1,i=1,2,...,n-1,    (24)

where

LεWij+1={εδ2Wij+1+(1Δt+αi)Wij+1, i=1(1)n/2,εδ2Wij+1+(1Δt+αi)Wij+1+(1τΔt)βiWin/2j+1, i=n/2+1(1)n1

and

ϑij+1={γij+1+WijΔt(1τΔt)βiψin/2j+1τΔtβiψin/2j,i=1(1)n/2,γij+1+WijΔtτΔtβiWin/2j, i=n/2+1(1)n1,

With δ2Wij+1=Wi+1j+1-2Wij+1+Wi-1j+1σi2, and the discrete initial and boundary conditions W0(si)=w0(si),Wj+1(si)=ψj+1(si), Wj+1(2)=φ(2,tj+1),si[0,2],tj+1[0,T].

Lemma 8. Let Zij+1 be a given mesh function satisfying Z0j+10 and Znj+10. If Lεn,mZij+10 for i = 1(1)n − 1, then we have Zij+10 for i = 0(1)n.

Proof. For some y ∈ {1, …, n − 1}, suppose that Zyj+1=mini=1,...,n-1Zij+1<0.

Case 1: For y=1,2,...,n2, we have Lε,1nZyj+1=-εδ2Zyj+1+(1Δt+αy)Zyj+1< 0.

Case 2: For y=n2+1,n2+2,...,n-1, we have Lε,2nZyj+1=-εδ2Zyj+1+(1Δt+αy)Zyj+1+(1-τΔt)βyZy-n/2j+1-εδ2Zyj+1+(1Δt+αy)Zyj+1+(1-τΔt)βyZy< 0. From the two cases, we see that the given hypothesis is contradicted, and hence our assumption fails. Thus, Zij+10, i = 0(1)n.

Lemma 9. The solution Wij+1 of the fully discrete scheme (24) is estimated as

|Wij+1|Δt||ϑ||1+λΔt+max{|Dn,m|},i=0(1)n.

Proof. Let us define πi,±j+1=Δt||ϑ||1+λΔt+max{|Dn,m|}±Wij+1. Then, we have

π0,±j+1=Δt||ϑ||1+λΔt+max{|Dn,m|}±W0j+1Δt||ϑ||1+λΔt0πn,±j+1=Δt||ϑ||1+λΔt+max{|Dn,m||}±Wnj+1Δt||ϑ||1+λΔt0

When i=1(1)1n, we obtain

Lε,1n,mπi,±j+1=-εδ2πi,±j+1+(1Δt+αi)πi,±j+1                 =(1Δt+αi)(Δt||ϑ||1+λΔt+max{|Dn,m|})±ϑij+1                    (1Δt+αi)max{|Dn,m|}0

When i=n2+1(1)n-1, we obtain

Lε,2n,mπi,±j+1=-εδ2πi,±j+1+(1Δt+αi)πi,±j+1+(1-τΔt)βiϕi-n/2,±j+1               =(1Δt+αi+(1-τΔt)βi)                    (Δt||ϑ||1+λΔt+max{|Dn,m|})±ϑij+1               (1+Δtαi+(1-τΔt)βi)max{|Dn,m|}0

Thus, applying Lemma 8, we have πi,±j+10, j = 0(1)n, which yields the stability estimate.

Lemma 10. For a given fixed mesh number n, we have

limε0{maxsi(0,1]exp(α(xi)εsi)+exp(α(si)ε(1si))εk/2=0,maxsi(1,2)exp(α(si)ε(si1))+exp(α(si)ε(2si))εk/2=0

for all i = 1(1)n − 1 and k = 1, 2, 3, ….

Proof. We refer to Lemma 3.3 of Woldaregay and Duressa [37].

Theorem 1. Let the solution of Equation (14) be W(s, tj+1) and the solution of (24) be Wij+1. Then, we have

maxi=0(1)n,j=0(1)m|Wij+1-W(si,tj+1)|Cn-2.

Proof. From the differential and difference equations, the truncation error is given by

|Lεn,m(Wj+1(si)-Wij+1)|=|-εd2Wj+1(si)ds2+εδ2Wij+1|=|-εd2Wj+1(si)ds2+εσi2(Wi+1j+1-2Wij+1+Wi-1j+1)|    (25)

By the Taylor series expansion of Wi+1j+1, Wi-1j+1, and 1σi truncated up to order five, we have

Wi+1j+1=Wij+1+hdWj+1ds+h22!d2Wj+1ds2+h33!d3Wj+1ds3+h44!d4Wj+1(ζi)ds4,ζi(si-1,si+1),
Wi-1j+1=Wij+1-hdWj+1ds+h22!d2Wj+1ds2-h33!d3Wj+1ds3+h44!d4Wj+1(ζi)ds4,ζi(si-1,si+1),1σi2=ωi24sinh2(ωih/2)=ωi24(4(ωih)2-13+(ωih)260).

Putting these expansions in Equation (25) yields

|Lεn,m(Wj+1(si)-Wij+1)|=|-εd2Wij+1ds2+εd2Wij+1ds2+(ε12d4Wj+1(ζi)ds4-εωi212d2Wij+1ds2)h2+(εωi4240d2Wij+1ds2-εωi2144d4Wj+1(ζi)ds4)h4+εωi42880d4Wj+1(ζi)ds4h6|.

Applying the bound of derivatives in Lemma 7 and then using Lemma 10 yields

|Lεn,m(Wj+1(si)-Wij+1)|C1h2+C2h4+C3h6Ch2.    (26)

Invoking Lemma 8, we have maxi=0(1)n,j=0(1)m|Wij+1-W(si,tj+1)|Cn-2, because h = 2n−1.

Theorem 2. For the solutions w(s, t) of Equation (6) and Wij+1 of Equation (24), the uniform error is estimated as

maxi=0(1)n,j=0(1)m|w(si,tj+1)-Wij+1|C(Δt+n-2).

Proof. By the triangular inequality, we can obtain that

|w(si,tj+1)-Wij+1|=|w(si,tj+1)-Wj+1(si)+Wj+1(si)-Wij+1||w(si,tj+1)-Wj+1(si)|+|Wj+1(si)-Wij+1|.

Then, the combination of Lemma 6 and Theorem 1 yields the required uniform error estimate.

4. Numerical examples, results, and discussions

To demonstrate the validity and applicability of the proposed numerical scheme, we solved examples of the problem under consideration. Since the exact solution for each example is not given, we use the variant of the double mesh principle [38] to determine the maximum absolute error as

eεn,Δt=max0(1)n,0(1)m|Wn,Δt(si,tj)-W2n,Δt/4(si,tj)|,

Where W2n,Δt/4(si,tj) is the approximate solution obtained by taking (2n, Δt/4) for fixed value of the transition parameter. The uniform absolute error is determined by en,Δt=maxε(eεn,Δt). The convergence rate of the method is computed by ρεn,Δt=log(eεn,Δt/eε2n,Δt/4)log2 and uniformly it is obtained as ρn,Δt=maxε(ρεn,Δt).

Example 1. Consider wt-ε2ws2+6w(s,t)-2w(s-1,t-τ)=1, subjected to w0(s) = 0, w(s, t) = 0, and w(2, t) = 0.

Example 2. Consider wt-ε2ws2+(s+4)w(s,t)-(s2+1)w(s-1,t-τ)=2, subjected to w0(s) = 0, w(s, t) = 0, and w(2, t) = 0.

Example 3. Consider wt-ε2ws2+5w(s,t)-2w(s-1,t-τ)=2, subjected to w0(s) = sin(πs), w(s, t) = 0, and w(2, t) = 0.

We treated each problem by applying the proposed numerical method with the help of MATLAB R2019a packages. Since the exact solutions of the examples are not given, we used a variant of the double mesh principle to determine the numerical results. Moreover, the obtained results are displayed in tables and graphs. The uniform convergence is shown by computing the maximum point-wise error and convergence rate as given in Tables 13 for each example, respectively. From these tables, we observe that for a fixed value of ε, increasing the mesh numbers minimizes the maximum absolute error. On the contrary, for a fixed number of meshes, decreasing ε yields stable point-wise errors after certain changes in the values of ε. This shows the ε-uniform convergence of the proposed scheme.

TABLE 1
www.frontiersin.org

Table 1. Maximum absolute errors and uniform convergence rates of Example 1 at τ = ε/4.

TABLE 2
www.frontiersin.org

Table 2. Maximum absolute errors and uniform convergence rates of Example 2 at τ = ε/4.

TABLE 3
www.frontiersin.org

Table 3. Maximum absolute errors and uniform convergence rates of Example 3 at τ = ε/4.

In Figure 1, we observe the effect of the perturbation parameter in layer resolving for each example. Surface plots are simulated in Figures 24 for Examples 1–3, respectively. From each figure, we observe the effect of ε, that is, decreasing the values of ε decreases the width of the layers. The robustness of the developed scheme is illustrated by plotting log–log figures as given in Figure 5 for the considered examples.

FIGURE 1
www.frontiersin.org

Figure 1. Effect of the perturbation parameter in resolving the boundary layers for (A) Example 1, (B) Example 2, and (C) Example 3 when τ = ε/4.

FIGURE 2
www.frontiersin.org

Figure 2. Surface plot simulations of Example 1, for various values of ε at τ = ε/4. (A) ε = 20 and (B) ε = 216.

FIGURE 3
www.frontiersin.org

Figure 3. Surface plot simulations of Example 2, for various values of ε at τ = ε/4. (A) ε = 20 and (B) ε = 216.

FIGURE 4
www.frontiersin.org

Figure 4. Surface plot simulations of Example 3, for various values of ε at τ = ε/4. (A) ε = 20 and (B) ε = 216.

FIGURE 5
www.frontiersin.org

Figure 5. Log–log plots of N vs maximum absolute errors for (A) Example 1, (B) Example 2, and (C) Example 3.

5. Conclusion

In this study, we proposed and analyzed a fitted numerical method for a singularly perturbed differential equation involving spatial and temporal delays in the reaction term. The solution varies abruptly in the layers due to the presence of the perturbation parameter. The rapidly changing behavior of the layers and the effect of the delays cause difficulties to find the analytical solution. To solve the problem, we proposed a fitted numerical method. The method is obtained by using the implicit Euler method in the temporal variable and the nonstandard finite difference method in the spatial variable on uniform meshes. The effect of the temporal delay is handled by applying Taylor's series approximation and the spatial delay is handled by choosing a special mesh so that the delay term lies on the mesh point xi = 1. We investigated and proved that the proposed method is stable and uniformly convergent. We considered and solved three model examples to test the validity and applicability of the proposed method. The solutions and accuracy of the results of the examples are shown in graphical and tabular forms. From the theoretical and numerical findings discussed in the article, we can conclude that the proposed numerical method is uniformly convergent of order one in time and of order two in space.

Data availability statement

The original contributions presented in the study are included in the article/supplementary material, further inquiries can be directed to the corresponding author.

Author contributions

GD initiated the plan of this study. AE proposed the numerical scheme and analysis of the results. MW and TD revised the procedures, analysis, and results of the study. All authors have equal contributions to the article and agreed on the submitted version.

Acknowledgments

The authors kindly thank the editor for handling the manuscript carefully and referees for their valuable contributions in modifying the original version of the manuscript.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

References

1. Stépán G, Haller G. Quasiperiodic oscillations in robot dynamics. Nonlinear Dyn. (1995) 8:513–28. doi: 10.1007/BF00045711

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Feng Z, Chicone C. A delay differential equation model for surface acoustic wave sensors. Sensors Actuators A Phys. (2003) 104:171–8. doi: 10.1016/S0924-4247(03)00052-9

CrossRef Full Text | Google Scholar

3. Epstein IR. Delay effects and differential delay equations in chemical kinetics. Int Reviews in Phy Chem. (1992) 11:135–60. doi: 10.1080/01442359209353268

CrossRef Full Text | Google Scholar

4. Barjau A, Gibiat V. Delayed models for simplified musical instruments. J Acoust Soc Am. (2003) 114:496–504. doi: 10.1121/1.1577558

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Nagatani T, Nakanishi K. Delay effect on phase transitions in traffic dynamics. Phys Rev E. (1998) 57:6415. doi: 10.1103/PhysRevE.57.6415

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Culshaw RV, Ruan S. A delay-differential equation model of HIV infection of CD4+ T-cells. Math Biosci. (2000) 165:27–39. doi: 10.1016/S0025-5564(00)00006-7

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Gopalsamy K. Stability and Oscillations in Delay Differential Equations of Population Dynamics. Vol. 74. Springer Science & Business Media (2013). doi: 10.1007/978-94-015-7920-9

CrossRef Full Text | Google Scholar

8. Szydłowski M, Krawiec A. The Kaldor-Kalecki model of business cycle as a two-dimensional dynamical system. J Nonlinear Math Phy. (2001) 8:266–71. doi: 10.2991/jnmp.2001.8.s.46

CrossRef Full Text | Google Scholar

9. Mahendran R, Subburayan V. Fitted finite difference method for third order singularly perturbed delay differential equations of convection diffusion type. Int J Comput Methods. (2019) 16:1840007. doi: 10.1142/S0219876218400078

CrossRef Full Text | Google Scholar

10. Longtin A, Milton JG. Complex oscillations in the human pupil light reflex with mixed and delayed feedback. Math Biosci. (1988) 90:183–99. doi: 10.1016/0025-5564(88)90064-8

CrossRef Full Text | Google Scholar

11. Yang Y, Ouyang J, Ma L, Tseng RH, Chu CW. Electrical switching and bistability in organic/polymeric thin films and memory devices. Adv Funct Mater. (2006) 16:1001–14. doi: 10.1002/adfm.200500429

CrossRef Full Text | Google Scholar

12. Saksena VR, O'reilly J, Kokotovic PV. Singular perturbations and time-scale methods in control theory: survey 1976-1983. Automatica. (1984) 20:273–93. doi: 10.1016/0005-1098(84)90044-X

CrossRef Full Text | Google Scholar

13. Bocharov G, Romanyukha A. Numerical treatment of the parameter identification problem for delay-differential systems arising in immune response modelling. Appl Numer Math. (1994) 15:307–26. doi: 10.1016/0168-9274(94)00007-7

CrossRef Full Text | Google Scholar

14. Kot M. Elements of Mathematical Ecology. Cambridge: Cambridge University Press (2001).

Google Scholar

15. Schneider IAN. Spatio-Temporal Feedback Control of Partial Differential Equations. Berlin: Fereie University (2016). doi: 10.17169/refubium-5712

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Batzel JJ, Kappel F. Time delay in physiological systems: analyzing and modeling its impact. Math Biosci. (2011) 234:61–74. doi: 10.1016/j.mbs.2011.08.006

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Sharma KK, Rai P, Patidar KC. A review on singularly perturbed differential equations with turning points and interior layers. Appl Math Comput. (2013) 219:10575–609. doi: 10.1016/j.amc.2013.04.049

CrossRef Full Text | Google Scholar

18. Clavero C, Jorge JC. Another uniform convergence analysis technique of some numerical methods for parabolic singularly perturbed problems. Comput Math Appl. (2015) 70:222–35. doi: 10.1016/j.camwa.2015.04.006

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Erdogan F, Cen Z. A uniformly almost second order convergent numerical method for singularly perturbed delay differential equations. J Comput Appl Math. (2018) 333:382–94. doi: 10.1016/j.cam.2017.11.017

CrossRef Full Text | Google Scholar

20. Alam MP, Kumar D, Khan A. Trigonometric quintic B-spline collocation method for singularly perturbed turning point boundary value problems. Int J Comput Math. (2021) 98:1029–48. doi: 10.1080/00207160.2020.1802016

CrossRef Full Text | Google Scholar

21. Gupta V, Kumar M, Kumar S. Higher order numerical approximation for time dependent singularly perturbed differential-difference convection-diffusion equations. Numer Methods Partial Diff Equat. (2018) 34:357–80. doi: 10.1002/num.22203

CrossRef Full Text | Google Scholar

22. Kehinde OO, Munyakazi JB, Appadu AR. A NSFD discretization of two-dimensional singularly perturbed semilinear convection-diffusion problems. Front Appl Math Stat. (2022) 8:861276. doi: 10.3389/fams.2022.861276

CrossRef Full Text | Google Scholar

23. Mbroh NA, Munyakazi JB. A fitted operator finite difference method of lines for singularly perturbed parabolic convection-diffusion problems. Math Comput Simulat. (2019) 165:156–71. doi: 10.1016/j.matcom.2019.03.007

CrossRef Full Text | Google Scholar

24. Ansari A, Bakr S, Shishkin G. A parameter-robust finite difference method for singularly perturbed delay parabolic partial differential equations. J Comput Appl Math. (2007) 205:552–66. doi: 10.1016/j.cam.2006.05.032

CrossRef Full Text | Google Scholar

25. Sahoo SK, Gupta V. Higher order robust numerical computation for singularly perturbed problem involving discontinuous convective and source term. Math Methods Appl Sci. (2022) 45:4876–98. doi: 10.1002/mma.8077

CrossRef Full Text | Google Scholar

26. Appadu AR, Tijani YO. 1D generalised burgers-huxley: proposed solutions revisited and numerical solution using FTCS and NSFD methods. Front Appl Math Stat. (2022) 7:773733. doi: 10.3389/fams.2021.773733

CrossRef Full Text | Google Scholar

27. Ejere AH, DuressaGF, Woldaregay MM, Dinka TG. An exponentially fitted numerical scheme via domain decomposition for solving singularly perturbed differential equations with large negative shift. J Math. (2022) 2022:1–13. doi: 10.1155/2022/7974134

CrossRef Full Text | Google Scholar

28. Kumari P, Kumar D. Parameter-uniform numerical treatment of singularly perturbed initial-boundary value problems with large delay. Appl Num Math. (2020) 163:412–29. doi: 10.1016/j.apnum.2020.02.021

CrossRef Full Text | Google Scholar

29. Bansal K, Sharma KK. Parameter-Robust numerical scheme for time-dependent singularly perturbed reaction-diffusion problem with large delay. Numer Funct Anal Optim. (2018) 39:127–54. doi: 10.1080/01630563.2016.1277742

CrossRef Full Text | Google Scholar

30. Ejere AH, Duressa GF, Woldaregay MM, Dinka TG. A uniformly convergent numerical scheme for solving singularly perturbed differential equations with large spatial delay. SN Appl Sci. (2022) 4:1–15. doi: 10.1007/s42452-022-05203-9

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Alam MP, Khan A. A new numerical algorithm for time-dependent singularly perturbed differential-difference convection-diffusion equation arising in computational neuroscience. Comput Appl Math. (2022) 41:402. doi: 10.1007/s40314-022-02102-y

CrossRef Full Text | Google Scholar

32. Kadalbajoo MK, Gupta V. A brief survey on numerical methods for solving singularly perturbed problems. Appl Math Comput. (2010) 217:3641–16. doi: 10.1016/j.amc.2010.09.059

CrossRef Full Text | Google Scholar

33. Roos HG, Stynes M, Tobiska L. Robust Numerical Methods for Singularly Perturbed Differential Equations: Convection-Diffusion-Reaction and Flow Problems. Vol. 24. Springer Science & Business Media (2008). doi: 10.1007/978-3-540-34467-4

CrossRef Full Text | Google Scholar

34. Kellogg RB, Tsan A. Analysis of some difference approximations for a singular perturbation problem without turning points. Math Comput. (1978) 32:1025–39. doi: 10.1090/S0025-5718-1978-0483484-9

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Clavero C, Gracia JL. On the uniform convergence of a finite difference scheme for time dependent singularly perturbed reaction-diffusion problems. Appl Math Computat. (2010) 216:1478–88. doi: 10.1016/j.amc.2010.02.050

CrossRef Full Text | Google Scholar

36. Ehrhardt M, Mickens RE. A nonstandard finite difference scheme for convection-diffusion equations having constant coefficients. Appl Math Comput. (2013) 219:6591–604. doi: 10.1016/j.amc.2012.12.068

CrossRef Full Text | Google Scholar

37. Woldaregay MM, Duressa GF. Uniformly convergent numerical method for singularly perturbed delay parabolic differential equations arising in computational neuroscience. Kragujevac J Math. (2022) 46:65–84. doi: 10.46793/KgJMat2201.065W

CrossRef Full Text | Google Scholar

38. Munyakazi JB, Patidar KC. A fitted numerical method for singularly perturbed parabolic reaction-diffusion problems. Comput Appl Math. (2013) 32:509–19. doi: 10.1007/s40314-013-0033-7

CrossRef Full Text | Google Scholar

Keywords: singularly perturbed problem, spatio-temporal delays, nonstandard finite difference, implicit Euler method, uniform convergence

Citation: Ejere AH, Duressa GF, Woldaregay MM and Dinka TG (2023) A robust numerical scheme for singularly perturbed differential equations with spatio-temporal delays. Front. Appl. Math. Stat. 9:1125347. doi: 10.3389/fams.2023.1125347

Received: 16 December 2022; Accepted: 03 March 2023;
Published: 24 March 2023.

Edited by:

Yulong Xing, The Ohio State University, United States

Reviewed by:

Vikas Gupta, LNM Institute of Information Technology, India
Appanah Rao Appadu, Nelson Mandela University, South Africa
Arshad Khan, Jamia Millia Islamia, India

Copyright © 2023 Ejere, Duressa, Woldaregay and Dinka. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Ababi Hailu Ejere, hailu_ababi@yahoo.com

ORCID: Gemechis File Duressa orcid.org/0000-0003-1889-4690
Mesfin Mekuria Woldaregay orcid.org/0000-0002-6555-7534

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.