# Uniformly Convergent Scheme for Singularly Perturbed Space Delay Parabolic Differential Equation with Discontinuous Convection Coefficient and Source Term

Authors: Mulunesh Amsalu Ayele; Awoke Andargie Tiruneh; Getachew Adamu Derese
Journal: Journal of Mathematics
 (2022)

---

## Abstract

A singularly perturbed delay parabolic problem of convection-diffusion type with a discontinuous convection coefficient and source term is examined. In the problem, strong interior layers and weak boundary layers are exhibited due to a large delay in the spatial variable and discontinuity of convection coefficient and source. The problem is discretized by a nonstandard finite difference scheme in the spatial variable and for the time derivative, we used the Crank–Nicolson scheme. To enhance the order of convergence of the spatial variable, the Richardson extrapolation technique is applied. The error analysis of the proposed scheme was carried out and proved that the scheme is uniformly convergent of second order in both spatial and temporal variables. Numerical experiments are performed to verify the theoretical estimates.


---

## Body

## 1. Introduction

In many real life situations, we encounter problems with having small parameters multiplying the highest order derivative terms, involving at least one shift term. We call these singularly perturbed delay differential equations. Such problems arise frequently in the mathematical modeling of various physical and biological phenomena. The delay terms in the models enable us to include some past behavior to get more practical models for the phenomena. For example, Stein’s model is a well-known space dependent model which represents a commonly-used description of spontaneous neuronal activity [1]. Many other examples can be found in [2, 3]. Extensive numerical methods have been developed for singularly perturbed delay differential equations, such as [4–7] and references therein.Among the recently conducted studies on time dependent large spatial delay differential equations, some to mention are [8–11] but still, all these are reaction-diffusion problems with smooth data.Nonetheless, there are numerical methods for singularly perturbed ordinary differential equations with nonsmooth data (discontinuous source term and/or convection coefficient) using special piecewise uniform meshes; see [12–16] and references therein.When we came to time dependent differential equations, in [17] the authors studied singularly perturbed parabolic delay differential equation with discontinuous coefficient and source term based on the upwind finite difference method on a specially generated mesh in the spatial direction with backward Euler method for the discretization of the time variable. Providing ε−uniform numerical method for singularly perturbed differential equations with discontinuous coefficients and source terms is not that easy. In this case, the situation is more complicated, especially when the delay is large.Motivated by the above-given works, we have designed a numerical scheme for singularly perturbed time dependent delay differential equation with discontinuous convection coefficient and source therm, using a nonstandard finite difference scheme in the spatial variable and Crank–Nicolson method in a temporal variable. To increase the order of convergence of the spatial variable, the Richardson extrapolation technique is applied. It is proved that the proposed scheme is uniformly convergent of orderh2+Δt2.The paper is organized as follows: in Section2, we formulate the problem. In Section 3, the bounds on the solution and its derivatives is discussed. Section 4 shows the derivation of the numerical scheme. In Section 5, the convergence of the full discrete scheme is analyzed. In Section 6, we describe the Richardson extrapolation technique and its convergence analysis. Numerical results are presented in Section 7, and lastly, conclusions are given in Section 8.


## 2. Problem Formulation

Consider the following singularly perturbed delay parabolic problems with a discontinuous convection coefficient and source term:(1)ε∂2u∂x2+αx∂u∂x−βxux,t−γxux−1,t−∂u∂t=fx,t,x,t∈Ω,subject to the following initial condition and interval boundary conditions(2)ux,0=ψ0x,x,t∈Γ0=x,t:0≤x≤2andt=0,ux,t=ψlx,t,x,t∈Γl=x,t:x,t∈−1,0×0,T,ux,t=ψrx,t,x,t∈Γr=x,t:x=2,and0≤t≤2,where 0<ε≪1 is the perturbation parameter, Ω=Ω1∪Ω2,Ω1=0,1×0,T,Ω2=1,2×0,T,Ω¯=0,2×0,T, the function βx and γx are sufficiently smooth functions such that βx>0,γx<0 and βx+γx≥0 for all x∈0,2. Moreover, assume that(3)αx=α1x,if0≤x≤1,α2x,if1≤x≤2,fx,t=f1x,t,ifx,t∈Ω¯1,f2x,t,ifx,t∈Ω2∗,α1x<−η1<−2η,α2x>η2>2η>0,α≤C,f≤C,where η=minη1,η2,Ω¯1=0,1×0,T and Ω2∗=1,2×0,T. The solution of (1) satisfies u=u1+,t−u1−,t=0 and ux=∂/∂x1+,t−∂/∂x1−,t=0 at x=1, here u1−,t and u1+,t are left and right side limit of u at x=1.The problem defined in (1) can be rewritten as follows:(4)ℒ˜εux,t=Fx,t,where(5a)ℒ˜εux,t≡ε∂2u∂x2+αx∂u∂x−βxux,t−∂u∂t,forx,t∈Ω1,ε∂2u∂x2+αx∂u∂x−βxux−γx,tux−1,t−∂u∂t,forx,t∈Ω2,(5b)Fx,t=fx,t−γx,tψx−1,t,forx,t∈Ω1,fx,t,forx,t∈Ω2,with initial and boundary condition(6)ux,0=ψ0x,x,t∈Γ0=x,t;0≤x≤2andt=0,ux,t=ψlx,t,x,t∈Γl=x,t;x,t∈−1,0×0,T,u1−,t=u1+,t,∂∂x1−,t=∂∂x1+,tu2,t=ψrx,t,x,t∈Γr=x,t;x=2and0≤t≤T,The solutionux,t of problem (1) exhibits a strong interior layer and weak boundary layer in the neighborhood of the point x=1 and x=2, respectively [17].


## 3. Bounds on the Solution and Its Derivatives

In this section, the analytical aspects of the solution of problem (1) and its derivatives are studied. The differential operator ℒ˜ε satisfies the following minimum principle.Lemma 1 (Minimum principle).
Supposeϑx,t∈C0,0Ω¯∩C1,0Ω∩C2,1Ω1∪Ω2 and assume that ϑ0,t≥0,ϑx,0≥0,ϑ2,t≥0, and ℒ˜εϑx,t≤0,∀x,t∈Ω and ϑx1,t=ϑx1+,t−ϑx1−,t≤0. Then, ϑx,t≥0,∀x,t∈Ω¯.Proof.
Define a test function(7)sx,t=32+x2,x,t∈0,1×0,2,3−x,x,t∈1,2×0,2.
Note that,sx,t>0,∀x,t∈Ω¯,ℒ˜εsx,t>0,∀x,t∈Ω1∪Ω2 and sx1,t<0. Let(8)μ=max−ϑx,tφx,t:x,t∈Ω¯.
Then, there existsx⋆,t⋆ such that ϑx⋆,t⋆+μsx⋆,t⋆=0 and ϑx,t+μsx,t≥0, ∀x,t∈Ω¯. Therefore, the function attains minimum at x⋆,t⋆. Suppose the theorem does not hold true, then μ>0.Case 1:x⋆,t⋆∈Ω1(9)0>ℒ˜εϑ+μsx⋆,t⋆=ε∂2ϑ+μs∂x2x⋆,t⋆+α1x⋆∂ϑ+μs∂xx⋆,t⋆−βx⋆ϑ+μsx⋆,t⋆−∂ϑ+μs∂tx⋆,t⋆≥0.Case 2:x⋆,t⋆=1,t⋆(10)0≤ϑ+μsx1,t⋆=ϑx1,t⋆+μsx1,t⋆<0.Case 3:x⋆,t⋆∈Ω2(11)0>ℒ˜εϑ+μsx⋆,t⋆=ε∂2ϑ+μs∂x2x⋆,t⋆+α2x⋆∂ϑ+μs∂xx⋆,t⋆−βx⋆ϑ+μsx⋆,t⋆−γx⋆ϑ+μsx⋆−1,t⋆−∂ϑ+μs∂tx⋆,t⋆=ε∂2ϑ+μs∂x2x⋆,t⋆+α2x⋆∂ϑ+μs∂xx⋆,t⋆−βx⋆+γx⋆ϑ+μsx⋆,t⋆−γx⋆ϑ+μsx⋆−1,t⋆−ϑ+μsx⋆−1,t⋆−∂ϑ+μs∂tx⋆,t⋆≥0.
In all the cases we reached at a contradiction. Hence, the required result follows.Lemma 2.
(Stability result). Ifux,t satisfies problem (1), then the bound(12)ux,tj+1Ω¯≤CmaxuΓl,uΓ0,ℒ˜εuΩ,uΓr.Proof.
Defining the following barrier functions(13)Φx,t±=CmaxuΓl,uΓ0,ℒ˜εuΩ,uΓrsx,t±ux,t,and using the minimum principle in Lemma 1 we can obtain the required estimate.


## 4. Description of the Numerical Scheme

To obtain the totally discrete scheme, we discretized the temporal variable and space variable separately, then formulate the fully discrete scheme.


### 4.1. Temporal Semidiscretization

On the time domain0,T, we use uniform mesh given by: ΩtM=tj:tj=jΔt,Δt=T/M, forj=0,1,…,M, where M is the number of mesh intervals in 0,T and Δt is the time step. Then, on ΩtM the continuous problem (1), is discretized by using the Crank–Nicolson method, defined by the following equation:(14)ε∂2∂x2uj+1/2+αx∂∂xuj+1/2−βxux,tj+1/2−γxux−1,tj+1/2−uj+1−ujΔt=fx,tj+1/2,where(15)ux,tj+1/2=ux,tj+1+ux,tj2,fx,tj+1/2=fx,tj+1+fx,tj2.After some rearrangement, the temporal discretized form of (1) is given by the following equation:(16)ℒεMux,tj+1=Qx,tj,where(17)ℒεMux,tj+1=ε∂2∂x2uj+1+αx∂∂xuj+1−βx+2Δtuj+1,forx∈0,1,ε∂2∂x2uj+1+αx∂∂xuj+1−βx+2Δtuj+1−γx,tux−1,tj+1,forx∈0,2,where(18)Qx,tj=−ε∂2∂x2uj−αx∂∂xuj+βx−2Δtuj+2γxux−1,tj+1/2+2fx,tj+1/2,x∈0,1,−ε∂2∂x2uj−αx∂∂xuj+βx−2Δtuj+γx,tux−1,tj+2fx,tj+1/2,x∈1,2.and uj+1x=ux,tj+1 is the semi-discrete approximation to the exact solution ux,t of (1) at the j+1th time level. The local truncation error of the semi-discrete method (16) is given by the following equation:(19)e˜j+1=ux,tj+1−u˜j+1x,where u˜jx is the solution evaluated after one step of the semi-discrete scheme (16) taking the exact value ux,tj instead of uj as the initial data.Lemma 3.
Letϑx,tj+1 be a smooth function such that ϑ0,tj+1≥0,ϑ2,tj+1≥0, ℒεMϑx,tj+1≤0, and ϑx1+,tj+1=ϑx1+,tj+1−ϑx1−,tj+1=0,∀x∈0,2. Then, ϑx,tj+1≥0,∀x∈0,2.Proof.
Apply the same procedure as Lemma 1.Lemma 4.
Ifux,tj+1 satisfies problem (16), then the bound(20)ux,tj+1≤maxu0,tj+1,Qη,u2,tj+1,for allx∈0,2.Proof.
Consider a barrier function(21)Φ±x,tj+1=maxu0,tj+1,Qη,u​2,tj+1±ux,tj+1,Clearly, Φ±x,tj+1≥0 and Φ±2,tj+1≥0. Also, for ℒεMΦ±0,tj+1 we have two cases.Case 1:x∈0,1(22)ℒεMΦ±x,tj+1=−β1˜xmaxux,tj+1,Qη,u​2,tj+1±ℒεMux,tj+1≤0,sinceβ1˜x=βx+2Δt>0,and using equation16Case 2:x∈0,2(23)ℒεMΦ±x,tj+1=−β2˜xmaxu0,tj+1,Qη,u​2,tj+1±ℒεMux,tj+1≤0,sinceβx+γ>0,β2˜x=βx+2Δt+γ>0and using equation16.
Therefore, from Lemma3 we get Φ±x,tj+1≥0 for all x∈0,2.Lemma 5.
Suppose that∂v/∂tvux,t≤C,x,t∈Ω¯,v=0,1,2,3, the local truncation error associated to scheme (16) satisfies:(24)e˜j+1∞≤C1Δt3,j=1,2,…,M.Proof.
Using Taylor’s series expansion, expandux,tj+1 and ux,tj centered at tj+1/2, we get the following equation:(25)ux,tj+1−ux,tjΔt=utx,tj+1/2+OΔt2.
Then, substitute (25) into (1), we obtain the following equation:(26)ux,tj+1−ux,tjΔt=εddx2ux,tj+1/2+αxddxux,tj+1/2−βxux,tj+1/2−γxux−1,tj+1/2−fx,tj+1/2+OΔt2,where(27)fx,tj+1/2=fx,tj+1+fx,tj2+OΔt2,ux,tj+1/2=ux,tj+1+ux,tj2+OΔt2.
From (26), the local truncation error e˜j+1 is the solution of the following BVP:(28)ℒεMe˜j+1=OΔt3e˜j+10=0,e˜j+12=0.
Next, using the maximum principle for the operatorℒεM proves the result, for further detail one can refer to [18].Theorem 1.
(Global error estimate). Under the hypothesis of Lemma5, the global error estimate Ej+1=ux,tj+1−uj+1x=∑k=1je˜k, associated with the Crank–Nicolson scheme in the time direction at j+1th time level is given by the following equation:(29)Ej+1∞≤CΔt2,forj=1,2,…,M.Proof.
Using the local error estimate given by Lemma 5, we obtain the following global error estimates atj+1th time level:(30)Ej+1=∑ξ=1je˜ξ≤e˜1+e˜2+⋯+e˜j≤C0jΔt3,by Lemma5≤C0jΔtΔt2,≤C0TΔt2,sincejΔt≤T≤CΔt2,C0T=C.Lemma 6.
The solution of the semidiscretized problem (16) satisfies the following equation:(31)dpux,tj+1dxp≤C1+ε−pexp−η1−xε,x∈0,1,1+ε−pexp−ηx−1ε+ε−p+1exp−η2−xε,x∈1,2,.forp=0,1,2,3,4.Proof.
See [14, 15].


### 4.2. Spatial Discretization

On the spatial domain [0, 2], is discretized uniformly as follows:(32)Ω¯xN=0=x0,x1,…,xN=2,where xi=x0+ih,h=2/Nfori=1,2,..,N−1,andN is the number of mesh intervals in the spatial variable. Thus, the discretized domain is defined as Ω¯N,M=Ω¯xN×ΩtM. In the spatial discretization, the problem (16) is further discretized using the nonstandard finite difference method. The main idea of the nonstandard finite difference method is to replace the denominator of the finite difference approximation of the derivatives by positive functions.


### 4.3. Nonstandard Finite Difference

For the construction of an exact finite difference scheme, consider the following constant coefficient homogeneous differential equation, corresponding to (16):(33)εddx2uj+1+α^ddxuj+1−β^uj+1=0,(34)εddx2uj+1+α^ddxuj+1=0,where αx>α^>0,βx+2/Δt>β^. Now, the homogeneous differential equation (33) possesses two linearly independent solutions, given as follows:(35)expλ1and expλ2,λ1,2=−α^±α^2+4εβ^2ε.The target is to find a difference equation which has the same general solution as the differential equations in (16) at the mesh point xi given by Ui,j+1=a1expλ1+a2λ2. Using the theory of difference equations in [19], by taking there consecutive points we can get the following set of equations:(36)Ui−1,j+1−uxi−1,tj+1=0,Ui,j+1−uxi,tj+1=0,Ui+1,j+1−uxi+1,tj+1=0,or,(37)Ui−1,j+1expλ1xi−1expλ2xi−1Ui,j+1expλ1xiexpλ2xiUi+1,j+1expλ1xi+1expλ2xi+1=0.Simplifying and substituting the values ofλ1 and λ2, we obtain the exact difference scheme(in the sense that it has the same general solution as the corresponding differential equation) for (33) as follows:(38)exp−α^h2εUi−1,j+1−2coshhα^2+4εβ^2εUi,j+1+expα^h2εUi+1,j+1=0.The main goal is finding a suitable denominator function for the second order derivative, the extraction of the denominator function from equation (38) is not straightforward. As a result, assume layer behaviors of the solution of problem (1) and that of the problem in the case when βx=0 are similar, so that in equation (38) use the approximation hα^2+4εβ^/2ε≈hα^/2ε and following steps in [20], we obtain the following equation:(39)εUi−1,j+1−2Ui,j+1+Ui+1,j+1hε/−α^exp−hα^/ε−1+α^Ui,j+1−Ui−1,j+1h=0,fori=1,2,…,N2εUi−1,j+1−2Ui,j+1+Ui+1,j+1hε/α^exphα^/ε−1+α^Ui+1,j+1−Ui,j+1h=0,fori=N2+1,N2+2,…,N−1.Therefore, the denominator function for the second order derivative approximation becomes(40)σ2ε,h=hε−α^exp−hα^ε−1,for0≤x≤1,hεα^exphα^ε−1,for1<x≤2.This denominator function can be modified to variable coefficient problem as follows:(41)σi2ε,h=hε−αxiexp−hαxiε−1,fori=1,2,…,N2,hεαxiexphαxiε−1,fori=N2+1,N2+2,…,N.Remark 1.
For the discretization of the first derivative of the spatial variable we use backward and forward finite differences depending on the convection coefficient term,(42)Dx−Uij+1=Uij+1−Ui−1jh,Dx+Uij+1=Ui+1j+1−Uij+1h.


### 4.4. Fully Discrete Scheme

Using the denominator functionσi2 in (41) into the discretized form of the scheme in (16), we obtain the following difference scheme:(43)ℒM,NUi,j+1=Qi,j+1,i=1,2,…,N−1,j=0,1,2,…,M−1,with the conditionsUi,0=ψ0xi,fori=0,1,2,…,N−1,Ui,j+1=ψlxi,tj+1,fori=−N2,−N2−1,…,0,UN,j+1=ψrtj+1,forj=0,1,2,…,M−1,where(44)ℒεN,MUi,j+1=εUi−1,j+1−2Ui,j+1+Ui+1,j+1σi2ε,h+α1xiUi,j+1−Ui−1,j+1h−β1xiUi,j+1,fori=1,2,…,N2,j=0,1,…,M−1,εUi−1,j+1−2Ui,j+1+Ui+1,j+1σi2ε,h+α2xiUi+1,j+1−Ui,j+1h−β1xiUi,j+1−γxiUi−N/2,j+1,fori=N2+1,N2+2,…,N−1,j=0,1,…,M−1.Qi,j=−εUi−1,j−2Ui,j+Ui+1,jσi2ε,h−α1xiUi,j−Ui−1,jh+β2xiUi,j+fxi,tj+1+fxi,tj+γxiUi−N/2,j+1+Ui−N/2,j,fori=1,2,…,N2,j=0,1,…,M−1,−εUi−1,j−2Ui,j+Ui+1,jσi2ε,h−α2xiUi+1,j−Ui,jh+β2xiUi,j+fxi,tj+1+fxi,tj+γxiUi−N/2,j,fori=N2+1,N2+2,…,N−1,j=0,1,…,M−1.where β1xi=βxi+2/Δt,β2xi=βxi−2/Δt.


## 4.1. Temporal Semidiscretization

On the time domain0,T, we use uniform mesh given by: ΩtM=tj:tj=jΔt,Δt=T/M, forj=0,1,…,M, where M is the number of mesh intervals in 0,T and Δt is the time step. Then, on ΩtM the continuous problem (1), is discretized by using the Crank–Nicolson method, defined by the following equation:(14)ε∂2∂x2uj+1/2+αx∂∂xuj+1/2−βxux,tj+1/2−γxux−1,tj+1/2−uj+1−ujΔt=fx,tj+1/2,where(15)ux,tj+1/2=ux,tj+1+ux,tj2,fx,tj+1/2=fx,tj+1+fx,tj2.After some rearrangement, the temporal discretized form of (1) is given by the following equation:(16)ℒεMux,tj+1=Qx,tj,where(17)ℒεMux,tj+1=ε∂2∂x2uj+1+αx∂∂xuj+1−βx+2Δtuj+1,forx∈0,1,ε∂2∂x2uj+1+αx∂∂xuj+1−βx+2Δtuj+1−γx,tux−1,tj+1,forx∈0,2,where(18)Qx,tj=−ε∂2∂x2uj−αx∂∂xuj+βx−2Δtuj+2γxux−1,tj+1/2+2fx,tj+1/2,x∈0,1,−ε∂2∂x2uj−αx∂∂xuj+βx−2Δtuj+γx,tux−1,tj+2fx,tj+1/2,x∈1,2.and uj+1x=ux,tj+1 is the semi-discrete approximation to the exact solution ux,t of (1) at the j+1th time level. The local truncation error of the semi-discrete method (16) is given by the following equation:(19)e˜j+1=ux,tj+1−u˜j+1x,where u˜jx is the solution evaluated after one step of the semi-discrete scheme (16) taking the exact value ux,tj instead of uj as the initial data.Lemma 3.
Letϑx,tj+1 be a smooth function such that ϑ0,tj+1≥0,ϑ2,tj+1≥0, ℒεMϑx,tj+1≤0, and ϑx1+,tj+1=ϑx1+,tj+1−ϑx1−,tj+1=0,∀x∈0,2. Then, ϑx,tj+1≥0,∀x∈0,2.Proof.
Apply the same procedure as Lemma 1.Lemma 4.
Ifux,tj+1 satisfies problem (16), then the bound(20)ux,tj+1≤maxu0,tj+1,Qη,u2,tj+1,for allx∈0,2.Proof.
Consider a barrier function(21)Φ±x,tj+1=maxu0,tj+1,Qη,u​2,tj+1±ux,tj+1,Clearly, Φ±x,tj+1≥0 and Φ±2,tj+1≥0. Also, for ℒεMΦ±0,tj+1 we have two cases.Case 1:x∈0,1(22)ℒεMΦ±x,tj+1=−β1˜xmaxux,tj+1,Qη,u​2,tj+1±ℒεMux,tj+1≤0,sinceβ1˜x=βx+2Δt>0,and using equation16Case 2:x∈0,2(23)ℒεMΦ±x,tj+1=−β2˜xmaxu0,tj+1,Qη,u​2,tj+1±ℒεMux,tj+1≤0,sinceβx+γ>0,β2˜x=βx+2Δt+γ>0and using equation16.
Therefore, from Lemma3 we get Φ±x,tj+1≥0 for all x∈0,2.Lemma 5.
Suppose that∂v/∂tvux,t≤C,x,t∈Ω¯,v=0,1,2,3, the local truncation error associated to scheme (16) satisfies:(24)e˜j+1∞≤C1Δt3,j=1,2,…,M.Proof.
Using Taylor’s series expansion, expandux,tj+1 and ux,tj centered at tj+1/2, we get the following equation:(25)ux,tj+1−ux,tjΔt=utx,tj+1/2+OΔt2.
Then, substitute (25) into (1), we obtain the following equation:(26)ux,tj+1−ux,tjΔt=εddx2ux,tj+1/2+αxddxux,tj+1/2−βxux,tj+1/2−γxux−1,tj+1/2−fx,tj+1/2+OΔt2,where(27)fx,tj+1/2=fx,tj+1+fx,tj2+OΔt2,ux,tj+1/2=ux,tj+1+ux,tj2+OΔt2.
From (26), the local truncation error e˜j+1 is the solution of the following BVP:(28)ℒεMe˜j+1=OΔt3e˜j+10=0,e˜j+12=0.
Next, using the maximum principle for the operatorℒεM proves the result, for further detail one can refer to [18].Theorem 1.
(Global error estimate). Under the hypothesis of Lemma5, the global error estimate Ej+1=ux,tj+1−uj+1x=∑k=1je˜k, associated with the Crank–Nicolson scheme in the time direction at j+1th time level is given by the following equation:(29)Ej+1∞≤CΔt2,forj=1,2,…,M.Proof.
Using the local error estimate given by Lemma 5, we obtain the following global error estimates atj+1th time level:(30)Ej+1=∑ξ=1je˜ξ≤e˜1+e˜2+⋯+e˜j≤C0jΔt3,by Lemma5≤C0jΔtΔt2,≤C0TΔt2,sincejΔt≤T≤CΔt2,C0T=C.Lemma 6.
The solution of the semidiscretized problem (16) satisfies the following equation:(31)dpux,tj+1dxp≤C1+ε−pexp−η1−xε,x∈0,1,1+ε−pexp−ηx−1ε+ε−p+1exp−η2−xε,x∈1,2,.forp=0,1,2,3,4.Proof.
See [14, 15].


## 4.2. Spatial Discretization

On the spatial domain [0, 2], is discretized uniformly as follows:(32)Ω¯xN=0=x0,x1,…,xN=2,where xi=x0+ih,h=2/Nfori=1,2,..,N−1,andN is the number of mesh intervals in the spatial variable. Thus, the discretized domain is defined as Ω¯N,M=Ω¯xN×ΩtM. In the spatial discretization, the problem (16) is further discretized using the nonstandard finite difference method. The main idea of the nonstandard finite difference method is to replace the denominator of the finite difference approximation of the derivatives by positive functions.


## 4.3. Nonstandard Finite Difference

For the construction of an exact finite difference scheme, consider the following constant coefficient homogeneous differential equation, corresponding to (16):(33)εddx2uj+1+α^ddxuj+1−β^uj+1=0,(34)εddx2uj+1+α^ddxuj+1=0,where αx>α^>0,βx+2/Δt>β^. Now, the homogeneous differential equation (33) possesses two linearly independent solutions, given as follows:(35)expλ1and expλ2,λ1,2=−α^±α^2+4εβ^2ε.The target is to find a difference equation which has the same general solution as the differential equations in (16) at the mesh point xi given by Ui,j+1=a1expλ1+a2λ2. Using the theory of difference equations in [19], by taking there consecutive points we can get the following set of equations:(36)Ui−1,j+1−uxi−1,tj+1=0,Ui,j+1−uxi,tj+1=0,Ui+1,j+1−uxi+1,tj+1=0,or,(37)Ui−1,j+1expλ1xi−1expλ2xi−1Ui,j+1expλ1xiexpλ2xiUi+1,j+1expλ1xi+1expλ2xi+1=0.Simplifying and substituting the values ofλ1 and λ2, we obtain the exact difference scheme(in the sense that it has the same general solution as the corresponding differential equation) for (33) as follows:(38)exp−α^h2εUi−1,j+1−2coshhα^2+4εβ^2εUi,j+1+expα^h2εUi+1,j+1=0.The main goal is finding a suitable denominator function for the second order derivative, the extraction of the denominator function from equation (38) is not straightforward. As a result, assume layer behaviors of the solution of problem (1) and that of the problem in the case when βx=0 are similar, so that in equation (38) use the approximation hα^2+4εβ^/2ε≈hα^/2ε and following steps in [20], we obtain the following equation:(39)εUi−1,j+1−2Ui,j+1+Ui+1,j+1hε/−α^exp−hα^/ε−1+α^Ui,j+1−Ui−1,j+1h=0,fori=1,2,…,N2εUi−1,j+1−2Ui,j+1+Ui+1,j+1hε/α^exphα^/ε−1+α^Ui+1,j+1−Ui,j+1h=0,fori=N2+1,N2+2,…,N−1.Therefore, the denominator function for the second order derivative approximation becomes(40)σ2ε,h=hε−α^exp−hα^ε−1,for0≤x≤1,hεα^exphα^ε−1,for1<x≤2.This denominator function can be modified to variable coefficient problem as follows:(41)σi2ε,h=hε−αxiexp−hαxiε−1,fori=1,2,…,N2,hεαxiexphαxiε−1,fori=N2+1,N2+2,…,N.Remark 1.
For the discretization of the first derivative of the spatial variable we use backward and forward finite differences depending on the convection coefficient term,(42)Dx−Uij+1=Uij+1−Ui−1jh,Dx+Uij+1=Ui+1j+1−Uij+1h.


## 4.4. Fully Discrete Scheme

Using the denominator functionσi2 in (41) into the discretized form of the scheme in (16), we obtain the following difference scheme:(43)ℒM,NUi,j+1=Qi,j+1,i=1,2,…,N−1,j=0,1,2,…,M−1,with the conditionsUi,0=ψ0xi,fori=0,1,2,…,N−1,Ui,j+1=ψlxi,tj+1,fori=−N2,−N2−1,…,0,UN,j+1=ψrtj+1,forj=0,1,2,…,M−1,where(44)ℒεN,MUi,j+1=εUi−1,j+1−2Ui,j+1+Ui+1,j+1σi2ε,h+α1xiUi,j+1−Ui−1,j+1h−β1xiUi,j+1,fori=1,2,…,N2,j=0,1,…,M−1,εUi−1,j+1−2Ui,j+1+Ui+1,j+1σi2ε,h+α2xiUi+1,j+1−Ui,j+1h−β1xiUi,j+1−γxiUi−N/2,j+1,fori=N2+1,N2+2,…,N−1,j=0,1,…,M−1.Qi,j=−εUi−1,j−2Ui,j+Ui+1,jσi2ε,h−α1xiUi,j−Ui−1,jh+β2xiUi,j+fxi,tj+1+fxi,tj+γxiUi−N/2,j+1+Ui−N/2,j,fori=1,2,…,N2,j=0,1,…,M−1,−εUi−1,j−2Ui,j+Ui+1,jσi2ε,h−α2xiUi+1,j−Ui,jh+β2xiUi,j+fxi,tj+1+fxi,tj+γxiUi−N/2,j,fori=N2+1,N2+2,…,N−1,j=0,1,…,M−1.where β1xi=βxi+2/Δt,β2xi=βxi−2/Δt.


## 5. Uniform Convergence of the Fully Discrete Scheme

Lemma 7 (Discrete Minimum principle).
Letϑi,j+1 be a mesh function such that ϑi,j+1≥0,∀i=0,1,…,N,ℒεN,Mϑi,j+1≤0 for all i∈0,1,…,N/N/2 and Dx+ϑN/2,j+1−Dx−ϑN/2,j+1=0. Then, ϑi,j+1≥0 for all i=0,1,…,N.Proof.
Define a test function(45)sxi,tj+1=32+xi2,xi,tj+1∈Ω1N,M,3−xi,xi,tj+1∈Ω2N,M,where Ω1N,M=0,1∩Ω¯xN×ΩtM,Ω2N,M=1,2∩Ω¯xN×ΩtM. Note that, sxi,tj+1>0,∀xi,tj+1∈Ω¯N,M,ℒεsxi,tj+1>0,∀xi,tj+1∈Ω∩Ω¯xN×ΩtM and sxN/2,tj+1<0. Let(46)μ=max−ϑxi,tj+1sxi,tj+1:xi,tj+1∈Ω¯N,M.
Then, there existsxw,tj+1∈Ω¯N,M such that ϑxw,tj+1+μsxw,tj+1=0 and ϑxi,tj+1+μsxi,tj+1≥0, ∀xi,tj+1∈Ω¯N,M. Therefore, the function attains a minimum at xw,tj+1. Suppose the theorem does not hold true, then μ>0.Case 1:w∈1,2,…N/2−1(47)0>ℒεN,Mϑ+μsw,j+1=εϑ+μsw−1,j+1−2ϑ+μsw,j+1+ϑ+μsw+1,j+1σw2ε,h+α1xwϑ+μsw,j+1−ϑ+μsw−1,j+1h−β1xwϑ+μsw,j+1=εϑ+μsw−1,j+1−ϑ+μsw,j+1σw2ε,h+ϑ+μsw+1,j+1−ϑ+μsw,j+1σw2ε,h+α1xwϑ+μsw,j+1−ϑ+μsw−1,j+1h−β1xwϑ+μsw,j+1>0,sinceα1xw<0.Case 2:xw,tj+1=xN/2,tj+1(48)0≤ϑ+μsxxN/2,tj+1=ϑxxN/2,tj+1+μsxxN/2,tj+1<0.Case 3:w∈N/2+1,N/2+2,…,N−1(49)0>ℒεN,Mϑ+μsw,j+1=εϑ+μsw−1,j+1−2ϑ+μsw,j+1+ϑ+μsw+1,j+1σw2ε,h+α2xwϑ+μsw,j+1−ϑ+μsw−1,j+1h−β1xwϑ+μsw,j+1−γxwϑ+μsw−N/2,j+1=εϑ+μsw−1,j+1−ϑ+μsw,j+1σw2ε,h+ϑ+μsw+1,j+1−ϑ+μsw,j+1σw2+α2xwϑ+μsw+1,j+1−ϑ+μsw,j+1h−β1xw+γxwϑ+μsw,j+1−γxwϑ+μsw−N/2,j+1−ϑ+μsw,j+1>0,sinceα2xw>0,β1xw+γxw>0andγxw<0,in all cases, the result contradicts with our assumption. Hence, the required result follows.
An immediate consequence of the discrete minimum principle is the following uniform stability property of the operatorℒεN,M.Lemma 8.
The solutionUi,j+1 of the discrete scheme in (25) satisfies the bound(50)Ui,j+1≤Cmaxψl0,tj+1,ψ00,tj+1,ℒεN,MUi,j+1,ψr2,tj+1,∀i=0,1,…N−1,j=1,2,…,M−1.Proof.
Define barrier function(51)Φi,j+1±=Cmaxψl0,tj+1,ℒεN,MUi,j+1sxi,tj+1,ψr1,tj+1±ϑi,j+1,then following the procedure applied for the proof of Lemma 7, we obtain Φ0,j+1±≥0,ΦN,j+1±≥0 and ℒεN,MΦi,j+1±≥0. The required result follows from Lemma 7.


### 5.1. Error Estimates

Theorem 2.
Letαx and Qx,tj+1 be sufficiently smooth functions so that Uj+1x∈C40,2. Then, the truncation error of the discrete scheme satisfies the bound.(52)ℒεN,MUj+1xi−Ui,j+1≤Ch1+maxxi∈0,1∩Ω¯xNexp−α1−xi/εε3,Ch1+maxxi∈1,2∩Ω¯xNexp−ηxi−1/εε3+maxxi∈1,2∩Ω¯xNexp−η2−xi/εε5.Proof.
The truncation error in the spatial discretization is(53)ℒεN,MUj+1xi−Ui,j+1=ℒεN,MUj+1xi−ℒεN,MUi,j+1≤Cεd2dx2−D+D−h2σi2ε,hUj+1xi+αiddx−Dx−Uj+1xi≤Cεd2dx2−Dx+Dx−Uj+1xi+Cεh2σi2ε,h−1Dx+Dx−Uj+1xi+Ch​d2dx2Uj+1xi≤Cεh2d4dx4Uj+1xi+Chd2dx2Uj+1xi,depending on the behavior of σi2ε,h, the estimate ε|h2/σε,h−1|≤Ch. Then apply Lemma 6 in to the truncation error bound in (53).Case 1: Wheni∈1,2,…N/2−1(54)ℒεN,MUj+1xi−Ui,j+1≤Cεh21+ε−4exp−η1−xiε+Ch1+ε−2exp−η1−xiε≤Ch2ε+ε−3exp−η1−xiε+Ch21+ε−2exp−η1−xiε≤Ch1+maxexp−η1−xi/εε3.Case 2: Wheni∈N/2+1,N/2+2,…,N−1, applying the same procedure like case 1, we will get the following equation:(55)ℒεN,MUj+1xi−Ui,j+1≤Ch1+maxxi∈Ω2Nexp−ηxi−1/εε3+maxxi∈Ω2Nexp−η2−xi/εε5.Lemma 9.
For a fixed mesh, asε⟶0, it holds(56)limε⟶0max1≤i≤N−1exp−ηxi/εεp=0,limε⟶0max1≤i≤N−1exp−η2−xi/εεp=0,∀p∈Z+,where xi=ih,h=2/N,∀i=1,2,…,N−1.Proof.
For the proof see [21].Theorem 3.
Letu and Ui,j+1 are solutions of problem (1) and problem (25), respectively. Then,(57)Sup0<ε≪1uxi,tj+1−Ui,j+1≤CN−1+Δt2.Proof.
Using Lemma 9 in Theorem 1 and also apply discrete minimum principle Lemma 8, we obtain the following equation:(58)Sup0<ε≪1uxi,tj−Ui,j+1≤Sup0<ε≪1uxi,tj+1−Uj+1xi+Uj+1xi−Ui,j+1≤CN−1+Δt2.


## 5.1. Error Estimates

Theorem 2.
Letαx and Qx,tj+1 be sufficiently smooth functions so that Uj+1x∈C40,2. Then, the truncation error of the discrete scheme satisfies the bound.(52)ℒεN,MUj+1xi−Ui,j+1≤Ch1+maxxi∈0,1∩Ω¯xNexp−α1−xi/εε3,Ch1+maxxi∈1,2∩Ω¯xNexp−ηxi−1/εε3+maxxi∈1,2∩Ω¯xNexp−η2−xi/εε5.Proof.
The truncation error in the spatial discretization is(53)ℒεN,MUj+1xi−Ui,j+1=ℒεN,MUj+1xi−ℒεN,MUi,j+1≤Cεd2dx2−D+D−h2σi2ε,hUj+1xi+αiddx−Dx−Uj+1xi≤Cεd2dx2−Dx+Dx−Uj+1xi+Cεh2σi2ε,h−1Dx+Dx−Uj+1xi+Ch​d2dx2Uj+1xi≤Cεh2d4dx4Uj+1xi+Chd2dx2Uj+1xi,depending on the behavior of σi2ε,h, the estimate ε|h2/σε,h−1|≤Ch. Then apply Lemma 6 in to the truncation error bound in (53).Case 1: Wheni∈1,2,…N/2−1(54)ℒεN,MUj+1xi−Ui,j+1≤Cεh21+ε−4exp−η1−xiε+Ch1+ε−2exp−η1−xiε≤Ch2ε+ε−3exp−η1−xiε+Ch21+ε−2exp−η1−xiε≤Ch1+maxexp−η1−xi/εε3.Case 2: Wheni∈N/2+1,N/2+2,…,N−1, applying the same procedure like case 1, we will get the following equation:(55)ℒεN,MUj+1xi−Ui,j+1≤Ch1+maxxi∈Ω2Nexp−ηxi−1/εε3+maxxi∈Ω2Nexp−η2−xi/εε5.Lemma 9.
For a fixed mesh, asε⟶0, it holds(56)limε⟶0max1≤i≤N−1exp−ηxi/εεp=0,limε⟶0max1≤i≤N−1exp−η2−xi/εεp=0,∀p∈Z+,where xi=ih,h=2/N,∀i=1,2,…,N−1.Proof.
For the proof see [21].Theorem 3.
Letu and Ui,j+1 are solutions of problem (1) and problem (25), respectively. Then,(57)Sup0<ε≪1uxi,tj+1−Ui,j+1≤CN−1+Δt2.Proof.
Using Lemma 9 in Theorem 1 and also apply discrete minimum principle Lemma 8, we obtain the following equation:(58)Sup0<ε≪1uxi,tj−Ui,j+1≤Sup0<ε≪1uxi,tj+1−Uj+1xi+Uj+1xi−Ui,j+1≤CN−1+Δt2.


## 6. Richardson Extrapolation Technique

To increase the accuracy of the numerical solution of the proposed scheme, we used the Richardson extrapolation technique in a spatial variable. LetUi,j+12N be the solution of the discrete problem (25) on the mesh Ω2N,M. From Theorem 2, we have the following equation:(59)Uj+1xi−Ui,j+1≤CN−1+ON−2≤CN−1+RN,where RN is the remainder term in the spatial direction. Then, (59) also holds for any h/2≠0, which is(60)Uj+1xi−Ui,j+12N≤CN−1/2+R2N−1.Next, to eliminate the termON−1 subtract twice of (60) from (59), we get the following equation:(61)Uj+1xi−Ui,j+1−2Uj+1xi−Ui,j+12N≤RN−2R2N,after simplifying, we obtain the following equation:(62)Uj+1xi−2Ui,j+12N−Ui,j+1≤ON−2.Therefore, we used the following extrapolation formula:(63)U¯i,j+1=2Ui,j+12N−Ui,j+1.The approximate solutionU¯i,j+1 is more accurate than either Ui,j+1N or Ui,j+12N. The truncation error of the spatial discretization in the approximation of (63) becomes(64)max∀xi∈0,2Uj+1xi−U¯i,j+1≤CN−2.For more detail see [22].Theorem 4.
Letuxi,tj+1 and U¯i,j+1 be the solution of problems in (1) and (63) respectively, then the proposed scheme satisfies the following error estimate(65)Sup0<ε≪1max0≤xi,tj+1≤2uxi,tj+1−U¯i,j+1≤CN−2+Δt2.Proof.
Combining the error in the temporal and spatial discretization gives the required bound.


## 7. Numerical Illustration

In this section, we test the performance of the proposed scheme through numerical experiments. The exact solution of the problem is not known, so to compute the maximum point-wise errors, we use the double mesh principle given by the formula,(i)
Before extrapolation(66)EεN,M=max0≤i,j≤N,MUN,Mxi,tj−U2N,2Mx2i,t2j.(ii)
After extrapolation(67)EεN,M=max0≤i,j≤N,MU¯N,Mxi,tj−U¯2N,2Mx2i,t2j,where UN,Mxi,tj and U¯xi,tj denote the numerical solutions obtained by N mesh intervals in the spatial direction and M mesh intervals in the time direction, such that M=T/Δt. The corresponding rate of convergence by the following formula is(68)ROCεN,M=log2EεN,MEε2N,2M.Also, theε−uniform maximum point-wise error EN,M is computed as follows:(69)EN,M=maxεEεN,M,and the corresponding ε−uniform rate of convergence ROCN,M is given by the following equation:(70)ROCN,M=lnEN,ME2N,2M.Example 1.
Consider the following singularly perturbed convection diffusion problem [17]:(71)ε∂2u∂x2+αx∂u∂x−5ux,t+2ux−1,t−∂u∂t=fx,t,x,t∈0,2×0,2,ux,0=0,x∈0,2,ux,t=0,x,t∈−1,0×0,2,u2,0=0,t∈0,2,where(72)αx=−4+x2,x∈0,1,8−x2,x∈1,2,fx,t=4xt2e−t,x,t∈0,1×0,2,42−xt2e−t,x,t∈1,2×0,2.Example 2.
Consider the following singularly perturbed convection diffusion problem:(73)ε∂2u∂x2+αx∂u∂x−x+3ux,t+ux−1,t−∂u∂t=fx,t,x,t∈0,2×0,2,ux,0=0,x∈0,2,ux,t=0,x,t∈−1,0×0,2,u2,0=0,t∈0,2,where(74)αx=exx−4,x∈0,1,ex+x,x∈1,2,fx,t=xt,x,t∈0,1×0,2,cosπx2t,x,t∈1,2×0,2.


## 8. Conclusion

Singularly perturbed time dependent convection-diffusion problem with discontinuous convection coefficient and source term is treated via nonstandard finite difference method for the spatial derivative and Crank–Nicolson method for the time derivative. In addition to this, to enhance the order of convergence of the spatial variable, the Richardson extrapolation technique is used. Due to the presence of large delay and discontinuity in coefficient and source term, the problem exhibits strong interior layers atx = 1 and weak boundary layer atx = 2. Figures 1(a), 1(b) and 2–4 of the surface plot for the numerical solution of the problem in Examples 1 and 2, respectively, demonstrate the existence of a strong boundary layer at x = 1. Figures 3(a) and 3(b) The maximum point-wise error of Examples 1 and 2 is plotted in Figures 3(b) and 4(b), respectively, in the log-log plot. From these figures, it is clear to observe that the maximum error is ε−uniform convergent of second order in both space and time variable. The error analysis of the proposed scheme is carried out and the scheme is ε-uniform convergent of order Oh2+Δt2. This result clearly shown in Tables 1–4 of the analogous computed errors and uniform convergence rates for the given examples, before and after the extrapolation technique. The results confirm that the theoretical error estimates are in agreement with the numerical results.Figure 1
Surface plot of numerical solution of Example1 at N = 32, M = 128. (a) ε=100, (b) ε=10−4.
(a)(b)Figure 2
Surface plot of numerical solution of Example2 at N = 32, M = 128. (a) ε=100, (b) ε=10−4.
(a)(b)Figure 3
One dimensional plot and log-log plot for Example1 (a) Numerical solution at t = 2, when N = 32, (b) Log-log plot.
(a)(b)Figure 4
One dimensional plot and log-log plot for Example2 (a) Numerical solution at t = 2, when N = 32, (b) Log-log plot.
(a)(b)Table 1
Computed maximum point-wise errors,ε-uniform errors EN,M, and the ε−uniform rate of convergence ROCN,M for Example 1 when M = N, before extrapolation.
ε↓N⟶326412825651210−32.3872e − 031.1497e − 035.6604e − 042.8110e − 041.4002e − 0410−42.3872e − 031.1497e − 035.6604e − 042.8110e − 041.4002e − 0410−52.3872e − 031.1497e − 035.6604e − 042.8110e − 041.4002e − 0410−62.3872e − 031.1497e − 035.6604e − 042.8110e − 041.4002e − 0410−72.3872e − 031.1497e − 035.6604e − 042.8110e − 041.4002e − 0410−82.3872e − 031.1497e − 035.6604e − 042.8110e − 041.4002e − 0410−92.3872e − 031.1497e − 035.6604e − 042.8110e − 041.4002e − 0410−102.3872e − 031.1497e − 035.6604e − 042.8110e − 041.4002e − 04EN,M2.3872e − 031.1497e − 035.6604e − 042.8110e − 041.4002e − 04ROCN,M1.05411.02231.00981.0055Table 2
Computed maximum point-wise errors,ε-uniform errors EN,M, and the ε−uniform rate of convergence ROCN,M for Example 1, after extrapolation.
ε↓N/M⟶16/6432/12864/256128/512256/102410−31.6250e − 034.6468e − 041.2399e − 043.2008e − 058.1246e − 0610−41.6250e − 034.6468e − 041.2399e − 043.2008e − 058.1304e − 0610−51.6250e − 034.6468e − 041.2399e − 043.2008e − 058.1304e − 0610−61.6250e − 034.6468e − 041.2399e − 043.2008e − 058.1304e − 0610−71.6250e − 034.6468e − 041.2399e − 043.2008e − 058.1304e − 0610−81.6250e − 034.6468e − 041.2399e − 043.2008e − 058.1304e − 0610−91.6250e − 034.6468e − 041.2399e − 043.2008e − 058.1304e − 0610−101.6250e − 034.6468e − 041.2399e − 043.2008e − 058.1304e − 06EN,M1.6250e − 034.6468e − 041.2399e − 043.2008e − 058.1304e − 06ROCN,M1.80621.90601.95371.9770Table 3
Computed maximum point-wise errors,ε-uniform errors EN,M and the ε−uniform rate of convergence ROCN,M for Example 2 when M = N, before extrapolation.
ε↓N⟶326412825651210−33.6703e − 031.9480e − 031.0030e − 035.0885e − 042.5616e − 0410−43.6703e − 031.9480e − 031.0030e − 035.0885e − 042.5627e − 0410−53.6703e − 031.9480e − 031.0030e − 035.0885e − 042.5627e − 0410−63.6703e − 031.9480e − 031.0030e − 035.0885e − 042.5627e − 0410−73.6703e − 031.9480e − 031.0030e − 035.0885e − 042.5627e − 0410−83.6703e − 031.9480e − 031.0030e − 035.0885e − 042.5627e − 0410−93.6703e − 031.9480e − 031.0030e − 035.0885e − 042.5627e − 0410−103.6703e − 031.9480e − 031.0030e − 035.0885e − 042.5627e − 04EN,M3.6703e − 031.9480e − 031.0030e − 035.0885e − 042.5627e − 04ROCN,M0.91390.95770.97900.9896Table 4
Computed maximum point-wise errors,ε-uniform errors EN,M, and the ε−uniform rate of convergence ROCN,M for Example 2, after extrapolation.
ε↓N/M⟶16/6432/12864/256128/512256/102410−31.3438e − 034.0562e − 041.1130e − 042.9142e − 057.4504e − 0610−41.3438e − 034.0562e − 041.1130e − 042.9142e − 057.4556e − 0610−51.3438e − 034.0562e − 041.1130e − 042.9142e − 057.4556e − 0610−61.3438e − 034.0562e − 041.1130e − 042.9142e − 057.4556e − 0610−71.3438e − 034.0562e − 041.1130e − 042.9142e − 057.4556e − 0610−81.3438e − 034.0562e − 041.1130e − 042.9142e − 057.4556e − 0610−91.3438e − 034.0562e − 041.1130e − 042.9142e − 057.4556e − 0610−101.3438e − 034.0562e − 041.1130e − 042.9142e − 057.4556e − 06EN,M1.3438e − 034.0562e − 031.1130e − 032.9142e − 047.4556e − 04ROCN,M1.72811.86581.93321.9677



---

*Source: 1874741-2022-10-22.xml*