Front. Appl. Math. Stat. Frontiers in Applied Mathematics and Statistics Front. Appl. Math. Stat. 2297-4687 Frontiers Media S.A. 10.3389/fams.2023.1244490 Applied Mathematics and Statistics Original Research Fitted computational method for singularly perturbed convection-diffusion equation with time delay Tesfaye Sisay Ketema 1 * † Duressa Gemechis File 2 † Woldaregay Mesfin Mekuria 1 † Dinka Tekle Gemechu 1 1Department of Applied Mathematics, Adama Science and Technology University, Adama, Ethiopia 2Department of Mathematics, Jimma University, Jimma, Ethiopia

Edited by: Vikas Gupta, LNM Institute of Information Technology, India

Reviewed by: Sunil Kumar, Indian Institute of Technology (BHU), India; Sanjay Ku Sahoo, LNM Institute of Information Technology, India

 *Correspondence: Sisay Ketema Tesfaye sisayk12@gmail.com

†ORCID: Sisay Ketema Tesfaye orcid.org/0009-0005-0602-3013

Gemechis File Duressa orcid.org/0000-0003-1889-4690

Mesfin Mekuria Woldaregay orcid.org/0000-0002-6555-7534

 31 08 2023 2023 9 1244490 22 06 2023 07 08 2023 Copyright © 2023 Tesfaye, Duressa, Woldaregay and Dinka. 2023 Tesfaye, 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.

A uniformly convergent numerical scheme is proposed to solve a singularly perturbed convection-diffusion problem with a large time delay. The diffusion term of the problem is multiplied by a perturbation parameter, ε. For a small ε, the problem exhibits a boundary layer, which makes it challenging to solve it analytically or using standard numerical methods. As a result, the backward Euler scheme is applied in the temporal direction. Non-symmetric finite difference schemes are applied for approximating the first-order derivative terms, and a higher-order finite difference method is applied for approximating the second-order derivative term. Furthermore, an exponential fitting factor is computed and induced in the difference scheme to handle the effect of the small parameter. Using the discrete maximum principle, the stability of the scheme is examined and analyzed. The developed scheme is parameter-uniform with a linear order of convergence in both space and time. To examine the accuracy of the method, two model examples are considered. Further, the boundary layer behavior of the solutions is given graphically.

 singularly perturbed delay differential equation exponentially fitted finite difference non-symmetric finite difference uniform convergence section-at-acceptance Numerical Analysis and Scientific Computation 1. Introduction

Delay differential equations (DDEs) are differential equations in which the evolution of the system is influenced by its past history. DDEs are called retarded types if the delay argument does not appear in the highest-order derivative term; otherwise, they are neutral types. DDEs play an important role in a variety of fields, including robotics, biosciences [1], economics, epidemiology and mechanics [2], fluid dynamics, reaction-diffusion equations [3], and population dynamics [4].

A singularly perturbed delay differential equation (SPDDE) is a delay differential equation in which its higher-order derivative term is multiplied by a small perturbation parameter (0 < ε ≪ 1) and contains at least one delay parameter on the term different from the highest derivative. In contrast to the magnitude of the delay parameter with the perturbation parameter, the delay is classified as a large delay or a small delay. If the magnitude of the delay parameter of the SPDDE is smaller than the perturbation parameter, then the equation is said to be a singularly perturbed delay differential equation with a small delay, whereas when the magnitude of the delay parameter is higher than the perturbation parameter, it is said to be a singularly perturbed delay differential equation with a large delay [5]. A singularly perturbed problem, which arises as a time delay, occurs in many application areas of science and engineering, for instance, in the simulation of oil extraction from underground reservoirs, chemical processes, fluid flows, water quality problems in river networks, and mechanical systems [6].

The presence of ε as a multiple of the higher-order derivative term causes a boundary layer. The boundary layer is an asymptotically narrow region located in the neighborhood of the endpoints of the domain, where the solution has a steep gradient as ε tends to zero [6]. With the rapidly changing behavior of the solution in the boundary layer, one encounters computational difficulties in treating a singularly perturbed problem using analytically or classical numerical schemes. On the contrary, classical numerical schemes lead to spurious non-physical oscillations in the numerical solution, unless an unacceptably large number of mesh points are considered, which leads to a massive computational cost [7]. In response to this, different authors have to look for sounding numerical schemes which converge uniformly regardless of ε.

Recently, the authors in [8], proposed the implicit-Euler scheme in the time direction and the central difference scheme in the space direction. The authors in [7, 9], proposed the implicit-Euler scheme in the time direction and the hybrid scheme by a proper combination of the midpoint upwind in the outer region and the central difference scheme in the inner region in the spatial direction on the Shishkin mesh. Moreover, this method is addressed in [10], for the two-parameter problem. In [11], the authors proposed the implicit-Euler scheme in the time direction and the hybrid scheme on a generalized Shishkin mesh in the spatial direction. Gowrisankar and Natisan in [12] developed the backward Euler scheme in time direction and the upwind finite difference scheme in the spatial direction using a piecewise uniform mesh. The implicit Euler scheme in the time direction and the upwind scheme in the spatial direction are considered in [13]. In [14], the implicit trapezoidal scheme in the time direction and the hybrid scheme by proper combination of the midpoint upwind in the outer region and the central difference in the inner region in the spatial direction are used.

The implicit Euler scheme in the time direction and the central difference scheme in the space direction are used in [4]. The extended cubic B-spline is considered in [15]. A domain decomposition method is considered in [16, 17]. The authors in [18, 19] proposed hybrid scheme on both Shishkin and Bakhvalov meshes. Podila and Kumar [20] proposed a new stable finite difference scheme on a uniform mesh and also on an adaptive mesh. The backward Euler scheme in the time direction and exponentially fitted difference method is considered in [21]. The Crank-Nicolson method in the time direction and a novel fitted finite difference scheme in spatial direction are proposed in [22]. The Crank-Nicolson method in the time direction and an exponentially fitted spline in the spatial direction are discussed in [23]. The implicit Euler scheme in the time direction and the non-standard finite difference method in the space direction are considered in [24]. In [25], the authors proposed Crank-Nicolson method in the time direction and the operator compact implicit (OCI) method on the Shishkin mesh in the space direction. The backward Euler in the time direction and method of line following Micken's type discretization for the space derivatives are used in [26]. Sahoo and Gupta [27] used higher-order difference with an identity expansion (HODIE) on a piecewise uniform mesh. A similar technique was also used in [28] for a coupled system of singularly perturbed problems. The authors [29, 30] proposed the numerical schemes that work for both cases when the delay term is large or small.

The main aim of this work was to develop a ε-uniform numerical scheme for the class of singularly perturbed convection-diffusion problem with a large time delay. The method comprises the backward Euler scheme in the time direction and an exponentially fitted higher-order finite difference scheme in the spatial direction. Error bound and uniform convergence of the developed scheme is investigated and proved. The proposed scheme gives more accurate, stable, and uniformly convergent results.

In this study, C has been considered as a generic positive constant, which does not depend on Δs, Δt, and ε. The maximum norm is denoted by ‖.‖, which is defined by ‖γ‖ = maxs,t∈Ω|γ(s, t)|.

 2. Continuous problem

Let Ω = Ωs × Ωt = (0, 1) × (0, 𝕋] for 𝕋 > 0, we consider SPDDE of the form {zt(s,t)+𝔏εz(s,t)=-κ(s,t)z(s,t-δ)+γ(s,t),(s,t)∈Ωz(s,t)=ψb(s,t),  (s,t)∈ηb=[0,1]×[-δ,0],z(0,t)=ψl(t),t∈ηl={(0,t):0≤t≤𝕋},z(1,t)=ψr(t),t∈ηr={(1,t):0≤t≤𝕋}, where 𝔏εz(s, t) = −εzss(s, t) + β(s, t)zs(s, t) + α(s, t)z(s, t).

Here, ε ∈ (0, 1] and δ > 0 are the perturbation parameter and the delay parameter, respectively. We pretended that the functions β(s, t), α(s, t), κ(s, t), γ(s, t) on Ω̄=[0,1]×[0,𝕋] and ψb(s, t), ψl(t), ψr(t) on η = ηl ∪ ηr ∪ ηb are smooth enough and bounded which meet α(s, t) ≥ ϖ > 0, κ(s, t) ≥ φ > 0, β(s, t) ≥ μ > 0 on Ω̄. These conditions assure that problem (1) has a boundary layer near s = 1 [7].

2.1. A priori bounds

Under the premises that the data are Hölder continuous and satisfy the following compatibility conditions at the corner points and the delay terms [31], we confirm the existence and uniqueness of the solution of (1) ψl(0)=ψb(0,0), ψr(0)=ψb(1,0), dψl(0)dt-ε∂2ψb(0,0)∂s2+β(0,0)∂ψb(0,0)∂s+              α(0,0)ψb(0,0)=-κ(0,0)ψb(0,-δ)                        +γ(0,0),dψr(0)dt-ε∂2ψb(1,0)∂s2+β(1,0)∂ψb(1,0)∂s+              α(1,0)ψb(1,0)=-κ(1,0)ψb(1,-δ)                        +γ(1,0). These assumptions and conditions are fulfilled. Then, the problem (1) admits a unique solution [31].

Setting ε = 0, the reduced problem of (1) is given as {∂z0(s,t)∂t+β(s,t)∂z0(s,t)∂s+α(s,t)z0(s,t)=-κ(s,t)z0(s,t-δ)+γ(s,t),z0(s,t)=ψb(s,t),(s,t)∈ηb,z0(0,t)=ψl(t),t∈Ω̄t, where z0(s, t) is the solution of the reduced problem.

Lemma 2.1. Let z(s, t) be the solution of (1). Then, we have |z(s,t)-ψb(s,0)|≤Ct,     (s,t)∈Ω̄, where C does not depends on ε.

Proof: The proof is considered in [7].

The operator 𝔏=(∂∂t+𝔏ε) in (1) satisfies the next lemma.

Lemma 2.2. (Maximum principle). Let ν(s, t) ∈ C2(Ω) ∪ C0(η) satisfies ν(s, t) ≥ 0 (s, t) ∈ η. If 𝔏ν(s, t) ≥ 0, (s, t) ∈ Ω, then ν(s,t)≥0,(s,t)∈Ω̄.

Proof: The proof is considered in [14].

Lemma 2.3. (Stability result). Let z(s, t) be the solution of (1). Then, we have |z(s,t)|≤ϖ-1‖γ‖+max{|ψl(t)|,|ψb(s,t)|,|ψr(t)|}, where ϖ ≤ α(s, t).

Proof: The proof is considered in [22].

Lemma 2.4. The derivative of the solution z(s, t) of (1) with respect to s and t satisfy |∂iz(s,t)∂si|≤C(1+ε-iexp(-μ(1-s)ε)), (s,t)∈Ω̄, i=0(1)4,|∂lz(s,t)∂tl|≤C, (s,t)∈Ω̄, l=0(1)2, where μ ≤ β(s, t).

Proof: The proof is considered in [7].

 3. Numerical scheme 3.1. Temporal semi-discretization

The time domain [0, 𝕋] is discretized uniformly with step size Δt as ΩtM={tm=mΔt,m=0,1,2,...,M,tM=𝕋,Δt=𝕋/M,} and Ωtj={tm=mΔt,m=0,1,2,...,j,tj=δ,Δt=δ/j,} with M + 1 mesh points in [0, 𝕋] and j + 1 mesh points in [−δ, 0]. We have 𝕋 = rδ for some positive integer r.

Applying the backward Euler scheme for time derivative, we get Zm(s)-Zm-1(s)Δt-εd2Zm(s)ds2+β(s,tm)dZm(s)ds+α(s,tm)Zm(s)=-κ(s,tm)Z(s,tm-δ)+γ(s,tm). Simplifying (8), we have 𝔏ΔtZm(s)={-κ(x,tm)ψb(s)+γ(s,tm)+Zm-1(s)Δt,for m=0,1,2,...,j,s∈Ωs,-κ(s,tm)Zm-j(s)+γ(s,tm)+Zm-1(s)Δt,for m=j+1,...,M-1,s∈Ωs, where 𝔏ΔtZm(s)=-εd2Zm(s)ds2+β(s,tm)dZm(s)ds+P(s,tm)Zm(s) and P(s,tm)=(1Δt+α(s,tm)) with the boundaries Zm(0)=ψl(tm),    Zm(1)=ψr(tm),    m=0(1)M. Now, (9) rewrite as 𝔏ΔtZ(s)={-κ(s,tm)ψb(s)+γ(s,tm)+Zm-1(s)Δt,for m=0,1,2,...,j,s∈Ωs,-κ(s,tm)Zm-j(s)+γ(s,tm)+Zm-1(s)Δt,for m=j+1,...,M-1,s∈Ωs, where 𝔏ΔtZ(s)=-εd2Z(s)ds2+β(s,tm)dZ(s)ds+P(s,tm)Z(s) and  Z(s)=Zm(s)≈z(s,tm) and Q(s)=Qm(s)=Q(s,tm).

The local truncation error in the time direction is given as em(s):=z(s,tm)-Zm(s),m=0(1)M.

Lemma 3.1. The local error em at tm satisfies the bound ‖em‖≤C(Δt)2.

Lemma 3.2. The global error Em at tm satisfies the bound ‖Em‖≤C(Δt), m=1(1)M-1. Proof: Using Lemma 3.1, the global error Em bound at mth time step is given as ‖Em‖=‖∑l=1mel‖≤‖e1‖+‖e2‖+‖e3‖+...+‖em‖        ≤C1T(Δt),   since(m)Δt≤T        =C(Δt),   whereC1T=C,

Lemma 3.3. For every m = 0(1)M − 1, the solution Zm(s) of (9)-(10) satisfies the estimate |diZm(s)dsi|≤C(1+ε-iexp(-με(1-s))),s∈Ω̄s,i=0(1)4. Proof: From (11), -εd2Z(s)ds2+β(s,tm)dZ(s)ds=g, where g = Q(s) − P(s, tm)Z(s).

Now, we integrate twice and we obtain Z(s)=ZP(s)+C1+C2∫s1exp(-ε-1(B(1)-B(y)))dy, where ZP(s)=-∫s1u(y)dy,u(s)=∫s1ε-1g(y)exp(-ε-1(B(y)-B(s)))dy,B(s)=∫0sβ(y)dy.

Using inequality exp(-ε-1(B(y)-B(s)))≤exp(-ε-1μ(y-s)),s≤y, and the bound |u(s)|≤Cε-1∫s1(exp(-ε-1μ(y-s))+Cε-1exp(-ε-1μ(1-s)))dy       ≤C(1+ε-2(1-s)exp(-ε-1μ(1-s))). Hence, |Zp(s)| ≤ C. Here, Z′(1)=-C2. The boundary condition Z(1) = 0 yields C1 = 0. Now, the constants C1 and C2 must satisfy C2∫01exp(-ε-1(B(1)-B(y)))dy=-Zp(0). Since B(s) is bounded on (0, 1), B(1) − B(y) ≤ C(1 − y). Then, ∫01exp(-ε-1(B(1)-B(y)))dy≥∫01exp(-ε-1-μ(1-y))dy≥Cε. It follows that |C2|≤Cε-1. Hence, |Z′(1)|=|C2|≤Cε-1. Finally, Z′(s)=u(s)-C2exp(-ε-1(B(1)-B(s))) implies that |dZ(s)ds|≤C(1+ε-1exp(-με(1-s))),s∈Ω̄s. The proof is done for i = 1. For i > 1 follows by induction and repeated differentiation. For the details, refer [32].

 3.2. Spatial discretization

We discretize the spatial domain [0, 1] into N equal number of sub-intervals with the length of h as 0 = s0, s1, ..., sN = 1, and sn = nh, n = 0(1)N. Consider a smooth function Z(s) in the interval [0, 1]. From Taylor's series approximation, we get Zn+1=Zn+hZn ′+h22!Zn ″+h33!Zn ‴+h44!Zn(4)+h55!Zn(5)+h66!Zn(6)+         h77!Zn(7)+h88!Zn(8)+O(h9),Zn-1=Zn-hZn ′+h22!Zn ″-h33!Zn ‴+h44!Zn(4)-h55!Zn(5)+h66!Zn(6)-         h77!Zn(7)+h88!Zn(8)-O(h9). Following a similar relation of (15), it holds Zn-1-2Zn+Zn+1=2h22!Zn ″+2h44!Zn(4)+2h66!Zn(6)+2h88!Zn(8)+O(h10),Zn-1 ″-2Zn ″+Zn+1 ″=2h22!Zn(4)+2h44!Zn(6)+2h66!Zn(8)+2h88!Zn(10)+O(h12). From (16), we have h412Zn(6)=Zn-1 ″-2Zn ″+Zn+1 ″-h2Zn(4)-2h66!Zn(8)-2h88!Zn(10)-O(h12). Substituting (17) into (16) and simplifying, we obtain Zn-1-2Zn+Zn+1=h230(Zn-1 ″+28Zn ″+Zn+1 ″)+𝔗, where 𝔗=h420Zn(4)-13h8302,400Zn(8)+O(h10).

From (11), we draw -εd2Z(s)ds2=-β(s,tm)dZ(s)ds- P(s,tm)Z(s)+Q(s), where Q(s)={-κ(s,tm)ψb(s)+γ(s,tm)+Zm-1(s)Δt,for m=0,1,2,...,j,s∈Ωs,-κ(s,tm)Zm-j(s)+γ(s,tm)+Zm-1(s)Δt,for m=j+1,...,M-1,s∈Ωs. Using (19), we have -εZn+1 ″=-β(sn+1,tm)Zn+1 ′-P(sn+1,tm)Zn+1+Qn+1,  -εZn ″=-β(sn,tm)Zn ′-P(sn,tm)Zn+Qn,-εZn-1 ″=-β(sn-1,tm)Zn-1 ′-P(sn-1,tm)Zn-1+Qn-1. From the Taylor series approximations of Zn-1 ′,Zn ′ and Zn+1 ′, we get    Zn ′=Zn+1-Zn-12h,Zn+1 ′=3Zn+1-4Zn+Zn-12h-hZn ″,Zn-1 ′=-Zn+1+4Zn-3Zn-12h+hZn ″. Substituting (21) into (20), we have -εZn+1 ″=-β(sn+1,tm)(3Zn+1-4Zn+Zn-12h-hZn ″)-P(sn+1,tm)Zn+1+Qn+1,-εZn ″=-β(sn,tm)Zn+1-Zn-12h-P(sn,tm)Zn+Qn,-εZn-1 ″=-β(sn-1,tm)(-Zn+1+4Zn-3Zn-12h+hZn ″)-P(sn-1,tm)Zn-1+Qn-1. From (18), we draw -ε(Zn-1-2Zn+Zn+1h2)=130(-εZn-1 ″-28εZn ″-εZn+1 ″)+𝔗. Substituting (22) into (23) and rearranging, we obtain -(ε-hβ(sn-1,tm)30+hβ(sn+1,tm)30)(Zn-1m-2Znm+Zn+1mh2)+β(sn-1,tm)60h(-3Zn-1m+4Znm-Zn+1m)+28β(sn,tm)60h(Zn+1m-Zn-1m)+β(sn+1,tm)60h(Zn-1m-4Znm+3Zn+1m)+P(sn-1,tm)30Zn-1m+28P(sn,tm)30Znm+P(sn+1,tm)30Zn+1m=130(Qn-1m+28Qnm+Qn+1m)+𝔗, where Qnm={-κ(sn,tm)ψb(sn)+γ(sn,tm)+Znm-1Δt,for m=0,1,2,...,j,s∈Ωs,-κ(sn,tm)Znm-j+γ(sn,tm)+Znm-1Δt,for m=j+1,...,M-1,s∈Ωs.

3.2.1. Computing the exponential fitting factor

We introduce the exponential fitting factor σ to handle the effect of ε in the layer. From the singular perturbation theory stated in [33], the zero order asymptotic solution of the problem of the form {-εZ ″(s)+β(s)Z ′(s)+α(s)Z(s)=q(s),s∈Ωs=(0,1),Z(0)=ωl,Z(1)=ωr is given as      Z(s)≈Z0(s)+β(0)β(s)(ωr-Z0(1))exp(-∫s1(β(s)ε-α(s)β(s))ds)+O(ε). From Taylor's series, approximation for β(s) and α(s) restricting to their first terms about s = 1 is given as Z(s)≈Z0(s)+(ωr-Z0(1))exp(-β(1)(1-s)ε), where Z0(s) is the solution of reduced problem. Taking h → 0 and solving (28) at sn−1, sn, and sn+1, we get    Z(sn)≈Z0(0)+(ωr-Z0(1))exp(-β(1)(1ε-nρ)),Z(sn-1)≈Z0(0)+(ωr-Z0(1))exp(-β(1)(1ε-(n-1)ρ)),Z(sn+1)≈Z0(0)+(ωr-Z0(1))exp(-β(1)(1ε-(n+1)ρ)), where ρ=hε. Multiplying (24) by h and the term containing ε by σ and evaluating the limit of the resulting equation as h → 0, we get -σρlimh→0(Zn-1m-2Znm+Zn+1m)+β(sn-1,tm)60limh→0(-3Zn-1m+4Znm-Zn+1m)+28β(sn,tm)60limh→0(Zn+1m-Zn-1m)+β(sn+1,tm)60limh→0(Zn-1m-4Znm+3Zn+1m)=0. From (29), we have limh→0(Z((n-1)h)-2Z(nh)+Z((n+1)h))=(ωr-Z0(1))e-β(1)(1-snε)(e-β(1)ρ+eβ(1)ρ-2),limh→0(-3Z((n-1)h)+4Z(nh)-Z((n+1)h))=(ωr-Z0(1))e-β(1)(1-snε)(-3e-β(1)ρ-eβ(1)ρ+4),limh→0(Z((n-1)h)-4Z(nh)+3Z((n+1)h))=(ωr-Z0(1))e-β(1)(1-snε)(e-β(1)ρ+3eβ(1)ρ-4),limh→0(Z((n+1)h)-Z((n-1)h))=(ωr-Z0(1))e-β(1)(1-snε)(eβ(1)ρ-e-β(1)ρ). Substituting (31) into (30) and simplifying yields σρ(eβ(1)ρ+e-β(1)ρ-2)=β(sn,tm)60(30eβ(1)ρ-30e-β(1)ρ). Then, we get the fitting factor σ σ=ρβ(sn,tm)2coth(ρβ(1)2). Therefore, the required scheme is taken as 𝔏Δt,hZnm=130(Qn-1m+28Qnm+Qn+1m), where 𝔏Δt,hZnm=-(εσ-hβ(sn-1,tm)30+hβ(sn+1,tm)30)(Zn-1m-2Znm+Zn+1mh2)+β(sn-1,tm)60h(-3Zn-1m+4Znm-Zn+1m)+28β(sn,tm)60h(Zn+1m-Zn-1m)+β(sn+1,tm)60h(Zn-1m-4Znm+3Zn+1m)+P(sn-1,tm)30Zn-1m+28P(sn,tm)30Znm+P(sn+1,tm)30Zn+1m. In the explicit form, it becomes Rn-Zn-1m+Rn0Znm+Rn+Zn+1m=Hnm, where Rn-=-1h2(εσ-hβ(sn-1,tm)30+hβ(sn+1,tm)30)-3β(sn-1,tm)60h+P(sn-1,tm)30-28β(sn,tm)60h+β(sn+1,tm)60h,Rn0=2h2(εσ-hβ(sn-1,tm)30+hβ(sn+1,tm)30)+4β(sn-1,tm)60h-4β(sn+1,tm)60h+28P(sn,tm)30,Rn+=-1h2(εσ-hβ(sn-1,tm)30+hβ(sn+1,tm)30)-β(sn-1,tm)60h+28β(sn,tm)60h+3β(sn+1,tm)60h+P(sn+1,tm)30,Hnm=130(Qn-1m+28Qnm+Qn+1m).

 3.3. Stability and uniform convergence analysis

Lemma 3.4. (Discrete maximum principle). Assume that Z0m≥0,ZNm≥0 and 𝔏Δt,hZnm≥0,∀n=1(1)N-1, then Znm≥0,∀n=0(1)N.

Proof: Assume that there is k ∈ {0, 1, 2, …, N}, such that Zkm=min0≤n≤NZnm<0. Assume that Zkm<0 and from the assumption, it is shown that k ∉ {0, 1}. So, we have Zk+1m-Zkm>0 and Zkm-Zk-1m<0. Then, we get 𝔏Δt,hZkm<0 for k = 1(1)N − 1. So, the assumption Znm<0∀n=0(1)N is wrong. Therefore, Znm≥0 and ∀n = 0(1)N.

Lemma 3.5. (Uniform stability result). Let Znm be the solution of (33), then we have |Znm|≤‖ 𝔏Δt,hZnm‖ζ+max{|ψl(tm)|,|ψr(tm)|}, where P(sn, tm) ≥ ζ > 0.

Proof: Let R=‖ 𝔏Δt,hZnm‖ζ+max{|ψl(tm)|,|ψr(tm)|} and define the barrier functions ϑn,m± by ϑn,m±=R±Znm. On the boundaries, we get ϑ0,m±=R±Z0m=‖ 𝔏Δt,hZnm‖ζ+max{|ψl(tm)|,|ψr(tm)|}± ψl(tm) ≥ 0 and ϑN,m±=R±ZNm=‖ 𝔏Δt,hZnm‖ζ+max{|ψl(tm)|,|ψr(tm)|}±ψr(tm)≥0. For sn, n = 1(1)N − 1, we obtain 𝔏Δt,hϑn,m=-(εσ-hβ(sn-1,tm)30+hβ(sn+1,tm)30)(R±Zn-1m-2(R±Znm)+R±Zn+1mh2)+β(sn-1,tm)60h(-3(R±Zn-1m)+4(R±Znm)-(R±Zn+1m))+28β(sn,tm)60h(R±Zn+1m-(R±Zn-1m))+β(sn+1,tm)60h(R±Zn-1m-4(R±Znm)+3(R±Zn+1m))+P(sn-1,tm)30(R±Zn-1m)+28P(sn,tm)30(R±Znm)+P(sn+1,tm)30(R±Zn+1m)=(P(sn-1,tm)30+28P(sn,tm)30+P(sn+1,tm)30)R±𝔏Δt,hZnm=(P(sn-1,tm)30+28P(sn,tm)30+P(sn+1,tm)30)(‖ 𝔏Δt,hZnm‖ζ+max{|ψl(tm)|,|ψr(tm)|})±130(Qn-1m+28Qnm+Qn+1m)≥0. By Lemma 3.4, we get ϑn,m±≥0, n=0(1)N. Therefore, the needed bound is obtained.

From Taylor's series expansion, we get |-(dds2-δs2)Zm(sn)|≤Ch2‖d4Zm(sn)ds4‖,|dZm(sn-1)ds-(-Zn+1m+4Znm-3Zn-1m2h+hd2Zm(sn)ds2)|≤Ch2‖d3Zm(sn)ds3‖,|dZm(sn+1)ds-(3Zn+1m-4Znm+Zn-1m2h-hd2Zm(sn)ds2)|≤Ch2‖d3Zm(sn)ds3‖,|(dds-δs0)Zm(sn)|≤Ch2‖d3Zm(sn)ds3‖,|δs2Zm(sn)|≤C‖d2Zm(sn)ds2‖, where ‖Z(k)(sn)‖=maxsn∈(s0,sN)|Zm(k)(sn)|,k=2,3,4.

The next theorem provides the truncation error estimate for the developed scheme.

Theorem 3.6. Let the coefficients α(s, tm), β(s, tm), and κ(s, tm) of (9)-(10) be sufficiently smooth such that Zm(s) ∈ C4[0, 1]. Then, the solution Znm of (33) satisfies the next bound |𝔏Δt,h(Zm(sn)-Znm)|≤Ch2h+ε(1+ε-3exp(-μ(1-sn)ε)). Proof: The error estimate in the spatial direction is given as |𝔏Δt,h(Zm(sn)-Znm)|=|-ε(dds2-σδs2)Zm(sn)+β(sn-1,tm)30(dZm(sn-1)ds-(-Zn+1m+4Znm-3Zn-1m2h+hd2Zm(sn)ds2))+28β(sn,tm)30(dds-δs0)Zm(sn)+β(sn+1,tm)30(dZm(sn+1)ds-(3Zn+1m-4Znm+Zn-1m2h-hd2Zm(sn)ds2))|≤|ε(β(sn,tm)ρ2coth(β(1)ρ2)-1)δs2Zm(sn)|+|ε(d2ds2-δs2)Zm(sn)|+|β(sn-1,tm)30(dZm(sn-1)ds-(-Zn+1m+4Znm-3Zn-1m2h+hd2Zm(sn)∂s2))|+|28β(sn,tm)30(dds-δs0)Zm(sn)|+|β(sn+1,tm)30(dZn(sn+1)ds-(3Zn+1m-4Znm+Zn-1m2h-hd2Zm(sn)ds2))|, where σ=β(sn,tm)ρ2coth(β(1)ρ2) and ρ=hε.

For the constants C1 and C2, we have |ρcoth(ρ)-1|≤C1ρ2 for ρ ≤ 1. For ρ → ∞, since limρ→∞coth(ρ)=1 which gives |ρ coth(ρ) − 1| ≤ C1ρ.

Generally, ∀ρ > 0, we express as C1ρ2ρ+1≤ρcoth(ρ)-1≤C2ρ2ρ+1 and we have ε[β(sn,tm)ρ2coth(β(1)ρ2)-1]≤ε(h/ε)2h/ε+1=h2h+ε. From the bounds in (38), (40), and (42), we have |𝔏Δt,h(Zm(sn)-Znm)|≤Ch2h+ε‖d2Zm(sn)ds2‖+Ch2‖d3Zm(sn)ds3‖                  +Cεh2‖d4Zm(sn)ds4‖. By Lemma 3.3, we have |𝔏Δt,h(Zm(sn)-Znm)|≤Ch2h+ε(1+ε-2exp(-μ(1-sn)ε))                    +Ch2[(1+ε-3exp(-μ(1-sn)ε))                    +(ε+ε-3exp(-μ(1-sn)ε))]. Obviously, ε−3 ≥ ε−2, then we draw |𝔏Δt,h(Zm(sn)-Znm)|≤Ch2h+ε(1+ε-3exp(-μ(1-sn)ε)) thus, it gives the wanted bound.

Lemma 3.7. For a fixed mesh and as ε → 0, it holds          limε→0maxiexp(-μsn/ε)εi=0,limε→0maxiexp(-μ(1-sn)/ε)εi=                    0, i=1,2,3,..., where sn = nh, 1 ≤ n ≤ N − 1.

Proof: The proof is in [22].

Theorem 3.8. Let Znm be the solution of (33), then we have the following uniform error bound supε∈(0,1]maxn|Zm(sn)-Znm|≤Ch,n=0(1)N. Proof: Substituting Lemma 3.7 into (39), we arrive at |𝔏Δt,h(Zm(sn)-Znm)|≤Ch2h+ε.

Hence, the result leads |Zm(sn)-Znm|≤Ch2h+ε. Using the sup over all ε ∈ (0, 1], we get supε∈(0,1]maxn|Zm(sn)-Znm|≤Ch,n=0(1)N. From (46), when ε > h, the obtained method uniformly converges uniformly with order two in the space direction. When ε ≪ h, the method converges uniformly with order one in the space direction.

Theorem 3.9. Let z and Z are the solutions of (1) and (33), respectively, then we have the following uniform error bound supε∈(0,1]|z-Z|≤C(h+(Δt)). Proof: The proof is considered by combining of Lemma 3.1 and Theorem 3.8.

 4. Numerical results

Considering two test examples we carry out some numerical inquiries to confirm the developed scheme is ε-uniform convergent. Since the exact solution of the examples are not known, we used a variant of double mesh principle is applied for the numerical inquiries. So, we calculate the maximum pointwise error by EεN,M=maxn,m|Zn,mN,M-Zn,m2N,2M|, the ε-uniform error by EN,M=maxn,m(EεN,M), the rate of convergence by rεN,M=log2(EεN,M/Eε2N,2M), and the ε-uniform rate of convergence by rN,M = log2(EN,M/E2N, 2M).

4.1. Example

Consider the problem [7] ∂z∂t-ε∂2z∂s2+(2-s2)∂z∂s+su(s,t)+u(s,t-δ)=            10t2exp(-t)s(1-s),

(s, t) ∈ (0, 1) × (0, 2] with interval condition z(s, t) = 0, on (s, t) ∈ [0, 1] × [−1, 0] and the boundary conditions z(0, t) = 0 and z(1, t) = 0, t ∈ [0, 2].

 4.2. Example

Consider the problem [13] ∂z∂t-ε∂2z∂s2+(2-s2)∂z∂s+(s+1)(t+1)z(s,t)+z(s,t-δ)=10t2exp(-t)s(1-s),

(s, t) ∈ (0, 1) × (0, 2] with interval condition z(s, t) = 0, on (s, t) ∈ [0, 1] × [−1, 0] and the boundary conditions z(0, t) = 0 and z(1, t) = 0, t ∈ [0, 2].

For distinguishable values of ε and N, the obtained results for the model Examples 4.1 and 4.2, respectively, EεN,M,EN,M, rεN,M, and rN,M of the developed scheme are delineated in Tables 1, 2. From these tables, one can observe that the maximum absolute error decreases as the step sizes decrease for every value of ε, and as ε approaches to zero, the maximum absolute error after getting large becomes constant, which displays ε-uniform convergence of the proposed scheme regardless of ε. On the other hand, the calculated EN,M and the corresponding rN,M using the proposed scheme are given in the last two rows, which confirms that the theoretical finding of the developed scheme is order one in both space and time direction.

EεN,M,EN,M,rεN,M, and rN,M for Example 4.1.

ε↓ Number of intervals N = M
16 32 64 128 256
20 2.1285e-04 1.4601e-05 2.7352e-05 2.0267e-05 1.1789e-05
3.8657 -0.9056 0.4325 0.7817 -
2−2 9.8466e-04 8.9095e-05 1.1378e-04 8.8529e-05 5.2297e-05
3.4662 -0.3528 0.3620 0.7594 -
2−4 3.6424e-03 5.2203e-04 2.1648e-04 1.6404e-04 9.7451e-05
2.8027 1.2699 0.4002 0.7513 -
2−6 1.2902e-02 3.8264e-03 7.1066e-04 3.0405e-04 1.5008e-04
1.7535 2.4288 1.2249 1.0186 -
2−8 1.5384e-02 7.4554e-03 3.1045e-03 9.0706e-04 2.8724e-04
1.0451 1.2639 1.7751 1.6589 -
2−10 1.5389e-02 7.6241e-03 3.7907e-03 1.8490e-03 7.7045e-04
1.0133 1.0081 1.0357 1.2630 -
2−12 1.5389e-02 7.6241e-03 3.7923e-03 1.8943e-03 9.4730e-04
1.0133 1.0075 1.0014 0.9998 -
2−14 1.5389e-02 7.6241e-03 3.7923e-03 1.8943e-03 9.4771e-04
1.0133 1.0075 1.0014 0.9992 -
2−16 1.5389e-02 7.6241e-03 3.7923e-03 1.8943e-03 9.4771e-04
1.0133 1.0075 1.0014 0.9992 -
2−18 1.5389e-02 7.6241e-03 3.7923e-03 1.8943e-03 9.4771e-04
1.0133 1.0075 1.0014 0.9992 -
2−20 1.5389e-02 7.6241e-03 3.7923e-03 1.8943e-03 9.4771e-04
1.0133 1.0075 1.0014 0.9992 -
EN,M 1.5389e-02 7.6241e-03 3.7923e-03 1.8943e-03 9.4771e-04
rN,M 1.0133 1.0075 1.0014 0.9992 −

EεN,M,EN,M,rεN,M, and rN,M for Example 4.2.

ε↓ Number of intervals N = M
16 32 64 128 256
20 1.3602e-04 1.4927e-05 1.4266e-05 1.0773e-05 6.3185e-06
3.1878 0.0653 0.4052 0.7698 -
2−2 3.8939e-04 6.6953e-05 6.2477e-05 3.9473e-05 2.1855e-05
2.5400 0.0998 0.6625 0.8529 -
2−4 9.8460e-04 2.5366e-04 1.3556e-04 6.9903e-05 3.5485e-05
1.9566 0.9040 0.9555 0.9782 -
2−6 3.5833e-03 1.2258e-03 4.0647e-04 1.4312e-04 5.5885e-05
1.5476 1.5925 1.5059 1.3567 -
2−8 4.8773e-03 2.9551e-03 1.2250e-03 4.0706e-04 1.2709e-04
0.7229 1.2704 1.5895 1.6794 -
2−10 4.8804e-03 3.0579e-03 1.6631e-03 8.3185e-04 3.3471e-04
0.6745 0.8787 0.9995 1.3134 -
2−12 4.8804e-03 3.0579e-03 1.6641e-03 8.6096e-04 4.3655e-04
0.6745 0.8778 0.9507 0.9798 -
2−14 4.8804e-03 3.0579e-03 1.6641e-03 8.6096e-04 4.3682e-04
0.6745 0.8778 0.9507 0.9789 -
2−16 4.8804e-03 3.0579e-03 1.6641e-03 8.6096e-04 4.3682e-04
0.6745 0.8778 0.9507 0.9789 -
2−18 4.8804e-03 3.0579e-03 1.6641e-03 8.6096e-04 4.3682e-04
0.6745 0.8778 0.9507 0.9789 -
2−20 4.8804e-03 3.0579e-03 1.6641e-03 8.6096e-04 4.3682e-04
0.6745 0.8778 0.9507 0.9789 -
EN,M 4.8804e-03 3.0579e-03 1.6641e-03 8.6096e-04 4.3682e-04
rN,M 0.6745 0.8778 0.9507 0.9789 −

In Figures 1, 2, the numerical solutions of the method for Examples 4.1 and 4.2 for different values of ε are given, respectively, for N = 80 and M = 40. Figure 3 displays the effect ε on the solutions profile of the developed scheme for Examples 4.1 and 4.2. From the figures, we see that a strong boundary layer is created on the right side of the spatial domain as ε close to zero. Furthermore, in Figure 4, the maximum point wise errors of the scheme is shown by the log-log scale plot. From these figures, one can observe that maximum absolute error decreases as the step sizes decrease for every values of ε, which confirm ε-uniform convergence of the proposed scheme.

Numerical solution of Example 4.1 for ε = 2−6 and ε = 2−20, respectively.

Numerical solution of Example 4.2 for ε = 2−6 and ε = 2−20, respectively.

Effect of ε on solution profiles for Examples 4.1 and 4.2, respectively.

Maximum point wise error in log-log scale plot for Examples 4.1 and 4.2, respectively.

In Table 3, the comparison with results of the developed method with the existing recently published studies of [23, 29] are given for Example 4.1. In Table 4, the comparison with results of the developed method with the existing number of recently published studies of [15, 24, 29, 30] are given for Example 4.2. As one follows, the developed scheme holds more accurate.

EN,M and rN,M for Example 4.1.

Number of intervals N = M
Schemes ↓ 16 32 64 128 256
Present method
EN,M 1.5389e-02 7.6241e-03 3.7923e-03 1.8943e-03 9.4771e-04
rN,M 1.0133 1.0075 1.0014 0.9992 −
Method in [23]
EN,M − 7.2307e-03 3.8523e-03 1.9892e-03 1.0107e-03
rN,M − 0.90842 1.0062 1.0155 0.98837
Method in [29]
EN,M 3.41e-02 1.84e-02 9.38e-03 4.67e-03 2.31e-03
rN,M 0.8901 0.9720 1.0062 1.0155 1.0063

EN,M and rN,M for Example 4.2.

Number of intervals N = M
Schemes ↓ 16 32 64 128 256
Present method
EN,M 4.8804e-03 3.0579e-03 1.6641e-03 8.6096e-04 4.3682e-04
rN,M 0.6745 0.8778 0.9507 0.9789 −
Method in [30]
EN,M 6.40e-03 3.43e-03 1.75e-03 8.85e-04 4.44e-04
rN,M 0.89986 0.97085 0.98361 0.99512 −
Method in [15]
EN,M 1.86e-2 1.00e-2 5.48e-3 2.86e-3 1.46e-3
rN,M 0.89 0.87 0.94 0.97 1.11
Method in [29]
EN,M 3.06e-02 1.72e-02 9.00e-03 4.58e-03 2.30e-03
rN,M 0.8311 0.9344 0.9746 0.9937 1.0000
Method in [24]
EN,M 4.9006e-02 2.8622e-02 1.5141e-02 7.7173e-03 3.8858e-03
rN,M 0.775836 0.918608 0.972338 0.989898 0.995894
5. Conclusion

We have developed a numerical method for solving singularly perturbed parabolic convection-diffusion equation with a large time delay. The solution of the problem exhibits a boundary layer on the right side of the domain. The solution has a steep gradient in the layer region due to the presence of ε. In the rapidly changing behavior of the solution in the layer region, one encounters computational difficulties to find the solution using analytically or using classical numerical methods. To handle this effect, we developed method comprises of the backward Euler scheme in the time direction and an exponentially fitted higher order finite difference scheme in the spatial direction. Using comparison principle, the stability of the discrete scheme is analyzed. The stability and uniformly convergent of the method are discussed theoretically. Numerical results are delineated by applying maximum point wise error, ε-uniform error and ε-uniform rate of convergence in tables which are in acceptable agreement with the theoretical analysis. The developed method contributes more accurate, stable, and ε-uniform with a linear order of convergence in the spatial and in the time direction. The proposed scheme can be extended for singularly perturbed turning point problems.

 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

ST and MW carried out the scheme development, algorithms writing, MATLAB code writing, the numerical simulations, and write final version of the manuscript. GD and TD planned the problem, design, wrote draft of the manuscript, and revised the manuscript. All authors read, commented, and approved the submitted version of the manuscript.

The authors would like to thank the referees for their constructive comments that improved the quality of this article.

 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 Baker CT Bocharov GA Rihan FA. A Report on the Use of Delay Differential Equations in Numerical Modelling in the Biosciences. Citeseer (1999). Gopalsamy K. Stability and Oscillations in Delay Differential Equations of Population Dynamics. vol. 74. Springer Science and Business Media (1992). Bestehorn M Grigorieva EV. Formation and propagation of localized states in extended systems. Ann Phys. (2004) 516:423–31. 10.1002/andp.200451607-80630070523 Govindarao L Mohapatra J Das A. A fourth-order numerical scheme for singularly perturbed delay parabolic problem arising in population dynamics. J Appl Math Comput. (2020) 63:171–95. 10.1007/s12190-019-01313-7 Tian H. Numerical Treatment of Singularly Perturbed Delay Differential Equations. Manchester, UK: The University of Manchester (2000).34526134 Ansari AR Bakr SA Shishkin GI. A parameter-robust finite difference method for singularly perturbed delay parabolic partial differential equations. J Comput and Appl Math. (2007) 205:552–66. 10.1016/j.cam.2006.05.032 Das A Natesan S. Uniformly convergent hybrid numerical scheme for singularly perturbed delay parabolic convection–diffusion problems on Shishkin mesh. Appl Math Comput. (2015) 271:168–86. 10.1016/j.amc.2015.08.137 Gowrisankar S Natesan S. A robust numerical scheme for singularly perturbed delay parabolic initial-boundary-value problems on equidistributed grids. Electron Trans Numer Anal. (2014) 41:376–95. 10.1016/j.cpc.2014.04.004 Gupta V Kumar M Kumar S. Higher order numerical approximation for time dependent singularly perturbed differential-difference convection-diffusion equations. Num Methods Partial Differ Equ. (2018) 34:357–80. 10.1002/num.22203 Sumit Kumar S Kuldeep Kumar M. A robust numerical method for a two-parameter singularly perturbed time delay parabolic problem. Comput Appl Math. (2020) 39:1–25. 10.1007/s40314-020-01236-1 Kumar S Kumar M. High order parameter-uniform discretization for singularly perturbed parabolic partial differential equations with time delay. Comput Math Appl. (2014) 68:1355–67. 10.1016/j.camwa.2014.09.004 Gowrisankar S Natesan S. ε-Uniformly convergent numerical scheme for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. (2017) 94:902–21. 10.1080/00207160.2016.1154948 Das A Natesan S. Second-order uniformly convergent numerical method for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. (2018) 95:490–510. 10.1080/00207160.2017.1290439 Govindarao L Mohapatra J. A second order numerical method for singularly perturbed delay parabolic partial differential equation. Eng Comput. (2018) 36:420–44. 10.1108/EC-08-2018-0337 Kumar D Kumari P. A parameter-uniform scheme for singularly perturbed partial differential equations with a time lag. Num Methods Partial Differ Equ. (2020) 36:868–86. 10.1002/num.22455 Singh J Kumar S Kumar M. A domain decomposition method for solving singularly perturbed parabolic reaction-diffusion problems with time delay. Num Methods Partial Differ Equ. (2018) 34:1849–66. 10.1002/num.22256 Rao SCS Kumar S. An almost fourth order uniformly convergent domain decomposition method for a coupled system of singularly perturbed reaction–diffusion equations. J Comput Appl Math. (2011) 235:3342–54. 10.1016/j.cam.2011.01.047 Kumar S Kumar M. A second order uniformly convergent numerical scheme for parameterized singularly perturbed delay differential problems. Num Algor. (2017) 76:349–60. 10.1007/s11075-016-0258-9 Kumar S Kumar M. Analysis of some numerical methods on layer adapted meshes for singularly perturbed quasilinear systems. Num Algor. (2016) 71:139–50. 10.1007/s11075-015-9989-2 Podila PC Kumar K. A new stable finite difference scheme and its convergence for time-delayed singularly perturbed parabolic PDEs. Comput Appl Math. (2020) 39:1–16. Tesfaye SK Woldaregay MM Dinka TG Duressa GF. Fitted computational method for solving singularly perturbed small time lag problem. BMC Res Notes. (2022) 15:318. 10.1186/s13104-022-06202-036221103 Woldaregay MM Aniley WT Duressa GF. Novel numerical scheme for singularly perturbed time delay convection-diffusion equation. Adv Math Phys. (2021) 2021:6641236. 10.1155/2021/6641236 Negero NT Duressa GF. An efficient numerical approach for singularly perturbed parabolic convection-diffusion problems with large time-lag. J Math Model. (2022) 10:173–90. 10.22124/jmm.2021.19608.1682 Babu G Bansal K. A high order robust numerical scheme for singularly perturbed delay parabolic convection diffusion problems. J Appl Math Comput. (2021) 68:363–89. 10.1007/s12190-021-01512-1 Salama AA Al-Amery DG. A higher order uniformly convergent method for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. (2017) 94:2520–46. 10.1080/00207160.2017.1284317 Negero NT Duressa GF. A method of line with improved accuracy for singularly perturbed parabolic convection–diffusion problems with large temporal lag. Results Appl Math. (2021) 11:100174. 10.1016/j.rinam.2021.100174 Sahoo SK Gupta V. Parameter robust higher-order finite difference method for convection-diffusion problem with time delay. Num Methods Partial Differ Equ. (in press). 10.1002/num.23039 Rao SCS Kumar S. Second order global uniformly convergent numerical method for a coupled system of singularly perturbed initial value problems. Appl Math Comput. (2012) 219:3740–53. 10.1016/j.amc.2012.09.075 Kumar D Kumari P. A parameter-uniform numerical scheme for the parabolic singularly perturbed initial boundary value problems with large time delay. J Appl Math Comput. (2019) 59:179–206. 10.1007/s12190-018-1174-z Negero NT Duressa GF. Uniform convergent solution of singularly perturbed parabolic differential equations with general temporal-lag. Iran J Sci Technol Trans A Sci. (2022) 46:507–24. 10.1007/s40995-021-01258-2 Protter M Weinberger H. Maximum Principles in Differential Equations. Englewood Cliffs, NJ: Prentice-Hall Inc. (1967). Kellogg RB Tsan A. Analysis of some difference approximations for a singular perturbation problem without turning points. Math Comp. (1978) 32:1025–39. O'Malley RE. Singular Perturbation Methods for Ordinary Differential Equations. vol. 89. New York, NY: Springer (1991).23077920