# Uniformly Convergent Hybrid Numerical Method for Singularly Perturbed Delay Convection-Diffusion Problems

Authors: Mesfin Mekuria Woldaregay; Gemechis File Duressa
Journal: International Journal of Differential Equations
 (2021)

---

## Abstract

This paper deals with numerical treatment of nonstationary singularly perturbed delay convection-diffusion problems. The solution of the considered problem exhibits boundary layer on the right side of the spatial domain. To approximate the term with the delay, Taylor’s series approximation is used. The resulting time-dependent singularly perturbed convection-diffusion problems are solved using Crank-Nicolson method for temporal discretization and hybrid method for spatial discretization. The hybrid method is designed using mid-point upwind in regular region with central finite difference in boundary layer region on piecewise uniform Shishkin mesh. Numerical examples are used to validate the theoretical findings and analysis of the proposed scheme. The present method gives accurate and nonoscillatory solutions in regular and boundary layer regions of the solution domain. The stability and the uniform convergence of the scheme are proved. The scheme converges uniformly with almost second-order rate of convergence.


---

## Body

## 1. Introduction

Singularly perturbed differential equations are differential equations in which the highest-order derivative term is multiplied by a small perturbation parameterε. Singularly perturbed problems appear in several research areas of applied mathematics [1–3], for example, in water quality problem in river networks, in simulation of oil extraction from underground reservoirs, in convective heat transport problem with large Peclet numbers, in atmospheric pollution, and in fluid flow at high Reynolds number.One of the common examples of singularly perturbed problems is the Navier-Stokes equation of fluid dynamics:(1)∂∂xu2+p+∂∂yuv+1Re∂2u∂x2+∂2u∂y2,with suitable initial and boundary conditions. In (1), p is the pressure and u and v are the velocity components in the x and y directions, respectively. Parameter Re is known as Reynolds number which is proportional to the length scale, velocity scale, and density and inversely proportional to the kinematic viscosity of the fluid [4].Singularly Perturbed Delay Differential Equations (SPDDEs) are differential equations in which the highest-order derivative term is multiplied by small perturbation parameterε and involves at least one delay term. SPDDEs model process for which the evaluation does not only depend on the current state of the system but also includes the past history. So, it is called problem with memory.In general, when the perturbation parameter tends to zero, the smoothness of the solution of SPDDEs deteriorates and it forms boundary layer [5]. As a result, standard numerical methods like FDM, FEM, and collocation method give accurate results only when mesh size is less than the perturbation parameter (i.e., h≤ε). But for the case where h>ε (or ε is very small), it leads to oscillations in the computed solution while using the standard numerical methods [5]. To overcome these oscillations, while using the standard numerical methods, a large number of mesh points are required, which is not practical due to round-off error, computer processing ability, and computer memory issue.To handle this drawback, recently, different authors have developed numerical schemes for solving singularly perturbed parabolic problems having delay on the spatial variable. Chakravarthy and Kumar [6] used Taylor’s approximation for the delay and discretize the domain using uniform mesh in the temporal direction and adaptive mesh is generated using the concept of entropy function for the spatial direction. The method is based on central finite difference scheme. Kumar and Kadalbajoo [7] first treated the delay using Taylor’s series approximation and developed numerical scheme using upwind for time derivative with B-spline collocation for spatial derivatives approximation on Shishkin mesh.In [8], Kumar used Taylor’s series approximation to tackle the shift terms and semidiscretized the time direction using the Crank-Nicolson method. The resulting set of ODEs is discretized using a mid-point upwind finite difference scheme on Shishkin mesh. Rai and Sharma [9] developed numerical scheme for turning layer problem. They first treated the delays using Taylor’s series approximation and used upwind FDM on Shishkin mesh for spatial derivative approximation. The scheme converges with rate of convergence of almost one. Ramesh and Kadalbajo [5] first treated the delays using Taylor’s series approximation and developed numerical scheme using upwind for time derivative discretization with upwind and mid-point upwind for the spatial derivatives discretization on Shishkin mesh. Rao and Chakravarthy [10] developed numerical scheme using upwind for time derivative with fitted operator FDM for spatial derivatives discretization. Kumar et al. [11] developed complete flux-finite volume method for treating the considered problem. Woldaregay and Duressa [12] developed uniformly convergent scheme using exponentially fitted operator finite difference method. They proved the uniform convergence of the scheme with linear order of convergence.Numerical treatment of singularly perturbed parabolic delay differential equations still needs improvement; developing higher-order uniformly convergent numerical scheme is an active research [13]. In this paper, we formulate almost second-order uniformly convergent numerical scheme for treating the considered problem. The proposed scheme consists of Crank-Nicolson method for temporal discretization with hybrid FDM (mid-point upwind in outer layer and central difference in boundary layer) on Shishkin mesh for spatial discretization.Notation 1.
The symbolsNandMare denoted for the number of mesh intervals in spatial and temporal discretization, respectively. SymbolC(in some cases indexed) denotes positive constant independent ofεandN. The norm.represents the maximum norm.


## 2. Statement of the Problem

We consider a class of singularly perturbed parabolic delay convection-diffusion problems containing delay on the reaction term of the spatial variable of the form(2)∂u∂t−ε2∂2u∂x2+ax∂u∂x+αxux−δ,t+βxux,t=fx,t,x,t∈D=Ω×Λ,ux,0=u0x,x∈D0=x,0:x∈Ω¯,ux,t=ϕx,t,x,t∈DL=x,t:x,t∈−δ,0×Λ,u1,t=ψ1,t,x,t∈DR=x,t:x=1,t∈Λ,where D=Ω×Λ=0,1×0,T for some fixed number T>0, with smooth boundary ∂D=D¯\D, and ε,0<ε≪1 is singular perturbation parameter and δ is delay parameter satisfying δ<ε. The functions ax,αx,βx,fx,t,u0x,ϕx,t, and ψ1,t are assumed to be sufficiently smooth, bounded, and independent of parameter ε to guarantee the existence of unique solution. The coefficient functions αx and βx are assumed to satisfy(3)αx+βx≥α∗+β∗>0,∀x∈Ω¯,where α∗ and β∗ are the lower bounds for αx and βx, respectively.


### 2.1. Properties of the Continuous Solution and Bounds

In case ofδ<ε, using Taylor’s series approximation for the terms containing delay is appropriate [14]. So, we approximate(4)ux−δ,t≈ux,t−δuxx,t+δ22uxxx,t.Using (3) into (2) gives(5)∂u∂t−ε2−δ22αx∂2u∂x2+ax−δαx∂u∂x+αx+βxux,t=fx,t,x,t∈D,ux,0=u0x,x∈Ω¯,u0,t=ϕ0,t,t∈Λ¯,u1,t=ψ1,t,t∈Λ¯.For small values ofδ, (2) and (5) are asymptotically equivalent, since the difference between the two is Oδ3. We have 0<ε2−δ2/2αx≤ε2−δ2/2α∗ since α∗ is lower bound of αx.The problem in (5) exhibits regular boundary layer and the position of the boundary layer depends on the following condition: If ax−δαx<0, left boundary layer exists and, for ax−δαx<0, right boundary layer exists. In case ax−δαx changes sign, inside layer will exist [15, 16]. For the case where δ=0, the problem reduces to singularly perturbed parabolic PDEs, in which for small ε it exhibits boundary layer. The layer is maintained for δ≠0 but it is sufficiently small. We assume also that ax−δαx≥a∗−δα∗>0 implies occurrence of right boundary layer.Lemma 1.
Let0<ε≪1; there exists a constantCindependent ofεsuch that the solutionux,tsatisfies(6)ux,t−u0x≤Ct,ux,t−ϕ0,t≤C1−x,x,t∈D¯.Remark 1.
There does not exist a constantCindependent ofεsuch thatux,t−ψ1,t≤Cx, since the boundary layer is on the right side of the domain.The problem obtained by settingε2−δ2/2α∗=0 in (5) is called reduced problem and is given as(7)∂u0∂t+ax−δαx∂u0∂x+αx+βxu0x,t=fx,t,x,t∈D,u0x,0=u0x,x∈Ω¯,u00,t=ϕ0,t,u01,t≠ψ1,t,x∈Λ¯.For small values ofε2−δ2/2α∗, the solution ux,t of (5) is very close to the solution u0x,t of (5).To show the bounds of the solutionux,t of (5), we assume without loss of generality that u0x=0 [17]. Since u0x is sufficiently smooth and using the property of norm, we prove the following lemma.Lemma 2.
The solutionux,tof (5) is bounded as(8)ux,t≤C,x,t∈D¯.Proof.
From the inequalityux,t−u0x≤Ct, we have ux,t−u0x≤ux,t−u0x≤Ct which implies that ux,t≤Ct+u0x,∀x,t∈D¯; since t∈0,T and u0x is bounded, this implies that there exists a constant C1 such that ux,t≤C1.The differential operatorL is denoted for the differential equations in (4) as(9)Lux,t=∂u∂t−ε2−δ22α∗∂2u∂x2+ax−δαx∂u∂x+αx+βxux,t.Lemma 3 (the maximum principle).
Letz be a sufficiently smooth function defined on D which satisfies zx,t≥0,x,t∈∂D. Then Lzx,t≥0,x,t∈D implies that zx,t≥0,∀x,t∈D¯.Proof.
Suppose that there existsx∗,t∗∈D¯ such that zx∗,t∗=minx,t∈D¯zx,t<0. It is clear that x∗,t∗∉∂D; that is, x∗,t∗∈D. Since zx∗,t∗=minx,t∈D¯zx,t which implies that zxx∗,t∗=0, ztx∗,t∗=0 and zxxx∗,t∗≥0. Giving that Lzx∗,t∗<0 which is contradiction to the assumption made Lzx∗,t∗≥0,x,t∈D. Therefore, zx,t≥0,∀x,t∈D¯.Lemma 4.
Suppose that Lemmas2 and 3 hold. Then, the bound on the derivatives of ux,t with respect to t is given by(10)∂iux,t∂ti≤C,x,t∈D,i=0,1,2.Proof.
See [18].Lemma 5 (stability estimate).
The solutionux,t of the continuous problem in (5) is bounded as(11)ux,t≤fα∗+β∗+maxu0x,ϕ0,t,ψ1,t.Proof.
Define two barrier functions asϑ±x,t=f/α∗+β∗+maxu0x,ϕ0,t,ψ1,t±ux,t. At the initial stage, we have(12)ϑ±x,0=fα∗+β∗+maxu0x,ϕ0,t,ψ1,t±u0x≥0.
On the boundaries, we have(13)ϑ±0,t=fα∗+β∗+maxu0x,ϕ0,t,|ψ1,t±u0,t≥0,ϑ±1,t=fα∗+β∗+maxu0x,ϕ0,t,ψ1,t±u1,t≥0.
On domainD, we have(14)Lϑ±x,t=ϑt±x,t−ε2−δ22α∗ϑxx±x,t+ax−δαxϑx±x,t+αx+βxϑ±x,t=0±utx,t−ε2−δ22α∗0±uxxx,t+ax−δαx0±uxx,t+αx+βxfα∗+β∗+maxu0x,ϕ0,t|,ψ1,t±ux,t=αx+βxfα∗+β∗+maxu0x,ϕ0,t,ψ1,t±fx,t≥0,sinceαx+βx≥α∗+β∗>0.
Sinceϑ±x,0≥0,ϑ±0,t≥0,ϑ±1,t≥0 and Lϑ±x,t≥0, by the hypothesis of the maximum principle, we obtain ϑ±x,t≥0,∀x,t∈D¯, which gives the required bound.Lemma 6.
The derivative of the solutionux,t of the problem in (5) satisfies the bound(15)∂k∂lux,t∂xk∂tl≤C1+ε2−δ22α∗−kexp−α∗+β∗ε2−δ2/2α∗1−x,x,t∈D¯,0≤k≤4,0≤l≤2.Proof.
See [18].


## 2.1. Properties of the Continuous Solution and Bounds

In case ofδ<ε, using Taylor’s series approximation for the terms containing delay is appropriate [14]. So, we approximate(4)ux−δ,t≈ux,t−δuxx,t+δ22uxxx,t.Using (3) into (2) gives(5)∂u∂t−ε2−δ22αx∂2u∂x2+ax−δαx∂u∂x+αx+βxux,t=fx,t,x,t∈D,ux,0=u0x,x∈Ω¯,u0,t=ϕ0,t,t∈Λ¯,u1,t=ψ1,t,t∈Λ¯.For small values ofδ, (2) and (5) are asymptotically equivalent, since the difference between the two is Oδ3. We have 0<ε2−δ2/2αx≤ε2−δ2/2α∗ since α∗ is lower bound of αx.The problem in (5) exhibits regular boundary layer and the position of the boundary layer depends on the following condition: If ax−δαx<0, left boundary layer exists and, for ax−δαx<0, right boundary layer exists. In case ax−δαx changes sign, inside layer will exist [15, 16]. For the case where δ=0, the problem reduces to singularly perturbed parabolic PDEs, in which for small ε it exhibits boundary layer. The layer is maintained for δ≠0 but it is sufficiently small. We assume also that ax−δαx≥a∗−δα∗>0 implies occurrence of right boundary layer.Lemma 1.
Let0<ε≪1; there exists a constantCindependent ofεsuch that the solutionux,tsatisfies(6)ux,t−u0x≤Ct,ux,t−ϕ0,t≤C1−x,x,t∈D¯.Remark 1.
There does not exist a constantCindependent ofεsuch thatux,t−ψ1,t≤Cx, since the boundary layer is on the right side of the domain.The problem obtained by settingε2−δ2/2α∗=0 in (5) is called reduced problem and is given as(7)∂u0∂t+ax−δαx∂u0∂x+αx+βxu0x,t=fx,t,x,t∈D,u0x,0=u0x,x∈Ω¯,u00,t=ϕ0,t,u01,t≠ψ1,t,x∈Λ¯.For small values ofε2−δ2/2α∗, the solution ux,t of (5) is very close to the solution u0x,t of (5).To show the bounds of the solutionux,t of (5), we assume without loss of generality that u0x=0 [17]. Since u0x is sufficiently smooth and using the property of norm, we prove the following lemma.Lemma 2.
The solutionux,tof (5) is bounded as(8)ux,t≤C,x,t∈D¯.Proof.
From the inequalityux,t−u0x≤Ct, we have ux,t−u0x≤ux,t−u0x≤Ct which implies that ux,t≤Ct+u0x,∀x,t∈D¯; since t∈0,T and u0x is bounded, this implies that there exists a constant C1 such that ux,t≤C1.The differential operatorL is denoted for the differential equations in (4) as(9)Lux,t=∂u∂t−ε2−δ22α∗∂2u∂x2+ax−δαx∂u∂x+αx+βxux,t.Lemma 3 (the maximum principle).
Letz be a sufficiently smooth function defined on D which satisfies zx,t≥0,x,t∈∂D. Then Lzx,t≥0,x,t∈D implies that zx,t≥0,∀x,t∈D¯.Proof.
Suppose that there existsx∗,t∗∈D¯ such that zx∗,t∗=minx,t∈D¯zx,t<0. It is clear that x∗,t∗∉∂D; that is, x∗,t∗∈D. Since zx∗,t∗=minx,t∈D¯zx,t which implies that zxx∗,t∗=0, ztx∗,t∗=0 and zxxx∗,t∗≥0. Giving that Lzx∗,t∗<0 which is contradiction to the assumption made Lzx∗,t∗≥0,x,t∈D. Therefore, zx,t≥0,∀x,t∈D¯.Lemma 4.
Suppose that Lemmas2 and 3 hold. Then, the bound on the derivatives of ux,t with respect to t is given by(10)∂iux,t∂ti≤C,x,t∈D,i=0,1,2.Proof.
See [18].Lemma 5 (stability estimate).
The solutionux,t of the continuous problem in (5) is bounded as(11)ux,t≤fα∗+β∗+maxu0x,ϕ0,t,ψ1,t.Proof.
Define two barrier functions asϑ±x,t=f/α∗+β∗+maxu0x,ϕ0,t,ψ1,t±ux,t. At the initial stage, we have(12)ϑ±x,0=fα∗+β∗+maxu0x,ϕ0,t,ψ1,t±u0x≥0.
On the boundaries, we have(13)ϑ±0,t=fα∗+β∗+maxu0x,ϕ0,t,|ψ1,t±u0,t≥0,ϑ±1,t=fα∗+β∗+maxu0x,ϕ0,t,ψ1,t±u1,t≥0.
On domainD, we have(14)Lϑ±x,t=ϑt±x,t−ε2−δ22α∗ϑxx±x,t+ax−δαxϑx±x,t+αx+βxϑ±x,t=0±utx,t−ε2−δ22α∗0±uxxx,t+ax−δαx0±uxx,t+αx+βxfα∗+β∗+maxu0x,ϕ0,t|,ψ1,t±ux,t=αx+βxfα∗+β∗+maxu0x,ϕ0,t,ψ1,t±fx,t≥0,sinceαx+βx≥α∗+β∗>0.
Sinceϑ±x,0≥0,ϑ±0,t≥0,ϑ±1,t≥0 and Lϑ±x,t≥0, by the hypothesis of the maximum principle, we obtain ϑ±x,t≥0,∀x,t∈D¯, which gives the required bound.Lemma 6.
The derivative of the solutionux,t of the problem in (5) satisfies the bound(15)∂k∂lux,t∂xk∂tl≤C1+ε2−δ22α∗−kexp−α∗+β∗ε2−δ2/2α∗1−x,x,t∈D¯,0≤k≤4,0≤l≤2.Proof.
See [18].


## 3. Numerical Scheme Formulation

### 3.1. Temporal Semidiscretization

Subdivide the time domain0,T into M intervals as t1=0,tj=jΔt,j=1,2,…,M, where Δt=T/M−1. Let Uj+1x denote the approximation of ux,tj+1 at the j+1th time level discretization. Using Crank-Nicolson method for semidiscretizing the problem in (4), we obtain a system of boundary value problems:(16)LΔtUj+1x=gj+1x,x∈Ω,j=1,2,…,M−1,U1x=u0x,x∈Ω¯,Uj+10=ϕ0,tj+1,Uj+11=ψ1,tj+1,where(17)LΔtUj+1x=−ε2−δ2/2α∗2d2dx2Uj+1x+ax−δαx2ddxUj+1x+dxUj+1x,gj+1x=ε2−δ2/2α∗2d2dx2Ujx−ax−δαx2ddxUjx−exUjx+fx,tj+1+fx,tj2,for dx=αx+βx/2+1/Δt and ex=αx+βx/2−1/Δt.Next, we state the discrete maximum principle for semidiscretized problem as follows.Lemma 7.
(semidiscrete maximum principle). For eachj=1,2,…,M−1, let zj+1 be sufficiently smooth function on domain Ω¯. If zj+10≥0,zj+11≥0 and LΔtzj+1x≥0,x∈Ω, then zj+1x≥0,∀x∈Ω¯.Proof.
Suppose thatx∗ is such that zj+1x∗=minx∈Ω¯zj+1x<0. From assumption, it is clear that x∗∉0,1 implies that x∗∈0,1. Applying the property of extrema values in calculus gives d/dxzj+1x∗=0, and d2/dx2zj+1x∗≥0, given that LΔtzj+1x∗<0 which is a contradiction to LΔtzj+1x∗≥0,x∈Ω. Therefore, we conclude that zj+1x≥0,∀x∈Ω¯. Hence, operator LΔt satisfies the semidiscrete maximum principle.The local error for each time step is denoted byej+1 and given by ej+1:=ux,tj+1−Uj+1x.Lemma 8 (local error estimate).
Considering the result in Lemma4, the local error in the temporal semidiscretization is given by(18)ej+1≤C2Δt3,j=1,2,…,M−1.Proof.
Using Taylor’s series expansion abouttj+1/2, we obtain(19)ux,tj+1=ux,tj+1/2+Δt2∂∂tux,tj+1/2+Δt28∂2∂t2ux,tj+1/2+OΔt3,ux,tj=ux,tj+1/2−Δt2∂∂tux,tj+1/2+Δt28∂2∂t2ux,tj+1/2+OΔt3.
Substituting (19) into (5), we obtain(20)ux,tj+1−ux,tjΔt=∂∂tux,tj+1/2+OΔt2=ε2−δ22α∗∂2∂x2ux,tj+1/2−ax−δαx∂∂xux,tj+1/2−αx+βxux,tj+1/2+fx+tj+1/2+OΔt2,where ux,tj+1/2=ux,tj+1+ux,tj/2+OΔt2. By simplifying, (20) can be rewritten as(21)LΔtux,tj+1=gx,tj+1+OΔt3,and, from (8), we have(22)LΔtUj+1x=gj+1x,x∈Ω.
From (21) and (22), we obtain(23)LΔtux,tj+1−Uj+1x=LΔtej+1x≤C2Δt3,with the boundary conditions u0,tj+1−Uj+10=ej+10=0 and u1,tj+1−Uj+11=ej+11=0. Hence, applying the maximum principle gives ej+1≤C2Δt3.Next, we need to show the bound for the global error of the temporal discretization. LetEj+1 denote the global error estimate up to the j+1th time step.Lemma 9 (Global error estimate).
The global error estimate up totj+1 time step is given by(24)Ej+1≤C2Δt2,j=1,2,…,M−1.Proof.
Using the local error estimate up to thej+1th time step given in Lemma 8, we obtain the global error estimate at the j+1th time step as(25)Ej+1=∑l=1j+1el≤e1+e2+⋯+ej+1≤C2j+1Δt3=j+1ΔtC2Δt2≤C2TΔt2, sincej+1Δt≤T=C1Δt2,C2T≔C1,where Ci is a constant independent of ε2−δ2/2α∗ and Δt.Lemma 10.
Forj=1,2,…,M−1, let Uj+1x be a solution of (16); then the solution Uj+1x satisfies the following bound:(26)Uj+1x≤LΔtUj+1d∗+maxϕ0,tj+1,ψ1,tj+1,where dx≥d∗>0.Proof.
Construct barrier functionsπj+1±x=LΔtUj+1/d∗+maxϕ0,tj+1,ψ1,tj+1±Uj+1x. We need to show that if πj+1±0≥0,πj+1±1≥0 and LΔtπj+1±x≥0,x∈Ω, this implies that πj+1±x≥0,∀x∈Ω¯. We easily obtain that πj+1±0≥0 and πj+1±1≥0, and(27)LΔtπj+1±x=−ε2−δ2/2α∗2d2dx2πj+1±x+ax−δαx2ddxπj+1±x+dxπj+1±x=−ε2−δ2/2α∗2d2dx2Uj+1x+ax−δαx2ddxUj+1x+dxπj+1±x=dxLΔtUj+1d∗+maxϕ0,tj+1,ψ1,tj+1±LΔtUj+1x≥0,sincedx≥d∗.
Using the semidiscrete maximum principle, we conclude thatπj+1±x≥0,∀x∈Ω¯. Hence, the required bound is satisfied.


### 3.2. Decomposition of Solution

Since the bounds on the derivatives of the solution are not sharp enough for the proof of uniform convergence, we need to derive more stronger error bounds. This can be achieved by decomposing the solution into regular and singular component and doing error bound separately for each. Let(28)Uj+1x=Vj+1x+Wj+1x,x∈Ω¯.The regular component is the solution to the nonhomogeneous problem:(29)LΔtVj+1x=gj+1x,x∈Ω,Vj+10=Uj+10,V0x=U0x,x∈Ω¯,and the singular component is the solution to the homogeneous problem:(30)LΔtWj+1x=0,x∈Ω,Wj+10=0,Wj+11=Uj+11−Vj+11,W0x=0,x∈Ω.Lemma 11.
Forj=1,2,…,M−1, the regular component of the solution and its derivatives satisfies the bound(31)dkdxkVj+1x≤C1+ε2−δ22α∗2−k,x∈Ω,k=0,1,2.
In general, the derivatives are bounded as(32)dkdxkVj+1x≤C,x∈Ω,k=0,1,2,3.Proof.
See [18].Lemma 12.
Forj=1,2,…,M−1, the singular component of the solution and its derivatives satisfies the bound:(33)dkdxkWj+1x≤Cε2−δ22α∗−kexp−a∗−δα∗ε2−δ2/2α∗1−x,x∈Ω,k=0,1,2,3.Proof.
See [18].The next lemma gives a bound for the derivatives of solution.Lemma 13.
The solutions of the boundary value problems in (16) satisfy the bound(34)dkdxkUj+1x≤C1+ε2−δ22α∗−kexp−a∗−δα∗ε2−δ2/2α∗1−x,x∈Ω¯,k=0,1,2,3.Proof.
The results from Lemmas11 and 12 give the required bound.


### 3.3. Spatial Discretization and the Hybrid FDM

Letxii=0N be the discretized domain of the spatial domain, where N be the number of grid points and must be even positive integer. For each i≥1, we define hi=xi−xi−1. Since the considered problem exhibits right boundary layer, we set a mesh transition parameter τ which divides the domain 0,1 into outer layer region Ω1=0,1−τ and boundary layer region Ω2=1−τ,1, where(35)τ=min12,σε2−δ2/2α∗d∗lnN,where N denotes the number of mesh points in spatial discretization and σ denotes a constant that represents the order of the scheme. We set a uniform mesh N/2 in Ω1 with mesh spacing H=21−τ/N and similarly a uniform mesh N/2 is placed in Ω2 with the mesh spacing h=2τ/N. So, the mesh point xi is given as(36)xi=21−τNi,i=0,1,2,…,N2,1−τ+2τNi−N/2,i=N/2,…,N.In a similar fashion, one can generate a Shishkin mesh for the left boundary layer problem.For the mesh functionsUj+1xi at the grid points xi, the backward and forward difference operators are given as(37)Dx−Uj+1xi=Ui,j+1−Ui−1,j+1hi,Dx+Uj+1xi=Ui+1,j+1−Ui,j+1hi+1,respectively. The central difference operator Dx0 is(38)Dx0Uj+1xi=Ui+1,j+1−Ui−1,j+1hi+hi+1.Similarly, the second-order derivativeDx+Dx− is given by(39)Dx−Dx+Uj+1xi=2hi+hi+1Ui+1,j+1−Ui,j+1hi+1−Ui,j+1−Ui−1,j+1hi.Remark 2.
The approximationDx0Ui,j+1is second-order consistent on any mesh, just as on equidistant meshes, but this is not the case forDx+Dx−Ui,j+1: a termhi+1−hid3/dx3Uj+1xi, which arises in the consistency error analysis, is only of first order on arbitrary mesh. For more details, the reader is referred to [3, 19]. Whenεis small, the central difference approximationDx0Uj+1xitod/dxUj+1xon equidistant grid leads to nonphysical oscillations of the computed solution which is due to the loss in stability, unless the grid size is very small.To overcome the challenge indicated in Remark2, hybrid numerical scheme containing central difference in the boundary layer region and mid-point upwind scheme on outer region is developed in order to get better accuracy and rate of convergence. Using the operators above, the discrete hybrid scheme can be formulated as(40)LΔt,NUi,j+1=gi,j+1,i=1,2,…,N−1,U0,j+1=ϕ0,tj+1,UN,j+1=ψ1,tj+1,where(41)LΔt,NUi,j+1=−ε2−δ2/2α∗2Dx+Dx−Ui,j+1+axi−1/2−δαxi−1/22Dx−Ui,j+1+dxi−1/2Ui,j+1,i=1,2,…,N/2,−ε2−δ2/2α∗2Dx+Dx−Ui,j+1+axi−δαxi2Dx0Ui,j+1+dxiUi,j+1,i=N/2+1,…,N−1,gi,j+1=ε2−δ2/2α∗2Dx+Dx−Ui,j−axi−1/2−δαxi−1/22Dx−Ui,j−exi−1/2Ui,j+12fi−1/2,j+fi−1/2,j+1,i=1,2,…,N/2,ε2−δ2/2α∗2Dx+Dx−Ui,j+1−axi−δαxi2Dx0Ui,j+1−exiUi,j+1+12fi,j+fi,j+1,i=N/2+1,…,N−1,and the notation dxi−1/2 stands for dxi+dxi−1/2 and similarly for others.


### 3.4. Properties of the Discrete Scheme and Uniform Convergence

Assume thatτ=min1/2,σε2−δ2/2α∗/d∗lnN. We assign each of the subintervals 0,1−τ and 1−τ,1 with N/2 number of equidistant grid points. Let H be the mesh width in the subinterval 0,1−τ and let h be the mesh width in 1−τ,1. These mesh widths satisfy(42)N−1≤H≤2N−1,h=4d∗ε2−δ22α∗N−1lnN.To attain stability condition while using central difference scheme on the singular component, we need the following condition [19]:(43)NlnN≥2maxdxd∗.Lemma 14 (discrete maximum principle).
The discrete schemeLΔt,NUi,j+1=gi,j+1,i=1,2,…,N−1and for givenU0,j+1andUN,j+1for fixedjhas a solution. IfLΔt,NUi,j+1≤LΔt,NBi,j+1fori=1,2,…,N−1and ifU0,j+1≤B0,j+1andUN,j+1≤BN,j+1, thenUi,j+1≤Bi,j+1,∀i=0,1,2,…,N.Proof.
The matrix associated with operatorLΔt,N is of size N+1×N+1 and satisfies the property of M-matrix with the given condition in (43). See the detailed proof in [20].Lemma 15.
LetZi,j+1=1+xi,i=0,1,…,N. Then, there exists a positive constant C such that LΔt,NZi,j+1≥C,i=1,2,…,N−1.Proof.
Using the discrete operator on the mesh functionZi,j+1, for 1≤i≤N/2, we have(44)LΔt,NZi,j+1=−ε2−δ2/2α∗2Dx+Dx−Zi,j+1+axi−1/2−δαxi−1/22Dx−Zi,j+1+dxi−1/2Zi,j+1,=−ε2−δ2/2α∗2Dx+Dx−1+xi+axi−1/2−δαxi−1/22Dx−1+xi+dxi−1/21+xi,=0+axi−1/2−δαxi−1/22Dx−1+xi+dxi−1/21+xi,≥C,sincedxi−1/2is bounded.For the case whereN/2+1≤i≤N−1, we have(45)LΔt,NZi,j+1=−ε2−δ2/2α∗2Dx+Dx−Zi,j+1+axi−δαxi2Dx0Zi,j+1+dxiZi,j+1,=−ε2−δ2/2α∗2Dx+Dx−1+xi+axi−δαxi2Dx01+xi+dxi1+xi,=0+axi−δαxi2Dx01+xi+dxi1+xi,≥C,sincedxiis bounded.Lemma 16.
Fori=0,1,…,Nand fixedj, define a mesh function(46)Si,j+1=∏k=1i1+a∗−δα∗hkε2−δ2/2α∗,with the convention that if i=0, then S0=1. Then, for i=1,2,…,N−1, we have(47)LΔt,NSi,j+1≥Cmaxε2−δ2/2α∗,hiSi,j+1,for some constant C.Proof.
SinceSi,j+1−Si−1,j+1/hi=a∗−δα∗/ε2−δ2/2α∗Si−1,j+1, we have Si−1,j+1=ε2−δ2/2α∗/ε2−δ2/2α∗+a∗−δα∗hiSi,j+1. Now, putting these on the discrete operator, we have(48)LΔt,NSi,j+1=−ε2−δ2/2α∗hi+hi+1Si+1,j+1−Si,j+1hi+1−Si,j+1−Si−1,j+1hi+axi−1/2−δαxi−1/22Si,j+1−Si−1,j+1hi+dxi−1/2Si,j+1,i=1,2,…,N/2,−ε2−δ2/2α∗hi+hi+1Si+1,j+1−Si,j+1hi+1−Si,j+1−Si−1,j+1hi+axi−δαxi2Si+1,j+1−Si−1,j+1hi+hi+1+dxiSi,j+1,i=N/2+1,…,N−1.For i=1,2,…,N/2, we have(49)LΔt,NSi,j+1=−ε2−δ2/2α∗hi+hi+1a∗−δα∗Si,j+1−Si−1,j+1ε2−δ2/2α∗+axi−1/2−δαxi−1/22a∗−δα∗Si−1,j+1ε2−δ2/2α∗+dxi−1/2Si,j+1=a∗−δα∗2ε2−δ2/2α∗Si−1,j+1axi−1/2−δαxi−1/2−2a∗−δα∗hihi+hi+1+dxi−1/2Si,j+1=a∗−δα∗2ε2−δ2/2α∗+a∗−δα∗hiSi,j+1axi−1/2−δαxi−1/2−2a∗−δα∗hihi+hi+1+dxi−1/2ε2−δ2/2α∗+a∗−δα∗hia∗−δα∗≥Si,j+1Cmaxε2−δ2/2α∗,hi.
Hence, the required result holds. For the case wherei=N/2+1,…,N−1,(50)LΔt,NSi,j+1=−ε2−δ2/2α∗hi+hi+1a∗−δα∗Si,j+1−Si−1,j+1ε2−δ2/2α∗+axi−δαxi2×1+a∗−δα∗hihi+1ε2−δ2/2α∗hi+hi+1a∗−δα∗Si−1,j+1ε2−δ2/2α∗+dxiSi,j+1=−a∗−δα∗2hiε2−δ2/2α∗2hi+hi+1Si−1,j+1+axi−δαxia∗−δα∗2ε2−δ2/2α∗×1+a∗−δα∗hihi+1ε2−δ2/2α∗hi+hi+1Si−1,j+1+dxiSi,j+1=−a∗−δα∗2hiε2−δ2/2α∗hi+hi+1ε2−δ2/2α∗+a∗−δα∗hiSi,j+1+axi−δαxia∗−δα∗2ε2−δ2/2α∗+a∗−δα∗hi×1+a∗−δα∗hihi+1ε2−δ2/2α∗hi+hi+1Si,j+1+dxiSi,j+1=a∗−δα∗Si,j+1ε2−δ2/2α∗+a∗−δα∗hi−a∗−δα∗hiε2−δ2/2α∗hi+hi+1+axi−δαxi1+a∗−δα∗hihi+1ε2−δ2/2α∗hi+hi+1+dxiε2−δ2/2α∗+a∗−δα∗hia∗−δα∗≥Si,j+1Cmaxε2−δ2/2α∗,hi.Lemma 17.
Letqxbe any smooth function defined on0,1. Fori=1,2,…,N−1, define truncation errorTEtoqxatxiwhich is to beLΔtqxi−LΔt,Nqi; then there exists a constantCsuch that(51)TE=LΔtqxi−LΔt,Nqi≤Cε2−δ22α∗∫xi−1xi+1d3qξdx3dξ+Chi∫xi−1xid3qξdx3dξ,for ξ∈xi−1,xi+1.Proof.
See Lemma 4.1 in [16].Lemma 18.
For eachi=1,2,…,N−1, we have(52)exp−a∗−δα∗1−xiε2−δ2/2α∗≤∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1.Proof.
For eachi, we have(53)exp−a∗−δα∗hiε2−δ2/2α∗=ea∗−δα∗hi/ε2−δ2/2α∗−1≤1+a∗−δα∗hiε2−δ2/2α∗−1.
Now multiplying fork=i+1,…,N, we obtain(54)∏k=i+1Nea∗−δα∗hi/ε2−δ2/2α∗−1≤∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1⇒e−a∗−δα∗1−xiε2−δ2/2α∗≤∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1.Lemma 19.
For Shishkin mesh defined above, there exists a constantC such that(55)∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1≤CN−41−i/N,i=N/2,…,N.Proof.
LetN/2≤i≤N by Lemma 1 in [20].(56)1+a∗−δα∗hkε2−δ2/2α∗−1=1−a∗−δα∗hiε2−δ2/2α∗+a∗−δα∗hi≤e−a∗−δα∗hi/ε2−δ2/2α∗+a∗−δα∗hi.
Therefore, we have(57)∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1≤∏k=i+1Ne−a∗−δα∗hi/ε2−δ2/2α∗+a∗−δα∗hi=e−a∗−δα∗1−xi/ε2−δ2/2α∗+a∗−δα∗hi=e−4N−iN−1lnN/1+4N−1lnN,using,h=4d∗ε2−δ22α∗N−1lnN=N−4N−iN−1/1+4N−1lnN=N−41−i/NN161−i/NN−1lnN/1+4N−1lnN≤CN−41−i/N.
SinceN161−i/NN−1lnN/1+4N−1lnN is bounded for N≥2, the desired result is obtained.


### 3.5. Decomposition of Numerical Solution

Here, we decompose the numerical solutionUi,j+1 into regular and singular components in a similar fashion to the continuous case:(58)Ui,j+1=Vi,j+1+Wi,j+1,where the regular part satisfies the nonhomogeneous equation:(59)LΔt,NVi,j+1=gi,j+1,Vj+10=Uj+10,and the singular component Wi,j+1 satisfies the homogeneous equation:(60)LΔt,NWi,j+1=0,Wj+10=0,Wj+11=Uj+11−Vj+11,and now the error in the numerical solution can also be decomposed as(61)Ui,j+1−uj+1xi=Vi,j+1−Vj+1xi+Wi,j+1−Wj+1xi.We estimate the error in the regular and singular solution separately.Theorem 1 (error in the regular component).
The error in the regular componentVi,j+1satisfies the following estimate at thej+1thtime level:(62)Vi,j+1−Vj+1xi≤CN−1ε2−δ22α∗+N−1,i=1,2,…,N−1.Proof.
For the error in the regular component,(63)LΔt,NVj+1xi−Vi,j+1=LΔt,NVj+1xi−LΔt,NVi,j+1≤Cε2−δ22α∗∫xi−1xi+1d3Vj+1ξdx3dξ+Chi∫xi−1xid3Vj+1ξdx3dξ.≤Cε2−δ22α∗hi+hi+1+hi2,using Lemma17.
LetH=maxihi; then, LΔt,NVj+1xi−Vi,j+1≤CHε2−δ2/2α∗+H. Using the barrier function Zj+1xi=CHε2−δ2/2α∗+H1+xi and applying Lemmas 14 and 15, we have(64)Vj+1xi−Vi,j+1≤CHε2−δ2/2α∗+H=CN−1ε2−δ2/2α∗+N−1.Theorem 2 (error in singular component).
The error in the singular component satisfies the following bound:(65)Wi,j+1−Wj+1xi≤CN−2,i=0,1,2,…,N/2,CN−2lnN2,i=N/2+1,…,N.Proof.
Let us analyse the error in singular component. We know that(66)W0,j+1=Wj+10≤Ce−a∗−δα∗/ε2−δ2/2α∗=C∏j=1Ne−a∗−δα∗hj/ε2−δ2/2α∗≤∏j=1N1+a∗−δα∗hjε2−δ2/2α∗−1,and W1,j+1=Wj+11≤C and LΔt,NWi,j+1=0 for i=1,2,…,N−1. Now, let us use Lemma 16 and considering the barrier function Zj+1xi as(67)Zj+1xi=C∏k=1N1+a∗−δα∗hkε2−δ2/2α∗−1Sj+1xi,for Wi,j+1 by discrete maximum principle, we get(68)Wi,j+1≤Zj+1xi=C∏k=1N1+a∗−δα∗hkε2−δ2/2α∗−1,∀i,i=0,1,…,N.
Therefore, using Lemmas18 and 12, we get(69)Wj+1xi−Wi,j+1=C∏k=1N1+a∗−δα∗hkε2−δ2/2α∗−1,∀i,i=0,1,…,N,on the coarse mesh part, we have(70)Wj+1xi−Wi,j+1≤CN−2,∀i,i=1,…,N/2,and, on the finer part, we compute the error as we did for the smooth component, that is, by using the consistency and barrier function, but we restrict our domain to 1−τ,1 rather than 0,1. Now, taking i=N/2 in (71), we get(71)Wj+1xN/2−WN/2,j+1≤CN−2,Wj+1xN−WN,j+1=0.
Now, modifying Lemma17 for N/2<i<N and using Lemma 12, we have(72)LΔt,NWj+1xi−Wi,j+1≤Chi+hi+1∫xi−1xi+1ε2−δ22α∗Wj+14x+Wj+13xdx≤Cε2−δ22α∗−3hi+hi+1∫xi−1xi+1e−a∗−δα∗1−x/ε2−δ2/2α∗dx=Cε2−δ22α∗−2hi+hi+1e−a∗−δα∗1−xi+1/ε2−δ2/2α∗−e−a∗−δα∗1−xi−1/ε2−δ2/2α∗=Cε2−δ22α∗−2hi+hi+1e−a∗−δα∗1−xi/ε2−δ2/2α∗×ea∗−δα∗hi/ε2−δ2/2α∗−e−a∗−δα∗hi/ε2−δ2/2α∗=Cε2−δ22α∗−2hi+hi+1e−a∗−δα∗1−xi/ε2−δ2/2α∗×sinha∗−δα∗hi/ε2−δ2/2α∗≤ε2−δ22α∗−3N−2e−a∗−δα∗1−xi/ε2−δ2/2α∗,≤ε2−δ22α∗−1N−2lnN2e−a∗−δα∗1−xi/ε2−δ2/2α∗,≤Cε2−δ22α∗−1N−2lnN∏j=i+1N1+a∗−δα∗hε2−δ2/2α∗−1.By Lemma 18 and using sinhx≤Cx for x∈0,1, now, consider a barrier function(73)Zi,j+1=CN−2+N−2lnN2∏k=1N1+a∗−δα∗hε2−δ2/2α∗−1Si,j+1,and, using the discrete maximum principle, we obtain(74)Wj+1xi−Wi,j+1≤N−2+N−2lnN2N−41−i/N≤CmaxN−2,N−2lnN2N−41−i/N.
Hence, the proof is complete.Theorem 3.
The error due to the spatial discretization of the computed solution satisfies the following estimate:(75)uj+1xi−Ui,j+1≤CN−1ε2−δ22α∗+N−1,i=0,1,2,…,N/2,CN−2lnN2,i=N/2+1,…,N.Proof.
Combining the error estimates given above for the regular and singular components in Theorems1 and 2 gives the following result.Theorem 4.
Letuxi,tj+1, Uj+1xi, and Ui,j+1 be the solutions of the problems in (5), (16), and (40), respectively; then the proposed scheme satisfies the following error estimate:(76)uj+1xi−Ui,j+1≤CN−1ε2−δ22α∗+N−1+Δt2,i=0,1,…,N/2,CN−2lnN2+Δt2,i=N/2+1,…,N.Proof.
Using the error for the temporal and spatial discretization in Lemma9 and Theorem 3 gives the required bound.


## 3.1. Temporal Semidiscretization

Subdivide the time domain0,T into M intervals as t1=0,tj=jΔt,j=1,2,…,M, where Δt=T/M−1. Let Uj+1x denote the approximation of ux,tj+1 at the j+1th time level discretization. Using Crank-Nicolson method for semidiscretizing the problem in (4), we obtain a system of boundary value problems:(16)LΔtUj+1x=gj+1x,x∈Ω,j=1,2,…,M−1,U1x=u0x,x∈Ω¯,Uj+10=ϕ0,tj+1,Uj+11=ψ1,tj+1,where(17)LΔtUj+1x=−ε2−δ2/2α∗2d2dx2Uj+1x+ax−δαx2ddxUj+1x+dxUj+1x,gj+1x=ε2−δ2/2α∗2d2dx2Ujx−ax−δαx2ddxUjx−exUjx+fx,tj+1+fx,tj2,for dx=αx+βx/2+1/Δt and ex=αx+βx/2−1/Δt.Next, we state the discrete maximum principle for semidiscretized problem as follows.Lemma 7.
(semidiscrete maximum principle). For eachj=1,2,…,M−1, let zj+1 be sufficiently smooth function on domain Ω¯. If zj+10≥0,zj+11≥0 and LΔtzj+1x≥0,x∈Ω, then zj+1x≥0,∀x∈Ω¯.Proof.
Suppose thatx∗ is such that zj+1x∗=minx∈Ω¯zj+1x<0. From assumption, it is clear that x∗∉0,1 implies that x∗∈0,1. Applying the property of extrema values in calculus gives d/dxzj+1x∗=0, and d2/dx2zj+1x∗≥0, given that LΔtzj+1x∗<0 which is a contradiction to LΔtzj+1x∗≥0,x∈Ω. Therefore, we conclude that zj+1x≥0,∀x∈Ω¯. Hence, operator LΔt satisfies the semidiscrete maximum principle.The local error for each time step is denoted byej+1 and given by ej+1:=ux,tj+1−Uj+1x.Lemma 8 (local error estimate).
Considering the result in Lemma4, the local error in the temporal semidiscretization is given by(18)ej+1≤C2Δt3,j=1,2,…,M−1.Proof.
Using Taylor’s series expansion abouttj+1/2, we obtain(19)ux,tj+1=ux,tj+1/2+Δt2∂∂tux,tj+1/2+Δt28∂2∂t2ux,tj+1/2+OΔt3,ux,tj=ux,tj+1/2−Δt2∂∂tux,tj+1/2+Δt28∂2∂t2ux,tj+1/2+OΔt3.
Substituting (19) into (5), we obtain(20)ux,tj+1−ux,tjΔt=∂∂tux,tj+1/2+OΔt2=ε2−δ22α∗∂2∂x2ux,tj+1/2−ax−δαx∂∂xux,tj+1/2−αx+βxux,tj+1/2+fx+tj+1/2+OΔt2,where ux,tj+1/2=ux,tj+1+ux,tj/2+OΔt2. By simplifying, (20) can be rewritten as(21)LΔtux,tj+1=gx,tj+1+OΔt3,and, from (8), we have(22)LΔtUj+1x=gj+1x,x∈Ω.
From (21) and (22), we obtain(23)LΔtux,tj+1−Uj+1x=LΔtej+1x≤C2Δt3,with the boundary conditions u0,tj+1−Uj+10=ej+10=0 and u1,tj+1−Uj+11=ej+11=0. Hence, applying the maximum principle gives ej+1≤C2Δt3.Next, we need to show the bound for the global error of the temporal discretization. LetEj+1 denote the global error estimate up to the j+1th time step.Lemma 9 (Global error estimate).
The global error estimate up totj+1 time step is given by(24)Ej+1≤C2Δt2,j=1,2,…,M−1.Proof.
Using the local error estimate up to thej+1th time step given in Lemma 8, we obtain the global error estimate at the j+1th time step as(25)Ej+1=∑l=1j+1el≤e1+e2+⋯+ej+1≤C2j+1Δt3=j+1ΔtC2Δt2≤C2TΔt2, sincej+1Δt≤T=C1Δt2,C2T≔C1,where Ci is a constant independent of ε2−δ2/2α∗ and Δt.Lemma 10.
Forj=1,2,…,M−1, let Uj+1x be a solution of (16); then the solution Uj+1x satisfies the following bound:(26)Uj+1x≤LΔtUj+1d∗+maxϕ0,tj+1,ψ1,tj+1,where dx≥d∗>0.Proof.
Construct barrier functionsπj+1±x=LΔtUj+1/d∗+maxϕ0,tj+1,ψ1,tj+1±Uj+1x. We need to show that if πj+1±0≥0,πj+1±1≥0 and LΔtπj+1±x≥0,x∈Ω, this implies that πj+1±x≥0,∀x∈Ω¯. We easily obtain that πj+1±0≥0 and πj+1±1≥0, and(27)LΔtπj+1±x=−ε2−δ2/2α∗2d2dx2πj+1±x+ax−δαx2ddxπj+1±x+dxπj+1±x=−ε2−δ2/2α∗2d2dx2Uj+1x+ax−δαx2ddxUj+1x+dxπj+1±x=dxLΔtUj+1d∗+maxϕ0,tj+1,ψ1,tj+1±LΔtUj+1x≥0,sincedx≥d∗.
Using the semidiscrete maximum principle, we conclude thatπj+1±x≥0,∀x∈Ω¯. Hence, the required bound is satisfied.


## 3.2. Decomposition of Solution

Since the bounds on the derivatives of the solution are not sharp enough for the proof of uniform convergence, we need to derive more stronger error bounds. This can be achieved by decomposing the solution into regular and singular component and doing error bound separately for each. Let(28)Uj+1x=Vj+1x+Wj+1x,x∈Ω¯.The regular component is the solution to the nonhomogeneous problem:(29)LΔtVj+1x=gj+1x,x∈Ω,Vj+10=Uj+10,V0x=U0x,x∈Ω¯,and the singular component is the solution to the homogeneous problem:(30)LΔtWj+1x=0,x∈Ω,Wj+10=0,Wj+11=Uj+11−Vj+11,W0x=0,x∈Ω.Lemma 11.
Forj=1,2,…,M−1, the regular component of the solution and its derivatives satisfies the bound(31)dkdxkVj+1x≤C1+ε2−δ22α∗2−k,x∈Ω,k=0,1,2.
In general, the derivatives are bounded as(32)dkdxkVj+1x≤C,x∈Ω,k=0,1,2,3.Proof.
See [18].Lemma 12.
Forj=1,2,…,M−1, the singular component of the solution and its derivatives satisfies the bound:(33)dkdxkWj+1x≤Cε2−δ22α∗−kexp−a∗−δα∗ε2−δ2/2α∗1−x,x∈Ω,k=0,1,2,3.Proof.
See [18].The next lemma gives a bound for the derivatives of solution.Lemma 13.
The solutions of the boundary value problems in (16) satisfy the bound(34)dkdxkUj+1x≤C1+ε2−δ22α∗−kexp−a∗−δα∗ε2−δ2/2α∗1−x,x∈Ω¯,k=0,1,2,3.Proof.
The results from Lemmas11 and 12 give the required bound.


## 3.3. Spatial Discretization and the Hybrid FDM

Letxii=0N be the discretized domain of the spatial domain, where N be the number of grid points and must be even positive integer. For each i≥1, we define hi=xi−xi−1. Since the considered problem exhibits right boundary layer, we set a mesh transition parameter τ which divides the domain 0,1 into outer layer region Ω1=0,1−τ and boundary layer region Ω2=1−τ,1, where(35)τ=min12,σε2−δ2/2α∗d∗lnN,where N denotes the number of mesh points in spatial discretization and σ denotes a constant that represents the order of the scheme. We set a uniform mesh N/2 in Ω1 with mesh spacing H=21−τ/N and similarly a uniform mesh N/2 is placed in Ω2 with the mesh spacing h=2τ/N. So, the mesh point xi is given as(36)xi=21−τNi,i=0,1,2,…,N2,1−τ+2τNi−N/2,i=N/2,…,N.In a similar fashion, one can generate a Shishkin mesh for the left boundary layer problem.For the mesh functionsUj+1xi at the grid points xi, the backward and forward difference operators are given as(37)Dx−Uj+1xi=Ui,j+1−Ui−1,j+1hi,Dx+Uj+1xi=Ui+1,j+1−Ui,j+1hi+1,respectively. The central difference operator Dx0 is(38)Dx0Uj+1xi=Ui+1,j+1−Ui−1,j+1hi+hi+1.Similarly, the second-order derivativeDx+Dx− is given by(39)Dx−Dx+Uj+1xi=2hi+hi+1Ui+1,j+1−Ui,j+1hi+1−Ui,j+1−Ui−1,j+1hi.Remark 2.
The approximationDx0Ui,j+1is second-order consistent on any mesh, just as on equidistant meshes, but this is not the case forDx+Dx−Ui,j+1: a termhi+1−hid3/dx3Uj+1xi, which arises in the consistency error analysis, is only of first order on arbitrary mesh. For more details, the reader is referred to [3, 19]. Whenεis small, the central difference approximationDx0Uj+1xitod/dxUj+1xon equidistant grid leads to nonphysical oscillations of the computed solution which is due to the loss in stability, unless the grid size is very small.To overcome the challenge indicated in Remark2, hybrid numerical scheme containing central difference in the boundary layer region and mid-point upwind scheme on outer region is developed in order to get better accuracy and rate of convergence. Using the operators above, the discrete hybrid scheme can be formulated as(40)LΔt,NUi,j+1=gi,j+1,i=1,2,…,N−1,U0,j+1=ϕ0,tj+1,UN,j+1=ψ1,tj+1,where(41)LΔt,NUi,j+1=−ε2−δ2/2α∗2Dx+Dx−Ui,j+1+axi−1/2−δαxi−1/22Dx−Ui,j+1+dxi−1/2Ui,j+1,i=1,2,…,N/2,−ε2−δ2/2α∗2Dx+Dx−Ui,j+1+axi−δαxi2Dx0Ui,j+1+dxiUi,j+1,i=N/2+1,…,N−1,gi,j+1=ε2−δ2/2α∗2Dx+Dx−Ui,j−axi−1/2−δαxi−1/22Dx−Ui,j−exi−1/2Ui,j+12fi−1/2,j+fi−1/2,j+1,i=1,2,…,N/2,ε2−δ2/2α∗2Dx+Dx−Ui,j+1−axi−δαxi2Dx0Ui,j+1−exiUi,j+1+12fi,j+fi,j+1,i=N/2+1,…,N−1,and the notation dxi−1/2 stands for dxi+dxi−1/2 and similarly for others.


## 3.4. Properties of the Discrete Scheme and Uniform Convergence

Assume thatτ=min1/2,σε2−δ2/2α∗/d∗lnN. We assign each of the subintervals 0,1−τ and 1−τ,1 with N/2 number of equidistant grid points. Let H be the mesh width in the subinterval 0,1−τ and let h be the mesh width in 1−τ,1. These mesh widths satisfy(42)N−1≤H≤2N−1,h=4d∗ε2−δ22α∗N−1lnN.To attain stability condition while using central difference scheme on the singular component, we need the following condition [19]:(43)NlnN≥2maxdxd∗.Lemma 14 (discrete maximum principle).
The discrete schemeLΔt,NUi,j+1=gi,j+1,i=1,2,…,N−1and for givenU0,j+1andUN,j+1for fixedjhas a solution. IfLΔt,NUi,j+1≤LΔt,NBi,j+1fori=1,2,…,N−1and ifU0,j+1≤B0,j+1andUN,j+1≤BN,j+1, thenUi,j+1≤Bi,j+1,∀i=0,1,2,…,N.Proof.
The matrix associated with operatorLΔt,N is of size N+1×N+1 and satisfies the property of M-matrix with the given condition in (43). See the detailed proof in [20].Lemma 15.
LetZi,j+1=1+xi,i=0,1,…,N. Then, there exists a positive constant C such that LΔt,NZi,j+1≥C,i=1,2,…,N−1.Proof.
Using the discrete operator on the mesh functionZi,j+1, for 1≤i≤N/2, we have(44)LΔt,NZi,j+1=−ε2−δ2/2α∗2Dx+Dx−Zi,j+1+axi−1/2−δαxi−1/22Dx−Zi,j+1+dxi−1/2Zi,j+1,=−ε2−δ2/2α∗2Dx+Dx−1+xi+axi−1/2−δαxi−1/22Dx−1+xi+dxi−1/21+xi,=0+axi−1/2−δαxi−1/22Dx−1+xi+dxi−1/21+xi,≥C,sincedxi−1/2is bounded.For the case whereN/2+1≤i≤N−1, we have(45)LΔt,NZi,j+1=−ε2−δ2/2α∗2Dx+Dx−Zi,j+1+axi−δαxi2Dx0Zi,j+1+dxiZi,j+1,=−ε2−δ2/2α∗2Dx+Dx−1+xi+axi−δαxi2Dx01+xi+dxi1+xi,=0+axi−δαxi2Dx01+xi+dxi1+xi,≥C,sincedxiis bounded.Lemma 16.
Fori=0,1,…,Nand fixedj, define a mesh function(46)Si,j+1=∏k=1i1+a∗−δα∗hkε2−δ2/2α∗,with the convention that if i=0, then S0=1. Then, for i=1,2,…,N−1, we have(47)LΔt,NSi,j+1≥Cmaxε2−δ2/2α∗,hiSi,j+1,for some constant C.Proof.
SinceSi,j+1−Si−1,j+1/hi=a∗−δα∗/ε2−δ2/2α∗Si−1,j+1, we have Si−1,j+1=ε2−δ2/2α∗/ε2−δ2/2α∗+a∗−δα∗hiSi,j+1. Now, putting these on the discrete operator, we have(48)LΔt,NSi,j+1=−ε2−δ2/2α∗hi+hi+1Si+1,j+1−Si,j+1hi+1−Si,j+1−Si−1,j+1hi+axi−1/2−δαxi−1/22Si,j+1−Si−1,j+1hi+dxi−1/2Si,j+1,i=1,2,…,N/2,−ε2−δ2/2α∗hi+hi+1Si+1,j+1−Si,j+1hi+1−Si,j+1−Si−1,j+1hi+axi−δαxi2Si+1,j+1−Si−1,j+1hi+hi+1+dxiSi,j+1,i=N/2+1,…,N−1.For i=1,2,…,N/2, we have(49)LΔt,NSi,j+1=−ε2−δ2/2α∗hi+hi+1a∗−δα∗Si,j+1−Si−1,j+1ε2−δ2/2α∗+axi−1/2−δαxi−1/22a∗−δα∗Si−1,j+1ε2−δ2/2α∗+dxi−1/2Si,j+1=a∗−δα∗2ε2−δ2/2α∗Si−1,j+1axi−1/2−δαxi−1/2−2a∗−δα∗hihi+hi+1+dxi−1/2Si,j+1=a∗−δα∗2ε2−δ2/2α∗+a∗−δα∗hiSi,j+1axi−1/2−δαxi−1/2−2a∗−δα∗hihi+hi+1+dxi−1/2ε2−δ2/2α∗+a∗−δα∗hia∗−δα∗≥Si,j+1Cmaxε2−δ2/2α∗,hi.
Hence, the required result holds. For the case wherei=N/2+1,…,N−1,(50)LΔt,NSi,j+1=−ε2−δ2/2α∗hi+hi+1a∗−δα∗Si,j+1−Si−1,j+1ε2−δ2/2α∗+axi−δαxi2×1+a∗−δα∗hihi+1ε2−δ2/2α∗hi+hi+1a∗−δα∗Si−1,j+1ε2−δ2/2α∗+dxiSi,j+1=−a∗−δα∗2hiε2−δ2/2α∗2hi+hi+1Si−1,j+1+axi−δαxia∗−δα∗2ε2−δ2/2α∗×1+a∗−δα∗hihi+1ε2−δ2/2α∗hi+hi+1Si−1,j+1+dxiSi,j+1=−a∗−δα∗2hiε2−δ2/2α∗hi+hi+1ε2−δ2/2α∗+a∗−δα∗hiSi,j+1+axi−δαxia∗−δα∗2ε2−δ2/2α∗+a∗−δα∗hi×1+a∗−δα∗hihi+1ε2−δ2/2α∗hi+hi+1Si,j+1+dxiSi,j+1=a∗−δα∗Si,j+1ε2−δ2/2α∗+a∗−δα∗hi−a∗−δα∗hiε2−δ2/2α∗hi+hi+1+axi−δαxi1+a∗−δα∗hihi+1ε2−δ2/2α∗hi+hi+1+dxiε2−δ2/2α∗+a∗−δα∗hia∗−δα∗≥Si,j+1Cmaxε2−δ2/2α∗,hi.Lemma 17.
Letqxbe any smooth function defined on0,1. Fori=1,2,…,N−1, define truncation errorTEtoqxatxiwhich is to beLΔtqxi−LΔt,Nqi; then there exists a constantCsuch that(51)TE=LΔtqxi−LΔt,Nqi≤Cε2−δ22α∗∫xi−1xi+1d3qξdx3dξ+Chi∫xi−1xid3qξdx3dξ,for ξ∈xi−1,xi+1.Proof.
See Lemma 4.1 in [16].Lemma 18.
For eachi=1,2,…,N−1, we have(52)exp−a∗−δα∗1−xiε2−δ2/2α∗≤∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1.Proof.
For eachi, we have(53)exp−a∗−δα∗hiε2−δ2/2α∗=ea∗−δα∗hi/ε2−δ2/2α∗−1≤1+a∗−δα∗hiε2−δ2/2α∗−1.
Now multiplying fork=i+1,…,N, we obtain(54)∏k=i+1Nea∗−δα∗hi/ε2−δ2/2α∗−1≤∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1⇒e−a∗−δα∗1−xiε2−δ2/2α∗≤∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1.Lemma 19.
For Shishkin mesh defined above, there exists a constantC such that(55)∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1≤CN−41−i/N,i=N/2,…,N.Proof.
LetN/2≤i≤N by Lemma 1 in [20].(56)1+a∗−δα∗hkε2−δ2/2α∗−1=1−a∗−δα∗hiε2−δ2/2α∗+a∗−δα∗hi≤e−a∗−δα∗hi/ε2−δ2/2α∗+a∗−δα∗hi.
Therefore, we have(57)∏k=i+1N1+a∗−δα∗hkε2−δ2/2α∗−1≤∏k=i+1Ne−a∗−δα∗hi/ε2−δ2/2α∗+a∗−δα∗hi=e−a∗−δα∗1−xi/ε2−δ2/2α∗+a∗−δα∗hi=e−4N−iN−1lnN/1+4N−1lnN,using,h=4d∗ε2−δ22α∗N−1lnN=N−4N−iN−1/1+4N−1lnN=N−41−i/NN161−i/NN−1lnN/1+4N−1lnN≤CN−41−i/N.
SinceN161−i/NN−1lnN/1+4N−1lnN is bounded for N≥2, the desired result is obtained.


## 3.5. Decomposition of Numerical Solution

Here, we decompose the numerical solutionUi,j+1 into regular and singular components in a similar fashion to the continuous case:(58)Ui,j+1=Vi,j+1+Wi,j+1,where the regular part satisfies the nonhomogeneous equation:(59)LΔt,NVi,j+1=gi,j+1,Vj+10=Uj+10,and the singular component Wi,j+1 satisfies the homogeneous equation:(60)LΔt,NWi,j+1=0,Wj+10=0,Wj+11=Uj+11−Vj+11,and now the error in the numerical solution can also be decomposed as(61)Ui,j+1−uj+1xi=Vi,j+1−Vj+1xi+Wi,j+1−Wj+1xi.We estimate the error in the regular and singular solution separately.Theorem 1 (error in the regular component).
The error in the regular componentVi,j+1satisfies the following estimate at thej+1thtime level:(62)Vi,j+1−Vj+1xi≤CN−1ε2−δ22α∗+N−1,i=1,2,…,N−1.Proof.
For the error in the regular component,(63)LΔt,NVj+1xi−Vi,j+1=LΔt,NVj+1xi−LΔt,NVi,j+1≤Cε2−δ22α∗∫xi−1xi+1d3Vj+1ξdx3dξ+Chi∫xi−1xid3Vj+1ξdx3dξ.≤Cε2−δ22α∗hi+hi+1+hi2,using Lemma17.
LetH=maxihi; then, LΔt,NVj+1xi−Vi,j+1≤CHε2−δ2/2α∗+H. Using the barrier function Zj+1xi=CHε2−δ2/2α∗+H1+xi and applying Lemmas 14 and 15, we have(64)Vj+1xi−Vi,j+1≤CHε2−δ2/2α∗+H=CN−1ε2−δ2/2α∗+N−1.Theorem 2 (error in singular component).
The error in the singular component satisfies the following bound:(65)Wi,j+1−Wj+1xi≤CN−2,i=0,1,2,…,N/2,CN−2lnN2,i=N/2+1,…,N.Proof.
Let us analyse the error in singular component. We know that(66)W0,j+1=Wj+10≤Ce−a∗−δα∗/ε2−δ2/2α∗=C∏j=1Ne−a∗−δα∗hj/ε2−δ2/2α∗≤∏j=1N1+a∗−δα∗hjε2−δ2/2α∗−1,and W1,j+1=Wj+11≤C and LΔt,NWi,j+1=0 for i=1,2,…,N−1. Now, let us use Lemma 16 and considering the barrier function Zj+1xi as(67)Zj+1xi=C∏k=1N1+a∗−δα∗hkε2−δ2/2α∗−1Sj+1xi,for Wi,j+1 by discrete maximum principle, we get(68)Wi,j+1≤Zj+1xi=C∏k=1N1+a∗−δα∗hkε2−δ2/2α∗−1,∀i,i=0,1,…,N.
Therefore, using Lemmas18 and 12, we get(69)Wj+1xi−Wi,j+1=C∏k=1N1+a∗−δα∗hkε2−δ2/2α∗−1,∀i,i=0,1,…,N,on the coarse mesh part, we have(70)Wj+1xi−Wi,j+1≤CN−2,∀i,i=1,…,N/2,and, on the finer part, we compute the error as we did for the smooth component, that is, by using the consistency and barrier function, but we restrict our domain to 1−τ,1 rather than 0,1. Now, taking i=N/2 in (71), we get(71)Wj+1xN/2−WN/2,j+1≤CN−2,Wj+1xN−WN,j+1=0.
Now, modifying Lemma17 for N/2<i<N and using Lemma 12, we have(72)LΔt,NWj+1xi−Wi,j+1≤Chi+hi+1∫xi−1xi+1ε2−δ22α∗Wj+14x+Wj+13xdx≤Cε2−δ22α∗−3hi+hi+1∫xi−1xi+1e−a∗−δα∗1−x/ε2−δ2/2α∗dx=Cε2−δ22α∗−2hi+hi+1e−a∗−δα∗1−xi+1/ε2−δ2/2α∗−e−a∗−δα∗1−xi−1/ε2−δ2/2α∗=Cε2−δ22α∗−2hi+hi+1e−a∗−δα∗1−xi/ε2−δ2/2α∗×ea∗−δα∗hi/ε2−δ2/2α∗−e−a∗−δα∗hi/ε2−δ2/2α∗=Cε2−δ22α∗−2hi+hi+1e−a∗−δα∗1−xi/ε2−δ2/2α∗×sinha∗−δα∗hi/ε2−δ2/2α∗≤ε2−δ22α∗−3N−2e−a∗−δα∗1−xi/ε2−δ2/2α∗,≤ε2−δ22α∗−1N−2lnN2e−a∗−δα∗1−xi/ε2−δ2/2α∗,≤Cε2−δ22α∗−1N−2lnN∏j=i+1N1+a∗−δα∗hε2−δ2/2α∗−1.By Lemma 18 and using sinhx≤Cx for x∈0,1, now, consider a barrier function(73)Zi,j+1=CN−2+N−2lnN2∏k=1N1+a∗−δα∗hε2−δ2/2α∗−1Si,j+1,and, using the discrete maximum principle, we obtain(74)Wj+1xi−Wi,j+1≤N−2+N−2lnN2N−41−i/N≤CmaxN−2,N−2lnN2N−41−i/N.
Hence, the proof is complete.Theorem 3.
The error due to the spatial discretization of the computed solution satisfies the following estimate:(75)uj+1xi−Ui,j+1≤CN−1ε2−δ22α∗+N−1,i=0,1,2,…,N/2,CN−2lnN2,i=N/2+1,…,N.Proof.
Combining the error estimates given above for the regular and singular components in Theorems1 and 2 gives the following result.Theorem 4.
Letuxi,tj+1, Uj+1xi, and Ui,j+1 be the solutions of the problems in (5), (16), and (40), respectively; then the proposed scheme satisfies the following error estimate:(76)uj+1xi−Ui,j+1≤CN−1ε2−δ22α∗+N−1+Δt2,i=0,1,…,N/2,CN−2lnN2+Δt2,i=N/2+1,…,N.Proof.
Using the error for the temporal and spatial discretization in Lemma9 and Theorem 3 gives the required bound.


## 4. Numerical Results and Discussion

In this section, we consider numerical examples to illustrate the theoretical findings of the developed hybrid numerical scheme.Example 1.
Consider the problem(77)∂u∂t−ε2∂2u∂x2+∂u∂x+2−x2ux−δ,t+1+x2+exux,t=50x1−x3.T=2 is subject to the initial condition ux,0=0,x∈0,1 and the interval-boundary conditions ϕx,t=0,−δ≤x≤0,ψ1,t=0,ont∈0,2.Example 2.
Consider the problem(78)∂u∂t−ε2∂2u∂x2+2+x+x2∂u∂x+1+x22ux−δ,t=sinπx1−x,with T=1 subject to the initial condition ux,0=0,x∈0,1 and the interval-boundary conditions ϕx,t=0,−δ≤x≤0,ψ1,t=0,ont∈0,1.In the considered examples, the exact solution of the problems is not known. We use the double mesh technique to calculate maximum absolute error. LetUi,jN,M denote the computed solution of the problem on N,M number of mesh points in x and t directions, respectively, and Ui,j2N,2M denotes the computed solution on double number of mesh points 2N,2M by including the mid-points xi+1/2=xi+1+xi/2 and tj+1/2=tj+1+tj/2 into the mesh points. The maximum point-wise absolute error is given by(79)Eε,δN,M=maxi,jUi,jN,M−Ui,j2N,2M.Theε-uniform error estimate is calculated using(80)EN,M=maxε,δEε,δN,M.The rate of convergence of the scheme is calculated using(81)rε,δN,M=log2Eε,δN,MEε,δ2N,2M,and the ε-uniform rate of convergence is calculated using(82)rN,M=log2EN,ME2N,2M.In Figure1, solution profiles of Examples 1 and 2 are plotted for different values of time level by taking ε=2−5 and δ=0.1ε. The solutions of the problems in the considered examples exhibit right boundary layer of thickness Oε2. One observes in Figures 2 and 3 that as the perturbation parameter ε goes small, the boundary layer formation becomes visible. In these figures, the effect of the singular perturbation parameter on the solution is shown. In addition, for small ε, there are a sufficient number of mesh points in the boundary layer region, which shows the layer resolving nature of the scheme. In Figure 4, the effect of delay parameter on the behaviour of the solution is shown by using different values for delay parameter for the test examples by keeping ε=2−5 at T=1. In Figure 5, the absolute errors of the proposed scheme for ε=2−20,N=27, and M=27 are given. As is seen on the figures, the absolute error is dominant in the boundary layer region.Figure 1
Solution profile for different time level on (a) Example1 and (b) Example 2 for ε=2−5.
(a)(b)Figure 2
3D view of solution of Example1, with boundary layer, on (a) ε=1 and (b) ε=2−20.
(a)(b)Figure 3
3D view of solution of Example2, with boundary layer, on (a) ε=1 and (b) ε=2−20.
(a)(b)Figure 4
Effect of the delay on solution on (a) Example1 and on (b) Example 2, for ε2=2−5 at T = 1.
(a)(b)Figure 5
Absolute error on (a) Example1 and on (b) Example 2 for ε=2−20,N=M=27.
(a)(b)In Tables1 and 2, the maximum point-wise absolute error, ε-uniform error, and the ε-uniform rate of convergence of the scheme are given. The results are given for the boundary layer and regular regions separately. In each column (or for each N and M), as ε⟶0, the maximum absolute error becomes uniform. This shows that the convergence of the scheme is independent of the perturbation parameter. In the last two rows of Tables 1 and 2, we have the ε-uniform error and the ε-uniform rate of convergence of the scheme. One observes in these tables that the scheme converges with order of convergence of almost two. In Tables 3 and 4, the maximum absolute error for different values of the delay parameter is given for ε=2−10. The results in these tables show that the convergence of the scheme is also independent of the delay parameter. In Table 5, the comparison of the ε-uniform error and the ε-uniform rate of convergence is given. The result in the proposed scheme is more accurate than those in the schemes in [5, 10, 21]. Based on the results given in Tables 1–5, the proposed scheme is accurate, stable, and ε-uniformly convergent with rate of convergence of almost two.Table 1
Maximum absolute error of Example1 using the proposed scheme at δ=0.3ε.
Inner layer region:N/2≤i≤NOuter layer region:1≤i≤N/2ε↓N⟶2526272825262728M⟶6012024048060120240480206.8812e 2013041.7201e −044.2719e −051.0658e −058.4699e −042.1210e −045.2856e −051.3200e −052−22.0498e −034.5521e −049.6120e −055.4087e −053.5913e −038.3333e −041.8490e −045.7215e −052−47.0311e −032.0223e −035.2379e −041.2986e −044.7288e −031.3098e −033.9806e −046.1801e −052−67.2256e −032.0602e −035.2983e −041.3318e −044.7755e −031.3222e −033.4259e −048.6541e −052−87.2299e −032.0616e −035.3015e −041.3325e −044.7640e −031.3219e −033.4265e −048.6341e −052−107.2283e −032.0615e −035.3016e −041.3349e −044.7596e −031.3212e −033.4248e −048.6462e −052−127.2278e −032.0614e −035.3016e −041.3349e −044.7584e −031.3210e −033.4243e −048.6462e −052−147.2276e −032.0614e −035.3016e −041.3349e −044.7581e −031.3209e −033.4241e −048.6463e −052−167.2276e −032.0614e −035.3016e −041.3349e −044.7580e −031.3209e −033.4241e −048.6467e −052−187.2276e −032.0614e −035.3016e −041.3349e −044.7580e −031.3209e −033.4241e −048.6467e −052−207.2276e −032.0614e −035.3016e −041.3349e −044.7580e −031.3209e −033.4241e −048.6467e −05EN,M7.2276e −032.0614e −035.3016e −041.3349e −044.7580e −031.3209e −033.4241e −048.6467e −05rN,M1.80991.95911.9896—1.84881.94771.9855—Table 2
Maximum absolute error of Example2 using the proposed scheme at δ=0.3ε.
Inner layer region:N/2≤i≤NOuter layer region:1≤i≤N/2ε↓N⟶2526272825262728M⟶6012024048060120240480201.1924e −033.0625e −047.7028e −051.9285e −051.1355e −033.0259e −047.6798e −051.9271e −052−23.9610e −039.4744e −042.4014e −045.7324e −053.4402e −038.9924e −042.3404e −046.1054e −052−44.3261e −031.2500e −033.7900e −041.0390e −044.5385e −031.5135e −034.5058e-041.2278e −042−64.3558e −031.2988e −034.2034e −041.04950e −044.6146e −031.5856e −034.9451e −041.3045e −042−84.3547e −031.3022e −034.2339e −041.0490e −044.6182e −031.5895e −034.9738e −041.3023e −042−104.3538e −031.3024e −034.2352e −041.6249e −044.6181e −031.5894e −034.9751e −041.5576e −042−124.3536e −031.3024e −034.2352e −041.6249e −044.6181e −031.5893e −034.9751e −041.5576e −042−144.3535e −031.3024e −034.2351e −041.6249e −044.6180e −031.5893e −034.9751e −041.5576e −042−164.3535e −031.3024e −034.2351e −041.6249e −044.6180e −031.5893e −034.9750e −041.5576e −042−184.3535e −031.3024e −034.2351e −041.6249e-044.6180e −031.5893e −034.9750e −041.5576e −042−204.3535e −031.3024e −034.2351e −041.6249e −044.6180e −031.5893e −034.9750e −041.5576e −04EN,M4.3535e −031.3024e −034.2351e −041.6249e −044.6180e −031.5893e −034.9750e −041.5576e −04rN,M1.74101.62071.3820—1.53891.67561.6754—Table 3
Maximum absolute error of Example1 using the proposed method for different values of delay parameter at ε=2−10.
δ↓, N=M⟶2526272807.1898e −032.0349e −035.2299e −041.3165e −040.1ε7.1902e −032.0349e −035.2300e −041.3165e −040.3ε7.1912e −032.0349e −035.2302e −041.3166e −040.5ε7.1921e −032.0349e −035.2304e −041.3166e −040.7ε7.1930e −032.0349e −035.2305e −041.3167e −04Table 4
Maximum absolute error of Example2 using the proposed method for different values of delay parameter at ε=2−10.
δ↓, N=M⟶2526272804.2349e −031.4965e −034.8385e −041.6247e −040.1ε4.2349e −031.4965e −034.8386e −041.6247e −040.3ε4.2350e −031.4966e −034.8387e −041.6248e −040.5ε4.2350e −031.4966e −034.8388e −041.6248e −040.7ε4.2351e −031.4966e −034.8389e −041.6249e −04Table 5
ε-Uniform error and ε-uniform rate of convergence of the proposed scheme and results in [5, 10, 21] of Example 2.
Schemes↓N⟶3264128256M⟶60120240480Proposed schemeEN,M4.3535e −031.3024e −034.2351e −041.6249e −04rN,M1.74101.62071.38201.5214Nonstandard FDM in [21]EN,M3.9288e −031.9657e −039.8350e −044.9194e −04rN,M0.99900.99910.99940.9997Fitted operator FDM in [10]EN,M5.9125e −033.6155e −032.1167e −031.1957e −03rN,M0.70960.77240.82400.8662Upwind in [5]EM,N1.1359e −027.2877e −034.4807e −032.6516e −03rM,N0.64030.70170.75690.8054Mid-point upwind in [5]EM,N8.3467e −035.3894e −033.3649e −032.0082e −03rM,N0.63110.67960.74470.7954


## 5. Conclusion

This paper deals with the numerical treatment of singularly perturbed parabolic delay convection-diffusion problems. The solution of the considered problem exhibits boundary layer behaviour. The delay term is approximated using Taylor’s series and resulted in asymptotically equivalent singularly perturbed parabolic partial differential equations. The developed numerical scheme is constituted of Crank-Nicolson method in temporal discretization with hybrid finite difference on Shishkin mesh for the spatial discretization. The stability of the scheme is investigated using construction of barrier function for the solution bound and maximum principle for the existence of unique discrete solution. Theε-uniform convergence of the scheme is proved. The applicability of the scheme is investigated by considering test examples exhibiting boundary layer. The numerical results are tabulated in terms of maximum absolute error and ε-uniform error and are observed to be in agreement with the theoretical results. Effect of the delay on the solution profile is depicted in graphs. The computational result shows uniform convergence of the scheme with order of convergence of almost two.



---

*Source: 6654495-2021-09-11.xml*