Skip to main content

ORIGINAL RESEARCH article

Front. Appl. Math. Stat., 31 August 2023
Sec. Numerical Analysis and Scientific Computation
Volume 9 - 2023 | https://doi.org/10.3389/fams.2023.1244490

Fitted computational method for singularly perturbed convection-diffusion equation with time delay

  • 1Department of Applied Mathematics, Adama Science and Technology University, Adama, Ethiopia
  • 2Department of Mathematics, Jimma University, Jimma, Ethiopia

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.

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):0t𝕋},z(1,t)=ψr(t),  tηr={(1,t):0t𝕋},    (1)

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),    (2)
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).    (3)

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,    (4)

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)Ω̄,    (5)

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)|},    (6)

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

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 𝕋 = 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).    (8)

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,    (9)

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.    (10)

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,    (11)

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

emC(Δt)2.    (12)

Lemma 3.2. The global error Em at tm satisfies the bound

EmC(Δt), m=1(1)M-1.    (13)

Proof: Using Lemma 3.1, the global error Em bound at mth time step is given as

Em=l=1mele1+e2+e3+...+em              C1T(Δt),   since   (m)ΔtT             =C(Δt),   where   C1T=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.    (14)

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+C2s1exp(-ε-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)),sy,

and the bound

|u(s)|Cε-1s1(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

C201exp(-ε-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)))dy01exp(-ε-1-μ(1-y))dyCε.

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).    (15)

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).    (16)

From (16), we have

h412Zn(6)=Zn-1-2Zn+Zn+1-h2Zn(4)-2h66!Zn(8)-2h88!Zn(10)-O(h12).    (17)

Substituting (17) into (16) and simplifying, we obtain

Zn-1-2Zn+Zn+1=h230(Zn-1+28Zn+Zn+1)+𝔗,    (18)

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),    (19)

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.    (20)

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.    (21)

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.    (22)

From (18), we draw

-ε(Zn-1-2Zn+Zn+1h2)=130(-εZn-1-28εZn-εZn+1)+𝔗.    (23)

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)+𝔗,    (24)

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.    (25)

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    (26)

is given as

           Z(s)Z0(s)+β(0)β(s)(ωr-Z0(1))exp(-s1(β(s)ε-α(s)β(s))ds)+O(ε).    (27)

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)ε),    (28)

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)ρ)),    (29)

where ρ=hε. Multiplying (24) by h and the term containing ε by σ and evaluating the limit of the resulting equation as h → 0, we get

-σρlimh0(Zn-1m-2Znm+Zn+1m)+β(sn-1,tm)60limh0(-3Zn-1m+4Znm-Zn+1m)+28β(sn,tm)60limh0(Zn+1m-Zn-1m)+β(sn+1,tm)60limh0(Zn-1m-4Znm+3Zn+1m)=0.    (30)

From (29), we have

limh0(Z((n-1)h)-2Z(nh)+Z((n+1)h))=(ωr-Z0(1))e-β(1)(1-snε)(e-β(1)ρ+eβ(1)ρ-2),limh0(-3Z((n-1)h)+4Z(nh)-Z((n+1)h))=(ωr-Z0(1))e-β(1)(1-snε)(-3e-β(1)ρ-eβ(1)ρ+4),limh0(Z((n-1)h)-4Z(nh)+3Z((n+1)h))=(ωr-Z0(1))e-β(1)(1-snε)(e-β(1)ρ+3eβ(1)ρ-4),limh0(Z((n+1)h)-Z((n-1)h))=(ωr-Z0(1))e-β(1)(1-snε)(eβ(1)ρ-e-β(1)ρ).    (31)

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).    (32)

Therefore, the required scheme is taken as

𝔏Δt,hZnm=130(Qn-1m+28Qnm+Qn+1m),    (33)

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.    (34)

In the explicit form, it becomes

Rn-Zn-1m+Rn0Znm+Rn+Zn+1m=Hnm,    (35)

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).    (36)

3.3. Stability and uniform convergence analysis

Lemma 3.4. (Discrete maximum principle). Assume that Z0m0,ZNm0 and 𝔏Δt,hZnm0,n=1(1)N-1, then Znm0,n=0(1)N.

Proof: Assume that there is k ∈ {0, 1, 2, …, N}, such that Zkm=min0nNZnm<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, Znm0 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.    (37)

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)|Ch2d4Zm(sn)ds4,|dZm(sn-1)ds-(-Zn+1m+4Znm-3Zn-1m2h+hd2Zm(sn)ds2)|Ch2d3Zm(sn)ds3,|dZm(sn+1)ds-(3Zn+1m-4Znm+Zn-1m2h-hd2Zm(sn)ds2)|Ch2d3Zm(sn)ds3,|(dds-δs0)Zm(sn)|Ch2d3Zm(sn)ds3,|δs2Zm(sn)|Cd2Zm(sn)ds2,    (38)

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)ε)).    (39)

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))|,    (40)

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(ρ)-1C2ρ2ρ+1    (41)

and we have

ε[β(sn,tm)ρ2coth(β(1)ρ2)-1]ε(h/ε)2h/ε+1=h2h+ε.    (42)

From the bounds in (38), (40), and (42), we have

|𝔏Δt,h(Zm(sn)-Znm)|Ch2h+εd2Zm(sn)ds2+Ch2d3Zm(sn)ds3                                        +Cεh2d4Zm(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)ε))].    (43)

Obviously, ε−3 ≥ ε−2, then we draw

|𝔏Δt,h(Zm(sn)-Znm)|Ch2h+ε(1+ε-3exp(-μ(1-sn)ε))    (44)

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 ≤ nN − 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.    (45)

Proof: Substituting Lemma 3.7 into (39), we arrive at

|𝔏Δt,h(Zm(sn)-Znm)|Ch2h+ε.    (46)

Hence, the result leads

|Zm(sn)-Znm|Ch2h+ε.    (47)

Using the sup over all ε ∈ (0, 1], we get

supε(0,1]maxn|Zm(sn)-Znm|Ch,n=0(1)N.    (48)

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)).    (49)

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]

zt-ε2zs2+(2-s2)zs+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]

zt-ε2zs2+(2-s2)zs+(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.

TABLE 1
www.frontiersin.org

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

TABLE 2
www.frontiersin.org

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

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.

FIGURE 1
www.frontiersin.org

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

FIGURE 2
www.frontiersin.org

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

FIGURE 3
www.frontiersin.org

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

FIGURE 4
www.frontiersin.org

Figure 4. 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.

TABLE 3
www.frontiersin.org

Table 3. EN,M and rN,M for Example 4.1.

TABLE 4
www.frontiersin.org

Table 4. EN,M and rN,M for Example 4.2.

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.

Acknowledgments

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

1. Baker CT, Bocharov GA, Rihan FA. A Report on the Use of Delay Differential Equations in Numerical Modelling in the Biosciences. Citeseer (1999).

Google Scholar

2. Gopalsamy K. Stability and Oscillations in Delay Differential Equations of Population Dynamics. vol. 74. Springer Science and Business Media (1992).

Google Scholar

3. Bestehorn M, Grigorieva EV. Formation and propagation of localized states in extended systems. Ann Phys. (2004) 516:423–31. doi: 10.1002/andp.200451607-806

PubMed Abstract | CrossRef Full Text | Google Scholar

4. 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. doi: 10.1007/s12190-019-01313-7

CrossRef Full Text | Google Scholar

5. Tian H. Numerical Treatment of Singularly Perturbed Delay Differential Equations. Manchester, UK: The University of Manchester (2000).

PubMed Abstract | Google Scholar

6. 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. doi: 10.1016/j.cam.2006.05.032

CrossRef Full Text | Google Scholar

7. 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. doi: 10.1016/j.amc.2015.08.137

CrossRef Full Text | Google Scholar

8. 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. doi: 10.1016/j.cpc.2014.04.004

CrossRef Full Text | Google Scholar

9. 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. doi: 10.1002/num.22203

CrossRef Full Text | Google Scholar

10. 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. doi: 10.1007/s40314-020-01236-1

CrossRef Full Text | Google Scholar

11. 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. doi: 10.1016/j.camwa.2014.09.004

CrossRef Full Text | Google Scholar

12. Gowrisankar S, Natesan S. ε-Uniformly convergent numerical scheme for singularly perturbed delay parabolic partial differential equations. Int J Comput Math. (2017) 94:902–21. doi: 10.1080/00207160.2016.1154948

CrossRef Full Text | Google Scholar

13. 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. doi: 10.1080/00207160.2017.1290439

CrossRef Full Text | Google Scholar

14. Govindarao L, Mohapatra J. A second order numerical method for singularly perturbed delay parabolic partial differential equation. Eng Comput. (2018) 36:420–44. doi: 10.1108/EC-08-2018-0337

CrossRef Full Text | Google Scholar

15. 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. doi: 10.1002/num.22455

CrossRef Full Text | Google Scholar

16. 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. doi: 10.1002/num.22256

CrossRef Full Text | Google Scholar

17. 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. doi: 10.1016/j.cam.2011.01.047

CrossRef Full Text | Google Scholar

18. Kumar S, Kumar M. A second order uniformly convergent numerical scheme for parameterized singularly perturbed delay differential problems. Num Algor. (2017) 76:349–60. doi: 10.1007/s11075-016-0258-9

CrossRef Full Text | Google Scholar

19. Kumar S, Kumar M. Analysis of some numerical methods on layer adapted meshes for singularly perturbed quasilinear systems. Num Algor. (2016) 71:139–50. doi: 10.1007/s11075-015-9989-2

CrossRef Full Text | Google Scholar

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

Google Scholar

21. 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. doi: 10.1186/s13104-022-06202-0

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Woldaregay MM, Aniley WT, Duressa GF. Novel numerical scheme for singularly perturbed time delay convection-diffusion equation. Adv Math Phys. (2021) 2021:6641236. doi: 10.1155/2021/6641236

CrossRef Full Text | Google Scholar

23. 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. doi: 10.22124/jmm.2021.19608.1682

CrossRef Full Text | Google Scholar

24. 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. doi: 10.1007/s12190-021-01512-1

CrossRef Full Text | Google Scholar

25. 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. doi: 10.1080/00207160.2017.1284317

CrossRef Full Text | Google Scholar

26. 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. doi: 10.1016/j.rinam.2021.100174

CrossRef Full Text | Google Scholar

27. 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). doi: 10.1002/num.23039

CrossRef Full Text | Google Scholar

28. 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. doi: 10.1016/j.amc.2012.09.075

CrossRef Full Text | Google Scholar

29. 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. doi: 10.1007/s12190-018-1174-z

CrossRef Full Text | Google Scholar

30. 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. doi: 10.1007/s40995-021-01258-2

CrossRef Full Text | Google Scholar

31. Protter M, Weinberger H. Maximum Principles in Differential Equations. Englewood Cliffs, NJ: Prentice-Hall Inc. (1967).

Google Scholar

32. Kellogg RB, Tsan A. Analysis of some difference approximations for a singular perturbation problem without turning points. Math Comp. (1978) 32:1025–39.

Google Scholar

33. O'Malley RE. Singular Perturbation Methods for Ordinary Differential Equations. vol. 89. New York, NY: Springer (1991).

PubMed Abstract | Google Scholar

Keywords: singularly perturbed, delay differential equation, exponentially fitted finite difference, non-symmetric finite difference, uniform convergence

Citation: Tesfaye SK, Duressa GF, Woldaregay MM and Dinka TG (2023) Fitted computational method for singularly perturbed convection-diffusion equation with time delay. Front. Appl. Math. Stat. 9:1244490. doi: 10.3389/fams.2023.1244490

Received: 22 June 2023; Accepted: 07 August 2023;
Published: 31 August 2023.

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

Copyright © 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.

*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

Download