# Uniformly Convergent Numerical Scheme for Singularly Perturbed Parabolic PDEs with Shift Parameters

Authors: Mesfin Mekuria Woldaregay; Gemechis File Duressa
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this article, singularly perturbed parabolic differential difference equations are considered. The solution of the equations exhibits a boundary layer on the right side of the spatial domain. The terms containing the advance and delay parameters are approximated using Taylor series approximation. The resulting singularly perturbed parabolic PDEs are solved using the Crank–Nicolson method in the temporal discretization and nonstandard finite difference method in the spatial discretization. The existence of a unique discrete solution is guaranteed using the discrete maximum principle. The uniform stability of the scheme is investigated using solution bound. The uniform convergence of the scheme is discussed and proved. The scheme converges uniformly with the order of convergenceON−1+Δt2, where N is number of subintervals in spatial discretization and Δt is mesh length in temporal discretization. Two test numerical examples are considered to validate the theoretical findings of the scheme.


---

## Body

## 1. Introduction

Differential equations play a prominent role in many disciplines for modelling real-world problems in which the evolution of the system depends on the present state of the system. Problems that depend on the past and current evolution of the system are modelled by the delay deferential equation [1].Singularly perturbed differential equations (SPDEs) are differential equations in which their highest order derivative term is multiplied by a small perturbation parameter. Such types of differential equations model various processes in applied mathematics, for instance, in fluid mechanics, fluid dynamics, elasticity, quantum mechanics, chemical-reactor theory, reaction-diffusion processes, and other domains of fluid motion [2, 3]. Singularly perturbed differential difference equations model physical problems for which the evaluation does not only depend on the present state of the system but also on the past history. In general, when the perturbation parameter tends to zero, the smoothness of the solution of the singularly perturbed differential difference equations deteriorates and it forms a boundary layer [4]. Such types of singularly perturbed differential difference equations have an application in the modelling of neuronal variability [5, 6].Musila and Lansky [7] give a generalization for the model in [5, 6] in the form of singularly perturbed parabolic differential difference equations (SPPDDEs):(1)utx,t=σ22uxxx,t+μD−xτuxx,t+λsux+as,t+ωsux+is,t−λs+ωsux,t,for studying time evolution of the trajectories of the membrane potential, where the first derivative term denotes the exponential decay between two consecutive jumps caused by the input processes. The membrane potential decays exponentially to the resting level with a membrane time constant τ. μD and σ are diffusion moments of the Wiener process characterizing the influence of dendritic synapses on the cell excitability. The excitatory input contributes to the membrane potential by an amplitude as with intensity λs, and similarly, the inhibitory input contributes by an amplitude is with intensity ωs.The mathematical model in (1) is in the form of SPPDDEs. There are no known analytical methods of solution for physiologically acceptable parameters [7]. So, to simulate the solution of this model, one has to develop suitable numerical methods. Numerical methods developed for regular problems become inapplicable for solving SPPDDEs when the perturbation parameter tends to zero [8]. Due to this, there is an interest in the development of numerical methods that converges uniformly (or converges independent of the perturbation parameter).In recent years, scholars have devoted for the numerical treatment of SPPDDEs. In a series of papers [8–11], Bansal and Sharma developed different numerical schemes for solving SPPDDEs with large shift arguments. Kumar and Kadalbajoo [12] applied the implicit Euler method in temporal discretization and B-spline collocation method on the Shishkin mesh in the spatial discretization, and they obtain a first-order uniformly convergent result. Ramesh and Kadalbajoo [4] used the implicit Euler method in temporal discretization and upwind and midpoint upwind FDM on Shishkin mesh for the spatial discretization. They obtain almost first-order uniformly convergent numerical result. Rao and Chakravarthy [13] applied the implicit Euler method in temporal discretization and exponentially fitted operator FDM in the spatial discretization. They obtain a first-order uniformly convergent result. Woldaregay and Duressa [14–16] developed numerical schemes using fitted operator methods on uniform mesh. They obtained first-order uniformly convergent result. Kumar et al. [17] developed stable FDM using the implicit Euler method for temporal discretization and a class of nonstandard FDM using infinite serious approximation for the derivative terms in spatial discretization, and they apply Richardson extrapolation technique to extend the rate of convergence of the scheme to order two. Shivhare et al. [18] developed a numerical scheme using the implicit Euler method in temporal discretization and quadratic B-spline collocation method for the spatial discretization on the exponentially graded mesh. Daba and Duressa [19] used the implicit Euler method in temporal discretization and extended cubic B-spline collocation method for the spatial discretization on the Shishkin mesh.In this article, we developed numerical scheme using the Crank–Nicolson method in the temporal discretization and nonstandard FDM in the spatial discretization. We used Rothe’s procedure for discretizing the problem (i.e., first we discretize the temporal direction and followed by spatial discretization). The proposed scheme satisfies the discrete maximum principle and uniform stability estimate. The developed scheme treats the considered problem accurately with parameter uniform convergence.


### 1.1. Organization of the Paper

The rest of the article is organized as follows. In Section2, the statement of the problem is given. In Section 3, the asymptotic approximation and its properties are discussed. In Section 4, the proposed numerical scheme and its uniform convergence analysis are given. In Section 5, numerical examples, results, and discussions are given. Finally, in Section 6, the conclusion of the article is included.


### 1.2. Notation

In this article,N and M are denoted for the number of mesh points in spatial and temporal discretization, respectively. The symbol C (in some cases indexed) denotes a positive constant independent of the perturbation parameter ε and N. The norm . denotes the supremum (or maximum) norm.


## 1.1. Organization of the Paper

The rest of the article is organized as follows. In Section2, the statement of the problem is given. In Section 3, the asymptotic approximation and its properties are discussed. In Section 4, the proposed numerical scheme and its uniform convergence analysis are given. In Section 5, numerical examples, results, and discussions are given. Finally, in Section 6, the conclusion of the article is included.


## 1.2. Notation

In this article,N and M are denoted for the number of mesh points in spatial and temporal discretization, respectively. The symbol C (in some cases indexed) denotes a positive constant independent of the perturbation parameter ε and N. The norm . denotes the supremum (or maximum) norm.


## 2. Statement of the Problem: Aim of the Research

On the domainD=Ω×Λ=0,1×0,T, for some fixed number T>0, with the boundary ∂D=D¯\D, a singularly perturbed parabolic differential difference equation is given by(2)∂∂t+Lε,δ,ηux,t=fx,t,x,t∈D,where(3)Lε,δ,ηux,t=−ε2uxxx,t+axuxx,t+αxux−δ,t+βxux,t+ωxux+η,t,ε,0<ε≪1 is a singular perturbation parameter and δ and η are small delay and advance parameters satisfying δ,η=oε, with the initial and interval conditions(4)ux,0=u0x,x∈Db=x,0:x∈Ω¯,ux,t=ϕx,t,x,t∈Dl=x,t:x,t∈−δ,0×Λ,ux,t=ψx,t,x,t∈Dr=x,t:x,t∈1,1+η×Λ.The problem in (2) and (4) is the simplified form of the model in (1) by including the source function f with specified initial and interval conditions. From the convection diffusion fluid flow point of view, the terms in (2) and (4) can be classified as −ε2∂2/∂x2ux,t is called the diffusion term; ∂/∂t+ax∂/∂xux,t is called the convection term and the nonderivative terms referred to as the reaction terms.For the existence of a unique solution, the functionsa,α, β,ω,f,u0,ϕ, and ψ are assumed to be sufficiently smooth, bounded, and independent of ε. The coefficient functions α,β, and ω are assumed to satisfy(5)αx+βx+ωx≥ζ>0,x∈Ω¯,for some constant ζ. For small values of ε, the solution of (2) and (4) exhibits left or right boundary layer depending on the sign of the coefficient function of the convection term ax. If ax<0, a regular boundary layer occurs in the neighbourhood of Dl, and for the case ax>0, the layer occurs in the neighbourhood of Dr [20].The main contribution of this article is, to develop uniformly convergent numerical scheme using nonstandard FDM for solving the problem in (2) and (4). The developed method is applicable on the uniform mesh.


## 3. Properties of the Analytical Solution

In this section, we discuss behaviours of the analytical solution of the problem in (2) and (4). When the shift parameters δ and η are smaller than ε, using Taylor series approximation for the terms with shifts is valid [21]. So, we approximate ux−δ,t and ux+η,t as(6)ux−δ,t≈ux,t−δuxx,t+δ22uxxx,t+Oδ3,ux+η,t≈ux,t+ηuxx,t+η22uxxx,t+Oη3.Using the approximations in (6) into (2), we obtain(7)∂∂t+Lcεux,t=fx,t,x,t∈D,with initial and boundary conditions(8)ux,0=u0x,x∈Db=x,0:x∈Ω¯,u0,t=ϕ0,t,x,t∈Dl=x,t:x=0,t∈Λ,u1,t=ψ1,t,x,t∈Dr=x,t:x=1,t∈Λ,where Lcεux,t=−cεxuxxx,t+pxuxx,t+qxux,t, cεx=ε2−αxδ2/2−ωxη2/2, px=ax−δαx+ηωx and qx=αx+βx+ωx.For small values of the shift parametersδ,η, (2) and (5) and (7) and (8) are asymptotically equivalent, since the difference between the two equations is Oδ3,η3. We assumed that 0<cεx≤ε2−δ2α/2−η2ω/2=cε, where α and β are the lower bounds for αx and ωx, respectively. Furthermore, we assumed px≥p∗>0, which implies the occurrence of the boundary layer in the neighbourhood of Dr. The other case px≤p∗<0 implies the occurrence of the boundary layer in the neighbourhood of Dl and can be treated similarly. The layer is maintained for sufficiently small shift parameters δ,η.Lemma 1.
To avoid the presence of a singularity in the numerical solution, it is required that the given data satisfy the compatibility conditionsu00=ϕ0,0,u01=ψ1,0, and(9)ϕt0,0−cεu0xx0+p0u0x0+q0u00=f0,0,ψt1,0−cεu0xx1+p1u0x1+q1u01=f1,0,so that the data match at the two corners0,0and1,0. Let the functions px,qx, and fx,t be continuous on D, then (7) and (8) have a unique solutionu∈C2D. In particular, when the compatibility conditions are not satisfied, a unique classical solution still exist but may not be differentiable on all of∂D [22].In this article, we considered the case, when boundary layer occurs in the neighbourhood ofDr. Using the compatibility conditions, we deduce that there exists a constant C independent of cε such that for all x,t∈D¯.(10)ux,t−ux,0=ux,t−u0x≤Ct,ux,t−u0,t=ux,t−ϕ0,t≤C1−x.For details, interested readers can refer [22].Remark 1.
Note that there does not exist a constantCindependent ofcεsuch thatux,t−ψ1,t≤Cx, since the layer occurs in the neighbourhood ofDr.The problem obtained by settingcε=0 in (7) and (8) is called reduced problem and given as(11)ut0x,t+pxux0x,t+qxu0x,t=fx,t,x,t∈D,with the initial conditions(12)u0x,0=u0x,x∈Ω¯,u00,t=ϕ0,t,t∈Λ¯.It is first-order hyperbolic PDEs with initial data given along sidest=0 and x=0 of the domain D¯. For small values of cε, the solution ux,t of the problem in (7) and (8) is very close to the solution u0x,t of (11) and (12).Lemma 2.
The maximum principle. Letz be any sufficiently smooth function defined on D, which satisfies zx,t≥0,x,t∈∂D. Then, ∂/∂t+Lcεzx,t>0,x,t∈D implies that zx,t≥0,∀x,t∈D¯.Proof.
Letx∗,t∗ be such that zx∗,t∗=minx,t∈D¯zx,t and suppose that zx∗,t∗<0. It is clear that x∗,t∗∉∂D implying that x∗,t∗∈D. Since zx∗,t∗=minx,t∈D¯zx,t from the theory of extrema of a function in calculus implies, zxx∗,t∗=0, ztx∗,t∗=0 and zxxx∗,t∗≥0 giving that ∂/∂t+Lcεzx∗,t∗<0, which is contradiction to the assumption. Therefore, ∂/∂t+Lcεzx∗,t∗>0,x,t∈D. Hence, we obtain zx,t≥0,∀x,t∈D¯.Lemma 3.
Uniform stability estimate. Letux,t be the solution of the problem in (7) and (8). Then, it satisfies the following bound:(13)ux,t≤ζ−1f+u∂D,where u∂D=maxu0x,ϕ0,t,ψ1,t.Proof.
The proof is by the construction of barrier functionsϑ±x,t as ϑ±x,t=ζ−1f+u∂D±ux,t, and applying the maximum principle, we obtain the required bound. For the detail of the proof, the interested reader can refer [16].Lemma 4.
Letux,t be solution of (7) and (8). Then, the derivatives of the solutions with respect to x and t satisfy following the bound:(14)∂k∂xkux,t≤C1+cε−kexp−p∗1−xcε,x,t∈D¯,0≤k≤4,∂l∂tlux,t≤C,x,t∈D¯,0≤l≤2,where p∗ is lower bound of px.Proof.
Refer in [23].


## 4. Formulation of the Numerical Scheme

### 4.1. Temporal Discretization

We discretize the temporal domain0,T into M subintervals using the grid points t0=0,tj=jΔt,j=0,1,2,…,M, where Δt=T/M. Let Uj+1x is denoted for the approximation of ux,tj+1 at the j+1th time level discretization. The continuous problem in (7) and (8) semidiscretized using the Crank–Nicolson method in the temporal direction resulting to a system of boundary value problems(15)LcεΔtUj+1x=gj+1x,x∈Ω¯,j=1,2,…,M−1,with the following discretized boundary conditions:(16)Uj+10=ϕ0,tj+1,Uj+11=ψ1,tj+1,where, LcεΔtUj+1x=−cε/2Uj+1xxx+px/2Uj+1xx+dxUj+1x and gj+1x=cε/2Ujxxx−px/2Ujxx−rxUjx+1/2fx,tj+1+fx,tj, for dx=1/Δt+qx/2>0 and rx=qx/2−1/Δt.Lemma 5.
Semidiscrete maximum principle. Letzj+1 be a sufficiently smooth function on the domain Ω¯. If zj+10≥0,zj+11≥0, and Lx,cεΔtzj+1x≥0,x∈Ω, then, zj+1x≥0,∀x∈Ω¯.Proof.
Letx∗∈Ω¯ be such that zj+1x∗=minx∈Ω¯zj+1x and suppose that zx∗<0. From the assumption, it is clear that x∗∉0,1, which implies that x∗∈0,1. Now, using the differential operator at the point x∗, we have(17)Lx,cεΔtzj+1x∗=−cε2zj+1x∗xx+px∗2zj+1x∗x+dx∗zj+1x∗.
Applying the property in calculus for minimum points, we obtain(18)zj+1x∗x=0,zj+1x∗xx≥0.
Using the estimates in (18) into (17), we obtain Lx,cεΔtzj+1x∗<0, which is contradiction to the assumption Lx,cεΔtzj+1x∗≥0,x∈Ω. Therefore, we conclude that zj+1x≥0,∀x∈Ω¯.
Next, let us analyse the bound of the error in temporal semidiscretization. Let us denoteux,tj+1 and Uj+1x be the exact and approximate solution of the problem in (7) and (8) using the above discretization. Let us denote the local error for each time step by ej+1:=ux,tj+1−Uj+1x.Lemma 6.
Local error estimate. Suppose that∂l/∂tlux,t≤C,0≤l≤2. The local error estimate in the temporal discretization is given by(19)ej+1≤C1Δt3.Proof.
Using Taylor’s series expansion forux,tj+1 and ux,tj centered at tj+1/2, we obtain(20)ux,tj+1=ux,tj+1/2+Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3,ux,tj=ux,tj+1/2−Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3.
Substituting (20) into (7), we obtain(21)ux,tj+1−ux,tjΔt=utx,tj+1/2+OΔt2=cεuxxx,tj+1/2−pxuxx,tj+1/2−qxux,tj+1/2+fx,tj+1/2+OΔt2,where, ux,tj+1/2=1/2ux,tj+1+ux,tj+OΔt2,fx,tj+1/2=1/2fx,tj+1+fx,tj+OΔt2.
Now, it can be seen that the local error is the solution of(22)Lx,cεΔtej+1=OΔt3,ej+10=0=ej+11.
Hence, applying the maximum principle on the operator gives(23)ej+1≤C1Δt3.
Next, we need to show the bound for the global error of the temporal semidiscretization. Let us denoteEj+1 for the global error estimate up to the j+1th time step.Lemma 7.
Global error estimate. The global error estimate at tj+1 is given by(24)Ej+1≤CΔt2,j=1,2,…,M−1.Proof.
Using the local error estimate up to thej+1th time step given in the above lemma, we obtain the global error estimate at the j+1th time step as(25)Ej+1=∑l=1j+1el≤e1+e2+⋯+ej+1≤C1jΔtΔt2≤C1TΔt2,sincej+1Δt≤T=CΔt2,C1T=C,where C is constant independent of cε and Δt.
Next, we set a bound for the derivatives of solution of the problem in (15).Lemma 8.
For each,j=0,1,…,M−1, the solution of the boundary value problems in (15) satisfies the following bound:(26)dkdxkUj+1x≤C1+cε−kexp−p∗1−xcε,x∈Ω¯,0≤k≤4.Proof.
See in [23].


### 4.2. Spatial Discretization

For the spatial discretization, we apply the nonstandard FDM. The construction of nonstandard FDM depends on the knowledge of the corresponding exact FDM, which is given in [24].


#### 4.2.1. Exact Finite Difference Method

For the problem of the form in (15), in order to construct exact finite difference scheme, we follow the procedures developed in [24] for BVPs. Consider the constant coefficient subequations from (15) as(27)−cεUj+1xxx+p∗Uj+1xx+2θUj+1x=0,(28)−cεUj+1xxx+p∗Uj+1xx=0,where px≥p∗ and dx≥2θ. Since equation (27) is constant coefficient second-order differential equations, it has two independent solutions, that is,(29)Uj+1x=A1expλ1x+A2expλ2x,where λ1=−p∗+p∗2+8cεθ/−2cεandλ2=−p∗−p∗2+8cεθ/−2cε.We discretize the spatial domain0,1, using uniform mesh length Δx=h such that ΩN=xi=x0+ih,i=1,2,…,N,x0=0,xN=1,h=1/N, where N is the number of mesh intervals in spatial discretization. We denote the approximate solution of Uj+1x at mesh point xi by Ui,j+1. The target is to calculate a difference equation, which has the same general solution as (27) at the mesh point xi given by(30)Ui,j+1=A1expλ1xi+A2expλ2xi.Using the theory of difference equations for second-order linear difference equations in [24], we obtain(31)Ui−1,j+1expλ1xi−1expλ2xi−1Ui,j+1expλ1xiexpλ2xiUi+1,j+1expλ1xi+1expλ2xi+1=0,and then, by simplifying and substituting the values of λ1 and λ2, we obtain(32)expp∗h2cεUi−1,j+1−2coshhp∗2+8cεθ2cεUi,j+1+exp−p∗h2cεUi+1,j+1=0,which is an exact difference scheme for (44). For ε⟶0, we use the approximation hp∗2+8cεθ/2cε≈p∗h/2cε in (32). Hence, multiplying both sides by expp∗h/2cε, simplifying we obtain(33)Ui−1,j+1−2Ui,j+1+Ui+1,j+1=expp∗hcε−1Ui,j+1−Ui−1,j+1.Rearranging, we obtain(34)−cεUi−1,j+1−2Ui,j+1+Ui+1,j+1hcε/p∗expp∗h/cε−1+p∗Ui,j+1−Ui−1,j+1h=0.The denominator function for the second-order derivative approximation in (34) is γ2:=hcε/p∗exphp∗/cε−1. For the variable coefficient problems, we modify the denominator function as(35)γi2=hcεpxiexphpxicε−1,where γi2 is a function of cε,pxi, and h.


### 4.3. The Discrete Scheme

Using the denominator function in (35) into the discretization of (15), we obtain the nonstandard finite difference schemes as(36)Lcεh,ΔtUi,j+1=gi,j+1,i=1,2,…,N−1,where Lcεh,ΔtUi,j+1=−cε/2Ui+1,j+1−2Ui,j+1+Ui−1,j+1/γi2+pxi/2Ui,j+1−Ui−1,j+1/h+dxiUi,j+1 and(37)gi,j+1=cε2Ui+1,j−2Ui,j+Ui−1,jγi2−pxi2Ui,j−Ui−1,jh−rxiUi,j+12fxi,tj+1+fxi,tj.


### 4.4. Stability and Uniform Convergence Analysis

In this subsection, we want to show that the discrete scheme in (24) satisfies the discrete maximum principle, uniform stability estimates. and uniform convergence.Lemma 9.
Discrete maximum principle. Let Ui,j+1 be any mesh function satisfying U0,J+1≥0,UN,j+1≥0. Then, Lcεh,ΔtUi,j+1≥0,i=1,2,…,N−1, which implies that Ui,j+1≥0,i=0,1,…,N.Proof.
Suppose that there existsk∈0,1,…,N such that Uk,j+1=min0≤i≤NUi,j+1<0. From the assumption, it implies k≠0,N. Also, we assume that Uk+1,j+1−Uk,j+1>0 and Uk,j+1−Uk−1,j+1<0. Now, we have(38)Lcεh,ΔtUk,j+1=−cε2Uk+1,j+1−2Uk,j+1+Uk−1,j+1γk2+pxk2Uk,j+1−Uk−1,j+1h+dxkUk,j+1,and using the assumptions made above, we obtain Lcεh,ΔtUk,j+1<0, for k=1,2,3,…,N−1. Thus, the supposition Ui,j+1<0, for i=0,1,…,N is wrong. Hence, we obtain. Ui,j+1≥0,i=0,1,…,N
Next, we want to prove that the scheme in (36) satisfies the uniform stability estimate.Lemma 10.
Uniform stability estimate. The solution Ui,j+1 of the discrete scheme in (24) satisfies the following bound:(39)Ui,j+1≤Lcεh,ΔtUi,j+1ζ∗+maxϕ0,tj+1,ψ1,tj+1,where ζ∗ is the lower bound of dx.Proof.
LetQ=Lcεh,ΔtUi,j+1/ζ∗+maxϕ0,tj+1,ψ1,tj+1 and define the barrier function ϑi,j+1± by ϑi,j+1±=Q±Ui,j+1. At the boundary points, we obtain(40)ϑ0,j+1±=Q±U0,j+1=Lcεh,ΔtUi,j+1ζ∗+maxϕ0,tj+1,ψ1,tj+1±ϕ0,tj+1≥0,ϑN,j+1±=Q±UN,j+1=Lcεh,ΔtUi,j+1ζ∗+maxϕ0,tj+1,ψ1,tj+1±ψ1,tj+1≥0.
On the discretized spatial domainxi,0<i<N, we obtain(41)Lcεh,Δtϑi,j+1±=−cε2Q±Ui+1,j+1−2Q±Ui,j+1+Q±Ui−1,j+1γi2+pxi2Q±Ui,j+1−Q±Ui−1,j+1h+dxiQ±Ui,j+1=dxiQ±Lcεh,ΔtUi,j+1=dxiLcεh,ΔtUi,j+1ζ∗+maxϕ0,tj+1,ψ1,tj+1±gi,j+1≥0,sincedxi≥ζ∗.
Using the maximum principle in Lemma9, we obtain ϑi,j+1±≥0,∀xi∈Ω¯N. Hence, the required bound is obtained.Lemma 11.
Letγi2 be the denominator function defined in (35), and then, the following bound holds(42)cεh2γi2−1≤C3h.Proof.
Let us defineρ=pxih/cε,ρ∈0,∞. Using the expression for γi2, we obtain(43)cεh2γi2−1=pxih1expρ−1−1ρ=hpxiexpρ−1−ρρexpρ−1≕hpxiRρ.
Next, let us set a bound forRρ=expρ−1−ρ/ρexpρ−1. Since(44)limρ⟶0Rρ=12,limρ⟶∞Rρ=0.
We obtain the bounds asRρ≤C3,∀ρ∈0,∞. Hence, from (43) and (44), we obtain(45)cεh2γi2−1≤C3h.
Let us define the forward and backward finite differences operators in spatial direction as(46)Dx+Zj+1xi=Zj+1xi+1−Zj+1xih,Dx−Zj+1xi=Zj+1xi−Zj+1xi−1h,respectively, and the second-order finite difference operator as(47)Dx+Dx−Zj+1xi=D+Zj+1xi−D−Zj+1xih=Zj+1xi+1−2Zj+1xi+Zj+1xih2.Theorem 1.
The discrete solutionUi,j+1 of the problem in (9) satisfies the truncation error bound:(48)Lcεh,ΔtUj+1xi−Ui,j+1≤Ch1+maxxi∈0,1exp−p∗1−xi/cεcε3.Proof.
We consider the truncation error(49)Lcεh,ΔtUj+1xi−Ui,j+1=−cεd2dx2−Dx+Dx−h2γi2Uj+1xi+pxiddx−Dx−Uj+1xi≤C1cεd2dx2−Dx+Dx−Uj+1xi+C1cεh2γi2−1Dx+Dx−Uj+1xi+C1hd2dx2Uj+1xi.
Using the estimate in Lemma11 for the bound of cεh2/γi2−1, the truncation error bound becomes(50)Lcεh,ΔtUj+1xi−Ui,j+1≤C3cεh2d4dx4Uj+1xi+C3hd2dx2Uj+1xi.
Using the bound of the derivatives of the solution in Lemma8 in (50) gives(51)Lcεh,ΔtUj+1xi−Ui,j+1≤C4h2cε+cε−3exp−p∗1−xicε+C4h1+cε−2exp−p∗1−xicε≤Ch1+maxxi∈0,1exp−p∗1−xi/cεcε3,since,cε3≤cε2.
This completes the proof.Lemma 12.
For a fixed mesh numberN and for ε⟶0, it holds(52)limcε⟶0max1≤i≤N−1exp−p∗1−xi/cεcεm=0,m=1,2,3,…,where xi=ih,h=1/N,∀i=1,2,…,N−1.Proof.
Refer in [25, 26].Theorem 2.
The discrete solutionUi,j+1 of the problem in (15) satisfies the error bound:(53)maxxi∈0,1Uj+1xi−Ui,j+1≤CN−1,where h=N−1.Proof.
By using the results of Lemma12 in Theorem 1 and applying the bound in Lemma 10 gives the required bound.Theorem 3.
Under the hypothesis of boundedness of discrete solution, the solution of (7) and (8) satisfies the parameter uniform bound:(54)sup0<cε≪1maxxi∈0,1uxi,tj+1−Ui,j+1≤CN−1+Δt2.Proof.
The parameter uniform convergence of the scheme follows from the bound in (15) and (31):(55)EN,M=maxε,δ,ηEε,δ,ηN,M.


## 4.1. Temporal Discretization

We discretize the temporal domain0,T into M subintervals using the grid points t0=0,tj=jΔt,j=0,1,2,…,M, where Δt=T/M. Let Uj+1x is denoted for the approximation of ux,tj+1 at the j+1th time level discretization. The continuous problem in (7) and (8) semidiscretized using the Crank–Nicolson method in the temporal direction resulting to a system of boundary value problems(15)LcεΔtUj+1x=gj+1x,x∈Ω¯,j=1,2,…,M−1,with the following discretized boundary conditions:(16)Uj+10=ϕ0,tj+1,Uj+11=ψ1,tj+1,where, LcεΔtUj+1x=−cε/2Uj+1xxx+px/2Uj+1xx+dxUj+1x and gj+1x=cε/2Ujxxx−px/2Ujxx−rxUjx+1/2fx,tj+1+fx,tj, for dx=1/Δt+qx/2>0 and rx=qx/2−1/Δt.Lemma 5.
Semidiscrete maximum principle. Letzj+1 be a sufficiently smooth function on the domain Ω¯. If zj+10≥0,zj+11≥0, and Lx,cεΔtzj+1x≥0,x∈Ω, then, zj+1x≥0,∀x∈Ω¯.Proof.
Letx∗∈Ω¯ be such that zj+1x∗=minx∈Ω¯zj+1x and suppose that zx∗<0. From the assumption, it is clear that x∗∉0,1, which implies that x∗∈0,1. Now, using the differential operator at the point x∗, we have(17)Lx,cεΔtzj+1x∗=−cε2zj+1x∗xx+px∗2zj+1x∗x+dx∗zj+1x∗.
Applying the property in calculus for minimum points, we obtain(18)zj+1x∗x=0,zj+1x∗xx≥0.
Using the estimates in (18) into (17), we obtain Lx,cεΔtzj+1x∗<0, which is contradiction to the assumption Lx,cεΔtzj+1x∗≥0,x∈Ω. Therefore, we conclude that zj+1x≥0,∀x∈Ω¯.
Next, let us analyse the bound of the error in temporal semidiscretization. Let us denoteux,tj+1 and Uj+1x be the exact and approximate solution of the problem in (7) and (8) using the above discretization. Let us denote the local error for each time step by ej+1:=ux,tj+1−Uj+1x.Lemma 6.
Local error estimate. Suppose that∂l/∂tlux,t≤C,0≤l≤2. The local error estimate in the temporal discretization is given by(19)ej+1≤C1Δt3.Proof.
Using Taylor’s series expansion forux,tj+1 and ux,tj centered at tj+1/2, we obtain(20)ux,tj+1=ux,tj+1/2+Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3,ux,tj=ux,tj+1/2−Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3.
Substituting (20) into (7), we obtain(21)ux,tj+1−ux,tjΔt=utx,tj+1/2+OΔt2=cεuxxx,tj+1/2−pxuxx,tj+1/2−qxux,tj+1/2+fx,tj+1/2+OΔt2,where, ux,tj+1/2=1/2ux,tj+1+ux,tj+OΔt2,fx,tj+1/2=1/2fx,tj+1+fx,tj+OΔt2.
Now, it can be seen that the local error is the solution of(22)Lx,cεΔtej+1=OΔt3,ej+10=0=ej+11.
Hence, applying the maximum principle on the operator gives(23)ej+1≤C1Δt3.
Next, we need to show the bound for the global error of the temporal semidiscretization. Let us denoteEj+1 for the global error estimate up to the j+1th time step.Lemma 7.
Global error estimate. The global error estimate at tj+1 is given by(24)Ej+1≤CΔt2,j=1,2,…,M−1.Proof.
Using the local error estimate up to thej+1th time step given in the above lemma, we obtain the global error estimate at the j+1th time step as(25)Ej+1=∑l=1j+1el≤e1+e2+⋯+ej+1≤C1jΔtΔt2≤C1TΔt2,sincej+1Δt≤T=CΔt2,C1T=C,where C is constant independent of cε and Δt.
Next, we set a bound for the derivatives of solution of the problem in (15).Lemma 8.
For each,j=0,1,…,M−1, the solution of the boundary value problems in (15) satisfies the following bound:(26)dkdxkUj+1x≤C1+cε−kexp−p∗1−xcε,x∈Ω¯,0≤k≤4.Proof.
See in [23].


## 4.2. Spatial Discretization

For the spatial discretization, we apply the nonstandard FDM. The construction of nonstandard FDM depends on the knowledge of the corresponding exact FDM, which is given in [24].


### 4.2.1. Exact Finite Difference Method

For the problem of the form in (15), in order to construct exact finite difference scheme, we follow the procedures developed in [24] for BVPs. Consider the constant coefficient subequations from (15) as(27)−cεUj+1xxx+p∗Uj+1xx+2θUj+1x=0,(28)−cεUj+1xxx+p∗Uj+1xx=0,where px≥p∗ and dx≥2θ. Since equation (27) is constant coefficient second-order differential equations, it has two independent solutions, that is,(29)Uj+1x=A1expλ1x+A2expλ2x,where λ1=−p∗+p∗2+8cεθ/−2cεandλ2=−p∗−p∗2+8cεθ/−2cε.We discretize the spatial domain0,1, using uniform mesh length Δx=h such that ΩN=xi=x0+ih,i=1,2,…,N,x0=0,xN=1,h=1/N, where N is the number of mesh intervals in spatial discretization. We denote the approximate solution of Uj+1x at mesh point xi by Ui,j+1. The target is to calculate a difference equation, which has the same general solution as (27) at the mesh point xi given by(30)Ui,j+1=A1expλ1xi+A2expλ2xi.Using the theory of difference equations for second-order linear difference equations in [24], we obtain(31)Ui−1,j+1expλ1xi−1expλ2xi−1Ui,j+1expλ1xiexpλ2xiUi+1,j+1expλ1xi+1expλ2xi+1=0,and then, by simplifying and substituting the values of λ1 and λ2, we obtain(32)expp∗h2cεUi−1,j+1−2coshhp∗2+8cεθ2cεUi,j+1+exp−p∗h2cεUi+1,j+1=0,which is an exact difference scheme for (44). For ε⟶0, we use the approximation hp∗2+8cεθ/2cε≈p∗h/2cε in (32). Hence, multiplying both sides by expp∗h/2cε, simplifying we obtain(33)Ui−1,j+1−2Ui,j+1+Ui+1,j+1=expp∗hcε−1Ui,j+1−Ui−1,j+1.Rearranging, we obtain(34)−cεUi−1,j+1−2Ui,j+1+Ui+1,j+1hcε/p∗expp∗h/cε−1+p∗Ui,j+1−Ui−1,j+1h=0.The denominator function for the second-order derivative approximation in (34) is γ2:=hcε/p∗exphp∗/cε−1. For the variable coefficient problems, we modify the denominator function as(35)γi2=hcεpxiexphpxicε−1,where γi2 is a function of cε,pxi, and h.


## 4.2.1. Exact Finite Difference Method

For the problem of the form in (15), in order to construct exact finite difference scheme, we follow the procedures developed in [24] for BVPs. Consider the constant coefficient subequations from (15) as(27)−cεUj+1xxx+p∗Uj+1xx+2θUj+1x=0,(28)−cεUj+1xxx+p∗Uj+1xx=0,where px≥p∗ and dx≥2θ. Since equation (27) is constant coefficient second-order differential equations, it has two independent solutions, that is,(29)Uj+1x=A1expλ1x+A2expλ2x,where λ1=−p∗+p∗2+8cεθ/−2cεandλ2=−p∗−p∗2+8cεθ/−2cε.We discretize the spatial domain0,1, using uniform mesh length Δx=h such that ΩN=xi=x0+ih,i=1,2,…,N,x0=0,xN=1,h=1/N, where N is the number of mesh intervals in spatial discretization. We denote the approximate solution of Uj+1x at mesh point xi by Ui,j+1. The target is to calculate a difference equation, which has the same general solution as (27) at the mesh point xi given by(30)Ui,j+1=A1expλ1xi+A2expλ2xi.Using the theory of difference equations for second-order linear difference equations in [24], we obtain(31)Ui−1,j+1expλ1xi−1expλ2xi−1Ui,j+1expλ1xiexpλ2xiUi+1,j+1expλ1xi+1expλ2xi+1=0,and then, by simplifying and substituting the values of λ1 and λ2, we obtain(32)expp∗h2cεUi−1,j+1−2coshhp∗2+8cεθ2cεUi,j+1+exp−p∗h2cεUi+1,j+1=0,which is an exact difference scheme for (44). For ε⟶0, we use the approximation hp∗2+8cεθ/2cε≈p∗h/2cε in (32). Hence, multiplying both sides by expp∗h/2cε, simplifying we obtain(33)Ui−1,j+1−2Ui,j+1+Ui+1,j+1=expp∗hcε−1Ui,j+1−Ui−1,j+1.Rearranging, we obtain(34)−cεUi−1,j+1−2Ui,j+1+Ui+1,j+1hcε/p∗expp∗h/cε−1+p∗Ui,j+1−Ui−1,j+1h=0.The denominator function for the second-order derivative approximation in (34) is γ2:=hcε/p∗exphp∗/cε−1. For the variable coefficient problems, we modify the denominator function as(35)γi2=hcεpxiexphpxicε−1,where γi2 is a function of cε,pxi, and h.


## 4.3. The Discrete Scheme

Using the denominator function in (35) into the discretization of (15), we obtain the nonstandard finite difference schemes as(36)Lcεh,ΔtUi,j+1=gi,j+1,i=1,2,…,N−1,where Lcεh,ΔtUi,j+1=−cε/2Ui+1,j+1−2Ui,j+1+Ui−1,j+1/γi2+pxi/2Ui,j+1−Ui−1,j+1/h+dxiUi,j+1 and(37)gi,j+1=cε2Ui+1,j−2Ui,j+Ui−1,jγi2−pxi2Ui,j−Ui−1,jh−rxiUi,j+12fxi,tj+1+fxi,tj.


## 4.4. Stability and Uniform Convergence Analysis

In this subsection, we want to show that the discrete scheme in (24) satisfies the discrete maximum principle, uniform stability estimates. and uniform convergence.Lemma 9.
Discrete maximum principle. Let Ui,j+1 be any mesh function satisfying U0,J+1≥0,UN,j+1≥0. Then, Lcεh,ΔtUi,j+1≥0,i=1,2,…,N−1, which implies that Ui,j+1≥0,i=0,1,…,N.Proof.
Suppose that there existsk∈0,1,…,N such that Uk,j+1=min0≤i≤NUi,j+1<0. From the assumption, it implies k≠0,N. Also, we assume that Uk+1,j+1−Uk,j+1>0 and Uk,j+1−Uk−1,j+1<0. Now, we have(38)Lcεh,ΔtUk,j+1=−cε2Uk+1,j+1−2Uk,j+1+Uk−1,j+1γk2+pxk2Uk,j+1−Uk−1,j+1h+dxkUk,j+1,and using the assumptions made above, we obtain Lcεh,ΔtUk,j+1<0, for k=1,2,3,…,N−1. Thus, the supposition Ui,j+1<0, for i=0,1,…,N is wrong. Hence, we obtain. Ui,j+1≥0,i=0,1,…,N
Next, we want to prove that the scheme in (36) satisfies the uniform stability estimate.Lemma 10.
Uniform stability estimate. The solution Ui,j+1 of the discrete scheme in (24) satisfies the following bound:(39)Ui,j+1≤Lcεh,ΔtUi,j+1ζ∗+maxϕ0,tj+1,ψ1,tj+1,where ζ∗ is the lower bound of dx.Proof.
LetQ=Lcεh,ΔtUi,j+1/ζ∗+maxϕ0,tj+1,ψ1,tj+1 and define the barrier function ϑi,j+1± by ϑi,j+1±=Q±Ui,j+1. At the boundary points, we obtain(40)ϑ0,j+1±=Q±U0,j+1=Lcεh,ΔtUi,j+1ζ∗+maxϕ0,tj+1,ψ1,tj+1±ϕ0,tj+1≥0,ϑN,j+1±=Q±UN,j+1=Lcεh,ΔtUi,j+1ζ∗+maxϕ0,tj+1,ψ1,tj+1±ψ1,tj+1≥0.
On the discretized spatial domainxi,0<i<N, we obtain(41)Lcεh,Δtϑi,j+1±=−cε2Q±Ui+1,j+1−2Q±Ui,j+1+Q±Ui−1,j+1γi2+pxi2Q±Ui,j+1−Q±Ui−1,j+1h+dxiQ±Ui,j+1=dxiQ±Lcεh,ΔtUi,j+1=dxiLcεh,ΔtUi,j+1ζ∗+maxϕ0,tj+1,ψ1,tj+1±gi,j+1≥0,sincedxi≥ζ∗.
Using the maximum principle in Lemma9, we obtain ϑi,j+1±≥0,∀xi∈Ω¯N. Hence, the required bound is obtained.Lemma 11.
Letγi2 be the denominator function defined in (35), and then, the following bound holds(42)cεh2γi2−1≤C3h.Proof.
Let us defineρ=pxih/cε,ρ∈0,∞. Using the expression for γi2, we obtain(43)cεh2γi2−1=pxih1expρ−1−1ρ=hpxiexpρ−1−ρρexpρ−1≕hpxiRρ.
Next, let us set a bound forRρ=expρ−1−ρ/ρexpρ−1. Since(44)limρ⟶0Rρ=12,limρ⟶∞Rρ=0.
We obtain the bounds asRρ≤C3,∀ρ∈0,∞. Hence, from (43) and (44), we obtain(45)cεh2γi2−1≤C3h.
Let us define the forward and backward finite differences operators in spatial direction as(46)Dx+Zj+1xi=Zj+1xi+1−Zj+1xih,Dx−Zj+1xi=Zj+1xi−Zj+1xi−1h,respectively, and the second-order finite difference operator as(47)Dx+Dx−Zj+1xi=D+Zj+1xi−D−Zj+1xih=Zj+1xi+1−2Zj+1xi+Zj+1xih2.Theorem 1.
The discrete solutionUi,j+1 of the problem in (9) satisfies the truncation error bound:(48)Lcεh,ΔtUj+1xi−Ui,j+1≤Ch1+maxxi∈0,1exp−p∗1−xi/cεcε3.Proof.
We consider the truncation error(49)Lcεh,ΔtUj+1xi−Ui,j+1=−cεd2dx2−Dx+Dx−h2γi2Uj+1xi+pxiddx−Dx−Uj+1xi≤C1cεd2dx2−Dx+Dx−Uj+1xi+C1cεh2γi2−1Dx+Dx−Uj+1xi+C1hd2dx2Uj+1xi.
Using the estimate in Lemma11 for the bound of cεh2/γi2−1, the truncation error bound becomes(50)Lcεh,ΔtUj+1xi−Ui,j+1≤C3cεh2d4dx4Uj+1xi+C3hd2dx2Uj+1xi.
Using the bound of the derivatives of the solution in Lemma8 in (50) gives(51)Lcεh,ΔtUj+1xi−Ui,j+1≤C4h2cε+cε−3exp−p∗1−xicε+C4h1+cε−2exp−p∗1−xicε≤Ch1+maxxi∈0,1exp−p∗1−xi/cεcε3,since,cε3≤cε2.
This completes the proof.Lemma 12.
For a fixed mesh numberN and for ε⟶0, it holds(52)limcε⟶0max1≤i≤N−1exp−p∗1−xi/cεcεm=0,m=1,2,3,…,where xi=ih,h=1/N,∀i=1,2,…,N−1.Proof.
Refer in [25, 26].Theorem 2.
The discrete solutionUi,j+1 of the problem in (15) satisfies the error bound:(53)maxxi∈0,1Uj+1xi−Ui,j+1≤CN−1,where h=N−1.Proof.
By using the results of Lemma12 in Theorem 1 and applying the bound in Lemma 10 gives the required bound.Theorem 3.
Under the hypothesis of boundedness of discrete solution, the solution of (7) and (8) satisfies the parameter uniform bound:(54)sup0<cε≪1maxxi∈0,1uxi,tj+1−Ui,j+1≤CN−1+Δt2.Proof.
The parameter uniform convergence of the scheme follows from the bound in (15) and (31):(55)EN,M=maxε,δ,ηEε,δ,ηN,M.


## 5. Numerical Results and Discussion

To validate the established theoretical results, we develop an algorithm and perform experiments using the proposed numerical scheme on the problem of the form given in (2) and (4).Example 1.
Consider the problem(56)∂u∂t−ε2∂2u∂x2+2−x2∂u∂x+2ux−δ,t+x−3ux,t+ux+η,t=10t2exp−tx1−x,with the initial condition u0x=0,0≤x≤1 and the interval conditions ϕx,t=0,−δ≤x≤0,ψx,t=0,1≤x≤1+η, and T=3.Example 2.
Consider the problem(57)∂u∂t−ε2∂2u∂x2+2−x2∂u∂x+1+xux−δ,t+1+x2+cosπxux,t+3ux+η,t=sinπx,with the initial condition u0x=0,0≤x≤1 and the interval conditions ϕx,t=0,−δ≤x≤0,ψx,t=0,1≤x≤1+η, and T=3.The exact solution of the considered examples is not available. So, the maximum absolute error,ε-uniform error and rate of convergence are calculated using the double mesh principle. Let Ui,jN,M be the computed solution of the problem using N,M numbers of mesh points and Ui,j2N,2M be the computed solution on the double number of mesh points 2N,2M by including the midpoints xi+1/2 and tj+1/2 into the mesh points. The maximum absolute error is calculated as(58)Eε,δ,ηN,M=maxi,jUi,jN,M−Ui,j2N,2M.The uniform error estimate is calculated using(59)EN,M=maxε,δ,ηEε,δ,ηN,M.The rate of convergence of the scheme is calculated using(60)rε,δ,ηN,M=log2Eε,δ,ηN,MEε,δ,η2N,2M,and the uniform rate of convergence is calculated using(61)rN,M=log2EN,ME2N,2M.The solution of the problems in example1 and example 2 exhibits a regular boundary layer on the right side of the spatial domain. In Figures 1 and 2, one can observe boundary layer formation as the perturbation parameter ε goes small, particularly at T=3, one can observe in Figures 3(d) and 4(d) . In Figures 3(a)–3(d) and 4(a)–4(d), the effect of the shift parameters and the perturbation parameter is shown on the profile of the solution by considering different values for the arameters.Figure 1
3D view of numerical solution of example1 on (a) ε=2−2 on (b) ε=2−20 for δ=0.6ε,η=0.5ε.
(a)(b)Figure 2
3D view of numerical solution of example2 exhibiting right boundary layer on (a) ε=2−2 on (b) ε=2−20 for δ=0.6ε,η=0.5ε.
(a)(b)Figure 3
Example1, on (a) effect of delay, (b) effect of advance, (c) effect of both advance and delay, and (d) effect of perturbation parameter on solution at T=3.
(a)(b)(c)(d)Figure 4
Example2, on (a) effect of delay, (b) effect of advance, (c) effect of both advance and delay, and (d) effect of perturbation parameter on solution at T=3.
(a)(b)(c)(d)The numerical result in Tables1 and 2 show the maximum absolute error and rate of convergence of the scheme. The maximum absolute errors presented in Tables 1 and 3 show the uniform convergence of the scheme. One can observe from these tables, as the value of the perturbation parameter goes small, the maximum absolute error becomes stable and constant (i.e., converges are independent of the influence of the perturbation parameters ε). Tables 3 and 4 show the maximum absolute error for different values of delay and advance parameters. As it is observed in these tables, the shift parameters do not influence the uniform convergence of the scheme. Table 5 shows the maximum absolute error and rate of convergence of the scheme for temporal discretization. As one observes in this table, the temporal discretization has the second order of convergence with correct agreement to the theoretical analysis. Tables 6 and 7 present the uniform error and uniform rate of convergence of the scheme and results from published papers in literature. The proposed scheme is more accurate than the upwind scheme in the Shishkin mesh, midpoint upwind scheme in the Shishkin mesh in [4], cubic B-spline scheme in the Shishkin mesh in [12], and the exponentially fitted scheme in [13].Table 1
Example1, maximum absolute error and rate of convergence for δ=0.6ε and η=0.5ε.
εN =  2526272829210↓M = 6012024048096019202−65.7729e −032.9297e −031.4758e −037.4067e −043.7103e −041.8557e −040.97850.98930.99460.99730.9996—2−85.7496e −032.9179e −031.4699e −037.3769e −043.6953e −041.8494e −040.97850.98920.99460.99730.9986—2−105.7439e −032.9150e −031.4684e −037.3694e −043.6916e −041.8475e −040.97850.98930.99460.99730.9987—2−125.7424e −032.9142e −031.4680e −037.3676e −043.6906e −041.8470e −040.97860.98920.99460.99730.9987—2−145.7421e −032.9141e −031.4679e −037.3671e −043.6904e −041.8469e −040.97850.98930.99460.99730.9987−2−165.7420e −032.9140e −031.4679e −037.3670e −043.6904e −041.8469e −040.97860.98920.99460.99730.9987—2−185.7420e −032.9140e −031.4679e −037.3670e −043.6903e −041.8469e −040.97860.98920.99460.99730.9986—2−205.7420e −032.9140e −031.4679e −037.3670e −043.6903e −041.8468e −040.97860.98920.99460.99730.9986—EN,M5.7729e −032.9297e −031.4758e −037.4067e −043.7103e −041.8557e −04rN,M0.97850.98920.99460.99730.9987—Table 2
Example2, maximum absolute error and rate of convergence for δ=0.6ε and η=0.5ε.
εN=242526272829↓M = 30601202404809602−63.9860e − 032.1291e − 031.0941e − 035.5455e − 042.7923e − 041.4010e − 040.90470.96050.98040.98990.9950—2−83.9858e − 032.1279e − 031.0931e − 035.5416e − 042.7902e − 041.3999e − 040.90540.96100.98010.98990.9950—2−103.9857e − 032.1276e − 031.0929e − 035.5406e − 042.7897e − 041.3997e − 040.90560.96110.98000.98990.9950—2−123.9857e − 032.1275e − 031.0928e − 035.5404e − 042.7896e − 041.3996e − 040.90570.96110.98000.98990.9950—2−143.9857e − 032.1275e − 031.0928e − 035.5403e − 042.7895e − 041.3996e − 040.90570.96110.98000.99000.9950—2−163.9857e − 032.1275e − 031.0928e − 035.5403e − 042.7895e − 041.3996e − 040.90570.96110.98000.99000.9950—2−183.9857e − 032.1274e − 031.0928e − 035.5403e − 042.7895e − 041.3996e − 040.90570.96110.98000.99000.9950—2−203.9857e − 032.1274e − 031.0928e − 035.5403e − 042.7895e − 041.3996e − 040.90570.96110.98000.99000.9950—EN,M3.9860e − 032.1291e − 031.0941e − 035.5455e − 042.7923e − 041.4010e − 04rN,M0.90470.960500.98000.98990.99500—Table 3
Example1, maximum absolute error for different values of δ and η for ε=2−10.
N=2526272829M= 60120240480960δ↓,η=0.5ε05.7406e − 032.9133e − 031.4676e − 037.3652e − 043.6895e − 040.1ε5.7411e − 032.9136e − 031.4677e − 037.3659e − 043.6898e − 040.3ε5.7422e − 032.9141e − 031.4680e − 037.3673e − 043.6905e − 040.5ε5.7433e − 032.9147e − 031.4683e − 037.3687e − 043.6912e − 040.7ε5.7444e − 032.9153e − 031.4685e − 037.3701e − 043.6919e − 04η↓,δ=0.6ε05.7452e − 032.9157e − 031.4688e − 037.3712e − 043.6925e − 040.1ε5.7450e − 032.9155e − 031.4687e − 037.3708e − 043.6923e − 040.3ε5.7444e − 032.9153e − 031.4685e − 037.3701e − 043.6919e − 040.5ε5.7439e − 032.9150e − 031.4684e − 037.3694e − 043.6916e − 040.7ε5.7433e − 032.9147e − 031.4683e − 037.3687e − 043.6912e − 04Table 4
Example2, maximum absolute error for different values of δ and η for ε=2−10.
N=2526272829M= 60120240480960δ↓,η=0.5ε02.1276e − 031.0930e − 035.5407e − 042.7898e − 041.3997e − 040.1ε2.1276e − 031.0929e − 035.5407e − 042.7897e − 041.3997e − 040.3ε2.1276e − 031.0929e − 035.5407e − 042.7897e − 041.3997e − 040.5ε2.1276e − 031.0929e − 035.5406e − 042.7897e − 041.3997e − 040.7ε2.1275e − 031.0929e − 035.5406e − 042.7897e − 041.3997e − 04η↓,δ=0.6ε02.1274e − 031.0928e − 035.5403e − 042.7895e − 041.3995e − 040.1ε2.1275e − 031.0928e − 035.5403e − 042.7895e − 041.3995e − 040.3ε2.1275e − 031.0929e − 035.5405e − 042.7896e − 041.3997e − 040.5ε2.1276e − 031.0929e − 035.5406e − 042.7897e − 041.3998e − 040.6ε2.1276e − 031.0929e − 035.5407e − 042.7899e − 041.3999e − 04Table 5
Maximum absolute error and rate of convergence forδ=0.6ε and η=0.5ε.
Example1εN =  252729252729↓M = 60120240601202402−105.7439e − 031.4684e − 033.6916e − 042.1276e − 035.5406e − 041.3997e − 041.96781.9919―1.94111.9849―2−125.7424e − 031.4680e − 033.6906e − 042.1274e − 035.5403e − 041.3996e − 041.96781.9919―1.94111.9850―2−145.7421e − 031.4679e − 033.6904e − 042.1274e − 035.5403e − 041.3996e − 041.96781.9919―1.94111.9850―2−165.7420e − 031.4679e − 033.6904e − 042.1274e − 035.5403e − 041.3996e − 041.96781.9919―1.94111.9850―2−185.7420e − 031.4679e − 033.6903e − 042.1274e − 035.5403e − 041.3996e − 041.96781.9919―1.94111.9850―2−205.7420e − 031.4679e − 033.6903e − 042.1274e − 035.5403e − 041.3996e − 041.96781.9919―1.94111.9850―Table 6
Example1, comparison of uniform absolute error and uniform rate of convergence.
SchemesN=26272829210↓M = 1202404809601920ProposedEN,M2.9297e − 031.4758e − 037.4067e − 043.7103e − 041.8557e − 04SchemerN,M0.98920.99460.99730.9987―Upwind in [4]EN,M9.2021e − 034.9863e − 032.6885e − 031.4245e − 037.5403e − 04rM,N0.88400.89120.91630.9178―Midpoint up in [4]EN,M6.5942e − 033.8199e − 032.1180e − 031.1399e − 036.0635e − 04rN,M0.78770.85080.89380.9107―B-splineEN,M4.4966e − 032.4450e − 031.2728e − 036.4609e − 043.2774e − 04Colloc. in [12]rN,M0.87910.94180.97150.9859―Expo.EN,M3.3107e − 031.7254e − 038.8049e − 044.4473e − 042.2349e − 04Fitted in [13]rN,M0.94020.97050.98540.9927―Table 7
Example2, comparison of uniform absolute error and uniform rate of convergence.
SchemesN=2526272829↓M = 60120240480960ProposedEN,M2.1291e − 031.0941e − 035.5455e − 042.7923e − 041.4010e − 04SchemerN,M0.960500.98000.98990.99500―B-splineEN,M7.6388e − 033.8384e − 031.9277e − 039.6627e − 044.8379e − 04Colloc. in [12]rN,M0.99250.99250.99640.9981―Figure5 shows the log-log plot of the maximum absolute error verses number of mesh points in space N of the scheme for different values of the perturbation parameter. As one observes from the figures, for all values of the perturbation parameter, the line graphs are overlapped and parallel to the reference line ON−1. This indicates that the convergence of the scheme is independent of the perturbation parameter ε.Figure 5
Log-log plot of maximum absolute error versesN on (a) example 1 on (b) example 2 for different values of ε.
(a)(b)


## 6. Conclusion

In this article, a parameter uniform numerical method is developed for solving a singularly perturbed parabolic differential equation involving both delay and advance on the spatial variable of reaction terms. The considered problems’ solution exhibits a boundary layer on the right side of the spatial domain. The developed scheme constitutes the Crank–Nicolson method in temporal discretization and nonstandard FDM in the spatial discretization. The uniform stability and parameter uniform convergence of the scheme are investigated and proved. The applicability of the scheme is validated by taking numerical test examples. The effect of the perturbation parameter, effect of the delay, and advance parameters on the solution are depicted using figures. The proposed scheme converges uniformly (i.e., converges independent of the perturbation parameter) with the order of convergence one in space and two in time directions. The performance of the scheme is investigated by comparing the results with some schemes in published papers. We found that the proposed scheme is more accurate, stable, and parameter uniformly convergent.



---

*Source: 6637661-2021-11-30.xml*