# Parameter-Uniform Numerical Scheme for Singularly Perturbed Delay Parabolic Reaction Diffusion Equations with Integral Boundary Condition

Authors: Wakjira Tolassa Gobena; Gemechis File Duressa
Journal: International Journal of Differential Equations
 (2021)

---

## Abstract

Numerical computation for the class of singularly perturbed delay parabolic reaction diffusion equations with integral boundary condition has been considered. A parameter-uniform numerical method is constructed via the nonstandard finite difference method for the spatial direction, and the backward Euler method for the resulting system of initial value problems in temporal direction is used. The integral boundary condition is treated using numerical integration techniques. Maximum absolute errors and the rate of convergence for different values of perturbation parameterε and mesh sizes are tabulated for two model examples. The proposed method is shown to be parameter-uniformly convergent.


---

## Body

## 1. Introduction

Singularly perturbed delay differential equations model physical problems for which the evaluation depends on the present state of the system and also on the past history. In general, when the perturbation parameter approaches zero, the smoothness of the solution of the singularly perturbed delay differential equations deteriorates and it forms a boundary layer(s) [1]. The time-dependent problem, i.e., singularly perturbed parabolic delay differential equations, has applications in the study of modeling of physiological kinetics, production of blood cell, and neuronal variability (see [2–4] and the references therein). But, for more applications, the theory and numerical treatment of delay differential equations can be found in [5].It is known that standard numerical methods for solving singular perturbation problems are unstable and fail to give accurate results when the perturbation parameter is small. The accuracy and convergence of the methods need attention, because the treatment of singular perturbation problems is not trivial, and the solution depends on perturbation parameter and mesh sizeh [6]. This suggests that the numerical treatment should be improved. The presence of the singular perturbation parameter leads to occurrences of oscillations or divergence in the computed solutions while using standard numerical methods [7]. To overcome these oscillations or divergence, a large number of mesh points are required as ε goes to zero. This is difficult and sometimes impossible to handle such cases. Therefore, it is necessary to develop suitable numerical methods which are uniformly convergent to solve this type of differential equations.Many researchers have proposed different numerical methods to solve singularly perturbed time delay parabolic differential equations; for instance, in [8], Ansari et al. proposed a parameter robust finite difference method for singularly perturbed delay parabolic partial differential equations of reaction diffusion type with Dirichlet boundary conditions on Shishkin mesh. They made use of the results of Miller et al. [9]. In [10], Bashier and Patidar improved the order of accuracy for the method suggested in [8]. Avudai Selvi and Ramanujam [11] studied singularly perturbed time delay parabolic reaction diffusion equation subject to mixed-type boundary conditions. A second-order fitted operator finite difference method for reaction diffusion type equation was suggested by Bashier and Patidar [12]. Singh et al. [13] developed an overlapping Schwarz domain decomposition method for parabolic reaction diffusion problems with time delay. An extended cubic B-spline collocation method on Shishkin mesh has been developed for singularly perturbed delay parabolic reaction diffusion problem with mixed-type boundary conditions in [14]. A numerical method is developed for solving a singularly perturbed parabolic delay differential equation whose solution exhibits a boundary layer behaviour in [15] and singularly perturbed parabolic differential difference equations arising in computational neuroscience in [16]. The authors in [17] gave a parameter robust numerical scheme to the solution of second-order singularly perturbed parabolic partial differential equations with large delay by using the finite difference scheme on the piecewise-uniform Shishkin mesh. In [18], Kumar and Kumari constructed a parameter-uniform implicit scheme for a class of parabolic singularly perturbed reaction diffusion initial-boundary value problems with large delay in the spatial direction, where the solution exhibits twin boundary layers and an interior layer on a piecewise-uniform mesh.But, in recent years, there has been a growing interest in numerical methods on singularly perturbed delay differential equations with integral boundary conditions. In [19, 20], the authors consider a first-order singularly perturbed quasilinear boundary value problem with integral boundary conditions on Shishkin mesh. The authors in [21, 22] studied singularly perturbed boundary value problem for second-order ordinary differential equations with nonlocal boundary condition on uniform mesh. In [23], Sekar and Tamilselvan have discussed singularly perturbed delay differential equations of convection-diffusion type with integral boundary condition, using the finite difference method on Shishkin mesh. The method has almost first-order convergence with respect to ε. Debela and Duressa [24] proposed first-order numerical methods using the exponentially fitted finite difference method on a uniform mesh. They also proposed the accelerated fitted finite difference method for singularly perturbed delay differential equations with nonlocal boundary condition in [25]. In [26], a uniformly convergent numerical method is constructed and a numerical integration method is used to treat the nonlocal boundary condition. However, up to the best of our knowledge, except the work in [27], not much work has been done to solve the problem under consideration.In the present paper, motivated by the studies of [27], we construct and analyze an accurate and uniformly convergent numerical scheme for solving singularly perturbed delay parabolic reaction diffusion equations with integral boundary condition. The proposed scheme uses the procedures of method of line (i.e., first discretizing in spatial direction followed by discretization in temporal direction) using the backward Euler method in temporal direction and nonstandard fitted operator FDM on spatial direction.The content of this paper is as follows: In Section2, definition of the problem and the properties of continuous solution are given. In Section 3, discretizing the spatial domain and techniques of nonstandard finite difference are discussed, and the ε−uniform convergence of the discrete problem is proved. Next, the backward Euler method used for the system of initial value problems resulted from spatial discretization and the convergence of the discrete scheme are discussed. In Section 4, numerical results and discussion are given to validate the theoretical analysis, and finally, in Section 5, the conclusion of the work done is presented.Notation 1.
The notation for jump in any function is at a point ‘p’ with χp=χp+−χp− and C is a positive constant independent of ε and N. The norm for a given function χ defined on the domain Ω is calculated by χΩ=supx,t∈Ωχx,t.


## 2. Statement of the Problem

Let us consider the following singularly perturbed delay parabolic differential equation with integral boundary condition onΩ=D×0,T≡0,2×0,T, where T is some fixed positive time and Γ=Γl∪Γb∪Γr, where Γl and Γr are the left and the right sides of the rectangular domain Ω corresponding to x=0 and x=2, respectively, and Γb is the base of the domain and given by Γb=D¯.(1)Lux,t≡∂ux,t∂t−ε∂2ux,t∂x2+ax,tux,t=−bx,tux−1,t+fx,t,forx,t∈Ω.Subject to initial conditions(2)ux,t=ϕbx,t,x,t∈Γb=x,0;0≤x≤2and boundary conditions(3)ux,t=ϕlx,t,x,t∈Γl=x,t;−1≤x≤0and0≤t≤T,κux,t=u2,t−ε∫02gxux,tdx=ϕrx,t,x,t∈Γr=2,t;0≤t≤T,where x,t∈Ω,Ω¯=D¯×0,T and ε∈0,1 is a given constant. ax,t,bx,t,fx,tonΩ¯ and ϕlx,t,ϕrx,t,ϕbx,tonΓ are sufficiently smooth, bounded functions that satisfy(4)ax,t≥α>0bx,t≤β<0α+β>0,onΩ¯.Furthermore,gx is a nonnegative, monotonic function and satisfies ∫02gxdx<1.Problems (1)–(3) can be rewritten as(5)Liux,t=Gix,t,i=1,2,where(6)L1ux,t=∂ux,t∂t−ε∂2ux,t∂x2+ax,tux,t=−bx,tϕlx−1,t+fx,t,x,t∈Ω1,(7)L2ux,t=∂ux,t∂t−ε∂2ux,t∂x2+ax,tux,t+bx,tux−1,t=fx,t,x,t∈Ω2,with boundary conditions(8)ux,t=ϕlx,t,x,t∈Γl=x,t;−1≤x≤0and0≤t≤T,u1−,t=u1+,t,∂u∂z1−,t=∂u∂z1+,t,κux,t≡u2,t−ε∫02gxux,tdx=ϕrx,t,x,t∈Γr=2,t;0≤t≤T,ux,t=ϕbx,t,x,t∈Γb=x,0;0≤x≤2,where Ω1=0,1×0,T,Ω2=1,2×0,T,andΩ∗=Ω1∪Ω2.


### 2.1. Properties of Continuous Problems

In the study of the numerical aspects of singularly perturbed problems, their analytical aspects play an important role. Setting the valueε=0, the reduced problem corresponding to (6)–(8) is(9)∂u0x,t∂t+ax,tu0x,t=−bx,tϕlx−1,t+fx,t,x,t∈Ω1,u0x,t=ϕbx,t,x,t∈Γb,∂u0x,t∂t+ax,tu0x,t+bx,tu0x−1,t=fx,t,x,t∈Ω2,u0x,t=ϕbx,t,x,t∈Γb.Asu0x,t need not satisfy u00,t=u0,t and u02,t=u2,t, the solution ux,t exhibits boundary layers at x=0andx=2. Furthermore, as u01−,t need not be equal to u01+,t, the solution ux,t exhibits interior layers at x=1. Problems (1)–(3) satisfy the following lemma.Lemma 1 (continuous maximum principle).
Letψx,t∈U∗=C0,0Ω¯∩C1,0Ω∩C2,1Ω1∪Ω2 such that ψ0,t≥0,ψx,0≥0, κψ2,t≥0,L1ψx,t≥0,∀x,t∈Ω1L2ψx,t≥0,∀x,t∈Ω2, and ψx1,t=ψx1+,t−ψx1−,t≤0. Then, ψx,t≥0,∀x,t∈Ω¯.Proof.
Define a test function(10)sx,t=18+x2,x,t∈0,1×0,T,38+x4,x,t∈1,2×0,T.
Note thatsx,t>0,∀x,t∈Ω¯,Lsx,t>0,∀x,t∈Ω1∪Ω2,s0,t>0,sx,0>0κs2,t>0, and sx1,t<0. Let(11)δ1=max−ψx,tsx,t:x,t∈D¯.
Then, there existsx0,t0∈Ω¯ such that ψx0,t0+δ1sx0,t0=0 and ψx,t+δ1sx,t≥0, ∀x,t∈Ω¯. Therefore, the function attains its minimum at x,t=x0,t0. Suppose the theorem does not hold true, then δ1>0:Case (i):x0,t0=0,t00<ψ+δ1s0,t0=ψ0,t0+δ1s0,t0=0. It is a contradiction.Case (ii):x0,t0∈Ω10<L1ψ+δ1sx0,t0=−εψ+δ1sxxx0,t0+ψ+δ1stx0,t0+ax0,t0ψ+δ1sx0,t0≤0. It is a contradiction.Case (iii):x0,t0=1,t00≤ψ+δ1s′1,t0=ψ′1,t0+δ1s′1,t0<0. It is a contradiction.Case (iv):x0,t0∈Ω20<L2ψ+δ1sx0,t0=−εψ+δ1sxxx0,t0+ψ+δ1stx0,t0+ax0,t0ψ+δ1sx0,t0+bx0,t0ψ+δ1sx0−1,t0≤0. It is a contradiction.Case (v):x0,t0=2,t00<κψ+δ1s2,t0=ψ+δ1s2,t0−ε∫02gxψ+δ1sx,tdx≤0. It is a contradiction.
Hence, the proof of the theorem is verified.Lemma 2 (stability result).
The solutionux,t of problems (6)–(8) satisfies the bound(12)uΩ¯≤CmaxuΓl,uΓb,κuΓr,LuΩ∗,x,t∈Ω¯.Proof.
It can be easily proved using the maximum principle (Lemma1) and the barrier functions(13)Θ±x,t=CMsx,t±ux,t,x,t∈Ω¯,where M=maxuΓl,uΓb,κuΓr,LuΩ∗ and sx,t are the test functions as in Lemma 1.
The existence and uniqueness of a solution for problems (6)–(8) can be established by assuming that the data are Holder¨ continuous and imposing appropriate compatibility conditions at the corner points (0, 0), (2, 0), (−1, 0), and (1, 0) (see [28]). Then, the required compatibility conditions are(14)ϕb0,0=ϕl0,0,ϕb2,0=ϕr2,0,(15)∂ϕl0,0∂t−ε∂2ϕb0,0∂x2+a0,0ϕb0,0=−b0,0ϕl−1,0+f0,0,∂ϕr2,0∂t−ε∂2ϕb2,0∂x2+a2,0ϕb2,0=−b2,0ϕb1,0+f2,0.
Note that,ϕlx,t,ϕrx,t,andϕbx,t are assumed to be smooth for (15) to make sense that ϕlx,t∈C10,T,ϕrx,t∈C10,Tandϕbx,t∈C2,1Γb.
The following theorem gives sufficient conditions for the existence of a unique solution of problems (6)–(8) (see [28]).Theorem 1.
Leta,b,f∈Cα1,α1/2Ω¯,ϕl,ϕr∈C1+α1/20,T,ϕb∈C2+α1,1+α1/2Γb,andα1∈0,1. Assume that compatibility conditions (14) and (15) are fulfilled. Then, problems (6)–(8) have a unique solution ux,t and u∈C2+α1,1+α1/2Ω¯.Proof.
For the proof, refer to [28].Theorem 2.
Leta,b,f∈C2+α1,1+α1/2Ω¯,ϕl,ϕr∈C2+α1/20,T,ϕb∈C4+α1,2+α1/2Γb,andα1∈0,1. Assume that compatibility conditions (14) and (15) are fulfilled. Then, problems (6)–(8) have a unique solution ux,t and u∈C4+α1,2+α1/2Ω¯. Furthermore, the derivatives of the solution ux,t satisfy(16)∂i+ju∂xi∂tjΩ¯≤Cε−i/2,∀i,j∈ℤ≥0such that0≤i+2j≤4,where the constant C is independent of ε.Proof.
For the proof, refer to [27].
The nonclassical bounds in singular and regular components and their derivatives are established in the following theorem since the bounds in Theorem2 are not adequate for the proof of parameter-uniform error estimate.Theorem 3.
Let the dataa,b,f∈C4+α1,2+α1/2Ω¯,ϕl,ϕr∈C3+α1/20,Tϕb∈C6+α1,3+α1/2Γb, and α1∈0,1. Assume that compatibility conditions (14) and (15) are satisfied. Then, we have(17)∂i+jv∂xi∂tjΩ¯≤C1+ε1−i/2,(18)∂i+jwlx,t∂xi∂tj≤Cε−i/2exp−xε,x,t∈Ω1,Cε−i/2exp−x−1ε,x,t∈Ω2,(19)∂i+jwrx,t∂xi∂tj≤Cε−i/2exp−1−xε,x,t∈Ω1,Cε−i/2exp−2−xε,x,t∈Ω2,where C is the constant independent of ε,x,t∈Ω¯,i,j≥0,0≤i+2j≤4.Proof.
One may refer [27] for the details.Theorem 4.
The partial derivative ofwx,t satisfies(20)∂i+jwx,t∂xi∂tj≤Cε−i/2exp−xε+exp−1−xε,x,t∈Ω1,Cε−i/2exp−x−1ε+exp−2−xε,x,t∈Ω2,∀i,j∈Z≥0 such that i+2j∈0,4.Proof.
The proof of the theorem is completed by using the estimates of (18) and (19) and the decomposition w=wl+wr.


## 2.1. Properties of Continuous Problems

In the study of the numerical aspects of singularly perturbed problems, their analytical aspects play an important role. Setting the valueε=0, the reduced problem corresponding to (6)–(8) is(9)∂u0x,t∂t+ax,tu0x,t=−bx,tϕlx−1,t+fx,t,x,t∈Ω1,u0x,t=ϕbx,t,x,t∈Γb,∂u0x,t∂t+ax,tu0x,t+bx,tu0x−1,t=fx,t,x,t∈Ω2,u0x,t=ϕbx,t,x,t∈Γb.Asu0x,t need not satisfy u00,t=u0,t and u02,t=u2,t, the solution ux,t exhibits boundary layers at x=0andx=2. Furthermore, as u01−,t need not be equal to u01+,t, the solution ux,t exhibits interior layers at x=1. Problems (1)–(3) satisfy the following lemma.Lemma 1 (continuous maximum principle).
Letψx,t∈U∗=C0,0Ω¯∩C1,0Ω∩C2,1Ω1∪Ω2 such that ψ0,t≥0,ψx,0≥0, κψ2,t≥0,L1ψx,t≥0,∀x,t∈Ω1L2ψx,t≥0,∀x,t∈Ω2, and ψx1,t=ψx1+,t−ψx1−,t≤0. Then, ψx,t≥0,∀x,t∈Ω¯.Proof.
Define a test function(10)sx,t=18+x2,x,t∈0,1×0,T,38+x4,x,t∈1,2×0,T.
Note thatsx,t>0,∀x,t∈Ω¯,Lsx,t>0,∀x,t∈Ω1∪Ω2,s0,t>0,sx,0>0κs2,t>0, and sx1,t<0. Let(11)δ1=max−ψx,tsx,t:x,t∈D¯.
Then, there existsx0,t0∈Ω¯ such that ψx0,t0+δ1sx0,t0=0 and ψx,t+δ1sx,t≥0, ∀x,t∈Ω¯. Therefore, the function attains its minimum at x,t=x0,t0. Suppose the theorem does not hold true, then δ1>0:Case (i):x0,t0=0,t00<ψ+δ1s0,t0=ψ0,t0+δ1s0,t0=0. It is a contradiction.Case (ii):x0,t0∈Ω10<L1ψ+δ1sx0,t0=−εψ+δ1sxxx0,t0+ψ+δ1stx0,t0+ax0,t0ψ+δ1sx0,t0≤0. It is a contradiction.Case (iii):x0,t0=1,t00≤ψ+δ1s′1,t0=ψ′1,t0+δ1s′1,t0<0. It is a contradiction.Case (iv):x0,t0∈Ω20<L2ψ+δ1sx0,t0=−εψ+δ1sxxx0,t0+ψ+δ1stx0,t0+ax0,t0ψ+δ1sx0,t0+bx0,t0ψ+δ1sx0−1,t0≤0. It is a contradiction.Case (v):x0,t0=2,t00<κψ+δ1s2,t0=ψ+δ1s2,t0−ε∫02gxψ+δ1sx,tdx≤0. It is a contradiction.
Hence, the proof of the theorem is verified.Lemma 2 (stability result).
The solutionux,t of problems (6)–(8) satisfies the bound(12)uΩ¯≤CmaxuΓl,uΓb,κuΓr,LuΩ∗,x,t∈Ω¯.Proof.
It can be easily proved using the maximum principle (Lemma1) and the barrier functions(13)Θ±x,t=CMsx,t±ux,t,x,t∈Ω¯,where M=maxuΓl,uΓb,κuΓr,LuΩ∗ and sx,t are the test functions as in Lemma 1.
The existence and uniqueness of a solution for problems (6)–(8) can be established by assuming that the data are Holder¨ continuous and imposing appropriate compatibility conditions at the corner points (0, 0), (2, 0), (−1, 0), and (1, 0) (see [28]). Then, the required compatibility conditions are(14)ϕb0,0=ϕl0,0,ϕb2,0=ϕr2,0,(15)∂ϕl0,0∂t−ε∂2ϕb0,0∂x2+a0,0ϕb0,0=−b0,0ϕl−1,0+f0,0,∂ϕr2,0∂t−ε∂2ϕb2,0∂x2+a2,0ϕb2,0=−b2,0ϕb1,0+f2,0.
Note that,ϕlx,t,ϕrx,t,andϕbx,t are assumed to be smooth for (15) to make sense that ϕlx,t∈C10,T,ϕrx,t∈C10,Tandϕbx,t∈C2,1Γb.
The following theorem gives sufficient conditions for the existence of a unique solution of problems (6)–(8) (see [28]).Theorem 1.
Leta,b,f∈Cα1,α1/2Ω¯,ϕl,ϕr∈C1+α1/20,T,ϕb∈C2+α1,1+α1/2Γb,andα1∈0,1. Assume that compatibility conditions (14) and (15) are fulfilled. Then, problems (6)–(8) have a unique solution ux,t and u∈C2+α1,1+α1/2Ω¯.Proof.
For the proof, refer to [28].Theorem 2.
Leta,b,f∈C2+α1,1+α1/2Ω¯,ϕl,ϕr∈C2+α1/20,T,ϕb∈C4+α1,2+α1/2Γb,andα1∈0,1. Assume that compatibility conditions (14) and (15) are fulfilled. Then, problems (6)–(8) have a unique solution ux,t and u∈C4+α1,2+α1/2Ω¯. Furthermore, the derivatives of the solution ux,t satisfy(16)∂i+ju∂xi∂tjΩ¯≤Cε−i/2,∀i,j∈ℤ≥0such that0≤i+2j≤4,where the constant C is independent of ε.Proof.
For the proof, refer to [27].
The nonclassical bounds in singular and regular components and their derivatives are established in the following theorem since the bounds in Theorem2 are not adequate for the proof of parameter-uniform error estimate.Theorem 3.
Let the dataa,b,f∈C4+α1,2+α1/2Ω¯,ϕl,ϕr∈C3+α1/20,Tϕb∈C6+α1,3+α1/2Γb, and α1∈0,1. Assume that compatibility conditions (14) and (15) are satisfied. Then, we have(17)∂i+jv∂xi∂tjΩ¯≤C1+ε1−i/2,(18)∂i+jwlx,t∂xi∂tj≤Cε−i/2exp−xε,x,t∈Ω1,Cε−i/2exp−x−1ε,x,t∈Ω2,(19)∂i+jwrx,t∂xi∂tj≤Cε−i/2exp−1−xε,x,t∈Ω1,Cε−i/2exp−2−xε,x,t∈Ω2,where C is the constant independent of ε,x,t∈Ω¯,i,j≥0,0≤i+2j≤4.Proof.
One may refer [27] for the details.Theorem 4.
The partial derivative ofwx,t satisfies(20)∂i+jwx,t∂xi∂tj≤Cε−i/2exp−xε+exp−1−xε,x,t∈Ω1,Cε−i/2exp−x−1ε+exp−2−xε,x,t∈Ω2,∀i,j∈Z≥0 such that i+2j∈0,4.Proof.
The proof of the theorem is completed by using the estimates of (18) and (19) and the decomposition w=wl+wr.


## 3. Description of the Numerical Method

In this section, the spatial and temporal semidiscretization was described and also their convergence was analyzed separately.


### 3.1. Spatial Semidiscretization

The theoretical basis of the nonstandard discrete modeling method is based on the development of the exact finite difference method. In [29], Mickens presented techniques and rules for developing nonstandard FDMs for different problem types. In Mickens rules, to develop a discrete scheme, a denominator function for the discrete derivatives must be expressed in terms of more complicated functions of step sizes than those used in the standard procedures. This complicated function constitutes a general property of the schemes, which is useful while designing reliable schemes for such problems. In line to this, we need to derive a nonstandard finite difference method which captures the layer behaviour of the problem on a uniform mesh. Mickens in [29, 30] gave a novel approach of nonstandard finite difference method, and the basic idea of this method is to replace the denominator function h2 of the second-order derivative with a suitable function in the differential equation, which comes under the category of the fitted operator finite difference method. On the spatial domain [0, 2], we introduce the equidistant meshes with uniform mesh length such that(21)ΩxN=x0=0,xi=ih,i=11N−1,xN/2=1,xN=2,h=2N,where h is the step size and N is the number of mesh points in the spatial direction.Moreover, we consider the domain0,1 which is discretized into N/2 equal number of subintervals, each of length h. For the problem of the form in (6)–(8), in order to construct the exact finite difference scheme, we follow the procedures used in [31]. Consider the constant coefficient subequations given in (22) by ignoring the time variable as(22)−εd2uxdx2+αux=0,where ax,t≥α>0. Thus, problem (22) has two independent solutions, namely, expσ1xandexpσ2x, with(23)σ1,2=±αε.We denote the approximate solution ofux at xi′s by Ui. Now, our objective is to calculate a difference equation which has the same general solution as problem (22) has at the grid point xi given by Ui=Aexpσ1xi+Bexpσ2xi. Using the procedures used in [31], we have(24)detUi−1expσ1xi−1expσ2xi−1Uiexpσ1xiexpσ2xiUi+1expσ1xi+1expσ2xi+1=0.Simplifying (24), we obtain that(25)Ui−1−2coshαεhUi+Ui+1=0is an exact difference scheme for (22).After doing the arithmetic manipulation and rearrangement on (25), we obtain(26)−εUi−1−2Ui+Ui+14ε/αsinh2α/εh/2+αUi=0.The denominator function becomesη2=4ε/αsinh2α/εh/2. Adopting this denominator function for the variable coefficient problem, we write it as(27)ηi2ε,h,t=4εaitsinh2aitεh2.Case 1: consider (6)–(8) on the domain Ω=0,2×0,T,fori=1,…,N−1.LetUit denote the approximation of uxi,t. By using the nonstandard finite difference approximation, at this stage, the problem in (6)–(8) reduces to the semidiscrete form as(28)LlNUit=Git,fori=1,2,…,N−1,l=1,2,where(29)L1NUit≡dUitdt−εUi+1t−2Uit+Ui−1tηi2+aitUit=−bitϕlxi−N/2,t+fit,xi∈Ω1N,(30)L2NUit≡dUitdt−εUi+1t−2Uit+Ui−1tηi2+aitUit+bitUxi−N/2,t=fix,t,xi∈Ω2N.Case 2: since the problem involves nonlocal boundary conditions, we considered the following to obtain an equation at the end boundary condition.Fori=N, the composite Simpson’s integration rule approximates the nonlocal condition in (1) as follows:(31)κux,t≡u2,t−ε∫02gxux,tdx=UNt−εh∑i=0NcigiUit=ϕrxi,t,0<t≤T,where c0=cN=1/3andc2k−1=4/3,fork=1,2,…,N/2,andc2k=2/3,fork=1,2,…,N/2−1,Now, using (31) and left boundary condition, we can obtain the following formula for UNt:(32)UNt=11−εhcNgNεhc0g0ϕl0,t+εh∑i=1N−1cigiUit+ϕrxi,t.Therefore, on the whole domainD¯=0,2, the basic schemes to solve (6)–(8) are the schemes given in (29), (30), and (32).The time domain is continuous at this stage, and the system of IVPs in (29), (30), and (32) can be written in the compact form as(33)dUitdt+EitUit=Git,i=1,2,…,N,where Eit is a matrix of size NxN and UitandGit are vectors of size N. The entries of the coefficient matrix Eit are, respectively, given by(34)Eii−1t=−εηi2,fori=2,3,…,N−1,Eiit=2εηi2+ait,fori=1,2,…,N−1,Eii+1t=−εηi2,fori=1,2,…,N−1,Eii−N/2t=bit,fori=N2+1,…,N−1,EN,N=−4εh3∑i=1Ng2i−1U2i−1t−2εh3∑i=1N−1g2iU2it+1−εh3gN,i=N,and the corresponding right-hand-side vector in the equation system has the entries(35)G1t=f1t−b1tϕlx1−N/2,t+εη12ϕl0,t,Git=fit−bitϕlxi−N/2,t,fori=2,3,…,N2,fit,fori=N2+1,…,N−1,GNt=ϕrxi,t+εh3g0U0t,fori=N,U=U1,U2,…,UNT.Before we proceed with the convergence analysis, we highlight some properties of the discrete problem to original (1) in the form of lemma which plays a major role in the said analysis.Fori=1,2,…,N/2,(36)L1NUit=fit−bitϕlxi−N/2,t.Fori=N/2+1,…,N−1,(37)L2NUit=fit.Subject to the boundary conditions,(38)Uit=ϕlxi,t,i=−N2−N2+1,…,0andt∈0,T,κNUxN,t=UxN,t−ε∑i=1Ngi−1Ui−1t+4giUit+gi+1Ui+1t3hi=ϕrxi,t,Dx−UxN/2,t=Dx+UxN/2,t,Uit=ϕbxi,t,i=1,2,…,N−1,where(39)L1NUit=Ui′t−εδ2Uit+aitUit,L2NUit=Ui′t−εδ2Uit+aitUit+bitUi−N/2t.Problems (29), (30), and (32) satisfy the following well-known semidiscrete maximum principle on Ω¯τN×0,T.Lemma 3 (semidiscrete maximum principle).
Assume that(40)∑i=1Ngi−1+4gi+gi+13hi=ς<1,and mesh function Zit satisfies Z0t≥0,Zit≥0,andκNZNt≥0 Then, L1NZit≥0,∀xi,t∈Ω1N×0,T, L2NZit≥0,∀xi,t∈Ω2N×0,T, and DxZxN/2,t=Dx+ZxN/2,t−Dx−ZxN/2,t≤0 imply that Zit≥0,∀xi,t∈Ω¯N×0,T.Proof.
Define a test functionSxi,t as(41)Sxi,t=18+xi2,xi,t∈Ω1N×0,T,38+xi4,xi,t∈Ω2N×0,T.
Note that(42)Sxi,t>0,∀xi,t∈Ω¯N×0,T,LNSxi,t>0,∀xi,t∈Ω1N∪Ω2N×0,T,Sx0,t>0,Sxi,t>0,kNSxN,t>0,DxSxN/2,t<0.
Letρ=max−xxi,t/Sxi,t:xi,t∈Ω¯N×0,T. Then, there exists x∗,t∈Ω¯N×0,T such that Zx∗,t+ρSx∗,t=0 and Zxi,t+ρSxi,t≥0,∀xi,t∈Ω¯N×0,T. Therefore, the function attains its minimum at x,t=x∗,t. Suppose the theorem does not hold true, then ρ>0:Case (i):x∗,t=x0,t0<x+ρSx0,t=Zx0,t+ρSx0,t=0. It is a contradiction.Case (ii):x∗,t∈Ω1N×0,T0<L1Nx+ρSx∗,t=−εδx2+Dt+ax+ρSx∗,t≤0. It is a contradiction.Case (iii):x∗,t=xN/2,t0≤Dxx+ρSxN/2,t<0. It is a contradiction.Case (iv):x∗,t∈Ω2N×0,T0<L2Nx+ρSx∗,t=−εδx2+Dt+ax+ρSx∗,t+bx+ρSx∗−1,t≤0. It is a contradiction.Case (v):x∗,t=xN,t0<kNx+ρSxN,t=x+ρSxN,t−ε∑i=1Ngi−1x+ρSi−1t+4gix+ρSit+gi+1x+ρSi+1t/3hi≤0. It is a contradiction.
Hence, the theorem is proved.Lemma 4 (semidiscrete stability result).
LetZit be any mesh function; then,(43)ZitΩ¯N×0,T≤CmaxZitΓlN×0,T,ZitΓbN×0,T,κNZitΓrN×0,T,maxxi,t∈Ω1N∪Ω2N×0,TLNZit.Proof.
Consider the barrier functions(44)ϑ±xi,t=CMSxi,t±Zxi,t,xi,t∈Ω¯N×0,T,where(45)M=maxZΓlN×0,T,ZΓbN×0,T,κNZΓrN×0,T,maxxi,t∈Ω1N∪Ω2N×0,TLNZ,and Sxi,t is the test function as in Lemma 3.
From (44), it is clear that ϑ±x0,t≥0,ϑ±xi,t≥0, and κNϑ±xN,t≥0L1Nϑ±xi,t≥0 and(46)∀xi,t∈Ω1N,L2Nϑ±xi,t≥0,∀xi,t∈Ω2N,Dx+ϑ±xN/2,t−Dx−ϑ±xN/2,t≤0.
Using Lemma3, ϑ±xi,t≥0,∀xi,t∈Ω¯N×0,T.


### 3.2. Convergence Analysis

Now, let us analyze the convergence of the spatial semidiscretization. We proved above the discrete problem satisfying the maximum principle and the uniform stability estimate. Note thatUit is denoted for the spatial discretization approximate solution to the exact solution ux,t at x=xi,i=0,1,…,N.Lemma 5.
For a fixed mesh and for∀q∈Z+, we have(47)limε⟶0max1<i<N/2ε−q/2exp−Cxiε=0,q=1,2,…,limε⟶0max1<i<N/2ε−q/2exp−C1−xiε=0,q=1,2,…,where xi=ih,h=2/N,for alli=1,2,…,N/2.Proof.
Consider the partition of the interval [0, 1], with pointsxi=ih,h=2/N,andi=01N/2. By using x1=hand1−xN/2=h, then we have(48)max1<i<N/2ε−q/2exp−Cxiε≤ε−q/2exp−Cx1ε=ε−q/2exp−Chε,max1<i<N/2ε−q/2exp−C1−xiε≤ε−q/2exp−C1−xN/2ε=ε−q/2exp−Chε.
Then, applying L’Hospital’s rule repeatedly results in(49)limε⟶0ε−q/2exp−Chε=limp=1/ε⟶∞pqexp−Chp=limp⟶∞q!Ch−qexp−Chp=0.
Hence, the proof is completed.
Now, the truncation error of scheme (29) is given by(50)L1NUit−uit=fit−bitϕlxi−N/2,t−L1Nuit=L1−L1Nuit,i=11N2=duitdt−εuxxti+aituit−duitdt−εui+1t−2uit+ui−1tηi2+aituit=−εuxxti+εηi2ui+1t−2uit+ui−1t,i=11N2.
Taylor series expansions of the termsui+1tandui−1t are given as follows:(51)ui+1t=uit+huxti+h22uxxti+h36uxxxti+h424uxxxxti+h5120uxxxxxti+…,ui−1t=uit−huxti+h22uxxti−h36uxxxti+h424uxxxxti−h5120uxxxxxti+….
Using the truncated Taylor series expansions of the termsui+1tandui−1t yields(52)L1NUit−uit=−εuxxti+εηi2h2uxxti+h412uxxxxtiγi,γi∈ui−1,ui+1.
Next, we use a truncated Taylor series expansion of the denominator function1/ηi2 of order five:(53)1ηi2=1h2−ζi212+ζi4h2240.
Now, substituting (53) into (52), we obtain(54)L1NUit−uit=ε12uxxxxtγi−ζi2uxxtih2+εζi248ζi25uxxti−13uxxxxtiγih4+εζi42,880uxxxxtiγiih6.
Using the bounds on the derivatives and Lemma5 gives(55)L1NUit−uit≤ε121−ζi2h2−εζi24813−ζi25h4+εζi42880h6≤CN−2,where we have used the relation N−2>N−4>N−6… The discrete problem satisfies the following bound:(56)L1NUit−uit≤CN−2.
Using Lemma5, we get the result:(57)Uit−uit≤CN−2.Remark 1.
A similar analysis for convergence may be carried out for the finite difference scheme on (30).
Next, the truncation error of scheme (31) at the point i=N is given by(58)κNUNt−uNt=κNUNt−κNuNt=ϕrxi,t−κNuNt=κuxN,t−κNuxN,t=uxN,t−∫x0xNgxux,tdx−uxN,t+∑i=1Ngi−1ui−1t+4giuit+gi+1ui+1t3hi=g0u0t+4g1u1t+g2u2t3h1+⋯+gN−1uN−1t+4gNuNt+gN+1uN+1t3hN−∫x0x1gxux,tdx−∫x1x2gxux,tdx−⋯−∫xN−1xNgxux,tdx.
Now, we take the norm ofκN:(59)κNUNt−uNt≤Cεh13u″γ1,t+⋯+hN3u″γN,t≤Cεh13+⋯+hN3≤CN−2,where γi∈xi−1,xi,for1≤i≤N.
Applying Lemma4, we have the bound U−uxN,t≤CN−2.
Hence,uNt−UNt≤CN−2Theorem 4.
The error associated with the spatial semidiscretization of problems (29), (30), and (32) satisfies(60)sup0<ε≤1max0<i<Nuxi,t−UitΩxN×Λt≤CN−2,where C is a constant independent of N−1 and ε.


### 3.3. Temporal Semidiscretization

IVPs (29), (30), and (32) are discretized by the backward Euler method on a uniform mesh. Now, we denote the approximation of uitj:=uij by Uij. We perform the time discretization as follows: let M be a positive integer and t0=0,tj=t0+jτ,τ=tj−tj−1, and j=11M−1,tM=T, where M denotes the number of mesh points in time direction; then,(61)Uj−Uj−1τ+EtjUj=Gtj,forj=1,2,…,M,with initial condition U0=ϕltj. Rearranging equation (61) gives(62)Uj=I+τEtj−1τGtj+Uj−1.Lemma 6.
The local truncation error associated with the time integration satisfies(63)ej≤Cτ2,where C is a constant independent of ε and M.Proof.
The local truncation error is defined as(64)ej=utj−Uj=utj−I+τEtj−1τGtj+utj−1.A Taylor series expansion of utj−1 takes the form(65)utj−1=utj−τ∂utj∂t+τ22∂2utj∂t2−τ36∂3utj∂t3+Oτ4,∂utj∂t=Gtj−Etjutj,utj−1=utj−τGtj−Etjutj+τ22∂2utj∂t2−τ36∂3utj∂t3+Oτ4.
The local truncation error by subtractingej yields(66)ej=utj−I+τEtj−1I+τEtjutj+τ22∂2utj∂t2−τ36∂3utj∂t3+Oτ4=I+τEtj−1τ22∂2utj∂t2−τ36∂3utj∂t3+Oτ4.
Since the matrixHt is invertible, using the relation τ3<τ2 for small τ and utj≤C, we obtain(67)ej=I+τEtj−1τ2≤Cτ2,which completes the proof.Lemma 7.
In this temporal direction, the global error estimates are given by(68)TEj=maxjUitj−UijΩxN×ΛtM≤Cτ,where TEj is the global error in the temporal direction at jth time level.Proof.
One may refer [15] for the details.
SinceC and τ are independent of the perturbation parameter ε, taking the supremum for all ε∈0,1, we obtain(69)sup0<ε≤1maxjUitj−UijΩxN×ΛtM≤Cτ.
This shows that the discretization in time direction is consistent and global error is bounded, with the error boundOτ.


### 3.4. Fully Discretized Problem

Now, using (60) and (69), the ε−uniform convergence of the fully discrete scheme was proved as(70)sup0<ε≤1maxi,juxi,tj−Uij=sup0<ε≤1maxi,juxi,tj−Uitj+Uitj−Uij≤sup0<ε≤1maxi,juxi,tj−Uitj+sup0<ε≤1maxi,jUitj−Uij.Hence, we obtain the required bound as follows:(71)sup0<ε≤1uxi,tj−Uij≤ΩxN×ΛtMCN−2+τ.Thus, the inequality in (71) shows the parameter-uniform convergence of the proposed scheme with order of convergence, second order in spatial direction and first order in temporal direction.


## 3.1. Spatial Semidiscretization

The theoretical basis of the nonstandard discrete modeling method is based on the development of the exact finite difference method. In [29], Mickens presented techniques and rules for developing nonstandard FDMs for different problem types. In Mickens rules, to develop a discrete scheme, a denominator function for the discrete derivatives must be expressed in terms of more complicated functions of step sizes than those used in the standard procedures. This complicated function constitutes a general property of the schemes, which is useful while designing reliable schemes for such problems. In line to this, we need to derive a nonstandard finite difference method which captures the layer behaviour of the problem on a uniform mesh. Mickens in [29, 30] gave a novel approach of nonstandard finite difference method, and the basic idea of this method is to replace the denominator function h2 of the second-order derivative with a suitable function in the differential equation, which comes under the category of the fitted operator finite difference method. On the spatial domain [0, 2], we introduce the equidistant meshes with uniform mesh length such that(21)ΩxN=x0=0,xi=ih,i=11N−1,xN/2=1,xN=2,h=2N,where h is the step size and N is the number of mesh points in the spatial direction.Moreover, we consider the domain0,1 which is discretized into N/2 equal number of subintervals, each of length h. For the problem of the form in (6)–(8), in order to construct the exact finite difference scheme, we follow the procedures used in [31]. Consider the constant coefficient subequations given in (22) by ignoring the time variable as(22)−εd2uxdx2+αux=0,where ax,t≥α>0. Thus, problem (22) has two independent solutions, namely, expσ1xandexpσ2x, with(23)σ1,2=±αε.We denote the approximate solution ofux at xi′s by Ui. Now, our objective is to calculate a difference equation which has the same general solution as problem (22) has at the grid point xi given by Ui=Aexpσ1xi+Bexpσ2xi. Using the procedures used in [31], we have(24)detUi−1expσ1xi−1expσ2xi−1Uiexpσ1xiexpσ2xiUi+1expσ1xi+1expσ2xi+1=0.Simplifying (24), we obtain that(25)Ui−1−2coshαεhUi+Ui+1=0is an exact difference scheme for (22).After doing the arithmetic manipulation and rearrangement on (25), we obtain(26)−εUi−1−2Ui+Ui+14ε/αsinh2α/εh/2+αUi=0.The denominator function becomesη2=4ε/αsinh2α/εh/2. Adopting this denominator function for the variable coefficient problem, we write it as(27)ηi2ε,h,t=4εaitsinh2aitεh2.Case 1: consider (6)–(8) on the domain Ω=0,2×0,T,fori=1,…,N−1.LetUit denote the approximation of uxi,t. By using the nonstandard finite difference approximation, at this stage, the problem in (6)–(8) reduces to the semidiscrete form as(28)LlNUit=Git,fori=1,2,…,N−1,l=1,2,where(29)L1NUit≡dUitdt−εUi+1t−2Uit+Ui−1tηi2+aitUit=−bitϕlxi−N/2,t+fit,xi∈Ω1N,(30)L2NUit≡dUitdt−εUi+1t−2Uit+Ui−1tηi2+aitUit+bitUxi−N/2,t=fix,t,xi∈Ω2N.Case 2: since the problem involves nonlocal boundary conditions, we considered the following to obtain an equation at the end boundary condition.Fori=N, the composite Simpson’s integration rule approximates the nonlocal condition in (1) as follows:(31)κux,t≡u2,t−ε∫02gxux,tdx=UNt−εh∑i=0NcigiUit=ϕrxi,t,0<t≤T,where c0=cN=1/3andc2k−1=4/3,fork=1,2,…,N/2,andc2k=2/3,fork=1,2,…,N/2−1,Now, using (31) and left boundary condition, we can obtain the following formula for UNt:(32)UNt=11−εhcNgNεhc0g0ϕl0,t+εh∑i=1N−1cigiUit+ϕrxi,t.Therefore, on the whole domainD¯=0,2, the basic schemes to solve (6)–(8) are the schemes given in (29), (30), and (32).The time domain is continuous at this stage, and the system of IVPs in (29), (30), and (32) can be written in the compact form as(33)dUitdt+EitUit=Git,i=1,2,…,N,where Eit is a matrix of size NxN and UitandGit are vectors of size N. The entries of the coefficient matrix Eit are, respectively, given by(34)Eii−1t=−εηi2,fori=2,3,…,N−1,Eiit=2εηi2+ait,fori=1,2,…,N−1,Eii+1t=−εηi2,fori=1,2,…,N−1,Eii−N/2t=bit,fori=N2+1,…,N−1,EN,N=−4εh3∑i=1Ng2i−1U2i−1t−2εh3∑i=1N−1g2iU2it+1−εh3gN,i=N,and the corresponding right-hand-side vector in the equation system has the entries(35)G1t=f1t−b1tϕlx1−N/2,t+εη12ϕl0,t,Git=fit−bitϕlxi−N/2,t,fori=2,3,…,N2,fit,fori=N2+1,…,N−1,GNt=ϕrxi,t+εh3g0U0t,fori=N,U=U1,U2,…,UNT.Before we proceed with the convergence analysis, we highlight some properties of the discrete problem to original (1) in the form of lemma which plays a major role in the said analysis.Fori=1,2,…,N/2,(36)L1NUit=fit−bitϕlxi−N/2,t.Fori=N/2+1,…,N−1,(37)L2NUit=fit.Subject to the boundary conditions,(38)Uit=ϕlxi,t,i=−N2−N2+1,…,0andt∈0,T,κNUxN,t=UxN,t−ε∑i=1Ngi−1Ui−1t+4giUit+gi+1Ui+1t3hi=ϕrxi,t,Dx−UxN/2,t=Dx+UxN/2,t,Uit=ϕbxi,t,i=1,2,…,N−1,where(39)L1NUit=Ui′t−εδ2Uit+aitUit,L2NUit=Ui′t−εδ2Uit+aitUit+bitUi−N/2t.Problems (29), (30), and (32) satisfy the following well-known semidiscrete maximum principle on Ω¯τN×0,T.Lemma 3 (semidiscrete maximum principle).
Assume that(40)∑i=1Ngi−1+4gi+gi+13hi=ς<1,and mesh function Zit satisfies Z0t≥0,Zit≥0,andκNZNt≥0 Then, L1NZit≥0,∀xi,t∈Ω1N×0,T, L2NZit≥0,∀xi,t∈Ω2N×0,T, and DxZxN/2,t=Dx+ZxN/2,t−Dx−ZxN/2,t≤0 imply that Zit≥0,∀xi,t∈Ω¯N×0,T.Proof.
Define a test functionSxi,t as(41)Sxi,t=18+xi2,xi,t∈Ω1N×0,T,38+xi4,xi,t∈Ω2N×0,T.
Note that(42)Sxi,t>0,∀xi,t∈Ω¯N×0,T,LNSxi,t>0,∀xi,t∈Ω1N∪Ω2N×0,T,Sx0,t>0,Sxi,t>0,kNSxN,t>0,DxSxN/2,t<0.
Letρ=max−xxi,t/Sxi,t:xi,t∈Ω¯N×0,T. Then, there exists x∗,t∈Ω¯N×0,T such that Zx∗,t+ρSx∗,t=0 and Zxi,t+ρSxi,t≥0,∀xi,t∈Ω¯N×0,T. Therefore, the function attains its minimum at x,t=x∗,t. Suppose the theorem does not hold true, then ρ>0:Case (i):x∗,t=x0,t0<x+ρSx0,t=Zx0,t+ρSx0,t=0. It is a contradiction.Case (ii):x∗,t∈Ω1N×0,T0<L1Nx+ρSx∗,t=−εδx2+Dt+ax+ρSx∗,t≤0. It is a contradiction.Case (iii):x∗,t=xN/2,t0≤Dxx+ρSxN/2,t<0. It is a contradiction.Case (iv):x∗,t∈Ω2N×0,T0<L2Nx+ρSx∗,t=−εδx2+Dt+ax+ρSx∗,t+bx+ρSx∗−1,t≤0. It is a contradiction.Case (v):x∗,t=xN,t0<kNx+ρSxN,t=x+ρSxN,t−ε∑i=1Ngi−1x+ρSi−1t+4gix+ρSit+gi+1x+ρSi+1t/3hi≤0. It is a contradiction.
Hence, the theorem is proved.Lemma 4 (semidiscrete stability result).
LetZit be any mesh function; then,(43)ZitΩ¯N×0,T≤CmaxZitΓlN×0,T,ZitΓbN×0,T,κNZitΓrN×0,T,maxxi,t∈Ω1N∪Ω2N×0,TLNZit.Proof.
Consider the barrier functions(44)ϑ±xi,t=CMSxi,t±Zxi,t,xi,t∈Ω¯N×0,T,where(45)M=maxZΓlN×0,T,ZΓbN×0,T,κNZΓrN×0,T,maxxi,t∈Ω1N∪Ω2N×0,TLNZ,and Sxi,t is the test function as in Lemma 3.
From (44), it is clear that ϑ±x0,t≥0,ϑ±xi,t≥0, and κNϑ±xN,t≥0L1Nϑ±xi,t≥0 and(46)∀xi,t∈Ω1N,L2Nϑ±xi,t≥0,∀xi,t∈Ω2N,Dx+ϑ±xN/2,t−Dx−ϑ±xN/2,t≤0.
Using Lemma3, ϑ±xi,t≥0,∀xi,t∈Ω¯N×0,T.


## 3.2. Convergence Analysis

Now, let us analyze the convergence of the spatial semidiscretization. We proved above the discrete problem satisfying the maximum principle and the uniform stability estimate. Note thatUit is denoted for the spatial discretization approximate solution to the exact solution ux,t at x=xi,i=0,1,…,N.Lemma 5.
For a fixed mesh and for∀q∈Z+, we have(47)limε⟶0max1<i<N/2ε−q/2exp−Cxiε=0,q=1,2,…,limε⟶0max1<i<N/2ε−q/2exp−C1−xiε=0,q=1,2,…,where xi=ih,h=2/N,for alli=1,2,…,N/2.Proof.
Consider the partition of the interval [0, 1], with pointsxi=ih,h=2/N,andi=01N/2. By using x1=hand1−xN/2=h, then we have(48)max1<i<N/2ε−q/2exp−Cxiε≤ε−q/2exp−Cx1ε=ε−q/2exp−Chε,max1<i<N/2ε−q/2exp−C1−xiε≤ε−q/2exp−C1−xN/2ε=ε−q/2exp−Chε.
Then, applying L’Hospital’s rule repeatedly results in(49)limε⟶0ε−q/2exp−Chε=limp=1/ε⟶∞pqexp−Chp=limp⟶∞q!Ch−qexp−Chp=0.
Hence, the proof is completed.
Now, the truncation error of scheme (29) is given by(50)L1NUit−uit=fit−bitϕlxi−N/2,t−L1Nuit=L1−L1Nuit,i=11N2=duitdt−εuxxti+aituit−duitdt−εui+1t−2uit+ui−1tηi2+aituit=−εuxxti+εηi2ui+1t−2uit+ui−1t,i=11N2.
Taylor series expansions of the termsui+1tandui−1t are given as follows:(51)ui+1t=uit+huxti+h22uxxti+h36uxxxti+h424uxxxxti+h5120uxxxxxti+…,ui−1t=uit−huxti+h22uxxti−h36uxxxti+h424uxxxxti−h5120uxxxxxti+….
Using the truncated Taylor series expansions of the termsui+1tandui−1t yields(52)L1NUit−uit=−εuxxti+εηi2h2uxxti+h412uxxxxtiγi,γi∈ui−1,ui+1.
Next, we use a truncated Taylor series expansion of the denominator function1/ηi2 of order five:(53)1ηi2=1h2−ζi212+ζi4h2240.
Now, substituting (53) into (52), we obtain(54)L1NUit−uit=ε12uxxxxtγi−ζi2uxxtih2+εζi248ζi25uxxti−13uxxxxtiγih4+εζi42,880uxxxxtiγiih6.
Using the bounds on the derivatives and Lemma5 gives(55)L1NUit−uit≤ε121−ζi2h2−εζi24813−ζi25h4+εζi42880h6≤CN−2,where we have used the relation N−2>N−4>N−6… The discrete problem satisfies the following bound:(56)L1NUit−uit≤CN−2.
Using Lemma5, we get the result:(57)Uit−uit≤CN−2.Remark 1.
A similar analysis for convergence may be carried out for the finite difference scheme on (30).
Next, the truncation error of scheme (31) at the point i=N is given by(58)κNUNt−uNt=κNUNt−κNuNt=ϕrxi,t−κNuNt=κuxN,t−κNuxN,t=uxN,t−∫x0xNgxux,tdx−uxN,t+∑i=1Ngi−1ui−1t+4giuit+gi+1ui+1t3hi=g0u0t+4g1u1t+g2u2t3h1+⋯+gN−1uN−1t+4gNuNt+gN+1uN+1t3hN−∫x0x1gxux,tdx−∫x1x2gxux,tdx−⋯−∫xN−1xNgxux,tdx.
Now, we take the norm ofκN:(59)κNUNt−uNt≤Cεh13u″γ1,t+⋯+hN3u″γN,t≤Cεh13+⋯+hN3≤CN−2,where γi∈xi−1,xi,for1≤i≤N.
Applying Lemma4, we have the bound U−uxN,t≤CN−2.
Hence,uNt−UNt≤CN−2Theorem 4.
The error associated with the spatial semidiscretization of problems (29), (30), and (32) satisfies(60)sup0<ε≤1max0<i<Nuxi,t−UitΩxN×Λt≤CN−2,where C is a constant independent of N−1 and ε.


## 3.3. Temporal Semidiscretization

IVPs (29), (30), and (32) are discretized by the backward Euler method on a uniform mesh. Now, we denote the approximation of uitj:=uij by Uij. We perform the time discretization as follows: let M be a positive integer and t0=0,tj=t0+jτ,τ=tj−tj−1, and j=11M−1,tM=T, where M denotes the number of mesh points in time direction; then,(61)Uj−Uj−1τ+EtjUj=Gtj,forj=1,2,…,M,with initial condition U0=ϕltj. Rearranging equation (61) gives(62)Uj=I+τEtj−1τGtj+Uj−1.Lemma 6.
The local truncation error associated with the time integration satisfies(63)ej≤Cτ2,where C is a constant independent of ε and M.Proof.
The local truncation error is defined as(64)ej=utj−Uj=utj−I+τEtj−1τGtj+utj−1.A Taylor series expansion of utj−1 takes the form(65)utj−1=utj−τ∂utj∂t+τ22∂2utj∂t2−τ36∂3utj∂t3+Oτ4,∂utj∂t=Gtj−Etjutj,utj−1=utj−τGtj−Etjutj+τ22∂2utj∂t2−τ36∂3utj∂t3+Oτ4.
The local truncation error by subtractingej yields(66)ej=utj−I+τEtj−1I+τEtjutj+τ22∂2utj∂t2−τ36∂3utj∂t3+Oτ4=I+τEtj−1τ22∂2utj∂t2−τ36∂3utj∂t3+Oτ4.
Since the matrixHt is invertible, using the relation τ3<τ2 for small τ and utj≤C, we obtain(67)ej=I+τEtj−1τ2≤Cτ2,which completes the proof.Lemma 7.
In this temporal direction, the global error estimates are given by(68)TEj=maxjUitj−UijΩxN×ΛtM≤Cτ,where TEj is the global error in the temporal direction at jth time level.Proof.
One may refer [15] for the details.
SinceC and τ are independent of the perturbation parameter ε, taking the supremum for all ε∈0,1, we obtain(69)sup0<ε≤1maxjUitj−UijΩxN×ΛtM≤Cτ.
This shows that the discretization in time direction is consistent and global error is bounded, with the error boundOτ.


## 3.4. Fully Discretized Problem

Now, using (60) and (69), the ε−uniform convergence of the fully discrete scheme was proved as(70)sup0<ε≤1maxi,juxi,tj−Uij=sup0<ε≤1maxi,juxi,tj−Uitj+Uitj−Uij≤sup0<ε≤1maxi,juxi,tj−Uitj+sup0<ε≤1maxi,jUitj−Uij.Hence, we obtain the required bound as follows:(71)sup0<ε≤1uxi,tj−Uij≤ΩxN×ΛtMCN−2+τ.Thus, the inequality in (71) shows the parameter-uniform convergence of the proposed scheme with order of convergence, second order in spatial direction and first order in temporal direction.


## 4. Numerical Results and Discussion

To validate the established theoretical results, we perform numerical experiments using the proposed numerical scheme on the problem given in (6)–(8). We consider two numerical examples to verify the parameter-uniform convergence of the proposed scheme.Exact solution is not available for these model problems; therefore, maximum nodal errors are calculated by using the double mesh technique as(72)EεN,τ=maxi,jUi,jN,τ−Ui,j2N,τ/2,where N is the number of mesh points in the x direction. For any value of the mesh points N and τ, the ε−uniform (parameter-uniform) error estimate is calculated by(73)EN,τ=maxεEεN,τ.The rate of convergence of the method is given by(74)rN,τ=log2EN,τE2N,τ/2=logEN,τ−logE2N,τ/2log2.The numerical results are presented for the value of the perturbation parameter:(75)ε∈10−6,10−7,…,10−10.Example 1.
Consider the following singularly perturbed problem [27]:(76)∂ux,t∂t−ε∂2ux,t∂x2+5ux,t=ux−1,t+e−x,x,t∈0,2×0,2.
Subject to the initial and boundary conditions,(77)ux,t=0,∀x,t∈Γl,κu2,t=u2,t−ε∫02x3ux,tdx=0,∀x,t∈Γr,ux,t=0,∀x,t∈Γb.Example 2.
Consider the following singularly perturbed problem:(78)∂ux,t∂t−ε∂2ux,t∂x2+5ux,t=xux−1,t+1,x,t∈0,2×0,2.
Subject to the initial and boundary conditions,(79)ux,t=0,∀x,t∈Γl,κu2,t=u2,t−ε∫0216ux,tdx=0,∀x,t∈Γr,ux,t=sinx,t,∀x,t∈Γb.
In Tables1 and 2 the numerical results of the proposed method are tabulated in terms of maximum absolute errors, numerical rate of convergence, and uniform errors. As one observes from tables, for each number of mesh intervals N as ε⟶0, the maximum absolute error becomes stable and uniform. This shows that convergence of the proposed method is independent of ε. The ε−uniform error and ε−uniform rate of convergence of the method are observed from the last row of each table. In Tables 3 and 4, the performance of the proposed scheme is investigated by comparing it with recently published paper in [27]. As one can see, the proposed method gives more accurate results.
In Figures1–4, we observe the behaviour of the numerical solution, and for very small values of ε, the problem under consideration exhibits strong parabolic boundary layer near x=0andx=2; furthermore, an interior layer at x=1 is created.Table 1
The maximum pointwiseEεN,τ for different values of ε for Example 1 with T = 2.
ε↓N=16,τ=0.1/232,0.1/2264,0.1/23128,0.1/24256,0.1/25512,0.1/2610−63.5045e − 032.0026e − 031.0723e − 035.5505e − 042.8241e − 041.4245e − 0410−73.5045e − 032.0026e − 031.0723e − 035.5505e − 042.8241e − 041.4245e − 0410−83.5045e − 032.0026e − 031.0723e − 035.5505e − 042.8241e − 041.4245e − 0410−93.5045e − 032.0026e − 031.0723e − 035.5505e − 042.8241e − 041.4245e − 0410−103.5045e − 032.0026e − 031.0723e − 035.5505e − 042.8241e − 041.4245e − 04EN,τ3.5045e − 032.0026e − 031.0723e − 035.5505e − 042.8241e − 041.4245e − 04rN,τ0.807330.901170.950020.974830.98734—Table 2
The maximum pointwiseEεN,τ for different values of ε for Example 2 with T = 2.
ε↓N=16,τ=0.1/232,0.1/2264,0.1/23128,0.1/24256,0.1/25512,0.1/2610−62.9520e − 021.5980e − 028.3421e − 034.2681e − 032.1581e − 031.0853e − 0310−72.9520e − 021.5980e − 028.3421e − 034.2681e − 032.1581e − 031.0853e − 0310−82.9520e − 021.5980e − 028.3421e − 034.2681e − 032.1581e − 031.0853e − 0310−92.9520e − 021.5980e − 028.3421e − 034.2681e − 032.1581e − 031.0853e − 0310−102.9520e − 021.5980e − 028.3421e − 034.2681e − 032.1581e − 031.0853e − 03EN,τ2.9520e − 021.5980e − 028.3421e − 034.2681e − 032.1581e − 031.0853e − 03rN,τ0.885430.937780.966820.983830.99180—Table 3
Comparison ofε−uniform values for the method and result in [27] for Example 1.
Scheme↓N=16,τ=0.1/232,0.1/2264,0.1/23128,0.1/24256,0.1/25512,0.1/26Our resultEN,τ3.5045e − 032.0026e − 031.0723e − 035.5505e − 042.8241e − 041.4245e − 04rN,τ0.807330.901170.950020.974830.98734—Result in [27]EN,τ2.0615e − 021.2534e − 026.9738e − 033.6873e − 031.8972e − 039.6241e − 04rN,τ0.717830.845840.919370.958730.97912—Table 4
Comparison ofε−uniform values for the method and result in [27] for Example 2.
Scheme↓N=16,τ=0.1/232,0.1/2264,0.1/23128,0.1/24256,0.1/25512,0.1/26Our resultEN,τ2.9520e − 021.5980e − 028.3421e − 034.2681e − 032.1581e − 031.0853e − 03rN,τ0.885430.937780.966820.983830.99180—Result in [27]EN,τ1.8765e − 011.4776e − 019.7571e − 025.7092e − 023.1057e − 021.6222e − 02rN,τ0.344790.598730.773160.878370.93697—Figure 1
Numerical solution of Example1 when ε=10−3 for N=M=32 on the left-hand side (a) and when ε=10−6 and N=M=64 on the right-hand side (b).
(a)(b)Figure 2
Line graph for Example1 when ε=10−10 for N=64andτ=0.1/23.Figure 3
Numerical solution of Example2 when ε=10−3 for N=M=32 on the left-hand side (a) and when ε=10−10 and N=M=64 on the right-hand side (b).
(a)(b)Figure 4
Line graph for Example2 when ε=10−10 for N=64andτ=0.1/23.


## 5. Conclusions

In this paper, a parameter-uniform numerical scheme has been developed for solving singularly perturbed delay parabolic differential equations with integral boundary condition exhibiting parabolic boundary layers and an interior layer for small values of the perturbation parameter. The behaviour of the analytical solution and its boundedness are discussed. The numerical scheme is developed based on the method of a line that constitutes the nonstandard finite difference for the spatial discretization, and a backward Euler method is used in the temporal direction for the system of initial value problems resulting from the spatial discretization. To treat the integral boundary condition, the numerical integration technique is applied. The stability and convergence of the proposed scheme are analyzed. Two model examples have been considered to validate the applicability of the method. The computational results are presented in terms of tables and figures. The proposed numerical scheme could be the product for new equations (see e.g. [32–34]), is shown to be first-order convergent in time and second-order convergent in space, and gives more accurate results than the previously developed methods existing in the literature.



---

*Source: 9993644-2021-12-10.xml*