# Novel Numerical Scheme for Singularly Perturbed Time Delay Convection-Diffusion Equation

Authors: Mesfin Mekuria Woldaregay; Worku Tilahun Aniley; Gemechis File Duressa
Journal: Advances in Mathematical Physics
 (2021)

---

## Abstract

This paper deals with numerical treatment of singularly perturbed parabolic differential equations having large time delay. The highest order derivative term in the equation is multiplied by a perturbation parameterε, taking arbitrary value in the interval 0,1. For small values of ε, solution of the problem exhibits an exponential boundary layer on the right side of the spatial domain. The properties and bounds of the solution and its derivatives are discussed. The considered singularly perturbed time delay problem is solved using the Crank-Nicolson method in temporal discretization and exponentially fitted operator finite difference method in spatial discretization. The stability of the scheme is investigated and analysed using comparison principle and solution bound. The uniform convergence of the scheme is discussed and proven. The formulated scheme converges uniformly with linear order of convergence. The theoretical analysis of the scheme is validated by considering numerical test examples for different values of ε.


---

## Body

## 1. Introduction

A large number of mathematical models appear in different areas of science and engineering such as in control theory, epidemiology, and laser optics that take into account not only the present state of a physical system but also its past history [1]. These models are described by certain classes of functional differential equations often called delay differential equations. Examples of delays include the time taken for a signal to travel to the controlled object, driver reaction time, the time for the body to produce red blood cells, and cell division time in the dynamics of viral exhaustion or persistence. In the life sciences, delays are often introduced to account for hidden variables and processes which, although not well understood, are known to cause a time lag [1]. Time delays are natural components of the dynamic processes of biology, ecology, physiology, economics, epidemiology, and mechanics [2], and to ignore them is to ignore reality [1].Singularly perturbed delay differential equations are differential equations in which its highest order derivative term is multiplied by small perturbation parameterε and having at least one delay term. A simplified mathematical model for the control system of a furnace used to process metal sheets is represented by singularly perturbed parabolic delay differential equation of the form [3]
(1)∂u∂t−ε∂2u∂x2=vgux,t−τ∂u∂x+cfux,t−τ−ux,t,where u is the temperature distribution in a metal sheet, moving at an instantaneous material strip velocity v and heated by distributed temperature source specified by the function f; both v and f are dynamically adapted by a controlling device monitoring the current temperature distribution. The finite speed of the controller, however, introduces a fixed delay of length τ. When τ=0, this problem becomes a thermal problem without time delay.The presence of the singular perturbation parameterε on the highest order derivative term leads to occurrences of oscillations in the computed solution while using central FDM, Galerkin FEM, and collocation methods on a uniform mesh [4]. While using these methods to avoid the oscillations, an unacceptably large number of mesh points are required when ε is very small. This is not practical due to the limited size of computer memory and rounding off error. Therefore, to overcome this drawback associated with classical numerical methods, authors use the fitted mesh technique which have a fine mesh in the boundary layer region or fitted operator technique. Recently, numerical treatment of singularly perturbed parabolic PDEs has gotten great attention, and different authors in [5–8] have developed uniformly convergent numerical schemes (scheme that converges independent of the influence of the singular perturbation parameter). Currently, authors in [9–13] have developed numerical schemes for solving singularly perturbed parabolic time delay reaction-diffusion equations. In papers [14, 15], the authors considered a singularly perturbed parabolic convection diffusion equation with time delay and degenerate coefficient. They treat the problems using fitted mesh techniques. Different authors in [16–23] developed numerical schemes using fitted mesh techniques for treating singularly perturbed parabolic time delay convection-diffusion equations. In [24], Podila and Kumar used nonstandard FDM for treating singularly perturbed parabolic time delay convection-diffusion equations. Their scheme gives linear order of convergence.The main motive of this paper is to formulate an accurate and uniformly convergent numerical scheme for the singularly perturbed parabolic time delay convection-diffusion-reaction equation using exponentially fitted operator FDM and to establish the stability and uniform convergence of the scheme. The proposed scheme uses the procedure of Roth (i.e., first discretizing in temporal direction followed by discretization in spatial direction) using the Crank-Nicolson method in temporal direction and exponentially fitted operator FDM on spatial direction. In this method, it is not required to have any restriction on the mesh generation.Notation 1.
Throughout this paper,N and M are denoted for the number of mesh intervals in space and time direction discretization, respectively; C is denoted for positive constant independent of ε,N,andM. The norm . denotes the maximum norm defined as g=maxx,tgx,t.


## 2. Continuous Problem

Consider a class of second order singularly perturbed parabolic convection-diffusion-reaction equations having a term with large time delay, on the domainD=Ωx×Ωt=0,1×0,T and on smooth boundary ∂D=Dl∪Db∪Dr is given by
(2)∂u∂t−ε∂2u∂x2+ax∂u∂x+bx,tux,t=−cx,tux,t−τ+fx,t,with the boundary conditions
(3)u0,t=ϕlt,t∈Ωt,u1,t=ϕrt,t∈Ωt,and the interval condition
(4)ux,t=ϕbx,t,x,t∈0,1×−τ,0,where ε∈0,1 is a singular perturbation parameter and τ>0 is the delay parameter. The terminal time T is assumed to satisfy the condition T=kτ for some positive integer k. The functions ax,bx,t,cx,t,andfx,tandϕlt,ϕbx,t,andϕrt are assumed to be sufficiently smooth and bounded and satisfy
(5)ax≥α>0,bx,t≥ς>0,cx,t≥γ>0,x,t∈D¯.This condition ensures that the solution of the problem in (2)–(4) exhibits a boundary layer of thickness Oε on the right side of the spatial domain.Our objective in this paper is to formulate an accurate and parameter uniformly convergent numerical scheme and to discuss the uniform stability and the parameter uniform convergence of the scheme for the considered problem in (2)–(4).


### 2.1. Bounds on the Solution and Its Derivatives

The existence and uniqueness of the solution of (2)–(4) can be established by assuming that the data is Holder continuous and imposing appropriate compatibility conditions at the corner points, using the assumptions of sufficiently smoothness ofϕlt,ϕbx,t, andϕrt.The required compatibility condition at the corner points and the delay term are(6)ϕl0=ϕb0,0,ϕr0=ϕb1,0,(7)dϕl0dt−ε∂2ϕb0,0∂x2=a0∂ϕb0,0∂x+b0,0ϕb0,0=−c0,0u0,−τ+f0,0,dϕr0dt−ε∂2ϕb1,0∂x2=a1∂ϕb1,0∂x+b1,0ϕb1,0=−c0,0u1,−τ+f1,0,so that the data matches at the corner points. For ax,bx,t,cx,t, and fx,t to be continuous on domain D, then (2)–(4) has the unique solution ux,t∈C2D∩C0D¯. In particular, when the compatibility conditions are not satisfied, a unique classical solution still exists but is not differentiable on all of ∂D.The reduced problem corresponding to the singularly perturbed parabolic delay PDE (2)–(4) is given as
(8)∂u0∂t+ax∂u0∂x+bx,tu0x,t=−cx,tu0x,t−τ+fx,t,x,t∈D,with the boundary conditions
(9)u00,t=ϕlt,t∈Ωt,u01,t≠ϕrt,t∈Ωt,and the interval condition
(10)u0x,t=ϕbx,t,x,t∈0,1×−τ,0.It is in the form of hyperbolic delay PDEs. The solutionux,t of (2)–(4) becomes very close to the solution u0x,t of (8)–(10) as ε⟶0.Lemma 2.
The solutionux,t of (2)–(4) satisfies the estimate
(11)ux,t−ϕbx,0≤Ct,x,t∈D¯,where C is a constant independent of ε.Proof.
The result follows from the compatibility condition. See the detailed proof in [16].LetL be a differential operator that denotes the differential equation in (2)–(4) as L=∂/∂t−ε∂2/∂x2+ax∂/∂x+bx,t.Lemma 3.
Suppose the functionzx,t∈C2D∩C0D¯, which satisfies zx,t≥0,x,t∈∂D, and Lzx,t≥0,x,t∈D implies that zx,t≥0,x,t∈D¯.Proof.
Assume there existsx∗,t∗∈D¯, such that
(12)zx∗,t∗=minx,t∈D¯zx,t<0.It is clear that the pointx∗,t∗∉∂D, which implies that x∗,t∗∈D. Since zx∗,t∗=minx,t∈D¯zx,t from extreme values in calculus which implies ∂/∂xzx∗,t∗=0,∂/∂tzx∗,t∗=0,and∂2/∂x2zx∗,t∗≥0 giving that Lzx∗,t∗<0, which is contradiction to the assumption that made Lzx∗,t∗≥0,∀x,t∈D. Therefore, zx,t≥0,∀x,t∈D¯.Lemma 4.
Letux,t be the solution of the continuous problem in (2)–(4). Then, we obtain the bound
(13)ux,t≤ζ−1f+maxϕlt,ϕbx,t,ϕrt,where ζ is the lower bound of bx,t.Proof.
By defining the barrier functionsϑ±x,t as ϑ±x,t=ζ−1f+maxϕlt,ϕbx,t,ϕrt±ux,t and applying the maximum principle, we obtain the required bound.Lemma 5.
The bound on the derivative of the solutionux,t of the problem in (2)–(4) with respect to x and t is given by
(14)∂kux,t∂xk≤C1+ε−kexp−α1−xε,x,t∈D¯,0≤k≤4,∂lux,t∂tl≤C,x,t∈D¯,l=0,1,2,,where α is lower bound of ax.Proof.
See on [5, 16].


## 2.1. Bounds on the Solution and Its Derivatives

The existence and uniqueness of the solution of (2)–(4) can be established by assuming that the data is Holder continuous and imposing appropriate compatibility conditions at the corner points, using the assumptions of sufficiently smoothness ofϕlt,ϕbx,t, andϕrt.The required compatibility condition at the corner points and the delay term are(6)ϕl0=ϕb0,0,ϕr0=ϕb1,0,(7)dϕl0dt−ε∂2ϕb0,0∂x2=a0∂ϕb0,0∂x+b0,0ϕb0,0=−c0,0u0,−τ+f0,0,dϕr0dt−ε∂2ϕb1,0∂x2=a1∂ϕb1,0∂x+b1,0ϕb1,0=−c0,0u1,−τ+f1,0,so that the data matches at the corner points. For ax,bx,t,cx,t, and fx,t to be continuous on domain D, then (2)–(4) has the unique solution ux,t∈C2D∩C0D¯. In particular, when the compatibility conditions are not satisfied, a unique classical solution still exists but is not differentiable on all of ∂D.The reduced problem corresponding to the singularly perturbed parabolic delay PDE (2)–(4) is given as
(8)∂u0∂t+ax∂u0∂x+bx,tu0x,t=−cx,tu0x,t−τ+fx,t,x,t∈D,with the boundary conditions
(9)u00,t=ϕlt,t∈Ωt,u01,t≠ϕrt,t∈Ωt,and the interval condition
(10)u0x,t=ϕbx,t,x,t∈0,1×−τ,0.It is in the form of hyperbolic delay PDEs. The solutionux,t of (2)–(4) becomes very close to the solution u0x,t of (8)–(10) as ε⟶0.Lemma 2.
The solutionux,t of (2)–(4) satisfies the estimate
(11)ux,t−ϕbx,0≤Ct,x,t∈D¯,where C is a constant independent of ε.Proof.
The result follows from the compatibility condition. See the detailed proof in [16].LetL be a differential operator that denotes the differential equation in (2)–(4) as L=∂/∂t−ε∂2/∂x2+ax∂/∂x+bx,t.Lemma 3.
Suppose the functionzx,t∈C2D∩C0D¯, which satisfies zx,t≥0,x,t∈∂D, and Lzx,t≥0,x,t∈D implies that zx,t≥0,x,t∈D¯.Proof.
Assume there existsx∗,t∗∈D¯, such that
(12)zx∗,t∗=minx,t∈D¯zx,t<0.It is clear that the pointx∗,t∗∉∂D, which implies that x∗,t∗∈D. Since zx∗,t∗=minx,t∈D¯zx,t from extreme values in calculus which implies ∂/∂xzx∗,t∗=0,∂/∂tzx∗,t∗=0,and∂2/∂x2zx∗,t∗≥0 giving that Lzx∗,t∗<0, which is contradiction to the assumption that made Lzx∗,t∗≥0,∀x,t∈D. Therefore, zx,t≥0,∀x,t∈D¯.Lemma 4.
Letux,t be the solution of the continuous problem in (2)–(4). Then, we obtain the bound
(13)ux,t≤ζ−1f+maxϕlt,ϕbx,t,ϕrt,where ζ is the lower bound of bx,t.Proof.
By defining the barrier functionsϑ±x,t as ϑ±x,t=ζ−1f+maxϕlt,ϕbx,t,ϕrt±ux,t and applying the maximum principle, we obtain the required bound.Lemma 5.
The bound on the derivative of the solutionux,t of the problem in (2)–(4) with respect to x and t is given by
(14)∂kux,t∂xk≤C1+ε−kexp−α1−xε,x,t∈D¯,0≤k≤4,∂lux,t∂tl≤C,x,t∈D¯,l=0,1,2,,where α is lower bound of ax.Proof.
See on [5, 16].


## 3. Numerical Scheme

In general, for singularly perturbed problems, there are two strategies for designing numerical methods which have a small error in the boundary layer region [25]. The first approach is the class of fitted mesh methods which uses fine mesh in the boundary layer region and coarse mesh in outer layer region. The stability and convergence analyses of this approach are well developed. The second approach is the fitted operator methods in which it uses uniform mesh and an exponentially fitting factor for stabilizing the term containing the singular perturbation parameter. In this approach, the difference schemes reflect the qualitative behaviour of the solution inside the boundary layer region. In this article, we formulate an exponentially fitted operator finite difference scheme to solve the problem in (2)–(4).


### 3.1. Temporal Semidiscretization

The time domain0,T is discretized using a uniform mesh with time step Δt as ΩtM=tj=jΔt,j=0,1,⋯,M,tM=T,Δt=T/M andΩtm=tj=jΔt,j=0,1,⋯,m,tm=τ,Δt=τ/m where M is the number of mesh points in time direction in the interval 0,T and m is the number of mesh points in −τ,0. Note that T=kτ for some positive integer k.For approximating the temporal derivative term of (2)–(4), we use the averaged Crank-Nicolson method, which gives a system of boundary value problems
(15)1+Δt2LΔtUj+1x=1−Δt2LΔtUjx−cx,tj+1/2ϕbx+fx,tj+1/2,forj=0,1,⋯,m,x∈Ωx1−Δt2LΔtUjx−cx,tj+1/2Uj+1/2−mx+fx,tj+1/2,forj=m+1,⋯,M−1,x∈Ωx,with the discretized boundary conditions
(16)Uj+10=ϕltj+1,Uj+11=ϕrtj+1,j=0,1,2,⋯,M,where
(17)LΔtUj+1x=−εd2dx2Uj+1x+axddxUj+1x+bx,tj+1/2Uj+1x.Here,Uj+1x is denoted for the approximation of ux,tj+1 at the j+1th time level andbx,tj+1/2=1/2bx,tj+bx,tj+1, similarly for cx,tj+1/2andfx,tj+1/2.The semidiscrete difference operator1+Δt/2LΔtUj+1x in (15)–(16) satisfies the maximum principle as follows.Lemma 6.
Semidiscrete maximum principle. LetZj+1x be a sufficiently smooth function on Ω¯x. If Zj+10≥0,Zj+11≥0and1+Δt/2LΔtZj+1x≥0,∀x∈Ωx, then Zj+1x≥0,∀x∈Ω¯x.Proof.
Assume that there existx∗∈Ω¯x=0,1 such that Zj+1x∗=minx∈Ω¯xZj+1x<0. From the above assumption, it is clear that x∗∉0,1impliesthatx∗∈0,1. Since Zj+1x∗=minx∈Ω¯xZj+1x, using the property in calculus, we have d/dxZj+1x∗=0 and d2/dx2Zj+1x∗≥0, then we obtain 1+Δt/2LΔtZj+1x∗<0 which is in contradiction to 1+Δt/2LΔtZj+1x∗≥0,∀x∈Ωx. Therefore, we conclude Zj+1x≥0,∀x∈Ω¯x..Next, let us analyse the truncation error for the temporal discretization made above.Let the local error at each time step be denoted byej+1x≔ux,tj+1−Uj+1x,j=0,1,2,⋯,M.Lemma 7.
Suppose that(18)∂kux,t∂tk≤C,x,t∈D¯,k=0,1,2.The local truncation error in the temporal direction is given by(19)ej+1≤C1Δt3.Proof.
Using Taylor’s series approximation forux,tj and ux,tj+1 centering at tj+1/2, we obtain
(20)ux,tj=ux,tj+1/2−Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3,ux,tj+1=ux,tj+1/2+Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3.From the approximation in (20), we obtain
(21)ux,tj+1−ux,tjΔt=utx,tj+1/2+OΔt2.Using the approximation into (2), we obtain
(22)ux,tj+1−ux,tjΔt=εuxxx,tj+1/2−axuxx,tj+1/2−bx,tj+1/2ux,tj+1/2+fx,tj+1/2+OΔt2+−cx,tj+1/2ϕbx,forj=0,1,⋯,m,−cx,tj+1/2ux,tj+1/2−m,forj=m+1,⋯,M−1,where
(23)ux,tj+1/2=ux,tj+1+ux,tj2+OΔt2,fx,tj+1/2=fx,tj+1+fx,tj2+OΔt2.Since the errorej+1x≕ux,tj+1−Uj+1x satisfies the semidiscrete differential operator
(24)1+Δt2LΔtej+1x=OΔt3,ej+10=0=ej+11.Hence, by applying the maximum principle, we obtain(25)ej+1≤C1Δt3.Next, we need to show the bound for the global error of the temporal discretization. Let us denoteTEj+1 as the global error up to the j+1th time step.Lemma 8.
The global error attj+1 is given by
(26)TE≤CΔt2,j=1,2,⋯,M−1.Proof.
Using the local error up to thej+1th time step given in the above lemma, we obtain the global error at the j+1th time step as
(27)TE=∑l=1j+1el≤e1+e2+⋯+ej+1≤C1TΔt2,sincej+1Δt≤T=CΔt2,C1T=C,where C is a constant independent of εandΔt. See the detailed proof in [26].Next, we set a bound for the derivatives of solution of (15)–(16).Lemma 9.
For eachj=0,1,⋯,M−1, the solution Uj+1x of the boundary value problems in (15)–(16) satisfies the bound
(28)dkUj+1xdxk≤C1+ε−kexp−α1−xε,x∈Ω¯x,0≤k≤4.Proof.
See the proof in [5].


### 3.2. Spatial Discretization

The spatial domain0,1 is discretized into N equal number of subintervals, each of length h. Let 0=x0,xN=1, and xi=ih,i=0,1,2,⋯,N, be the mesh points. For spatial discretization, we apply an exponentially fitted operator finite difference method which helps us to hinder the influence of the singular perturbation parameter.First, let us find the exponential fitting factor for anonymous BVPs and then apply discretization in the spatial direction.


#### 3.2.1. Computing the Exponential Fitting Factor

To develop the numerical method for (15)–(16), we use the technique designed in the theory of asymptotic method for solving singularly perturbed BVPs. In the considered case, the boundary layer occurs on the right side of the domain. From the theory of singular perturbation in [27], the zero order asymptotic solution of the singularly perturbed boundary value problems of the form
(29)−εu″x+axu′x+bxux=fx,x∈Ωx=0,1,u0=al,u1=ar,is given by
(30)ux=u0x+a1axar−u01exp−∫x1axε−bxaxdx+Oε.Using Taylor’s series expansion forax and bx about x=1 and restriction to their first terms and simplifying gives
(31)ux=u0x+ar−u01exp−a11−xε+Oε,where u0 is the solution of the reduced problems (obtained by setting ε=0) in (29). Considering h is reasonably small and evaluating the result in (31) at xi gives
(32)uih=u00+ar−u01exp−a11ε−iρ.Similarly, we have(33)ui−1=u00+ar−u01exp−a11ε−i−1ρ,ui+1=u00+ar−u01exp−a11ε−i+1ρ,where ρ=h/ε.Consider a uniform gridΩ¯xN=xii=0N and denote h=xi+1−xi. For any mesh function vi, define the following difference operators:
(34)D+vi=vi+1−vih,D−vi=vi−vi−1h,D0vi=vi+1−vi−12h,D+D−vi=vi+1−2vi+vi−1h2.To handle the effect of the perturbation parameter, we multiply artificial viscosity (exponentially fitting factorσρ) on the diffusive part of the problem. Introducing an exponentially fitting factor σρ in (29) and applying the central finite difference scheme gives
(35)−εσρD+D−uxi+axiD0uxi+bxiuxi=fxi.Multiplying (35) by h and considering h is small and truncating the term hfxi−bxiuxi gives
(36)σρρui−1−2ui+ui+1+axi2ui+1−ui−1=0.From (32), we have
(37)ui±1=u00+ar−u01exp−a11ε−i±1ρ.Substituting (37) and (33) into (36) and simplifying, the exponential fitting factor is obtained as
(38)σρ=ρaxi2cothρa12.


#### 3.2.2. The Discrete Scheme

Using the central finite difference method for the spatial discretization of (15)–(16) and applying the exponential fitting factor in (38), for i=1,2,⋯,N−1, the fully discrete scheme becomes
(39)1+Δt2LΔt,hUi,j+1=1−Δt2LΔt,hUi,j−Δtcxi,tj+1/2ϕbxi+Δtfxi,tj+1/2,forj=0,1,⋯,m,1−Δt2LΔt,hUi,j−Δtcxi,tj+1/2Ui,j+1/2−m+Δtfxi,tj+1/2,forj=m+1,⋯,M−1,where
(40)LΔt,hUi,j+1=−σρεD+D−Ui,j+1+axiD0Ui,j+1+bxi,tj+1/2Ui,j+1.In explicit form, the scheme is rewritten as(41)ri−Ui−1,j+1+ricUi,j+1+ri+Ui+1,j+1=si−Ui−1,j+sicUi,j+si+Ui+1,j+−Δtcxi,tj+1/2ϕbxi+Δtfxi,tj+1/2,forj=0,1,⋯,m,−Δtcxi,tj+1/2Ui,j+1/2−m+Δtfxi,tj+1/2,forj=m+1,⋯,M−1,where
(42)ri−=−Δt2εσρh2+axi2h,ric=Δtεσρh2+12bxi,tj+1/2+1,ri+=−Δt2εσρh2−axi2h,(43)si−=Δt2εσρh2+axi2h,sic=1−Δtεσρh2+12bxi,tj+1/2,si+=Δt2εσρh2−axi2h.


### 3.3. Stability and Uniform Convergence Analysis

First, we need to prove the discrete comparison principle for the discrete scheme in (39).Lemma 10.
Discrete comparison principle. There exist a comparison functionZi,j+1 such that 1+Δt/2LΔt,hUi,j+1≤1+Δt/2LΔt,hZi,j+1,i=1,2,⋯,N−1 and if U0,j+1≤Z0,j+1 and UN,j+1≤ZN,j+1, then Ui,j+1≤Zi,j+1,i=0,1,2,⋯,N.Proof.
The matrix associated with operator1+Δt/2LΔt,hUi,j+1 is of size N+1×N+1 with its entries for i=1,2,⋯,N−1 are
(44)ai,i−1=−Δt4haxicothha12ε+1<0,ai,i+1=−Δt4haxicothha12ε−1≤0,sincecothha12ε≥1,ai,i=Δt2haxicothha12ε+bx,tj+1/2+1>0.So, the coefficient matrix satisfies the property ofM matrix. So, the inverse matrix exists and it is nonnegative. This guarantees the existence and uniqueness of the discrete solution. See the detailed proof in [28].Lemma 11.
LetZi,j+1=1+xi, for, 0≤i≤N. Then, there exist a positive constant C such that 1+Δt/2LΔt,hZi,j+1≥C,for1≤i≤N−1.Proof.
The proof is a simple computation, enables one to give a bound, that is uniform inε for the norm of the inverse of 1+Δt/2LΔt,h.
(45)1+Δt2LΔt,hZi,j+1=1−Δt2σρεD+D−Zi,j+1+Δt2axiD0Zi,j+1+Δt2bxi,tj+1/2Zi,j+1=1−Δt2σρεD+D−1+xi+Δt2axiD01+xi+Δt2bxi,tj+1/21+xi=1+xi−Δt2h2σρε1+xi−1−21+xi+1+xi+1+Δt4haxi1+xi+1−1+xi−1+Δt2bxi,tj+1/21+xi=1+xi+Δt2axi+Δt2bxi,tj+1/21+xi≥C,sinceax≥α>0,bxi,tj+1/2≥ζ>0.Lemma 12.
(uniform stability estimate). The solutionUi,j+1 of the discrete scheme in (39) satisfies the bound
(46)Ui,j+1≤1+Δt/2LΔt,hUi,j+11+Δtζ/2+Cmaxϕltj+1,ϕrtj+1,
wherebxi,tj+1/2≥ζ>0.Proof.
Let us construct a barrier functionπi,j+1± as πi,j+1±=1+Δt/2LΔt,hUi,j+1/1+Δtζ/2+Cmaxϕltj+1,ϕrtj+1±Ui,j+1. We can easily show thatπ0,j+1±≥0andπN,j+1±≥0. Then,
(47)1+Δt2LΔt,hπi,j+1±=πi,j+1±+Δt2−εσρD+D−πi,j+1±+axiD0πi,j+1±+bxi,tj+1/2πi,j+1±=1+Δt2bxi,tj+1/21+Δt/2LΔt,hUi,j+11+Δt/2ζ+Cmaxϕltj+1,ϕrtj+1±Ui,j+1±Δt2−εσρD+D−Ui,j+1+axiD0Ui,j+1+bxi,tj+1/2Ui,j+1=1+Δt2bxi,tj+1/21+Δt/2LΔt,hUi,j+11+Δt/2ζ+Cmaxϕltj+1,ϕrtj+1±1+Δt2LΔt,hUi,j+1≥0.
Using the discrete comparison principle, we obtainπi,j+1±≥0,i=0,1,2,⋯,N.Lemma 13.
IfZi,j+1 be any mesh function such that Z0,j+1=ZN,j+1=0.Then,
(48)Zi,j+1≤1ζmax1≤k≤N−1LΔt,hZk,j+1.Proof.
Consider two barrier functions of the formπi,j+1±=P±Zi,j+1, where P=1/ζmax1≤k≤N−1LΔt,hZk,j+1. It is easily shown that π0,j+1±≥0,πN,j+1±≥0. Next, we show that LΔt,hπi,j+1±≥0.
(49)LΔt,hπi,j+1±=Lh,ΔtP±Zi,j+1=±LΔt,hZi,j+1+bxi,tj+1/2ζmax1≤k≤N−1LΔt,hZk,j+1≥0.Hence, using the discrete comparison principle gives(50)Zi,j+1≤1ζmax1≤k≤N−1LΔt,hZk,j+1.Next, we consider the semidiscrete problem in (15)–(16) and the fully discrete scheme in (39) to find the truncation error of the spatial direction discretization.Theorem 14.
Let the coefficient functionsax and bx,tj+1/2 in (15)–(16) be sufficiently smooth functions so that Uj+1x∈C40,1. Then, the computed solution Ui,j+1 of the problem in (39) satisfies the bound
(51)LΔt,hUj+1xi−Ui,j+1≤Ch2h+ε1+ε−3exp−α1−ηε,where η∈x0,xN such that exp−α1−η/ε=maxxiexp−α1−xi/ε.Proof.
The local truncation error in space discretization is given as(52)LΔt,hUj+1xi−Ui,j+1=−εσρd2dx2−D+D−Uj+1xi+axiddx−D0Uj+1xi≤εaxiρ2cotha1ρ2−1D+D−Uj+1xi+εd2dx2−D+D−Uj+1xi+axiddx−D0Uj+1xi,where σρ=axiρ/2cotha1ρ/2andρ=h/ε.Now forρ>0, C1 and C2 are constants and we have ρcothρ−1≤C1ρ2,forρ≤1. For ρ⟶∞, since limρ→∞cothρ=1 gives ρcothρ−1≤C1ρ. In general, for all ρ>0, we write
(53)C1ρ2ρ+1≤ρcothρ−1≤C2ρ2ρ+1,giving that
(54)εaxiρ2cotha1ρ2−1≤εh/ε2h/ε+1=h2h+ε.Using Taylor series expansion, we obtain the bound for(55)D+D−Uj+1xi≤CUj+1′′η,ddx−D0Uj+1xi≤Ch2Uj+1′′′η,d2dx2−D+D−Uj+1xi≤Ch2Uj+14η,where Uj+1kη=maxη∈x0,xNUj+1kxi.Using the bounds for the differences of the derivatives in (54) and (55), we obtainLΔt,hUj+1xi−Ui,j+1≤Ch2/h+εUj+1′′η+Cεh2Uj+14η+Ch2Uj+1′′′η≤Ch2/h+εUj+1′′η+Ch2εUj+14η+Uj+1′′′η.Using the bounds for the derivatives of the solution in Lemma 9 givesLΔt,hUj+1xi−Ui,j+1≤Ch2/h+ε1+ε−2exp−α1−η/ε+Ch2ε+ε−3exp−α1−η/ε+1+ε−3exp−α1−η/ε. Since ε−3≥ε−2, we obtain
(56)LΔt,hUj+1xi−Ui,j+1≤Ch2h+ε1+ε−3exp−α1−ηε,which gives the required bound.Lemma 15.
For a fixed number of mesh numbersN, and for ε⟶0, it holds
(57)limε→0max1≤i≤N−1exp−αxi/εεm=0,limε→0max1≤i≤N−1exp−α1−xi/εεm=0,m=1,2,3,⋯,where xi=ih,∀i=1,2,⋯,N−1..Proof.
On the discrete domainxi0N, for the interior grid points, we have
(58)max1≤i≤N−1exp−αxi/εεm≤exp−αx1/εεm=exp−αh/εεm,max1≤i≤N−1exp−α1−xi/εεm≤exp−α1−xN−1/εεm=exp−αh/εεm,Sincex1=h,1−xN−1=h. Repeatedly using L’Hospital’s rule gives
(59)limε→0exp−αh/εεm=limζ=1/ε→∞ζmexpαhζ=limζ=1/ε→∞m!αhmexpαhζ=0.This completes the proof.Theorem 16.
The numerical solutionUi,j+1 of the problem in (39) satisfies the following uniform error bound
(60)sup0<ε≤1max0≤i≤NUj+1xi−Ui,j+1≤Ch.Proof.
Using Lemma15 into (56) gives
(61)LΔt,hUj+1xi−Ui,j+1≤Ch2h+ε.Hence, using the result in Lemma13, we obtain
(62)Uj+1xi−Ui,j+1≤Ch2h+ε.Using the supremum over allε∈0,1, we obtain
(63)sup0<ε≤1max0≤i≤NUj+1xi−Ui,j+1≤Ch.Remark 17.
For the caseε>h, the developed scheme gives second order convergence. For small values of ε≪h, the scheme is first order uniformly convergent in spatial direction.Theorem 18.
LetuandU be the solution of (2)–(4) and (39), respectively. Then, the following uniform error bound holds:
(64)sup0<ε≤1u−U≤Ch+Δt2.Proof.
The combination of temporal and spatial error bounds gives the required result.


## 3.1. Temporal Semidiscretization

The time domain0,T is discretized using a uniform mesh with time step Δt as ΩtM=tj=jΔt,j=0,1,⋯,M,tM=T,Δt=T/M andΩtm=tj=jΔt,j=0,1,⋯,m,tm=τ,Δt=τ/m where M is the number of mesh points in time direction in the interval 0,T and m is the number of mesh points in −τ,0. Note that T=kτ for some positive integer k.For approximating the temporal derivative term of (2)–(4), we use the averaged Crank-Nicolson method, which gives a system of boundary value problems
(15)1+Δt2LΔtUj+1x=1−Δt2LΔtUjx−cx,tj+1/2ϕbx+fx,tj+1/2,forj=0,1,⋯,m,x∈Ωx1−Δt2LΔtUjx−cx,tj+1/2Uj+1/2−mx+fx,tj+1/2,forj=m+1,⋯,M−1,x∈Ωx,with the discretized boundary conditions
(16)Uj+10=ϕltj+1,Uj+11=ϕrtj+1,j=0,1,2,⋯,M,where
(17)LΔtUj+1x=−εd2dx2Uj+1x+axddxUj+1x+bx,tj+1/2Uj+1x.Here,Uj+1x is denoted for the approximation of ux,tj+1 at the j+1th time level andbx,tj+1/2=1/2bx,tj+bx,tj+1, similarly for cx,tj+1/2andfx,tj+1/2.The semidiscrete difference operator1+Δt/2LΔtUj+1x in (15)–(16) satisfies the maximum principle as follows.Lemma 6.
Semidiscrete maximum principle. LetZj+1x be a sufficiently smooth function on Ω¯x. If Zj+10≥0,Zj+11≥0and1+Δt/2LΔtZj+1x≥0,∀x∈Ωx, then Zj+1x≥0,∀x∈Ω¯x.Proof.
Assume that there existx∗∈Ω¯x=0,1 such that Zj+1x∗=minx∈Ω¯xZj+1x<0. From the above assumption, it is clear that x∗∉0,1impliesthatx∗∈0,1. Since Zj+1x∗=minx∈Ω¯xZj+1x, using the property in calculus, we have d/dxZj+1x∗=0 and d2/dx2Zj+1x∗≥0, then we obtain 1+Δt/2LΔtZj+1x∗<0 which is in contradiction to 1+Δt/2LΔtZj+1x∗≥0,∀x∈Ωx. Therefore, we conclude Zj+1x≥0,∀x∈Ω¯x..Next, let us analyse the truncation error for the temporal discretization made above.Let the local error at each time step be denoted byej+1x≔ux,tj+1−Uj+1x,j=0,1,2,⋯,M.Lemma 7.
Suppose that(18)∂kux,t∂tk≤C,x,t∈D¯,k=0,1,2.The local truncation error in the temporal direction is given by(19)ej+1≤C1Δt3.Proof.
Using Taylor’s series approximation forux,tj and ux,tj+1 centering at tj+1/2, we obtain
(20)ux,tj=ux,tj+1/2−Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3,ux,tj+1=ux,tj+1/2+Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3.From the approximation in (20), we obtain
(21)ux,tj+1−ux,tjΔt=utx,tj+1/2+OΔt2.Using the approximation into (2), we obtain
(22)ux,tj+1−ux,tjΔt=εuxxx,tj+1/2−axuxx,tj+1/2−bx,tj+1/2ux,tj+1/2+fx,tj+1/2+OΔt2+−cx,tj+1/2ϕbx,forj=0,1,⋯,m,−cx,tj+1/2ux,tj+1/2−m,forj=m+1,⋯,M−1,where
(23)ux,tj+1/2=ux,tj+1+ux,tj2+OΔt2,fx,tj+1/2=fx,tj+1+fx,tj2+OΔt2.Since the errorej+1x≕ux,tj+1−Uj+1x satisfies the semidiscrete differential operator
(24)1+Δt2LΔtej+1x=OΔt3,ej+10=0=ej+11.Hence, by applying the maximum principle, we obtain(25)ej+1≤C1Δt3.Next, we need to show the bound for the global error of the temporal discretization. Let us denoteTEj+1 as the global error up to the j+1th time step.Lemma 8.
The global error attj+1 is given by
(26)TE≤CΔt2,j=1,2,⋯,M−1.Proof.
Using the local error up to thej+1th time step given in the above lemma, we obtain the global error at the j+1th time step as
(27)TE=∑l=1j+1el≤e1+e2+⋯+ej+1≤C1TΔt2,sincej+1Δt≤T=CΔt2,C1T=C,where C is a constant independent of εandΔt. See the detailed proof in [26].Next, we set a bound for the derivatives of solution of (15)–(16).Lemma 9.
For eachj=0,1,⋯,M−1, the solution Uj+1x of the boundary value problems in (15)–(16) satisfies the bound
(28)dkUj+1xdxk≤C1+ε−kexp−α1−xε,x∈Ω¯x,0≤k≤4.Proof.
See the proof in [5].


## 3.2. Spatial Discretization

The spatial domain0,1 is discretized into N equal number of subintervals, each of length h. Let 0=x0,xN=1, and xi=ih,i=0,1,2,⋯,N, be the mesh points. For spatial discretization, we apply an exponentially fitted operator finite difference method which helps us to hinder the influence of the singular perturbation parameter.First, let us find the exponential fitting factor for anonymous BVPs and then apply discretization in the spatial direction.


### 3.2.1. Computing the Exponential Fitting Factor

To develop the numerical method for (15)–(16), we use the technique designed in the theory of asymptotic method for solving singularly perturbed BVPs. In the considered case, the boundary layer occurs on the right side of the domain. From the theory of singular perturbation in [27], the zero order asymptotic solution of the singularly perturbed boundary value problems of the form
(29)−εu″x+axu′x+bxux=fx,x∈Ωx=0,1,u0=al,u1=ar,is given by
(30)ux=u0x+a1axar−u01exp−∫x1axε−bxaxdx+Oε.Using Taylor’s series expansion forax and bx about x=1 and restriction to their first terms and simplifying gives
(31)ux=u0x+ar−u01exp−a11−xε+Oε,where u0 is the solution of the reduced problems (obtained by setting ε=0) in (29). Considering h is reasonably small and evaluating the result in (31) at xi gives
(32)uih=u00+ar−u01exp−a11ε−iρ.Similarly, we have(33)ui−1=u00+ar−u01exp−a11ε−i−1ρ,ui+1=u00+ar−u01exp−a11ε−i+1ρ,where ρ=h/ε.Consider a uniform gridΩ¯xN=xii=0N and denote h=xi+1−xi. For any mesh function vi, define the following difference operators:
(34)D+vi=vi+1−vih,D−vi=vi−vi−1h,D0vi=vi+1−vi−12h,D+D−vi=vi+1−2vi+vi−1h2.To handle the effect of the perturbation parameter, we multiply artificial viscosity (exponentially fitting factorσρ) on the diffusive part of the problem. Introducing an exponentially fitting factor σρ in (29) and applying the central finite difference scheme gives
(35)−εσρD+D−uxi+axiD0uxi+bxiuxi=fxi.Multiplying (35) by h and considering h is small and truncating the term hfxi−bxiuxi gives
(36)σρρui−1−2ui+ui+1+axi2ui+1−ui−1=0.From (32), we have
(37)ui±1=u00+ar−u01exp−a11ε−i±1ρ.Substituting (37) and (33) into (36) and simplifying, the exponential fitting factor is obtained as
(38)σρ=ρaxi2cothρa12.


### 3.2.2. The Discrete Scheme

Using the central finite difference method for the spatial discretization of (15)–(16) and applying the exponential fitting factor in (38), for i=1,2,⋯,N−1, the fully discrete scheme becomes
(39)1+Δt2LΔt,hUi,j+1=1−Δt2LΔt,hUi,j−Δtcxi,tj+1/2ϕbxi+Δtfxi,tj+1/2,forj=0,1,⋯,m,1−Δt2LΔt,hUi,j−Δtcxi,tj+1/2Ui,j+1/2−m+Δtfxi,tj+1/2,forj=m+1,⋯,M−1,where
(40)LΔt,hUi,j+1=−σρεD+D−Ui,j+1+axiD0Ui,j+1+bxi,tj+1/2Ui,j+1.In explicit form, the scheme is rewritten as(41)ri−Ui−1,j+1+ricUi,j+1+ri+Ui+1,j+1=si−Ui−1,j+sicUi,j+si+Ui+1,j+−Δtcxi,tj+1/2ϕbxi+Δtfxi,tj+1/2,forj=0,1,⋯,m,−Δtcxi,tj+1/2Ui,j+1/2−m+Δtfxi,tj+1/2,forj=m+1,⋯,M−1,where
(42)ri−=−Δt2εσρh2+axi2h,ric=Δtεσρh2+12bxi,tj+1/2+1,ri+=−Δt2εσρh2−axi2h,(43)si−=Δt2εσρh2+axi2h,sic=1−Δtεσρh2+12bxi,tj+1/2,si+=Δt2εσρh2−axi2h.


## 3.2.1. Computing the Exponential Fitting Factor

To develop the numerical method for (15)–(16), we use the technique designed in the theory of asymptotic method for solving singularly perturbed BVPs. In the considered case, the boundary layer occurs on the right side of the domain. From the theory of singular perturbation in [27], the zero order asymptotic solution of the singularly perturbed boundary value problems of the form
(29)−εu″x+axu′x+bxux=fx,x∈Ωx=0,1,u0=al,u1=ar,is given by
(30)ux=u0x+a1axar−u01exp−∫x1axε−bxaxdx+Oε.Using Taylor’s series expansion forax and bx about x=1 and restriction to their first terms and simplifying gives
(31)ux=u0x+ar−u01exp−a11−xε+Oε,where u0 is the solution of the reduced problems (obtained by setting ε=0) in (29). Considering h is reasonably small and evaluating the result in (31) at xi gives
(32)uih=u00+ar−u01exp−a11ε−iρ.Similarly, we have(33)ui−1=u00+ar−u01exp−a11ε−i−1ρ,ui+1=u00+ar−u01exp−a11ε−i+1ρ,where ρ=h/ε.Consider a uniform gridΩ¯xN=xii=0N and denote h=xi+1−xi. For any mesh function vi, define the following difference operators:
(34)D+vi=vi+1−vih,D−vi=vi−vi−1h,D0vi=vi+1−vi−12h,D+D−vi=vi+1−2vi+vi−1h2.To handle the effect of the perturbation parameter, we multiply artificial viscosity (exponentially fitting factorσρ) on the diffusive part of the problem. Introducing an exponentially fitting factor σρ in (29) and applying the central finite difference scheme gives
(35)−εσρD+D−uxi+axiD0uxi+bxiuxi=fxi.Multiplying (35) by h and considering h is small and truncating the term hfxi−bxiuxi gives
(36)σρρui−1−2ui+ui+1+axi2ui+1−ui−1=0.From (32), we have
(37)ui±1=u00+ar−u01exp−a11ε−i±1ρ.Substituting (37) and (33) into (36) and simplifying, the exponential fitting factor is obtained as
(38)σρ=ρaxi2cothρa12.


## 3.2.2. The Discrete Scheme

Using the central finite difference method for the spatial discretization of (15)–(16) and applying the exponential fitting factor in (38), for i=1,2,⋯,N−1, the fully discrete scheme becomes
(39)1+Δt2LΔt,hUi,j+1=1−Δt2LΔt,hUi,j−Δtcxi,tj+1/2ϕbxi+Δtfxi,tj+1/2,forj=0,1,⋯,m,1−Δt2LΔt,hUi,j−Δtcxi,tj+1/2Ui,j+1/2−m+Δtfxi,tj+1/2,forj=m+1,⋯,M−1,where
(40)LΔt,hUi,j+1=−σρεD+D−Ui,j+1+axiD0Ui,j+1+bxi,tj+1/2Ui,j+1.In explicit form, the scheme is rewritten as(41)ri−Ui−1,j+1+ricUi,j+1+ri+Ui+1,j+1=si−Ui−1,j+sicUi,j+si+Ui+1,j+−Δtcxi,tj+1/2ϕbxi+Δtfxi,tj+1/2,forj=0,1,⋯,m,−Δtcxi,tj+1/2Ui,j+1/2−m+Δtfxi,tj+1/2,forj=m+1,⋯,M−1,where
(42)ri−=−Δt2εσρh2+axi2h,ric=Δtεσρh2+12bxi,tj+1/2+1,ri+=−Δt2εσρh2−axi2h,(43)si−=Δt2εσρh2+axi2h,sic=1−Δtεσρh2+12bxi,tj+1/2,si+=Δt2εσρh2−axi2h.


## 3.3. Stability and Uniform Convergence Analysis

First, we need to prove the discrete comparison principle for the discrete scheme in (39).Lemma 10.
Discrete comparison principle. There exist a comparison functionZi,j+1 such that 1+Δt/2LΔt,hUi,j+1≤1+Δt/2LΔt,hZi,j+1,i=1,2,⋯,N−1 and if U0,j+1≤Z0,j+1 and UN,j+1≤ZN,j+1, then Ui,j+1≤Zi,j+1,i=0,1,2,⋯,N.Proof.
The matrix associated with operator1+Δt/2LΔt,hUi,j+1 is of size N+1×N+1 with its entries for i=1,2,⋯,N−1 are
(44)ai,i−1=−Δt4haxicothha12ε+1<0,ai,i+1=−Δt4haxicothha12ε−1≤0,sincecothha12ε≥1,ai,i=Δt2haxicothha12ε+bx,tj+1/2+1>0.So, the coefficient matrix satisfies the property ofM matrix. So, the inverse matrix exists and it is nonnegative. This guarantees the existence and uniqueness of the discrete solution. See the detailed proof in [28].Lemma 11.
LetZi,j+1=1+xi, for, 0≤i≤N. Then, there exist a positive constant C such that 1+Δt/2LΔt,hZi,j+1≥C,for1≤i≤N−1.Proof.
The proof is a simple computation, enables one to give a bound, that is uniform inε for the norm of the inverse of 1+Δt/2LΔt,h.
(45)1+Δt2LΔt,hZi,j+1=1−Δt2σρεD+D−Zi,j+1+Δt2axiD0Zi,j+1+Δt2bxi,tj+1/2Zi,j+1=1−Δt2σρεD+D−1+xi+Δt2axiD01+xi+Δt2bxi,tj+1/21+xi=1+xi−Δt2h2σρε1+xi−1−21+xi+1+xi+1+Δt4haxi1+xi+1−1+xi−1+Δt2bxi,tj+1/21+xi=1+xi+Δt2axi+Δt2bxi,tj+1/21+xi≥C,sinceax≥α>0,bxi,tj+1/2≥ζ>0.Lemma 12.
(uniform stability estimate). The solutionUi,j+1 of the discrete scheme in (39) satisfies the bound
(46)Ui,j+1≤1+Δt/2LΔt,hUi,j+11+Δtζ/2+Cmaxϕltj+1,ϕrtj+1,
wherebxi,tj+1/2≥ζ>0.Proof.
Let us construct a barrier functionπi,j+1± as πi,j+1±=1+Δt/2LΔt,hUi,j+1/1+Δtζ/2+Cmaxϕltj+1,ϕrtj+1±Ui,j+1. We can easily show thatπ0,j+1±≥0andπN,j+1±≥0. Then,
(47)1+Δt2LΔt,hπi,j+1±=πi,j+1±+Δt2−εσρD+D−πi,j+1±+axiD0πi,j+1±+bxi,tj+1/2πi,j+1±=1+Δt2bxi,tj+1/21+Δt/2LΔt,hUi,j+11+Δt/2ζ+Cmaxϕltj+1,ϕrtj+1±Ui,j+1±Δt2−εσρD+D−Ui,j+1+axiD0Ui,j+1+bxi,tj+1/2Ui,j+1=1+Δt2bxi,tj+1/21+Δt/2LΔt,hUi,j+11+Δt/2ζ+Cmaxϕltj+1,ϕrtj+1±1+Δt2LΔt,hUi,j+1≥0.
Using the discrete comparison principle, we obtainπi,j+1±≥0,i=0,1,2,⋯,N.Lemma 13.
IfZi,j+1 be any mesh function such that Z0,j+1=ZN,j+1=0.Then,
(48)Zi,j+1≤1ζmax1≤k≤N−1LΔt,hZk,j+1.Proof.
Consider two barrier functions of the formπi,j+1±=P±Zi,j+1, where P=1/ζmax1≤k≤N−1LΔt,hZk,j+1. It is easily shown that π0,j+1±≥0,πN,j+1±≥0. Next, we show that LΔt,hπi,j+1±≥0.
(49)LΔt,hπi,j+1±=Lh,ΔtP±Zi,j+1=±LΔt,hZi,j+1+bxi,tj+1/2ζmax1≤k≤N−1LΔt,hZk,j+1≥0.Hence, using the discrete comparison principle gives(50)Zi,j+1≤1ζmax1≤k≤N−1LΔt,hZk,j+1.Next, we consider the semidiscrete problem in (15)–(16) and the fully discrete scheme in (39) to find the truncation error of the spatial direction discretization.Theorem 14.
Let the coefficient functionsax and bx,tj+1/2 in (15)–(16) be sufficiently smooth functions so that Uj+1x∈C40,1. Then, the computed solution Ui,j+1 of the problem in (39) satisfies the bound
(51)LΔt,hUj+1xi−Ui,j+1≤Ch2h+ε1+ε−3exp−α1−ηε,where η∈x0,xN such that exp−α1−η/ε=maxxiexp−α1−xi/ε.Proof.
The local truncation error in space discretization is given as(52)LΔt,hUj+1xi−Ui,j+1=−εσρd2dx2−D+D−Uj+1xi+axiddx−D0Uj+1xi≤εaxiρ2cotha1ρ2−1D+D−Uj+1xi+εd2dx2−D+D−Uj+1xi+axiddx−D0Uj+1xi,where σρ=axiρ/2cotha1ρ/2andρ=h/ε.Now forρ>0, C1 and C2 are constants and we have ρcothρ−1≤C1ρ2,forρ≤1. For ρ⟶∞, since limρ→∞cothρ=1 gives ρcothρ−1≤C1ρ. In general, for all ρ>0, we write
(53)C1ρ2ρ+1≤ρcothρ−1≤C2ρ2ρ+1,giving that
(54)εaxiρ2cotha1ρ2−1≤εh/ε2h/ε+1=h2h+ε.Using Taylor series expansion, we obtain the bound for(55)D+D−Uj+1xi≤CUj+1′′η,ddx−D0Uj+1xi≤Ch2Uj+1′′′η,d2dx2−D+D−Uj+1xi≤Ch2Uj+14η,where Uj+1kη=maxη∈x0,xNUj+1kxi.Using the bounds for the differences of the derivatives in (54) and (55), we obtainLΔt,hUj+1xi−Ui,j+1≤Ch2/h+εUj+1′′η+Cεh2Uj+14η+Ch2Uj+1′′′η≤Ch2/h+εUj+1′′η+Ch2εUj+14η+Uj+1′′′η.Using the bounds for the derivatives of the solution in Lemma 9 givesLΔt,hUj+1xi−Ui,j+1≤Ch2/h+ε1+ε−2exp−α1−η/ε+Ch2ε+ε−3exp−α1−η/ε+1+ε−3exp−α1−η/ε. Since ε−3≥ε−2, we obtain
(56)LΔt,hUj+1xi−Ui,j+1≤Ch2h+ε1+ε−3exp−α1−ηε,which gives the required bound.Lemma 15.
For a fixed number of mesh numbersN, and for ε⟶0, it holds
(57)limε→0max1≤i≤N−1exp−αxi/εεm=0,limε→0max1≤i≤N−1exp−α1−xi/εεm=0,m=1,2,3,⋯,where xi=ih,∀i=1,2,⋯,N−1..Proof.
On the discrete domainxi0N, for the interior grid points, we have
(58)max1≤i≤N−1exp−αxi/εεm≤exp−αx1/εεm=exp−αh/εεm,max1≤i≤N−1exp−α1−xi/εεm≤exp−α1−xN−1/εεm=exp−αh/εεm,Sincex1=h,1−xN−1=h. Repeatedly using L’Hospital’s rule gives
(59)limε→0exp−αh/εεm=limζ=1/ε→∞ζmexpαhζ=limζ=1/ε→∞m!αhmexpαhζ=0.This completes the proof.Theorem 16.
The numerical solutionUi,j+1 of the problem in (39) satisfies the following uniform error bound
(60)sup0<ε≤1max0≤i≤NUj+1xi−Ui,j+1≤Ch.Proof.
Using Lemma15 into (56) gives
(61)LΔt,hUj+1xi−Ui,j+1≤Ch2h+ε.Hence, using the result in Lemma13, we obtain
(62)Uj+1xi−Ui,j+1≤Ch2h+ε.Using the supremum over allε∈0,1, we obtain
(63)sup0<ε≤1max0≤i≤NUj+1xi−Ui,j+1≤Ch.Remark 17.
For the caseε>h, the developed scheme gives second order convergence. For small values of ε≪h, the scheme is first order uniformly convergent in spatial direction.Theorem 18.
LetuandU be the solution of (2)–(4) and (39), respectively. Then, the following uniform error bound holds:
(64)sup0<ε≤1u−U≤Ch+Δt2.Proof.
The combination of temporal and spatial error bounds gives the required result.


## 4. Numerical Examples, Results, and Discussions

Here, we illustrate the proposed scheme using model examples. The exact solutions of the considered examples are not known. We investigate the theoretical results by performing experiments using the proposed scheme.Example 19.
Consider singularly perturbed time delay parabolic PDEs:(65)∂u∂t−ε∂2u∂x2+2−x2∂u∂x+xux,t+ux,t−τ=10t2exp−tx1−x,with interval condition ux,t=0,onx,t∈0,1×−1,0 and the boundary conditions
(66)u0,t=ϕlt=0,u1,t=ϕrt=0,fort∈0,2.Example 20.
Consider singularly perturbed time delay parabolic PDEs:(67)∂u∂t−ε∂2u∂x2+2−x2∂u∂x+x+1t+1ux,t+ux,t−τ=10t2exp−tx1−x,with interval condition ux,t=0,onx,t∈0,1×−1,0 and the boundary conditions
(68)u0,t=ϕlt=0,u1,t=ϕrt=0,fort∈0,2.Since the exact solutions of the considered examples are not known, with the help of double mesh techniques, we compute the maximum absolute errorEεN,M of the scheme. Let Ui,jN,M be the computed solution of the problem using N,M mesh points and let Ui,j2N,2M be the computed solution on double number of mesh points 2N,2M by including the midpoints xi+1/2=xi+1+xi/2 and tj+1/2=tj+1+tj/2 into the mesh points. We calculate the maximum absolute error as
(69)EεN,M=maxi,jUi,jN,M−Ui,j2N,2M,and the ε-uniform (parameter uniform) error estimate by using
(70)EN,M=maxεEεN,M.We calculate the rate of convergence of the proposed scheme using(71)rεN,M=log2EεN,M−log2Eε2N,2M.and the ε-uniform rate of convergence using
(72)rN,M=log2EN,M−log2E2N,2M.In Figures1 and 2, the numerical solution of Examples 19 and 20 is given, respectively, for different values of perturbation parameterε. As we observe on these figures, a strong boundary layer is maintained on the right side of the spatial domain as ε goes small starting from ε=2−5to2−20. In these figures, we observe that the proposed scheme approximates the solution without creating oscillations or divergence. Note that time delay does not have an effect on position and size of the boundary layer since the layer occurs on the spatial domain direction. If the delay was on the spatial variable, we may have an interior layer which is created because of the delay parameter. In Tables 1 and 2, the maximum absolute error of Examples 19 and 20 is given, respectively, for different values of perturbation parameter ε and mesh numbers N,M. As we observed the results in these tables for each number of mesh N,M as ε becomes small, the maximum absolute error after showing growth becomes stable and identical. This indicates the stability and uniform convergence of the scheme. In the last two rows of these tables, the ε-uniform error and ε-uniform rate of convergence of the scheme are given. As one observes, the scheme gives first order uniform convergence. The second order convergence of the temporal discretization is depicted in Table 3. In Tables 4 and 5, comparison of the ε-uniform error and ε-uniform rate of convergence of the proposed scheme with results of some recently published papers is given. As we observed, the proposed scheme gives a more accurate result than the results in [16, 18–20] and [24]. The proposed scheme has a limitation for solving nonlinear singularly perturbed problems.Figure 1
Solution of Example19, with boundary layer formation.Figure 2
Solution of Example20, with boundary layer formation.Table 1
Example19, maximum absolute errors of the scheme.
ε↓N=M=16326412825651211.0742e-046.6552e-054.4167e-052.7392e-051.5041e-057.8477e-062−23.8959e-043.3429e-042.5726e-041.5142e-048.1424e-054.2142e-052−41.6240e-035.6521e-044.3969e-042.6289e-041.4270e-047.4226e-052−68.3743e-032.0823e-038.0261e-043.8986e-041.9220e-049.5980e-052−81.0404e-024.8511e-031.9751e-037.6437e-043.0100e-041.2812e-042−101.0409e-024.9874e-032.6147e-031.3838e-035.8980e-042.2292e-042−121.0409e-024.9874e-032.6159e-031.4219e-037.3950e-043.6693e-042−141.0409e-024.9874e-032.6159e-031.4219e-037.3985e-043.7720e-042−161.0409e-024.9874e-032.6159e-031.4219e-037.3985e-043.7720e-042−181.0409e-024.9874e-032.6159e-031.4219e-037.3985e-043.7720e-042−201.0409e-024.9874e-032.6159e-031.4219e-037.3985e-043.7720e-04EN,M1.0409e-024.9874e-032.6159e-031.4219e-037.3985e-043.7720e-04rN,M1.06150.93100.87950.94250.9719—Table 2
Example20, maximum absolute errors of the scheme.
ε↓N=163264128256512M=1020408016032019.0612e-057.8840e-055.0913e-052.8361e-051.4905e-057.6151e-062−25.0971e-043.4843e-041.9886e-041.0569e-045.4425e-052.7604e-052−41.3108e-036.8037e-043.4323e-041.7164e-048.5788e-054.2882e-052−64.2800e-031.6597e-036.2601e-042.5497e-041.1234e-045.2343e-052−85.4803e-033.2949e-031.4175e-035.1091e-041.8211e-047.0533e-052−105.4832e-033.3921e-031.8283e-039.1599e-043.8114e-041.3511e-042−125.4832e-033.3921e-031.8292e-039.4367e-044.7807e-042.3354e-042−145.4832e-033.3921e-031.8292e-039.4367e-044.7832e-042.4067e-042−165.4832e-033.3921e-031.8292e-039.4367e-044.7832e-042.4067e-042−185.4832e-033.3921e-031.8292e-039.4367e-044.7832e-042.4067e-042−205.4832e-033.3921e-031.8292e-039.4367e-044.7832e-042.4067e-04EN,M5.4832e-033.3921e-031.8292e-039.4367e-044.7832e-042.4067e-04rN,M0.69280.89100.95490.98030.9909—Table 3
Maximum absolute error and rate of convergence of the scheme.
Example19Example20ε↓N=1664256N=1664256M=163264M=1020402−101.0409e-022.6147e-035.8980e-045.4832e-031.8283e-033.8004e-041.99312.1483—1.58452.0621—2−121.0409e-022.6159e-037.3950e-045.4832e-031.8292e-034.7807e-041.99251.8227—1.58381.9359—2−141.0409e-022.6159e-037.3985e-045.4832e-031.8292e-034.7832e-041.99251.8220—1.58381.9352—2−161.0409e-022.6159e-037.3985e-045.4832e-031.8292e-034.7832e-041.99251.8220—1.58381.9352—2−181.0409e-022.6159e-037.3985e-045.4832e-031.8292e-034.7832e-041.99251.8220—1.58381.9352—2−201.0409e-022.6159e-037.3985e-045.4832e-031.8292e-034.7832e-041.99251.8220—1.58381.9352—Table 4
Comparison of uniform error and uniform rate of convergence of Example19.
Schemes↓N=M=163264128256Proposed schemeEN,M1.0409e-024.9874e-032.6159e-031.4219e-037.3985e-04rN,M1.06150.93100.87950.94250.9719Scheme in [19]EN,M3.41e-021.84e-029.38e-034.67e-032.31e-03rN,M0.89010.97201.00621.01551.0063Scheme in [16]EN,M4.9485e-023.3203e-022.1165e-021.3320e-027.9345e-03rN,M0.57570.64960.66810.74740.7908Table 5
Comparison of uniform error and uniform rate of convergence of Example20.
Schemes↓N=163264128256M=10204080160Proposed schemeEN,M5.4832e-033.3921e-031.8292e-039.4367e-044.7832e-04rN,M0.69280.89100.95490.98030.9909Scheme in [20]EN,M1.86e-021.00e-025.48e-032.86e-031.46e-03rN,M0.890.870.940.971.11Scheme in [18]EN,M1.6119e-029.9504e-035.8541e-033.3439e-031.8650e-03rN,M0.69600.76530.80790.84240.8660Scheme in [24]EN,M7.4252e-034.0993e-032.1528e-031.1033e-035.5845e-04rN,M0.85700.92910.96440.9822—


## 5. Conclusion

In this paper, singularly perturbed parabolic convection-diffusion-reaction equation with large time delay is considered. The solution of the considered problem exhibits boundary layer of thicknessε on the right side of the spatial domain. The bounds and properties of the analytical solution are discussed. To stabilize the influence of the singular perturbation on the discrete solution, and exponential fitting parameter is developed using the derivation from the asymptotic solution. The numerical scheme is developed using the Crank-Nicolson method in temporal discretization and an exponentially fitted central finite difference method on the spatial discretization. The existence of the discrete solution is discussed using the comparison principle. The stability and uniform convergence of the scheme are investigated well theoretically. Numerical results are depicted using maximum absolute error, ε-uniform error, and ε-uniform rate of convergence in tables which are in good agreement with the theoretical analysis. The developed scheme gives stable and uniformly convergent result with linear order of convergence.



---

*Source: 6641236-2021-02-28.xml*