# Boundary Layer Resolving Exact Difference Scheme for Solving Singularly Perturbed Convection-Diffusion-Reaction Equation

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

---

## Abstract

This paper considers the numerical treatment of singularly perturbed time-dependent convection-diffusion-reaction equation. The diffusion term of the equation is multiplied by a small perturbation parameter (ε), which takes an arbitrary value in the interval (0, 1]. For small values of ε, the solution of the equation exhibits an exponential boundary layer which makes it difficult to solve it analytically or using classical numerical methods. We proposed numerical schemes using the Crank–Nicolson method in time derivative discretization and the nonstandard finite difference method (exact finite difference method) in space derivative discretization on a uniform and piecewise uniform Shishkin mesh. The existence of unique discrete solutions and the stability of the schemes are discussed and proved. Uniform convergence of the schemes is proved. The formulated schemes converge uniformly with linear order of convergence. The method on Shishkin mesh possesses boundary layer resolving property. We validated the methods by considering two numerical examples for different values of ε and mesh length.


---

## Body

## 1. Introduction

Differential equations are commonly used for modeling processes in various applied mathematics areas. It relates one or more functions and their derivatives; the functions represent physical quantities, the derivatives represent their rates of change, and the differential equations define the relationship between the two. Many authors are working on the analytical and numerical solutions of different types of differential equations that arise from different application areas[1, 2].The convection-diffusion-reaction equations provide a very useful and important mathematical model in a wide range of disciplines. The application areas include fluid flow with high Reynolds numbers, convective heat transport problems with large Peclet numbers, turbulence models, simulation of oil extraction from underground reservoirs, financial modeling of option pricing, the water quality problem in river networks, drift diffusion equations of semiconductor device modeling, atmospheric pollution; and so on [3]. In many of these applications, the unknown variables in the governing equation represent physical quantities that cannot take negative values, such as concentrations of chemical compounds, pollutants, and population [4]. In the listed scientific areas, the models are most often described by singularly perturbed differential equations.Convection dominated reaction-diffusion equations or singularly perturbed convection-diffusion-equations are differential equations in which the diffusion term (or the highest order derivative term) of the equation is multiplied by a small perturbation parameterε. In general, the solutions of such equations exhibit boundary layer behaviour. Thus, the solutions of singularly perturbed equations possess a multiscale property. The presence of both slow and fast phenomena makes the problem stiff from the standpoint of an approximate solution.To obtain an approximate solution for singularly perturbed problems, the methods are mainly classified into two classes: (i) the asymptotic (or analytical) approach and (ii) the numerical approach. The numerical approach provides quantitative information about the solution of a particular problem, while the asymptotic approach provides the qualitative behavior of the solution of a family of equations and semiquantitative information about the solution of any particular member of the family [5].There are two strategies for designing numerical schemes that give small errors inside the boundary layer region. The first approach is the class of fitted mesh methods, which consists of choosing a fine mesh in the boundary layer region and a coarse mesh in the outer layer region. The second approach is the fitted operator methods, in which the mesh remains uniform and the difference schemes reflect the qualitative behaviour of the solution inside the boundary layer region. A nice discussion of using one or both of the above strategies can be found in [3, 6–9].The fitted operator methods can be categorized into two classes based on the techniques of numerical scheme formulation. That is, the exponentially fitted operator FDM and the nonstandard FDM. In nonstandard FDM the discretized mesh is uniform and a complicated denominator function is developed for the diffusion part of the equation using the exact finite difference technique to handle the effect of the perturbation parameter. This method can be extended to solve higher dimensional and nonlinear singularly perturbed problems. Different authors have published papers on the designing and implementing of nonstandard FDM on uniform mesh for different classes of model equations [10]. In the exponentially fitted operator FDM, the discretized mesh is uniform and the exponential fitting parameter is induced using the asymptotic solution behaviour to control the rapid growth or decay of the solution in the boundary layer region. The main drawback of the exponentially fitted operator FDM is that it does not work for nonlinear equations and higher dimensional equations.In the case of the fitted mesh methods, one uses nonuniform grids, which are fine in the boundary layer region and coarse outside the layer. In designing appropriate numerical methods using the fitted mesh method, one must have some a priori information about sensitive regions where the solution undergoes some rapid changes. Different mesh generation techniques have been developed so far. The well-known meshes are the exponentially graded Bakhvalov mesh and the piecewise uniform Shishkin mesh [6].Standard numerical methods on uniform mesh fail to give accurate results and create oscillations in the computed solution due to the presence of a small perturbation parameterε on the diffusion term [11, 12]. To avoid the oscillations and obtain an accurate solution using the standard numerical methods on uniform mesh, it requires an unacceptably large number of mesh points when ε is very small. So, this becomes unpractical due to round-off error and the limited size of the computer’s memory. To obtain an accurate and uniformly convergent numerical solution on a small number of mesh points, authors use the fitted mesh methods or the fitted operator methods.Nowadays, the numerical treatment of singularly perturbed convection-diffusion-equations has attracted great attention [4]. The authors in [12–16] used the fitted mesh methods for solving different classes of singularly perturbed problems, and the authors in [17–19] used the fitted operator methods for developing uniformly convergent numerical schemes.The nonstandard FDM on Shishkin mesh has better accuracy and order of convergence than the equivalent standard FDM on Shishkin mesh [20]. In addition, the uniform convergence analysis of the nonstandard FDMs were restricted to uniform mesh [7]. In [21], He and Wang developed a new form of the nonstandard FDM for solving a convection dominated diffusion equation. In the study of [20], Woldaregay developed numerical schemes using the nonstandard FDM on uniform mesh and Shishkin mesh for solving singularly perturbed boundary value problem. Motivated by the results in [20, 21], in this paper we proposed a uniformly convergent numerical scheme for solving the problem in (1). The proposed scheme uses the Crank–Nicolson method for time derivative discretization and the midpoint upwind nonstandard FDM on the uniform mesh and the Shishkin mesh for the spatial derivative discretization. Furthermore, we discuss the convergence analysis of the proposed scheme. The proposed nonstandard FDM on the Shishkin mesh gives a boundary layer resolving result.


### 1.1. Organization of the Paper

The remaining parts of the paper are organised as follows: in Section2, the problem under consideration is defined; in Section 3, the numerical schemes and their convergence analysis are given. In Section 4, the numerical result and discussion are included; in Section 5, the conclusion is given.


### 1.2. Notations

Nx and Nt denotes the number of mesh intervals in space and time discretization, respectively; C is denoted for positive constant independent of ε and Nx. The norm . used in this article is the maximum or supremum norm.


## 1.1. Organization of the Paper

The remaining parts of the paper are organised as follows: in Section2, the problem under consideration is defined; in Section 3, the numerical schemes and their convergence analysis are given. In Section 4, the numerical result and discussion are included; in Section 5, the conclusion is given.


## 1.2. Notations

Nx and Nt denotes the number of mesh intervals in space and time discretization, respectively; C is denoted for positive constant independent of ε and Nx. The norm . used in this article is the maximum or supremum norm.


## 2. Continuous Problem

On the domainD=Ωx×Ωt=0,1×0,T with the boundary ∂D=D¯−D where D¯=Ω¯x×Ω¯t=0,1×0,T, we considered a singularly perturbed time dependent convection-diffusion-reaction equation.(1)∂∂t+Lεux,t=fx,t,x,t∈D,ux,t=u0x,x∈Ωx,u0,t=ϕlt,t∈Ωt,u1,t=ϕrt,t∈Ωt,where Lεux,t=−ε∂2u/∂x2+ax∂u/∂x+bxux,t and ε∈0,1 is the perturbation parameter.In (1), the term −ε∂2/∂x2ux,t is called the diffusion term, ∂/∂t+ax∂/∂xux,t is called the convection term and bxux,t is called the reaction term [9]. The functions ax,bx, fx,t and ϕlt,u0x,ϕrt are assumed to be sufficiently smooth and bounded. We assume that the coefficient functions ax and bx satisfies(2)ax≥α>0,bx≥β>0,x∈Ω¯x.where α and β are constants, for the existence of a boundary layer on the right side of the spatial domain.The existence and uniqueness of the solution of (1) can be established by assuming that the data is Holder continuous and by imposing appropriate compatibility conditions at the corner points. Using the assumptions of sufficiently smoothness of u0x,ϕlt and ϕrt. The required compatibility condition at the corner points are ϕl0=u00,ϕr0=u01 and(3)dϕl0dt−εd2u00dx2+a0du00dx+b0u00=f0,0,dϕr0dt−εd2u01dx2+a1du01dx+b1u01=f1,0.So that the data matches at the corner points0,0 and 1,0. For ax,bx and fx,t be continuous on domain D, then (1) has unique solution ux,t∈C2D∩C0D¯.Lemma 1.
The solutionux,t of (1) satisfies the following estimate:(4)ux,t−u0x≤Ct,ux,t−ϕlt≤C1−x,x,t∈D¯.where C is constant independent of ε.Proof.
The result follows directly from the compatibility condition. See the detailed proof in [22].
The reduced problem obtained by settingε=0 in (1) is given by the following equation:(5)∂u0∂t+ax∂u0∂x+bxu0x,t=fx,t,x,t∈D,u0x,t=u0x,x∈Ωx,u00,t=ϕlt,t∈Ωt,u01,t≠ϕrt,t∈Ωt.
For small values ofε, the solution of (1) is very close to the solution of (5).Lemma 2 (The maximum principle).
Letu be the solution of (1) which satisfies ux,t≥0,x,t∈∂D and ∂/∂t+Lεux,t≥0,x,t∈D implies that ux,t≥0,x,t∈D¯.Proof.
The proof is by contradiction. Suppose that there existx∗,t∗∈D¯ such that ux∗,t∗=minx,t∈D¯ux,t<0. It is clear that x∗,t∗∉∂D i.e. x∗,t∗∈D. Since ux∗,t∗=minx,t∈D¯ux,t, which implies uxx∗,t∗=0, utx∗,t∗=0 and uxxx∗,t∗≥0. Giving that ∂/∂t+Lεux∗,t∗<0, which is a contradiction to the assumption made ∂/∂t+Lεux∗,t∗≥0,∀x,t∈D. Therefore, ux,t≥0,∀x,t∈D¯.Lemma 3 (Stability estimate).
The solutionux,t of (1) satisfies the bounded(6)ux,t≤β−1f+u∂D,where β is lower bound of bx.Proof.
By defining the barrier functionsϑ±x,t as ϑ±x,t=β−1f+u∂D±ux,t, on the boundary lines we obtain ϑ±x,0=β−1f+u∂D±u0x≥0,;ϑ±0,t=β−1f+u∂D±u0,t≥0, and ϑ±1,t=β−1f+u∂D±u1,t≥0.
On the domainD=Ωx×Ωt, we have the following equation:(7)∂∂t+Lεϑ±x,t=ϑt±x,t−εϑxx±x,t+axϑx±x,t+bxϑ±x,t=∂∂tβ−1f+u∂D±ux,t−ε∂2∂x2β−1f+u∂D±ux,t+ax∂∂xβ−1f+u∂D±ux,t+bxβ−1f+u∂D±ux,t=bxβ−1f+u∂D±fx,t≥0,sincebx≥β>0.
With the hypothesis of the maximum principle, we obtainϑ±x,t≥0,∀x,t∈D¯. Hence, the proof is completed.
Using Shishkin decomposition in [12, 23], the solution ux,t of (1) can be decomposed into regular and singular component as follows:(8)ux,t=vx,t+wx,t.where vx,t denotes the regular component and wx,t denotes the singular component. The regular component satisfies the following nonhomogeneous equation:(9)∂∂t+Lεvx,t=fx,t,x,t∈D,vx,0=u0x,x∈Ω¯x,v0,t=u0,t,t∈Ω¯t,v1,t≠u1,t,t∈Ω¯t,and the singular component satisfies the following homogeneous equation:(10)∂∂t+Lεwx,t=0,x,t∈D,wx,0=0,x∈Ω¯x,w0,t=0,t∈Ω¯t,w1,t=u1,t−v1,t,t∈Ω¯t.Lemma 4 (see [14]).
The derivatives of the regular and singular component solutions respectively satisfies the bound(11)∂kvx,t∂xk≤C,k=0,1,2,3,∂kwx,t∂xk≤Cε−kexp−α1−xε,k=0,1,2,3,and the derivatives with respect to x and t of the solution ux,t of (1) satisfies the bounded(12)∂k∂lux,t∂xk∂tl≤C1+ε−kexp−α1−xε,0≤k≤4,0≤l≤2,where α is lower bound of ax.


## 3. The Numerical Schemes

### 3.1. Time Semidiscretization

We divide the time domainΩt=0,T into Nt subintervals using the grid points t0=0,tj=jht,j=1,2,…,Nt, where ht=T/Nt. We denote uj+1x for the approximation of ux,t in temporal discretization at j+1 grid point. Using the Crank–Nicolson finite difference method for semidiscretizing (1), we obtain a system of Nt−1 boundary value problems(13)1+ht2Lεhtuj+1x=1−ht2Lεhtujx+ht2fx,tj+1+fx,tj,x∈Ωx,uj+10=ϕltj+1,uj+11=ϕrtj+1,forj=1,2,…,Nt−1, where, Lεhtuj+1x=−εd2uj+1x/dx2+axduj+1x/dx+bxuj+1x. The semidiscrete scheme in (13) satisfies the maximum principle which is stated as follows.Lemma 5 (Semidiscrete maximum principle).
Letuj+1 be smooth function on Ω¯x. If uj+10≥0,uj+11≥0 and 1+ht/2Lεhtuj+1x≥0,x∈Ωx, then uj+1x≥0,x∈Ω¯x.Proof.
Similar techniques as the proof in Lemma2 gives the required result.
Next, we find a bound for the truncation error in temporal semidiscretization. Let us denote the local truncation error for each time step byej+1x:=ux,tj+1−uj+1x.Lemma 6 (Local error estimate).
The local error in the temporal semidiscretization is bounded as follows:(14)ej+1x≤Cht3.Proof.
Using the Taylor series expansion we obtain the local truncation error of the semidiscrete scheme as follows:(15)1+ht2Lεhtux,tj+1−1+ht2Lεhtuj+1x=1+ht2Lεhtej+1x≤Cht3.
Sinceej+10=0 and ej+11=0 and 1+ht/2Lεhtej+1x≥0. Using the semidiscrete maximum principle we obtain the required bound.
The global error estimate up totj+1 time step is bounded as follows:(16)Ej+1≤Cht2,j=1,2,…,Nt−1.


### 3.2. Spatial Discretization

In this subsection, we approximate the spatial derivatives using nonstandard FDM on uniform mesh and on the Shishkin mesh. Moreover, we discuss the convergence analysis of the schemes.


#### 3.2.1. Exact Finite Difference

For constructing the exact finite difference scheme for convection diffusion equation, we follow the techniques developed in [4, 7]. Consider a constant coefficient homogeneous convection diffusion equations(17)−εu″x+αu′x+βux=0,(18)−εu″x+αu′x=0,where ax≥α and bx≥β. The characteristic equation of (17) is given by the following equation:(19)−ελ2+αλ+β=0,and its solution is given as λ1,2=−α±α2+4εβ/−2ε. The constant coefficient equation in (17) has two independent solutions which is given as follows:(20)ux=A1eλ1x+A2eλ2x.Consider a uniform meshΩNx=xi=x0+ihxi=1Nx, x0=0,xNx=1,hx=1/Nx where Nx is the number of mesh intervals in spatial discretization. We want to develop a difference equation which has the same general solution as the differential equation in (18) has at the mesh point xi, which is given by the following equation:(21)Ui=A1eλ1xi+A2eλ2xi.From the theory of difference equations for second order linear difference equations, we obtain the following equation:(22)Ui−1eλ1xi−1eλ2xi−1Uieλ1xieλ2xiUi+1eλ1xi+1eλ2xi+1=0.rearranging gives the following equation:(23)eλ2hx−eλ1hxUi−1−eλ2−λ1hx−eλ1−λ2hxUi+e−λ1hx−e−λ2hxUi+1=0.Simplifying and substituting the values ofλ1,2, we obtain the following equation:(24)expαhx2εUi−1−2coshhx2εα2+4εβUi+exp−αhx2εUi+1=0,which is an exact difference scheme for (18). Rearranging, we obtain the following equation:(25)−εUi−1−2Ui+Ui+1hxε/αexpαhx/ε−1+αUi−Ui−1hx=0.The denominator of the second derivative approximation isγ2:=hxε/αexpαhx/ε−1 and the denominator of the first derivative approximation is hx.


#### 3.2.2. Spatial Discretization via Nonstandard FDM on Uniform Mesh

We modifyγ2 for the variable coefficient differential equations of the form in (13) as the following form:(26)γi2=εhxaiexphxaiε−1.where ai is denoted for axi. Using the denominator function γi2 and applying the midpoint upwind finite difference discretization the required difference schemes becomes the following equation:(27)1+ht2Lε,muhx,htUi,j+1=1−ht2Lε,muhx,htUi,j+ht2fxi−1/2,tj+fxi−1/2,tj+1,U0,j+1=ϕltj+1,UNx,j+1=ϕrtj+1,fori=1,2,…,Nx−1, where(28)Lε,muhx,htUi,j+1=−εUi+1,j+1−2Ui,j+1+Ui−1,j+1γi2+ai−1/2Ui,j+1−Ui−1,j+1hx+bi−1/2Ui,j+1.and ai−1/2,bi−1/2 and fxi−1/2,tj denoted for axi+axi−1/2, bxi+bxi−1/2 and fxi,tj+fxi−1,tj/2, respectively.


#### 3.2.3. Stability and Convergence Analysis on Uniform Mesh

Here, we want to show that the discrete scheme in (27) satisfies the maximum principle, stability estimates, and uniform convergence.Lemma 7.
LetUi,j+1 be any mesh function satisfying U0,j+1≥0,UNx,j+1≥0. Then, 1+ht/2Lε,muhx,htUi,j+1≥0,i=1,2,…,Nx−1 implies that Ui,j+1≥0,i=0,1,…,Nx.Proof.
Suppose there existk∈0,1,…,N such that Uk,j+1=min0≤i≤NUi,j+1. Suppose that Uk,j+1<0 which implies k≠0,N. So, we have that Uk+1,j+1−Uk,j+1>0 and Uk,j+1−Uk−1,j+1<0. So, we obtain 1+ht/2Lε,muhx,htUk,j+1<0 for k=1,2,3,…,Nx−1. Thus, the supposition Ui,j+1<0, for i=0,1,…,Nx is wrong. Hence, we obtain Ui,j+1≥0,i=0,1,…,Nx.
Using the discrete maximum principle, we want to prove the discrete scheme satisfies the uniform stability result.Lemma 8.
The solutionUi,j+1 of the discrete scheme in (27) satisfy the bound.(29)Ui,j+1≤maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1.Proof.
LetQ=maxi1+ht/2Lε,muhx,htUi,j+1/1+ht/2β+maxϕltj+1,ϕrtj+1 and define barrier functions ϑi,j+1±by ϑi,j+1±=Q±Ui,j+1. At the boundary points, we have the following equation:(30)ϑ0,j+1±=Q±U0,j+1=maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕltj+1≥0,ϑN,j+1±=Q±UN,j+1=maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕrtj+1≥0.
On the discretized domainxi,0<i<Nx, we obtain the following equation:(31)1+ht2Lε,muhx,htϑi,j+1±=Q±Ui,j+1−εht2Q±Ui+1,j+1−2Q±Ui,j+1+Q±Ui−1,j+1γi2+ai−1/2ht2Q±Ui,j+1−Q±Ui−1,j+1h+bi−1/2ht2Q±Ui,j+1=1+ht2bi−1/2maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±1+ht2Lε,muhx,htUi,j+1≥0,
By the discrete maximum principle in Lemma7, we obtain ϑi,j+1±≥0,xi∈Ω¯Nx. Hence the required bound is satisfied.Lemma 9.
The discrete solutionUi,j+1 of the scheme in (27) satisfy the bound(32)Ui,j+1≤β−1Lε,muhx,htUi,j+1+maxϕltj+1,ϕrtj+1.Proof.
The proof is similar to the above Lemma8 proof by construction a barrier functions ϑi,j+1±=β−1Lε,muhx,htUi,j+1+maxϕltj+1,ϕrtj+1±Ui,j+1.
Let us denote the forward and backward finite differences operators in spatial variable as follows:(33)D+zj+1xi=zj+1xi+1−zj+1xihx,D−zj+1xi=zj+1xi−zj+1xi−1hx.respectively and the second order finite difference operator as follows:(34)D+D−γzj+1xi=zi−1,j+1−2zi,j+1+zi+1,j+1γi2=hx2D+D−zj+1xiγi2.Lemma 10.
For a fixed meshNx and for ε⟶0, it holds(35)limε⟶0max1≤i≤Nx−1ε−mexp−αxiε=0,limε⟶0max1≤i≤Nx−1ε−mexp−α1−xiε=0,m=1,2,3,…where xi=ihx,hx=1/Nx,i=1,2,…,Nx−1.Lemma 11 (see [24]).
For allε∈0,1, the following estimate holds(36)εhx2γi−1≤Chx.Theorem 1.
The truncation error of nonstandard FDM on uniform mesh satisfies the bound(37)Lε,muhx,htuj+1xi−Ui,j+1≤Chx1+maxxiexp−α1−xi/εε3.Proof.
The spatial discretization truncation error is given by the following equation:(38)Lε,muhx,htuj+1xi−Ui,j+1≤εd2uj+1xidx2−Dx+Dx−hx2γi2uj+1xi+ai−1/2ddx−Dx−uj+1xi=Cεd2dx2−Dx+Dx−uj+1xi+hx2γi2−1Dx+Dx−uj+1xi+Chxd2dx2uj+1xi≤Cεd2dx2−Dx+Dx−uj+1xi+Cεhx2γi2−1Dx+Dx−uj+1xi+Chxd2dx2uj+1xi.
Using the result in Lemma11, the truncation error bound becomes the following equation:(39)Lε,muhx,htuj+1xi−Ui,j+1≤Cεhx2d4dx4uj+1xi+Chxd2dx2uj+1xi.
Using the bound in Lemma4 into (39), we obtain the following equation:(40)Lε,muhx,htuj+1xi−Ui,j+1≤Cεhx21+ε−4exp−α1−xicε+Chx1+ε−2exp−α1−xiε≤Chx2ε+ε−3exp−α1−xiε+Chx1+ε−2exp−α1−xiε≤Chx1+maxxiε−3exp−α1−xiε,since,ε−3≥ε−2.Theorem 2.
The nonstandard FDM on uniform mesh satisfies the error bound(41)uj+1xi−Ui,j+1≤Chx.Proof.
Using the results in Lemma10 and applying the bound in Lemma 9, we obtain the error bound.Theorem 3.
Letu and U be the solution of the problem in (1) and (27), respectively, then the nonstandard FDM on uniform mesh satisfies the error bound(42)supε∈0,1u−U≤Chx+ht2.Proof.
Combining the error bound for the temporal and spatial discretization in (16) and Theorem 2 gives the required bound.


### 3.3. Spatial Discretization via Nonstandard FDM on Shishkin Mesh

The development of nonstandard FDM and the uniform convergence of the schemes were restricted on a uniform mesh. In this paper, we extend the scheme and uniform convergence of the scheme on piecewise-uniform Shishkin mesh. Letxii=0Nx be the discretization of the spatial domain Ω¯x where Nx be the number of grid points be an even positive integer. For each i≥1, we define hi=xi−xi−1. Since the considered problem exhibits the right boundary layer, we set a mesh transition parameter τ as follows:(43)τ=min12,σεβlnNx.Which divides the domain [0, 1] in to outer layer regionΩx1=0,1−τ and boundary layer region Ωx2=1−τ,1, where σ is a constant that represents the order of the scheme and β is the lower bound of bx. Here, we set a uniform mesh Nx/2 in Ωx1 with mesh spacing H=21−τ/Nx and similarly a uniform mesh Nx/2 is placed in Ωx2 with the mesh spacing h=2τ/Nx. Hence, the mesh point xi is given as follows:(44)xi=21−τNxi,i=0,1,2,…,Nx2,1−τ+2τNxi−Nx2,i=Nx2+1,…,Nx.The nonstandard FDM are commonly applied on uniform mesh, so using it on layer adapted mesh requires the following change on the denominator function:(45)γ1,i=εaiexphiaiε−1,γ2,i=εaiexphi+1aiε−1.For the mesh functionuj+1x at the grid points xi, we denote the second derivatives term approximation as follows:(46)−εd2dx2uj+1xi≊−2εhi+hi+1Ui+1,j+1−Ui,j+1γ2,i−Ui,j+1−Ui−1,j+1γ1,i=−εDx+Dx−γuj+1xi.Fori=1,2,…,Nx−1, the midpoint upwind nonstandard FDM on Shishkin mesh is given as follows:(47)1+ht2Lε,smNx,htUi,j+1=1−ht2Lε,smNx,htUi,j+ht2fxi−1/2,tj+fxi−1/2,tj+1,U0,j+1=ϕltj+1,UNx,j+1=ϕrtj+1,where(48)Lε,smNx,htUi,j+1=−2εUi+1,j+1−Ui,j+1γ2,ihi+hi+1−Ui,j+1−Ui−1,j+1γ1,ihi+hi+1+ai−1/2Ui,j+1−Ui−1,j+1hi+bi−1/2Ui,j+1.


#### 3.3.1. Stability and Convergence Analysis on the Shishkin Mesh

LetH be the mesh width in the outer layer region 0,1−τ and h be the mesh width in layer region 1−τ,1. Then, the mesh widths satisfies the following conditions:(49)Nx−1≤H≤2Nx−1,h=σεβNx−1lnNx≤Nx−1.Lemma 12 (The discrete maximum principle).
LetUi,j+1 be any mesh function satisfying U0,j+1≥0,UNx,j+1≥0. Then, 1+ht/2Lε,smNx,htUi,j+1≥0,i=1,2,…,Nx−1 implies that Ui,j+1≥0,∀i=0,1,…,Nx.Proof.
The proof is similar to the proof of Lemma7.Lemma 13.
The solutionUi,j+1 of the discrete scheme in (47) satisfy the bound.(50)Ui,j+1≤maxi1+ht/2Lε,smNx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1.Proof.
LetQ=maxi1+ht/2Lε,smhx,htUi,j+1/1+ht/2β+maxϕltj+1,ϕrtj+1 and define barrier functions ϑi,j+1± by ϑi,j+1±=Q±Ui,j+1. At the boundary points, we have the following equation:(51)ϑ0,j+1±=Q±U0,j+1=maxi1+ht/2Lε,smNx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕltj+1≥0,ϑN,j+1±=Q±UNx,j+1=maxi1+ht/2Lε,smNx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕrtj+1≥0.
On the discretized spatial domainxi,i=1,2,…,Nx−1, we have the following equation:(52)1+ht2Lε,smNx,htϑi,j+1±=Q±Ui,j+1−2εht2Q±Ui+1,j+1−Q±Ui,j+1γ2,ihi+hi+1−Q±Ui,j+1−Q±Ui−1,j+1γ1,ihi+hi+1+ai−1/2ht2Q±Ui,j+1−Q±Ui−1,j+1hxi+bi−1/2ht2Q±Ui,j+1=1+ht2bi−1/2maxi1+ht/2Lε,smhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±1+ht2Lε,smNx,htUi,j+1≥0,sincebi−1/2≥β,
Using the maximum principle in Lemma12, we obtain ϑi,j+1±≥0,∀xi∈Ω¯Nx. Hence, the required bound is satisfied.Lemma 14.
The discrete solutionUi,j+1 of the scheme in (47) satisfy the bound(53)Ui,j+1≤β−1Lε,smNx,htUi,j+1+maxϕltj+1,ϕrtj+1.Proof.
The proof is similar to the above-given lemma proof by the construction of barrier functionsϑi,j+1±=β−1Lε,smNx,htUi,j+1+maxϕltj+1,ϕrtj+1±Ui,j+1.


#### 3.3.2. Decomposition of the Discrete Solution

Here, we decompose the numerical solutionUi,j+1 into regular and singular components using the Shishkin decomposition in [12] as follows:(54)Ui,j+1=Vi,j+1+Wi,j+1,where the regular part satisfied the following nonhomogeneous equation:(55)1+ht2Lε,smNx,htVi,j+1=1−ht2Lε,smNx,htVi,j+ht2fxi−1/2,tj+fxi−1/2,tj+1,Vj+10=vj+10,Vj+11=vj+11,and the singular component satisfies the following homogeneous equation:(56)1+ht2Lε,smNx,htWi,j+1=0,Wj+10=wj+10,Wj+11=wj+11,The error in the computed solution can also be decomposed as follows:(57)Ui,j+1−uj+1xi=Vi,j+1−vj+1xi+Wi,j+1−wj+1xi.First, let us consider the caseτ=1/2 which is a uniform mesh case say the mesh size is h. The truncation error becomes the following equation:(58)Lε,smNx,htuj+1xi−Ui,j+1=εd2dx2−Dx+Dx−γuj+1xi+ai−1/2ddx−Dx−uj+1xi.Since(59)εDx+Dx−γuj+1xi=2εUi+1,j+1−Ui,j+1γ2,ihi+hi+1−Ui,j+1−Ui−1,j+1γ1,ihi+hi+1=εUi+1,j+1−2Ui,j+1+Ui−1,j+1hε/aiexphai/ε−1.From the truncation error in uniform mesh, we obtain the following equation:(60)Lε,smNx,htuj+1xi−Ui,j+1≤Cεh2d4dx4uj+1xi+Chd2dx2uj+1xi≤Ch.Using Lemma14, we obtain uj+1xi−Ui,j+1≤Ch.For the caseτ=σε/βlnNx, we estimate the error in the regular and singular solutions separately. Here, σε/βlnNx≤0.5 implies that 1/ε≤ClnNx.Theorem 4.
The error in the regular component satisfies the bound(61)vj+1xi−Vi,j+1≤Ch.Proof.
On the outer layer region0,1−τ, it is clear that the mesh is uniform i.e. hi=hi+1=H,i=1,2,…,Nx/2 . So the truncation error in the regular component is given as follows:(62)Lε,smNx,htvj+1xi−Vi,j+1=−εd2dx2−Dx+Dx−γvj+1xi+ai−1/2ddx−Dx−vj+1xi≤CεNx−2d4dx4vj+1ξ+CNx−1d2dx2vj+1ξ,since,Nx−1≤H≤2Nx−1≤CεNx−2+CNx−1≤CNx−2lnNx+CNx−1≤CNx−1.
In the boundary layer region1−τ,1, the mesh is uniform i.e., hi=hi+1=h≤Nx−1,i=Nx/2+1,…,Nx−1. So, the truncation error is given as follows:(63)Lε,smNx,htvj+1xi−Vi,j+1=−εd2dx2−Dx+Dx−γvj+1xi+ai−1/2ddx−Dx−vj+1xi≤CεNx−2d4dx4vj+1ξ+CNx−1d2dx2vj+1ξ,since,h≤Nx−1≤CεNx−2+CNx−1≤CNx−2lnNx+CNx−1≤CNx−1.
Using Lemma14, we obtain vj+1xi−Vi,j+1≤CNx−1.Theorem 5.
The error in the singular component satisfies the bound(64)wj+1xi−Wi,j+1≤CNx−1lnNx2.Proof.
In outer layer region0,1−τ, we have hi=hi+1=H,i=1,2,…,Nx/2. So the truncation error becomes(65)Lε,smNx,htwj+1xi−Wi,j+1=−εd2dx2−Dx+Dx−γwj+1xi+ai−1/2ddx−Dx−wj+1xi≤CεNx−2d4dx4wj+1xi+CNx−1d2dx2wj+1xi,sinceNx−1≤H≤2Nx−1≤CεNx−2ε−4exp−α1−xiε+CNx−1ε−2exp−α1−xiε≤CNx−2ε−3+Nx−1ε−2exp−α1−xiε≤CNx−2ε−3+Nx−1ε−2Nx−1,sinceτ=σεβlnNx,exp−ατε≤CNx−1,≤CNx−2lnNx2,since1ε≤ClnNx.
In the layer region1−τ,1, the mesh is also uniform i.e. hi=hi+1=h≤Nx−1,i=Nx/2+1,…,Nx−1. Hence, we have the following equation:(66)Lε,smNx,htwj+1xi−Wi,j+1=−εd2dx2−Dx+Dx−γwj+1xi+ai−1/2ddx−Dx−wj+1xi≤CεNx−2d4dx4wj+1xi+CNx−1d2dx2wj+1xi,sinceh≤Nx−1≤CεNx−2ε−4exp−α1−xiε+CNx−1ε−2exp−α1−xiε≤CNx−2ε−3+Nx−1ε−2exp−α1−xiε,≤CNx−1ε−2,sincemaxiexp−α1−xiε≤1,≤CNx−1lnNx2,since1ε≤ClnNx.
Hence, we haveLε,smNx,htwj+1xi−Wi,j+1≤maxCNx−2lnNx2,CNx−1lnNx2. Using Lemma 14, we obtain wj+1xi−Wi,j+1≤CNx−1lnNx2.Theorem 6.
The nonstandard FDM on the Shishkin mesh satisfies the error bound(67)uj+1xi−Ui,j+1≤CNx−1lnNx2.Proof.
Combining the error estimate in the regular and singular components in Theorem4 and 5 gives the required bound.Theorem 7.
Letu and U be the solution of problem in (1) and (47), respectively, then the nonstandard FDM on the Shishkin mesh satisfies the error bound(68)supε∈0,1u−U≤CNx−1lnNx2+ht2.Proof.
Using the error bound for the temporal and spatial discretization gives the required bound.


## 3.1. Time Semidiscretization

We divide the time domainΩt=0,T into Nt subintervals using the grid points t0=0,tj=jht,j=1,2,…,Nt, where ht=T/Nt. We denote uj+1x for the approximation of ux,t in temporal discretization at j+1 grid point. Using the Crank–Nicolson finite difference method for semidiscretizing (1), we obtain a system of Nt−1 boundary value problems(13)1+ht2Lεhtuj+1x=1−ht2Lεhtujx+ht2fx,tj+1+fx,tj,x∈Ωx,uj+10=ϕltj+1,uj+11=ϕrtj+1,forj=1,2,…,Nt−1, where, Lεhtuj+1x=−εd2uj+1x/dx2+axduj+1x/dx+bxuj+1x. The semidiscrete scheme in (13) satisfies the maximum principle which is stated as follows.Lemma 5 (Semidiscrete maximum principle).
Letuj+1 be smooth function on Ω¯x. If uj+10≥0,uj+11≥0 and 1+ht/2Lεhtuj+1x≥0,x∈Ωx, then uj+1x≥0,x∈Ω¯x.Proof.
Similar techniques as the proof in Lemma2 gives the required result.
Next, we find a bound for the truncation error in temporal semidiscretization. Let us denote the local truncation error for each time step byej+1x:=ux,tj+1−uj+1x.Lemma 6 (Local error estimate).
The local error in the temporal semidiscretization is bounded as follows:(14)ej+1x≤Cht3.Proof.
Using the Taylor series expansion we obtain the local truncation error of the semidiscrete scheme as follows:(15)1+ht2Lεhtux,tj+1−1+ht2Lεhtuj+1x=1+ht2Lεhtej+1x≤Cht3.
Sinceej+10=0 and ej+11=0 and 1+ht/2Lεhtej+1x≥0. Using the semidiscrete maximum principle we obtain the required bound.
The global error estimate up totj+1 time step is bounded as follows:(16)Ej+1≤Cht2,j=1,2,…,Nt−1.


## 3.2. Spatial Discretization

In this subsection, we approximate the spatial derivatives using nonstandard FDM on uniform mesh and on the Shishkin mesh. Moreover, we discuss the convergence analysis of the schemes.


### 3.2.1. Exact Finite Difference

For constructing the exact finite difference scheme for convection diffusion equation, we follow the techniques developed in [4, 7]. Consider a constant coefficient homogeneous convection diffusion equations(17)−εu″x+αu′x+βux=0,(18)−εu″x+αu′x=0,where ax≥α and bx≥β. The characteristic equation of (17) is given by the following equation:(19)−ελ2+αλ+β=0,and its solution is given as λ1,2=−α±α2+4εβ/−2ε. The constant coefficient equation in (17) has two independent solutions which is given as follows:(20)ux=A1eλ1x+A2eλ2x.Consider a uniform meshΩNx=xi=x0+ihxi=1Nx, x0=0,xNx=1,hx=1/Nx where Nx is the number of mesh intervals in spatial discretization. We want to develop a difference equation which has the same general solution as the differential equation in (18) has at the mesh point xi, which is given by the following equation:(21)Ui=A1eλ1xi+A2eλ2xi.From the theory of difference equations for second order linear difference equations, we obtain the following equation:(22)Ui−1eλ1xi−1eλ2xi−1Uieλ1xieλ2xiUi+1eλ1xi+1eλ2xi+1=0.rearranging gives the following equation:(23)eλ2hx−eλ1hxUi−1−eλ2−λ1hx−eλ1−λ2hxUi+e−λ1hx−e−λ2hxUi+1=0.Simplifying and substituting the values ofλ1,2, we obtain the following equation:(24)expαhx2εUi−1−2coshhx2εα2+4εβUi+exp−αhx2εUi+1=0,which is an exact difference scheme for (18). Rearranging, we obtain the following equation:(25)−εUi−1−2Ui+Ui+1hxε/αexpαhx/ε−1+αUi−Ui−1hx=0.The denominator of the second derivative approximation isγ2:=hxε/αexpαhx/ε−1 and the denominator of the first derivative approximation is hx.


### 3.2.2. Spatial Discretization via Nonstandard FDM on Uniform Mesh

We modifyγ2 for the variable coefficient differential equations of the form in (13) as the following form:(26)γi2=εhxaiexphxaiε−1.where ai is denoted for axi. Using the denominator function γi2 and applying the midpoint upwind finite difference discretization the required difference schemes becomes the following equation:(27)1+ht2Lε,muhx,htUi,j+1=1−ht2Lε,muhx,htUi,j+ht2fxi−1/2,tj+fxi−1/2,tj+1,U0,j+1=ϕltj+1,UNx,j+1=ϕrtj+1,fori=1,2,…,Nx−1, where(28)Lε,muhx,htUi,j+1=−εUi+1,j+1−2Ui,j+1+Ui−1,j+1γi2+ai−1/2Ui,j+1−Ui−1,j+1hx+bi−1/2Ui,j+1.and ai−1/2,bi−1/2 and fxi−1/2,tj denoted for axi+axi−1/2, bxi+bxi−1/2 and fxi,tj+fxi−1,tj/2, respectively.


### 3.2.3. Stability and Convergence Analysis on Uniform Mesh

Here, we want to show that the discrete scheme in (27) satisfies the maximum principle, stability estimates, and uniform convergence.Lemma 7.
LetUi,j+1 be any mesh function satisfying U0,j+1≥0,UNx,j+1≥0. Then, 1+ht/2Lε,muhx,htUi,j+1≥0,i=1,2,…,Nx−1 implies that Ui,j+1≥0,i=0,1,…,Nx.Proof.
Suppose there existk∈0,1,…,N such that Uk,j+1=min0≤i≤NUi,j+1. Suppose that Uk,j+1<0 which implies k≠0,N. So, we have that Uk+1,j+1−Uk,j+1>0 and Uk,j+1−Uk−1,j+1<0. So, we obtain 1+ht/2Lε,muhx,htUk,j+1<0 for k=1,2,3,…,Nx−1. Thus, the supposition Ui,j+1<0, for i=0,1,…,Nx is wrong. Hence, we obtain Ui,j+1≥0,i=0,1,…,Nx.
Using the discrete maximum principle, we want to prove the discrete scheme satisfies the uniform stability result.Lemma 8.
The solutionUi,j+1 of the discrete scheme in (27) satisfy the bound.(29)Ui,j+1≤maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1.Proof.
LetQ=maxi1+ht/2Lε,muhx,htUi,j+1/1+ht/2β+maxϕltj+1,ϕrtj+1 and define barrier functions ϑi,j+1±by ϑi,j+1±=Q±Ui,j+1. At the boundary points, we have the following equation:(30)ϑ0,j+1±=Q±U0,j+1=maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕltj+1≥0,ϑN,j+1±=Q±UN,j+1=maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕrtj+1≥0.
On the discretized domainxi,0<i<Nx, we obtain the following equation:(31)1+ht2Lε,muhx,htϑi,j+1±=Q±Ui,j+1−εht2Q±Ui+1,j+1−2Q±Ui,j+1+Q±Ui−1,j+1γi2+ai−1/2ht2Q±Ui,j+1−Q±Ui−1,j+1h+bi−1/2ht2Q±Ui,j+1=1+ht2bi−1/2maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±1+ht2Lε,muhx,htUi,j+1≥0,
By the discrete maximum principle in Lemma7, we obtain ϑi,j+1±≥0,xi∈Ω¯Nx. Hence the required bound is satisfied.Lemma 9.
The discrete solutionUi,j+1 of the scheme in (27) satisfy the bound(32)Ui,j+1≤β−1Lε,muhx,htUi,j+1+maxϕltj+1,ϕrtj+1.Proof.
The proof is similar to the above Lemma8 proof by construction a barrier functions ϑi,j+1±=β−1Lε,muhx,htUi,j+1+maxϕltj+1,ϕrtj+1±Ui,j+1.
Let us denote the forward and backward finite differences operators in spatial variable as follows:(33)D+zj+1xi=zj+1xi+1−zj+1xihx,D−zj+1xi=zj+1xi−zj+1xi−1hx.respectively and the second order finite difference operator as follows:(34)D+D−γzj+1xi=zi−1,j+1−2zi,j+1+zi+1,j+1γi2=hx2D+D−zj+1xiγi2.Lemma 10.
For a fixed meshNx and for ε⟶0, it holds(35)limε⟶0max1≤i≤Nx−1ε−mexp−αxiε=0,limε⟶0max1≤i≤Nx−1ε−mexp−α1−xiε=0,m=1,2,3,…where xi=ihx,hx=1/Nx,i=1,2,…,Nx−1.Lemma 11 (see [24]).
For allε∈0,1, the following estimate holds(36)εhx2γi−1≤Chx.Theorem 1.
The truncation error of nonstandard FDM on uniform mesh satisfies the bound(37)Lε,muhx,htuj+1xi−Ui,j+1≤Chx1+maxxiexp−α1−xi/εε3.Proof.
The spatial discretization truncation error is given by the following equation:(38)Lε,muhx,htuj+1xi−Ui,j+1≤εd2uj+1xidx2−Dx+Dx−hx2γi2uj+1xi+ai−1/2ddx−Dx−uj+1xi=Cεd2dx2−Dx+Dx−uj+1xi+hx2γi2−1Dx+Dx−uj+1xi+Chxd2dx2uj+1xi≤Cεd2dx2−Dx+Dx−uj+1xi+Cεhx2γi2−1Dx+Dx−uj+1xi+Chxd2dx2uj+1xi.
Using the result in Lemma11, the truncation error bound becomes the following equation:(39)Lε,muhx,htuj+1xi−Ui,j+1≤Cεhx2d4dx4uj+1xi+Chxd2dx2uj+1xi.
Using the bound in Lemma4 into (39), we obtain the following equation:(40)Lε,muhx,htuj+1xi−Ui,j+1≤Cεhx21+ε−4exp−α1−xicε+Chx1+ε−2exp−α1−xiε≤Chx2ε+ε−3exp−α1−xiε+Chx1+ε−2exp−α1−xiε≤Chx1+maxxiε−3exp−α1−xiε,since,ε−3≥ε−2.Theorem 2.
The nonstandard FDM on uniform mesh satisfies the error bound(41)uj+1xi−Ui,j+1≤Chx.Proof.
Using the results in Lemma10 and applying the bound in Lemma 9, we obtain the error bound.Theorem 3.
Letu and U be the solution of the problem in (1) and (27), respectively, then the nonstandard FDM on uniform mesh satisfies the error bound(42)supε∈0,1u−U≤Chx+ht2.Proof.
Combining the error bound for the temporal and spatial discretization in (16) and Theorem 2 gives the required bound.


## 3.2.1. Exact Finite Difference

For constructing the exact finite difference scheme for convection diffusion equation, we follow the techniques developed in [4, 7]. Consider a constant coefficient homogeneous convection diffusion equations(17)−εu″x+αu′x+βux=0,(18)−εu″x+αu′x=0,where ax≥α and bx≥β. The characteristic equation of (17) is given by the following equation:(19)−ελ2+αλ+β=0,and its solution is given as λ1,2=−α±α2+4εβ/−2ε. The constant coefficient equation in (17) has two independent solutions which is given as follows:(20)ux=A1eλ1x+A2eλ2x.Consider a uniform meshΩNx=xi=x0+ihxi=1Nx, x0=0,xNx=1,hx=1/Nx where Nx is the number of mesh intervals in spatial discretization. We want to develop a difference equation which has the same general solution as the differential equation in (18) has at the mesh point xi, which is given by the following equation:(21)Ui=A1eλ1xi+A2eλ2xi.From the theory of difference equations for second order linear difference equations, we obtain the following equation:(22)Ui−1eλ1xi−1eλ2xi−1Uieλ1xieλ2xiUi+1eλ1xi+1eλ2xi+1=0.rearranging gives the following equation:(23)eλ2hx−eλ1hxUi−1−eλ2−λ1hx−eλ1−λ2hxUi+e−λ1hx−e−λ2hxUi+1=0.Simplifying and substituting the values ofλ1,2, we obtain the following equation:(24)expαhx2εUi−1−2coshhx2εα2+4εβUi+exp−αhx2εUi+1=0,which is an exact difference scheme for (18). Rearranging, we obtain the following equation:(25)−εUi−1−2Ui+Ui+1hxε/αexpαhx/ε−1+αUi−Ui−1hx=0.The denominator of the second derivative approximation isγ2:=hxε/αexpαhx/ε−1 and the denominator of the first derivative approximation is hx.


## 3.2.2. Spatial Discretization via Nonstandard FDM on Uniform Mesh

We modifyγ2 for the variable coefficient differential equations of the form in (13) as the following form:(26)γi2=εhxaiexphxaiε−1.where ai is denoted for axi. Using the denominator function γi2 and applying the midpoint upwind finite difference discretization the required difference schemes becomes the following equation:(27)1+ht2Lε,muhx,htUi,j+1=1−ht2Lε,muhx,htUi,j+ht2fxi−1/2,tj+fxi−1/2,tj+1,U0,j+1=ϕltj+1,UNx,j+1=ϕrtj+1,fori=1,2,…,Nx−1, where(28)Lε,muhx,htUi,j+1=−εUi+1,j+1−2Ui,j+1+Ui−1,j+1γi2+ai−1/2Ui,j+1−Ui−1,j+1hx+bi−1/2Ui,j+1.and ai−1/2,bi−1/2 and fxi−1/2,tj denoted for axi+axi−1/2, bxi+bxi−1/2 and fxi,tj+fxi−1,tj/2, respectively.


## 3.2.3. Stability and Convergence Analysis on Uniform Mesh

Here, we want to show that the discrete scheme in (27) satisfies the maximum principle, stability estimates, and uniform convergence.Lemma 7.
LetUi,j+1 be any mesh function satisfying U0,j+1≥0,UNx,j+1≥0. Then, 1+ht/2Lε,muhx,htUi,j+1≥0,i=1,2,…,Nx−1 implies that Ui,j+1≥0,i=0,1,…,Nx.Proof.
Suppose there existk∈0,1,…,N such that Uk,j+1=min0≤i≤NUi,j+1. Suppose that Uk,j+1<0 which implies k≠0,N. So, we have that Uk+1,j+1−Uk,j+1>0 and Uk,j+1−Uk−1,j+1<0. So, we obtain 1+ht/2Lε,muhx,htUk,j+1<0 for k=1,2,3,…,Nx−1. Thus, the supposition Ui,j+1<0, for i=0,1,…,Nx is wrong. Hence, we obtain Ui,j+1≥0,i=0,1,…,Nx.
Using the discrete maximum principle, we want to prove the discrete scheme satisfies the uniform stability result.Lemma 8.
The solutionUi,j+1 of the discrete scheme in (27) satisfy the bound.(29)Ui,j+1≤maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1.Proof.
LetQ=maxi1+ht/2Lε,muhx,htUi,j+1/1+ht/2β+maxϕltj+1,ϕrtj+1 and define barrier functions ϑi,j+1±by ϑi,j+1±=Q±Ui,j+1. At the boundary points, we have the following equation:(30)ϑ0,j+1±=Q±U0,j+1=maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕltj+1≥0,ϑN,j+1±=Q±UN,j+1=maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕrtj+1≥0.
On the discretized domainxi,0<i<Nx, we obtain the following equation:(31)1+ht2Lε,muhx,htϑi,j+1±=Q±Ui,j+1−εht2Q±Ui+1,j+1−2Q±Ui,j+1+Q±Ui−1,j+1γi2+ai−1/2ht2Q±Ui,j+1−Q±Ui−1,j+1h+bi−1/2ht2Q±Ui,j+1=1+ht2bi−1/2maxi1+ht/2Lε,muhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±1+ht2Lε,muhx,htUi,j+1≥0,
By the discrete maximum principle in Lemma7, we obtain ϑi,j+1±≥0,xi∈Ω¯Nx. Hence the required bound is satisfied.Lemma 9.
The discrete solutionUi,j+1 of the scheme in (27) satisfy the bound(32)Ui,j+1≤β−1Lε,muhx,htUi,j+1+maxϕltj+1,ϕrtj+1.Proof.
The proof is similar to the above Lemma8 proof by construction a barrier functions ϑi,j+1±=β−1Lε,muhx,htUi,j+1+maxϕltj+1,ϕrtj+1±Ui,j+1.
Let us denote the forward and backward finite differences operators in spatial variable as follows:(33)D+zj+1xi=zj+1xi+1−zj+1xihx,D−zj+1xi=zj+1xi−zj+1xi−1hx.respectively and the second order finite difference operator as follows:(34)D+D−γzj+1xi=zi−1,j+1−2zi,j+1+zi+1,j+1γi2=hx2D+D−zj+1xiγi2.Lemma 10.
For a fixed meshNx and for ε⟶0, it holds(35)limε⟶0max1≤i≤Nx−1ε−mexp−αxiε=0,limε⟶0max1≤i≤Nx−1ε−mexp−α1−xiε=0,m=1,2,3,…where xi=ihx,hx=1/Nx,i=1,2,…,Nx−1.Lemma 11 (see [24]).
For allε∈0,1, the following estimate holds(36)εhx2γi−1≤Chx.Theorem 1.
The truncation error of nonstandard FDM on uniform mesh satisfies the bound(37)Lε,muhx,htuj+1xi−Ui,j+1≤Chx1+maxxiexp−α1−xi/εε3.Proof.
The spatial discretization truncation error is given by the following equation:(38)Lε,muhx,htuj+1xi−Ui,j+1≤εd2uj+1xidx2−Dx+Dx−hx2γi2uj+1xi+ai−1/2ddx−Dx−uj+1xi=Cεd2dx2−Dx+Dx−uj+1xi+hx2γi2−1Dx+Dx−uj+1xi+Chxd2dx2uj+1xi≤Cεd2dx2−Dx+Dx−uj+1xi+Cεhx2γi2−1Dx+Dx−uj+1xi+Chxd2dx2uj+1xi.
Using the result in Lemma11, the truncation error bound becomes the following equation:(39)Lε,muhx,htuj+1xi−Ui,j+1≤Cεhx2d4dx4uj+1xi+Chxd2dx2uj+1xi.
Using the bound in Lemma4 into (39), we obtain the following equation:(40)Lε,muhx,htuj+1xi−Ui,j+1≤Cεhx21+ε−4exp−α1−xicε+Chx1+ε−2exp−α1−xiε≤Chx2ε+ε−3exp−α1−xiε+Chx1+ε−2exp−α1−xiε≤Chx1+maxxiε−3exp−α1−xiε,since,ε−3≥ε−2.Theorem 2.
The nonstandard FDM on uniform mesh satisfies the error bound(41)uj+1xi−Ui,j+1≤Chx.Proof.
Using the results in Lemma10 and applying the bound in Lemma 9, we obtain the error bound.Theorem 3.
Letu and U be the solution of the problem in (1) and (27), respectively, then the nonstandard FDM on uniform mesh satisfies the error bound(42)supε∈0,1u−U≤Chx+ht2.Proof.
Combining the error bound for the temporal and spatial discretization in (16) and Theorem 2 gives the required bound.


## 3.3. Spatial Discretization via Nonstandard FDM on Shishkin Mesh

The development of nonstandard FDM and the uniform convergence of the schemes were restricted on a uniform mesh. In this paper, we extend the scheme and uniform convergence of the scheme on piecewise-uniform Shishkin mesh. Letxii=0Nx be the discretization of the spatial domain Ω¯x where Nx be the number of grid points be an even positive integer. For each i≥1, we define hi=xi−xi−1. Since the considered problem exhibits the right boundary layer, we set a mesh transition parameter τ as follows:(43)τ=min12,σεβlnNx.Which divides the domain [0, 1] in to outer layer regionΩx1=0,1−τ and boundary layer region Ωx2=1−τ,1, where σ is a constant that represents the order of the scheme and β is the lower bound of bx. Here, we set a uniform mesh Nx/2 in Ωx1 with mesh spacing H=21−τ/Nx and similarly a uniform mesh Nx/2 is placed in Ωx2 with the mesh spacing h=2τ/Nx. Hence, the mesh point xi is given as follows:(44)xi=21−τNxi,i=0,1,2,…,Nx2,1−τ+2τNxi−Nx2,i=Nx2+1,…,Nx.The nonstandard FDM are commonly applied on uniform mesh, so using it on layer adapted mesh requires the following change on the denominator function:(45)γ1,i=εaiexphiaiε−1,γ2,i=εaiexphi+1aiε−1.For the mesh functionuj+1x at the grid points xi, we denote the second derivatives term approximation as follows:(46)−εd2dx2uj+1xi≊−2εhi+hi+1Ui+1,j+1−Ui,j+1γ2,i−Ui,j+1−Ui−1,j+1γ1,i=−εDx+Dx−γuj+1xi.Fori=1,2,…,Nx−1, the midpoint upwind nonstandard FDM on Shishkin mesh is given as follows:(47)1+ht2Lε,smNx,htUi,j+1=1−ht2Lε,smNx,htUi,j+ht2fxi−1/2,tj+fxi−1/2,tj+1,U0,j+1=ϕltj+1,UNx,j+1=ϕrtj+1,where(48)Lε,smNx,htUi,j+1=−2εUi+1,j+1−Ui,j+1γ2,ihi+hi+1−Ui,j+1−Ui−1,j+1γ1,ihi+hi+1+ai−1/2Ui,j+1−Ui−1,j+1hi+bi−1/2Ui,j+1.


### 3.3.1. Stability and Convergence Analysis on the Shishkin Mesh

LetH be the mesh width in the outer layer region 0,1−τ and h be the mesh width in layer region 1−τ,1. Then, the mesh widths satisfies the following conditions:(49)Nx−1≤H≤2Nx−1,h=σεβNx−1lnNx≤Nx−1.Lemma 12 (The discrete maximum principle).
LetUi,j+1 be any mesh function satisfying U0,j+1≥0,UNx,j+1≥0. Then, 1+ht/2Lε,smNx,htUi,j+1≥0,i=1,2,…,Nx−1 implies that Ui,j+1≥0,∀i=0,1,…,Nx.Proof.
The proof is similar to the proof of Lemma7.Lemma 13.
The solutionUi,j+1 of the discrete scheme in (47) satisfy the bound.(50)Ui,j+1≤maxi1+ht/2Lε,smNx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1.Proof.
LetQ=maxi1+ht/2Lε,smhx,htUi,j+1/1+ht/2β+maxϕltj+1,ϕrtj+1 and define barrier functions ϑi,j+1± by ϑi,j+1±=Q±Ui,j+1. At the boundary points, we have the following equation:(51)ϑ0,j+1±=Q±U0,j+1=maxi1+ht/2Lε,smNx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕltj+1≥0,ϑN,j+1±=Q±UNx,j+1=maxi1+ht/2Lε,smNx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕrtj+1≥0.
On the discretized spatial domainxi,i=1,2,…,Nx−1, we have the following equation:(52)1+ht2Lε,smNx,htϑi,j+1±=Q±Ui,j+1−2εht2Q±Ui+1,j+1−Q±Ui,j+1γ2,ihi+hi+1−Q±Ui,j+1−Q±Ui−1,j+1γ1,ihi+hi+1+ai−1/2ht2Q±Ui,j+1−Q±Ui−1,j+1hxi+bi−1/2ht2Q±Ui,j+1=1+ht2bi−1/2maxi1+ht/2Lε,smhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±1+ht2Lε,smNx,htUi,j+1≥0,sincebi−1/2≥β,
Using the maximum principle in Lemma12, we obtain ϑi,j+1±≥0,∀xi∈Ω¯Nx. Hence, the required bound is satisfied.Lemma 14.
The discrete solutionUi,j+1 of the scheme in (47) satisfy the bound(53)Ui,j+1≤β−1Lε,smNx,htUi,j+1+maxϕltj+1,ϕrtj+1.Proof.
The proof is similar to the above-given lemma proof by the construction of barrier functionsϑi,j+1±=β−1Lε,smNx,htUi,j+1+maxϕltj+1,ϕrtj+1±Ui,j+1.


### 3.3.2. Decomposition of the Discrete Solution

Here, we decompose the numerical solutionUi,j+1 into regular and singular components using the Shishkin decomposition in [12] as follows:(54)Ui,j+1=Vi,j+1+Wi,j+1,where the regular part satisfied the following nonhomogeneous equation:(55)1+ht2Lε,smNx,htVi,j+1=1−ht2Lε,smNx,htVi,j+ht2fxi−1/2,tj+fxi−1/2,tj+1,Vj+10=vj+10,Vj+11=vj+11,and the singular component satisfies the following homogeneous equation:(56)1+ht2Lε,smNx,htWi,j+1=0,Wj+10=wj+10,Wj+11=wj+11,The error in the computed solution can also be decomposed as follows:(57)Ui,j+1−uj+1xi=Vi,j+1−vj+1xi+Wi,j+1−wj+1xi.First, let us consider the caseτ=1/2 which is a uniform mesh case say the mesh size is h. The truncation error becomes the following equation:(58)Lε,smNx,htuj+1xi−Ui,j+1=εd2dx2−Dx+Dx−γuj+1xi+ai−1/2ddx−Dx−uj+1xi.Since(59)εDx+Dx−γuj+1xi=2εUi+1,j+1−Ui,j+1γ2,ihi+hi+1−Ui,j+1−Ui−1,j+1γ1,ihi+hi+1=εUi+1,j+1−2Ui,j+1+Ui−1,j+1hε/aiexphai/ε−1.From the truncation error in uniform mesh, we obtain the following equation:(60)Lε,smNx,htuj+1xi−Ui,j+1≤Cεh2d4dx4uj+1xi+Chd2dx2uj+1xi≤Ch.Using Lemma14, we obtain uj+1xi−Ui,j+1≤Ch.For the caseτ=σε/βlnNx, we estimate the error in the regular and singular solutions separately. Here, σε/βlnNx≤0.5 implies that 1/ε≤ClnNx.Theorem 4.
The error in the regular component satisfies the bound(61)vj+1xi−Vi,j+1≤Ch.Proof.
On the outer layer region0,1−τ, it is clear that the mesh is uniform i.e. hi=hi+1=H,i=1,2,…,Nx/2 . So the truncation error in the regular component is given as follows:(62)Lε,smNx,htvj+1xi−Vi,j+1=−εd2dx2−Dx+Dx−γvj+1xi+ai−1/2ddx−Dx−vj+1xi≤CεNx−2d4dx4vj+1ξ+CNx−1d2dx2vj+1ξ,since,Nx−1≤H≤2Nx−1≤CεNx−2+CNx−1≤CNx−2lnNx+CNx−1≤CNx−1.
In the boundary layer region1−τ,1, the mesh is uniform i.e., hi=hi+1=h≤Nx−1,i=Nx/2+1,…,Nx−1. So, the truncation error is given as follows:(63)Lε,smNx,htvj+1xi−Vi,j+1=−εd2dx2−Dx+Dx−γvj+1xi+ai−1/2ddx−Dx−vj+1xi≤CεNx−2d4dx4vj+1ξ+CNx−1d2dx2vj+1ξ,since,h≤Nx−1≤CεNx−2+CNx−1≤CNx−2lnNx+CNx−1≤CNx−1.
Using Lemma14, we obtain vj+1xi−Vi,j+1≤CNx−1.Theorem 5.
The error in the singular component satisfies the bound(64)wj+1xi−Wi,j+1≤CNx−1lnNx2.Proof.
In outer layer region0,1−τ, we have hi=hi+1=H,i=1,2,…,Nx/2. So the truncation error becomes(65)Lε,smNx,htwj+1xi−Wi,j+1=−εd2dx2−Dx+Dx−γwj+1xi+ai−1/2ddx−Dx−wj+1xi≤CεNx−2d4dx4wj+1xi+CNx−1d2dx2wj+1xi,sinceNx−1≤H≤2Nx−1≤CεNx−2ε−4exp−α1−xiε+CNx−1ε−2exp−α1−xiε≤CNx−2ε−3+Nx−1ε−2exp−α1−xiε≤CNx−2ε−3+Nx−1ε−2Nx−1,sinceτ=σεβlnNx,exp−ατε≤CNx−1,≤CNx−2lnNx2,since1ε≤ClnNx.
In the layer region1−τ,1, the mesh is also uniform i.e. hi=hi+1=h≤Nx−1,i=Nx/2+1,…,Nx−1. Hence, we have the following equation:(66)Lε,smNx,htwj+1xi−Wi,j+1=−εd2dx2−Dx+Dx−γwj+1xi+ai−1/2ddx−Dx−wj+1xi≤CεNx−2d4dx4wj+1xi+CNx−1d2dx2wj+1xi,sinceh≤Nx−1≤CεNx−2ε−4exp−α1−xiε+CNx−1ε−2exp−α1−xiε≤CNx−2ε−3+Nx−1ε−2exp−α1−xiε,≤CNx−1ε−2,sincemaxiexp−α1−xiε≤1,≤CNx−1lnNx2,since1ε≤ClnNx.
Hence, we haveLε,smNx,htwj+1xi−Wi,j+1≤maxCNx−2lnNx2,CNx−1lnNx2. Using Lemma 14, we obtain wj+1xi−Wi,j+1≤CNx−1lnNx2.Theorem 6.
The nonstandard FDM on the Shishkin mesh satisfies the error bound(67)uj+1xi−Ui,j+1≤CNx−1lnNx2.Proof.
Combining the error estimate in the regular and singular components in Theorem4 and 5 gives the required bound.Theorem 7.
Letu and U be the solution of problem in (1) and (47), respectively, then the nonstandard FDM on the Shishkin mesh satisfies the error bound(68)supε∈0,1u−U≤CNx−1lnNx2+ht2.Proof.
Using the error bound for the temporal and spatial discretization gives the required bound.


## 3.3.1. Stability and Convergence Analysis on the Shishkin Mesh

LetH be the mesh width in the outer layer region 0,1−τ and h be the mesh width in layer region 1−τ,1. Then, the mesh widths satisfies the following conditions:(49)Nx−1≤H≤2Nx−1,h=σεβNx−1lnNx≤Nx−1.Lemma 12 (The discrete maximum principle).
LetUi,j+1 be any mesh function satisfying U0,j+1≥0,UNx,j+1≥0. Then, 1+ht/2Lε,smNx,htUi,j+1≥0,i=1,2,…,Nx−1 implies that Ui,j+1≥0,∀i=0,1,…,Nx.Proof.
The proof is similar to the proof of Lemma7.Lemma 13.
The solutionUi,j+1 of the discrete scheme in (47) satisfy the bound.(50)Ui,j+1≤maxi1+ht/2Lε,smNx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1.Proof.
LetQ=maxi1+ht/2Lε,smhx,htUi,j+1/1+ht/2β+maxϕltj+1,ϕrtj+1 and define barrier functions ϑi,j+1± by ϑi,j+1±=Q±Ui,j+1. At the boundary points, we have the following equation:(51)ϑ0,j+1±=Q±U0,j+1=maxi1+ht/2Lε,smNx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕltj+1≥0,ϑN,j+1±=Q±UNx,j+1=maxi1+ht/2Lε,smNx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±ϕrtj+1≥0.
On the discretized spatial domainxi,i=1,2,…,Nx−1, we have the following equation:(52)1+ht2Lε,smNx,htϑi,j+1±=Q±Ui,j+1−2εht2Q±Ui+1,j+1−Q±Ui,j+1γ2,ihi+hi+1−Q±Ui,j+1−Q±Ui−1,j+1γ1,ihi+hi+1+ai−1/2ht2Q±Ui,j+1−Q±Ui−1,j+1hxi+bi−1/2ht2Q±Ui,j+1=1+ht2bi−1/2maxi1+ht/2Lε,smhx,htUi,j+11+ht/2β+maxϕltj+1,ϕrtj+1±1+ht2Lε,smNx,htUi,j+1≥0,sincebi−1/2≥β,
Using the maximum principle in Lemma12, we obtain ϑi,j+1±≥0,∀xi∈Ω¯Nx. Hence, the required bound is satisfied.Lemma 14.
The discrete solutionUi,j+1 of the scheme in (47) satisfy the bound(53)Ui,j+1≤β−1Lε,smNx,htUi,j+1+maxϕltj+1,ϕrtj+1.Proof.
The proof is similar to the above-given lemma proof by the construction of barrier functionsϑi,j+1±=β−1Lε,smNx,htUi,j+1+maxϕltj+1,ϕrtj+1±Ui,j+1.


## 3.3.2. Decomposition of the Discrete Solution

Here, we decompose the numerical solutionUi,j+1 into regular and singular components using the Shishkin decomposition in [12] as follows:(54)Ui,j+1=Vi,j+1+Wi,j+1,where the regular part satisfied the following nonhomogeneous equation:(55)1+ht2Lε,smNx,htVi,j+1=1−ht2Lε,smNx,htVi,j+ht2fxi−1/2,tj+fxi−1/2,tj+1,Vj+10=vj+10,Vj+11=vj+11,and the singular component satisfies the following homogeneous equation:(56)1+ht2Lε,smNx,htWi,j+1=0,Wj+10=wj+10,Wj+11=wj+11,The error in the computed solution can also be decomposed as follows:(57)Ui,j+1−uj+1xi=Vi,j+1−vj+1xi+Wi,j+1−wj+1xi.First, let us consider the caseτ=1/2 which is a uniform mesh case say the mesh size is h. The truncation error becomes the following equation:(58)Lε,smNx,htuj+1xi−Ui,j+1=εd2dx2−Dx+Dx−γuj+1xi+ai−1/2ddx−Dx−uj+1xi.Since(59)εDx+Dx−γuj+1xi=2εUi+1,j+1−Ui,j+1γ2,ihi+hi+1−Ui,j+1−Ui−1,j+1γ1,ihi+hi+1=εUi+1,j+1−2Ui,j+1+Ui−1,j+1hε/aiexphai/ε−1.From the truncation error in uniform mesh, we obtain the following equation:(60)Lε,smNx,htuj+1xi−Ui,j+1≤Cεh2d4dx4uj+1xi+Chd2dx2uj+1xi≤Ch.Using Lemma14, we obtain uj+1xi−Ui,j+1≤Ch.For the caseτ=σε/βlnNx, we estimate the error in the regular and singular solutions separately. Here, σε/βlnNx≤0.5 implies that 1/ε≤ClnNx.Theorem 4.
The error in the regular component satisfies the bound(61)vj+1xi−Vi,j+1≤Ch.Proof.
On the outer layer region0,1−τ, it is clear that the mesh is uniform i.e. hi=hi+1=H,i=1,2,…,Nx/2 . So the truncation error in the regular component is given as follows:(62)Lε,smNx,htvj+1xi−Vi,j+1=−εd2dx2−Dx+Dx−γvj+1xi+ai−1/2ddx−Dx−vj+1xi≤CεNx−2d4dx4vj+1ξ+CNx−1d2dx2vj+1ξ,since,Nx−1≤H≤2Nx−1≤CεNx−2+CNx−1≤CNx−2lnNx+CNx−1≤CNx−1.
In the boundary layer region1−τ,1, the mesh is uniform i.e., hi=hi+1=h≤Nx−1,i=Nx/2+1,…,Nx−1. So, the truncation error is given as follows:(63)Lε,smNx,htvj+1xi−Vi,j+1=−εd2dx2−Dx+Dx−γvj+1xi+ai−1/2ddx−Dx−vj+1xi≤CεNx−2d4dx4vj+1ξ+CNx−1d2dx2vj+1ξ,since,h≤Nx−1≤CεNx−2+CNx−1≤CNx−2lnNx+CNx−1≤CNx−1.
Using Lemma14, we obtain vj+1xi−Vi,j+1≤CNx−1.Theorem 5.
The error in the singular component satisfies the bound(64)wj+1xi−Wi,j+1≤CNx−1lnNx2.Proof.
In outer layer region0,1−τ, we have hi=hi+1=H,i=1,2,…,Nx/2. So the truncation error becomes(65)Lε,smNx,htwj+1xi−Wi,j+1=−εd2dx2−Dx+Dx−γwj+1xi+ai−1/2ddx−Dx−wj+1xi≤CεNx−2d4dx4wj+1xi+CNx−1d2dx2wj+1xi,sinceNx−1≤H≤2Nx−1≤CεNx−2ε−4exp−α1−xiε+CNx−1ε−2exp−α1−xiε≤CNx−2ε−3+Nx−1ε−2exp−α1−xiε≤CNx−2ε−3+Nx−1ε−2Nx−1,sinceτ=σεβlnNx,exp−ατε≤CNx−1,≤CNx−2lnNx2,since1ε≤ClnNx.
In the layer region1−τ,1, the mesh is also uniform i.e. hi=hi+1=h≤Nx−1,i=Nx/2+1,…,Nx−1. Hence, we have the following equation:(66)Lε,smNx,htwj+1xi−Wi,j+1=−εd2dx2−Dx+Dx−γwj+1xi+ai−1/2ddx−Dx−wj+1xi≤CεNx−2d4dx4wj+1xi+CNx−1d2dx2wj+1xi,sinceh≤Nx−1≤CεNx−2ε−4exp−α1−xiε+CNx−1ε−2exp−α1−xiε≤CNx−2ε−3+Nx−1ε−2exp−α1−xiε,≤CNx−1ε−2,sincemaxiexp−α1−xiε≤1,≤CNx−1lnNx2,since1ε≤ClnNx.
Hence, we haveLε,smNx,htwj+1xi−Wi,j+1≤maxCNx−2lnNx2,CNx−1lnNx2. Using Lemma 14, we obtain wj+1xi−Wi,j+1≤CNx−1lnNx2.Theorem 6.
The nonstandard FDM on the Shishkin mesh satisfies the error bound(67)uj+1xi−Ui,j+1≤CNx−1lnNx2.Proof.
Combining the error estimate in the regular and singular components in Theorem4 and 5 gives the required bound.Theorem 7.
Letu and U be the solution of problem in (1) and (47), respectively, then the nonstandard FDM on the Shishkin mesh satisfies the error bound(68)supε∈0,1u−U≤CNx−1lnNx2+ht2.Proof.
Using the error bound for the temporal and spatial discretization gives the required bound.


## 4. Numerical Examples, Results, and Discussions

Two numerical examples are considered for illustrating the theoretical analysis and findings of the proposed schemes. In the considered examples, the exact solution of the equations is not known.Example 1.
We take this example from [4, 15], were ax=2−x2,bx=x, fx,t=10t2exp−tx1−x with T=1 subject to the initial condition u0x=0 and the boundary conditions ϕlt=0=ϕrt,t∈0,T.Example 2.
We take this example from [16], were ax=2−x2,bx=1+x2+cosπx,fx,t=10t2exp−tx1−x with T=1 subject to the initial condition u0x=0 and the boundary conditions.
Since the exact solutions are not given, we use the double mesh technique to calculate the maximum absolute error. LetUi,jNx,Nt be the computed solution of the equation on Nx and Nt number of mesh points and Ui,j2Nx,2Nt be the computed solution on the double number of mesh points 2Nx and 2Nt. The maximum absolute error is calculated as follows:(69)EεNx,Nt=maxi,jUi,jNx,Nt−Ui,j2Nx,2Nt.
The uniform error is given by the following equation:(70)ENx,Nt=max∀ε∈0,1EεNx,Nt.
The rate of convergence is calculated as follows:(71)rateεNx,Nt=logEεNx,Nt−logEε2Nx,2Ntlog2,
and the uniform rate of convergence is calculated as follows:(72)rateNx,Nt=logENx,Nt−logE2Nx,2Ntlog2.
As one observes in Figures1 and 2, the solution of the considered Example 1 and 2 exhibits an exponential boundary layer on the right side of the spatial domain when ε goes small. In Figures 3(a) and 4(a), the computed solution of Example 1 and 2 using the midpoint upwind nonstandard FDM on uniform mesh for ε=10−6 and T=1 are given. In these figures, we observe that there is no mesh point in the boundary layer region. In Figures 3(b) and 4(b), the computed solution using the midpoint upwind nonstandard FDM on Shishkin mesh is depicted for ε=10−6. In these figures, we observe that there is a sufficient number of mesh points in the boundary layer region. Such property is called the boundary layer resolving property of the scheme.
In Tables1 and 2, the maximum absolute error, ε-uniform error and the ε-uniform rate of convergence of the schemes are given on uniform and Shishkin mesh, for different values of perturbation parameter ranging from ε=10−4 to ε=10−10. As it is observed in the tables, the schemes gives accurate and stable numerical results for different values of ε. We observe that the scheme on uniform mesh have better accuracy than the result on Shishkin mesh. In Table 3, a comparison of the uniform error of the proposed schemes with the results in [4, 15] are given. One observes that the proposed schemes gives more accurate result than the methods given in the literature. In Table 4, the comparison of the maximum absolute error of the proposed midpoint upwind nonstandard FDM on Shishkin mesh with the result in [16] is given for the boundary layer region. As it is observed in this table, the proposed scheme gives a more accurate result in the boundary layer region.Figure 1
The computed solution of Example1 for different values of ε, using the Shishkin mesh (a) ε=10−1 and (b) ε=10−6.
(a)(b)Figure 2
The computed solution of Example2 for different values of ε, using the Shishkin mesh (a) ε=10−1 and (b) ε=10−6.
(a)(b)Figure 3
Example1, boundary layer resolving property: (a) Uniform mesh and (b) the Shishkin mesh for N=64, ε=10−6 and T=1.
(a)(b)Figure 4
Example2, boundary layer resolving property: (a) Uniform mesh and (b) Shishkin mesh for N = 64, ε=10−6 and T=1.
(a)(b)Table 1
Example1, the maximum absolute error, ε-uniform error ENx,Nt and ε-uniform rate of convergence rNx,Nt.
εNx=2526272829↓ht=1/101/201/401/801/160On uniform mesh10−42.8825e −031.0224e −034.8019e −042.3423e −041.1576e −0410−52.8825e −031.0224e −034.8019e −042.3423e −041.1576e −0410−62.8825e −031.0224e −034.8019e −042.3423e −041.1576e −0410−82.8825e −031.0224e −034.8019e −042.3423e −041.1576e −0410−102.8825e −031.0224e −034.8019e −042.3423e −041.1576e −04ENx,Nt2.8825e −031.0224e −034.8019e −042.3423e −041.1576e −04rNx,Nt1.49541.09031.03571.0168—On Shishkin mesh10−43.3989e −031.7021e −038.5853e −044.3194e −042.1685e −0410−53.4008e −031.7033e −038.5927e −044.3237e −042.1708e −0410−63.4009e −031.7034e −038.5934e −044.3242e −042.1711e −0410−83.4010e −031.7034e −038.5935e −044.3242e −042.1711e −0410−103.4010e −031.7034e −038.5935e −044.3242e −042.1711e −04ENx,Nt3.4010e −031.7034e −038.5935e −044.3242e −042.1711e −04rNx,Nt0.99750.98710.99080.9940—Table 2
Example2, the maximum absolute error, ε-uniform error ENx,Nt and ε-uniform rate of convergence rNx,Nt.
εNx=2526272829↓ht=1/101/201/401/801/160On uniform mesh10−41.8404e −039.2256e −044.6117e −042.3048e −041.1521e −0410−51.8404e −039.2256e −044.6117e −042.3048e −041.1521e −0410−61.8404e −039.2256e −044.6117e −042.3048e −041.1521e −0410−81.8404e −039.2256e −044.6117e −042.3048e −041.1521e −0410−101.8404e −039.2256e −044.6117e −042.3048e −041.1521e −04ENx,Nt1.8404e −039.2256e −044.6117e −042.3048e −041.1521e −04rNx,Nt0.99631.00031.00071.0004—On Shishkin mesh10−42.5298e −031.2933e −036.6114e −043.3536e −041.6903e −0410−52.5302e −031.2939e −036.6166e −043.3570e −041.6924e −0410−62.5303e −031.2940e −036.6171e −043.3574e −041.6926e −0410−82.5303e −031.2940e −036.6172e −043.3574e −041.6926e −0410−102.5303e −031.2940e −036.6172e −043.3574e −041.6926e −04ENx,Nt2.5303e −031.2940e −036.6172e −043.3574e −041.6926e −04rNx,Nt0.96750.96750.97890.9881—Table 3
Example1, the comparison of ε-uniform error of the proposed schemes and results in literature.
SchemesNx=2526272829↓ht=1/101/201/401/801/160Scheme on uniform meshEN,M2.8825e −031.0224e −034.8019e −042.3423e −041.1576e −04Scheme on Shishkin meshEN,M3.4010e −031.7034e −038.5935e −044.3242e −042.1711e −04Method in [4]EN,M6.1336e −033.5748e −031.9245e −039.9750e −045.0763e −04Method in [15]EN,M1.4524e −026.7667e −033.6247e −031.7939e −038.9428e −04Table 4
Example2, the comparison of the maximum absolute error of the proposed scheme on the Shishkin mesh and the result in [16] for boundary layer region.
Nxε=10−6,Nx/2≤i≤Nxε=10−12,Nx/2≤i≤NxResult in [16]Proposed schemeResult in [16]Proposed scheme82.733708e −022.2928e −032.733715e −022.2928e −03161.420798e −021.6804e −031.420822e −021.6804e −03329.267409e −031.5971e −039.267692e −031.5969e −03645.863876e −031.0624e −035.864240e −031.0626e −031283.497379e −036.0492e −043.497815e −036.0503e −042561.969163e −033.2172e −041.969650e −033.2173e −04


## 5. Conclusion

Numerical schemes are developed using the nonstandard FDM on uniform and piecewise uniform Shishkin mesh for solving nonstationary or time dependent singularly perturbed convection-diffusion-reaction equation. The solution of the considered equation exhibit a boundary layer of widthOε. The developed numerical schemes uses the Crank–Nicolson method in temporal discretization and the midpoint upwind nonstandard FDM on uniform and Shishkin mesh for the spatial discretization. The stability of the schemes are investigated and discussed. The uniform convergence of the schemes are proved. The applicability of the schemes are investigated by considering test examples. The developed schemes are shown to be uniformly convergent with linear order of convergence. The proposed schemes gives accurate, stable, and uniformly convergent numerical results. Furthermore, the scheme on piecewise uniform Shishkin mesh gives boundary layer resolving result. The proposed schemes can be extended for solving nonlinear and higher dimensional singularly perturbed equations.



---

*Source: 2043323-2022-12-05.xml*