# Numerical Treatment on Parabolic Singularly Perturbed Differential Difference Equation via Fitted Operator Scheme

Authors: Dagnachew Mengstie Tefera; Awoke Andargie Tiruneh; Getachew Adamu Derese
Journal: Abstract and Applied Analysis
 (2021)

---

## Abstract

This paper proposes a new fitted operator strategy for solving singularly perturbed parabolic partial differential equation with delay on the spatial variable. We decomposed the problem into three piecewise equations. The delay term in the equation is expanded by Taylor series, the time variable is discretized by implicit Euler method, and the space variable is discretized by central difference methods. After developing the fitting operator method, we accelerate the order of convergence of the time direction using Richardson extrapolation scheme and obtainedOh2+k2 uniform order of convergence. Finally, three examples are given to illustrate the effectiveness of the method. The result shows the proposed method is more accurate than some of the methods that exist in the literature.


---

## Body

## 1. Introduction

The spatial delay parabolic singularly perturbed differential equation is a differential equation in which the perturbation parameter multiply the highest order derivative, and it has at least one retarded term on the spatial variable. Some mathematical problems can be treated as singularly perturbed problems such as Navier–Stokes equations, atmospheric pollution, turbulent transport, groundwater flow and solute transport, and Black–Scholes model, as presented in the survey paper by Kadalbajoo and Gupta [1]. There are two principle approaches for solving singular perturbation problems: numerical approach and asymptotic approach by Sharma et al. [2]. The numerical solution of the model problem is challenging due to the existence of a boundary layer. Furthermore, when using the classical finite difference method, we are unable to achieve an accurate solution, and the system becomes unstable. To tackle this issue, we need to create a fitted method for uniform or nonuniform meshes. As a result, ε-uniformly convergent numerical methods were constructed, with the order of convergence and error constant being independent of the perturbation parameter ε. For solving singular perturbation problems, some ε-uniform numerical schemes have been developed in the literature.Most of the scholars studied on time delayed singularly perturbed problems. Mbroh et al. [3] proposed parameter uniform method for solving a time delay nonautonomous singularly perturbed parabolic differential equation. Clavero and Gracia [4] studied singularly perturbed time-dependent problem of reaction–diffusion type using Richardson extrapolation technique. Woldaregay and Duressa [5] considered a numerical method for both small time delay and large time delay singularly perturbed boundary value problem. Kumar and Kumari [6] show that the influence of a small delay on the solution is extremely sensitive and that a small change in the delay can have a significant impact on the solution. Chakravarthy et al. [7] using the fitted technique, singularly perturbed differential equations with large delays were investigated.Nowadays, a few scholars have examined numerical solution for spatial delay singularly perturbed parabolic partial differential equations. Bansal and Sharma [8] numerical solutions for a large delay reaction-diffusion problem has been developed. Gupta et al. [9] examined spatial delay parabolic singularly perturbed partial differential equations and its solution using higher order fitted mesh method. Exponentially fitted operator method was developed to solve differential-difference singularly perturbed problem by Woldaregay and Duressa [10]. Das and Natesan [11] presented the solution of delay singularly perturbed partial differential equation using second-order convergent method. Singh and Srinivasan [12] developed Richardson extrapolation method for solving convection-diffusion equations with retarded term. Chahravarthy and Kumar [13] presented adaptive grid method for solving singularly perturbed convection-diffusion problems with spatial delay. Bansal et al. [14] designed numerical scheme for solving general shift singularly perturbed parabolic convectional diffusion problems.In this study, a singularly perturbed delayed partial differential equation with small spatial shift right boundary layer problem is decomposed into three piecewise equations which are treated using fitted operator difference methods. To accelerate the order of accuracy in the time variable, the Richardson extrapolation method is applied. The order of convergence of the present method is shown to be second order in both time and spatial variable, whereas the rate of convergence is two. Furthermore, the numerical results of the examples considered shows that the present method has better accuracy compared to some results that appear in the literature.


## 2. Statement of the Problem

Consider the parabolic singularly perturbed second-order differential equation with small spatial delay, on the domainΩ=D1×D2=0,1×0,T and Ω=Ωl∪Ωb∪Ωr, where Ωl=x,t:−γ≤x≤0and0≤t≤T, Ωb=x,t:0≤x≤1and0≤t≤T, and Ωr=x,t:1≤x≤1+μand0≤t≤T of the form:
(1)ut−εuxx+Axux+Bxux,t+Cxux−γ,t+Dxux+μ,t=Fx,t,x,t∈Ω,with
(2)ux,t=φlx,t,∀x,t∈Ωl,ux,0=φbx,x∈0,1,ux,t=φrx,t,∀x,t∈Ωr,where 0<ε≪1 and γ,μ denote the small delay parameters. The functions Ax,Bx,Cx,Dx,φlx,t,φbx,andφrx,t are supposed to be bounded and smooth functions on Ω¯, that satisfy the conditions Bx+Cx+Dx≥β>0 on Ωb=0,1. When γ=μ=0, the above problem (1) reduced to singularly perturbed parabolic partial differential problem. If Ax≥α>0, Cx<0 and Dx<0, ∀x∈Ωb=0,1, then the solution has boundary layer on the right side, i.e., at x=1, where α,β are some constants.Equation (1) together with initial boundary condition Equation (2) can be rewrite as
(3)Lεux,t=ut−εuxx+Axux+Bxu+Dxux+μ,t=Fx,t−Cxφlx−γ,t,if0<x≤γ,0<t≤T,ut−εuxx+Axux+Bxu+Cxux−γ,t+Dxux+μ,t=Fx,t,ifγ<x<1−μ,0<t≤T,ut−εuxx+Axux+Bxu+Cxux−γ,t=Fx,t−Dxφrx+μ,t,if1−μ≤x<1,0<t≤T.Under the assumption that the data are uniformly continuous and also meet relevant compatibility conditions at the corner points0,0,1,0,−γ,0,and1+μ,0, it is possible to establish the uniqueness of a solution to (1) [15]. The compatibility conditions are given as follows:
(4)φl0,0=φb0,φr1,0=φb1,∂φl0,0∂t−ε∂2φb0∂x2+A0∂φb0∂x+B0φb0+C0φl−γ,0+D0φbμ=F0,0,∂φr1,0∂t−ε∂2φb1∂x2+A1∂φb1∂x+B1φb1+C1φb1−γ+D0φr1+μ,0=F1,0.Using the above compatibility conditions, there exists a constantM independent of the perturbation parameter ε∀x,t∈Ω¯; we have
(5)∣ux,t−ux,0∣=∣ux,t−φb0∣≤Mt,ux,t−u0,t=ux,t−φl0,t≤M1−x,for the proof of Equation (5) see [16].Lemma 1 (Continuous Maximum Principle).
Let the functionψx,t∈C2,1Ω¯. If ψx,t≥0, for all x,t∈∂Ω and Lεψx,t≥0, for all x,t∈Ω, then ψx,t≥0, for all x,t∈Ω¯.Proof.
Assumexy,tz∈Ω¯−∂Ω such that ψxy,tz=minx,t∈Ω¯ψx,t and suppose ψxy,tz<0. This gives that ψtxy,tz=0, ψxxy,tz=0 and ψxxxy,tz≥0.☐Now, in order to showLεψx,t≤0, consider the following cases.Case 1.
0<xy≤γ and 0≤tz≤T(6)Lεuxy,tz=utxy,tz−εuxxxy,tz+Axyuxxy,tz+Bxyuxy,tz+Dxyuxy+μ,tz=−εuxxxy,tz+Bxyuxy,tz+Dxyuxy+μ,tz=−εuxxxy,tz+Bxy+Dxyuxy,tz+Dxyuxy+μ,tz−uxy,tz≤0.Case 2.
1−μ<xy<1 and 0≤tz≤T(7)Lεuxy,tz=utxy,tz−εuxxxy,tz+Axyuxxy,tz+Bxyuxy,tz+Cxyuxy−γ,tz=−εuxxxy,tz+Bxyuxy,tz+Cxyuxy−γ,tz=−εuxxxy,tz+Bxy+Cxyuxy,tz+Cxyuxy−γ,tz−uxy,tz≤0.Case 3.
γ<xy≤1−μ and 0≤tz≤T(8)Lεuxy,tz=utxy,tz−εuxxxy,tz+Axyuxxy,tz+Bxyuxy,tz+Cxyuxy−γ,tz+Dxyuxy+μ,tz=−εuxxxy,tz+Bxyuxy,tz+Cxyuxy−γ,tz+Dxyuxy+μ,tz≤0.SinceBxy+Cxy+Dxy≥α>0 and using the above assumption, we have Lεψx,t≤0 which contradicts our assumption.Thus,ψxy,tz≥0 which leads to ψx,t≥0∀x,y∈Ω¯.Lemma 2.
Letux,t be the solution of the continuous problem for Equations (1) and (2). Then, we get the bound
(9)∣ux,t∣≤∥f∥θ+max∣φlx,t∣,∣φbx∣,∣φrx,t∣,where Bx≥θ>0,∀x∈0,1.Proof.
Let the barrier functionψ± and defined as
(10)ψ±=∥f∥θ+maxφlx,t,φbx,φrx,t±ux,t,applying Lemma 1, and we get the above bound.☐Lemma 3.
The derivatives of the exact solutionux,t of the Equation (1) fulfill the following bound for v=0,1,and2.
(11)∂vux,t∂tv≤M,∀x,t∈Ω¯,where M∈ℝ is independent of ε.Proof.
The proof of this lemma can be found in [7].☐Theorem 4.
The analytical solution of Equation (1) satisfies
(12)∣diux,tdxi∣≤1+ε−iexp−α1−xε,0≤i≤4.Proof.
The proof for the bounds of its derivatives is given in [4].☐


## 3. Numerical Method

In this section, we develop an exponentially fitted operator difference method to solve Equation (1).


### 3.1. Time Discretization

A uniform mesh with a time step ofk is used to discretize the time domain 0,T as follows:
(13)tj=jk,forj=0,1,2,⋯,n,where k=T/n and in the interval 0,T, n is the number of subintervals in the time direction.We utilize the implicit Euler’s approach to approximate the time derivative term of Equation (1), which results in a system of boundary value problems.
(14)uj+1−ujk−εuxxj+1+Axuxj+1+Bxuj+1+Cxuj+1x−γ,t+Dxuj+1x+μ,t=Fx,t.Using Equations (3) and (14), we have
(15)Lεuj+1x=−εuxxj+1+Axuxj+1+Exuj+1+Dxuj+1x+μ=Fx,t−Cxφlx−γ,tj+1+ujk,if0<x≤γ,0<t≤T,−εuxxj+1+Axuxj+1+Exuj+1+Cxuj+1x−γ+Dxuj+1x+μ=Fx,t+ujk,ifγ<x<1−μ,0<t≤T,−εuxxj+1+Axuxj+1+Exuj+1+Cxuj+1x−γ=Fx,t−Dxφrx+μ,tj+1+ujk,if1−μ≤x≤1,0<t≤T,where Ex=Bx+1/k.


### 3.2. Spatial Discretization

The spatial domain0,1 is subdivided as follows using a uniform mesh with a step length of h:
(16)xi=ih,fori=0,1,2,⋯,m,where h=1/m and m is the number of subintervals in spatial direction in the interval 0,1. Using expansion of Taylor series expansion, we have
(17)uj+1x+μ=uj+1x+μuxj+1+μ22uj+1x+Oμ3,uj+1x−γ=uj+1x−γuxj+1+γ22uj+1x+Oγ3.Substituting Equation (17) into Equation (15) and rearranging gives
(18)Lεuj+1x=−εuxxj+1+p1xuxj+1+q1xuj+1=g1x,t,if0<x≤γ,0<t≤T,−εuxxj+1+p2xuxj+1+q2xuj+1=g2x,t,ifγ<x<1−μ,0<t≤T,−εuxxj+1+p3xuxj+1+q3xuj+1=g3x,t,if1−μ≤x≤1,0<t≤T,with the boundary conditions:
(19)u0,t=φl0,t,∀x,t∈Ωl,u1,t=φr1,t,∀x,t∈Ωr,and initial condition
(20)ux,0=φbx,x∈0,1, where
(21)p1x=εAx+μDxε−μ2/2Dx,p2x=εAx−γCx+μDxε−γ2/2Cx−μ2/2Dx,p3x=εAx−γCxε−γ2/2Cx,q1x=εDx+Exε−μ2/2Dx,q2x=εCx+Dx+Exε−γ2/2Cx−μ2/2Dx,q3x=εCx+Exε−γ2/2Cx,g1x,t=εFx,tj+1−Cxφlx−γ,tj+1+uj/kε−μ2/2Dx,g2x,t=εFx,tj+1+uj/kε−γ2/2Cx−μ2/2Dx,g3x,t=εFx,tj+1−Dxφrx+μ,tj+1+uj/kε−γ2/2Cx.Applying central difference approximation for spatial variable of Equation (18), we obtain
(22)−εui−1j+1−2uij+1+ui+1j+1h2+p1xiui+1j+1−ui−1j+12h+q1xiuij+1=g1xi,t,if0<x≤γ,0<t≤T,−εui−1j+1−2uij+1+ui+1j+1h2+p2xiui+1j+1−ui−1j+12h+q2xiuij+1=g2xi,t,ifγ<x<1−μ,0<t≤T,−εui−1j+1−2uij+1+ui+1j+1h2+p3xiui+1j+1−ui−1j+12h+q3xiuij+1=g3xi,t,if1−μ≤x≤1,0<t≤T.To obtain more accurate andε-uniform solution for Equation (22), introducing fitting factor σ as follows:
(23)−εσ1ui−1j+1−2uij+1+ui+1j+1h2+p1xiui+1j+1−ui−1j+12h+q1xiuij+1=g1xi,t,for0<x≤γ,0<t≤T,−εσ2ui−1j+1−2uij+1+ui+1j+1h2+p2xiui+1j+1−ui−1j+12h+q2xiuij+1=g2xi,t,forγ<x<1−μ,0<t≤T,−εσ3ui−1j+1−2uij+1+ui+1j+1h2+p3xiui+1j+1−ui−1j+12h+q3xiuij+1=g3xi,t,for1−μ≤x≤1,0<t≤T.To obtain the value of the fitting factor, multiply both sides of the equations in (23) by h and evaluating limit as h⟶0.
(24)limh⟶0−σ1ρui−1j+1−2uij+1+ui+1j+1+limh⟶0p1xiui+1j+1−ui−1j+12=0,where ρ=h/ε.The solution of Equation (18) comes from the theory of singular perturbation [16] given as
(25)ux,t=u0x,t+ae−pv11−x/ε,v=1,2,3,where u0x,t is outer solution.This implies the approximate solution atxi,tj is given as:
(26)limh⟶0uij=uoj0+ae−pv1/εepv1iρ.Using Equation (26), we have
(27a)limh⟶0ui+1j+1−ui−1j+1=ae−pv1/εepv1iρepv1ρ−e−pv1ρ,and
(27b)limh⟶0ui−1j+1−2uij+1+ui+1j+1=ae−pv1/εepv1iρepv1ρ−2+e−pv1ρ.Thus, using Equation (24), (27a), and (27b), we get
(28)σv=pviρ2cothpv1ρ2,v=1,2,3.Substituting Equation (28) into Equation (23), we get the following tri-diagonal system of equations that can be solved using Thomas algorithm.
(29)Ei,vj+1ui−1j+1+Fi,vj+1uij+1+Gi,vj+1ui+1j+1=Hi,vj+1,for i=1,2,⋯,m−1,j=0,1,⋯,n−1v=1,2,3, where
(30)Ei,vj+1=−εσvh2−pvxi2hFi,vj+1=2εσvh2+qvxi,Gi,vj+1=−εσvh2+pvxi2hHi,vj+1=gvxi,t.Sinceqvxi is nonnegative, then the system of Equation (29) becomes diagonally dominant.
(31)i.e.Fi,vj+1>Ei,vj+1+Gi,vj+1.Thus, the present method have convergent solution.


### 3.3. Convergence Analysis

Lemma 5 (Discrete Maximum Principle).
Suppose that the mesh functionwj+1xi satisfies wj+1x0≥0 and wj+1xm≥0. If Lεwj+1xi≥0 for 1≤i≤m−1, then wj+1xi>0, for all i,0≤i≤m.Proof.
Letwj+1xy=min1≤i≤m−1wj+1xi and suppose that wj+1xy<0, then
(32)Lεwj+1xy=−εσvwj+1xy−1−2wj+1xy+wj+1xy+1h2+pvxywj+1xy+1−wj+1xy−12h+qvwj+1xy<0,which contradicts our assumption.☐Thus,wj+1xi>0, for all i,0≤i≤m.Lemma 6 (Stability Estimate).
The solutionUij+1 of the discrete method satisfies the bound
(33)Uij+1≤θ−1maxLεUij+1+maxφl0,tj+1,φr1,tj+1,where qxi≥θ>0.Proof.
Define the barrier functionψi,j+1± as
(34)ψi,j+1±=θ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±Uij+1.
In the barrier function at the boundary condition, we obtain(35)ψj+1±0=θ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±φl0,tj+1≥0,ψj+1±1=θ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±φr1,tj+1≥0.
The barrier function at spatial domainxi, 1≤i≤m−1; we obtain
(36)Lεψi,j+1±=Lεθ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±LεUij+1=qvxiθ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±LεUij+1≥0.
Using Lemma5, we have that ψi,j+1±≥0 for all xi,tj+1∈Ω¯. Thus, the required bound in Equation (33) is satisfied.Consistency of the method: Consider the Taylor series expansion of
(37)ui+1j+1=uij+1+hux,ij+1+h22uxx,ij+1+h33!uxxx,ij+1+⋯ui−1j+1=uij+1−hux,ij+1+h22uxx,ij+1−h33!uxxx,ij+1+⋯uij=uij+1−kut,ij+1+k22utt,ij+1−k33!uttt,ij+1+⋯.Using the above expansion of Equation (22) becomes
(38)1kuij+1−uij+1−kut,ij+1+k22utt,ij+1−⋯−εh2uij+1+hux,ij+1+h22uxx,ij+1+h33!uxxx,ij+1+⋯−2uij+1+uij+1−hux,ij+1+h22uxx,ij+1−h33!uxxx,ij+1+⋯+Axi2huij+1+hux,ij+1+h22uxx,ij+1+h33!uxxx,ij+1+⋯−uij+1−hux,ij+1+h22uxx,ij+1−h33!uxxx,ij+1+⋯+Bxiuij+1+Cxiuxi−γ,tj+1+Dxiuxi+μ,tj+1−Fxi,tj+1.After simplifying the terms, the truncation error becomes(39)T.E.=ui,tj+1−εuxx,ij+1+Axiux,ij+1+Bxiuij+1+Cxiuxi−γ,tj+1+Dxiuxi+μ,tj+1uij+1−Fxi,tj+1−k2utt,ij+1−εh212uxxxx,ij+1+Axih26uxxx,ij+1,which gives
(40)T.E.=−k2utt,ij+1−εh212uxxxx,ij+1+Axih26uxxx,ij+1.Thus,limh,k⟶0,0T.E.=0, which shows that the present method is consistent. Since the method is consistent and stable, then using Lax equivalence theorem, the present method is convergent.Theorem 7.
Let the exact solution and numerical solution of Equation (1), respectively, are u and U. Then,
(41)sup0<ε≤1maxi,j∣uxi,tj−Uij∣≤Mh2+k.Proof.
For the proof, one can refer [17].☐


## 3.1. Time Discretization

A uniform mesh with a time step ofk is used to discretize the time domain 0,T as follows:
(13)tj=jk,forj=0,1,2,⋯,n,where k=T/n and in the interval 0,T, n is the number of subintervals in the time direction.We utilize the implicit Euler’s approach to approximate the time derivative term of Equation (1), which results in a system of boundary value problems.
(14)uj+1−ujk−εuxxj+1+Axuxj+1+Bxuj+1+Cxuj+1x−γ,t+Dxuj+1x+μ,t=Fx,t.Using Equations (3) and (14), we have
(15)Lεuj+1x=−εuxxj+1+Axuxj+1+Exuj+1+Dxuj+1x+μ=Fx,t−Cxφlx−γ,tj+1+ujk,if0<x≤γ,0<t≤T,−εuxxj+1+Axuxj+1+Exuj+1+Cxuj+1x−γ+Dxuj+1x+μ=Fx,t+ujk,ifγ<x<1−μ,0<t≤T,−εuxxj+1+Axuxj+1+Exuj+1+Cxuj+1x−γ=Fx,t−Dxφrx+μ,tj+1+ujk,if1−μ≤x≤1,0<t≤T,where Ex=Bx+1/k.


## 3.2. Spatial Discretization

The spatial domain0,1 is subdivided as follows using a uniform mesh with a step length of h:
(16)xi=ih,fori=0,1,2,⋯,m,where h=1/m and m is the number of subintervals in spatial direction in the interval 0,1. Using expansion of Taylor series expansion, we have
(17)uj+1x+μ=uj+1x+μuxj+1+μ22uj+1x+Oμ3,uj+1x−γ=uj+1x−γuxj+1+γ22uj+1x+Oγ3.Substituting Equation (17) into Equation (15) and rearranging gives
(18)Lεuj+1x=−εuxxj+1+p1xuxj+1+q1xuj+1=g1x,t,if0<x≤γ,0<t≤T,−εuxxj+1+p2xuxj+1+q2xuj+1=g2x,t,ifγ<x<1−μ,0<t≤T,−εuxxj+1+p3xuxj+1+q3xuj+1=g3x,t,if1−μ≤x≤1,0<t≤T,with the boundary conditions:
(19)u0,t=φl0,t,∀x,t∈Ωl,u1,t=φr1,t,∀x,t∈Ωr,and initial condition
(20)ux,0=φbx,x∈0,1, where
(21)p1x=εAx+μDxε−μ2/2Dx,p2x=εAx−γCx+μDxε−γ2/2Cx−μ2/2Dx,p3x=εAx−γCxε−γ2/2Cx,q1x=εDx+Exε−μ2/2Dx,q2x=εCx+Dx+Exε−γ2/2Cx−μ2/2Dx,q3x=εCx+Exε−γ2/2Cx,g1x,t=εFx,tj+1−Cxφlx−γ,tj+1+uj/kε−μ2/2Dx,g2x,t=εFx,tj+1+uj/kε−γ2/2Cx−μ2/2Dx,g3x,t=εFx,tj+1−Dxφrx+μ,tj+1+uj/kε−γ2/2Cx.Applying central difference approximation for spatial variable of Equation (18), we obtain
(22)−εui−1j+1−2uij+1+ui+1j+1h2+p1xiui+1j+1−ui−1j+12h+q1xiuij+1=g1xi,t,if0<x≤γ,0<t≤T,−εui−1j+1−2uij+1+ui+1j+1h2+p2xiui+1j+1−ui−1j+12h+q2xiuij+1=g2xi,t,ifγ<x<1−μ,0<t≤T,−εui−1j+1−2uij+1+ui+1j+1h2+p3xiui+1j+1−ui−1j+12h+q3xiuij+1=g3xi,t,if1−μ≤x≤1,0<t≤T.To obtain more accurate andε-uniform solution for Equation (22), introducing fitting factor σ as follows:
(23)−εσ1ui−1j+1−2uij+1+ui+1j+1h2+p1xiui+1j+1−ui−1j+12h+q1xiuij+1=g1xi,t,for0<x≤γ,0<t≤T,−εσ2ui−1j+1−2uij+1+ui+1j+1h2+p2xiui+1j+1−ui−1j+12h+q2xiuij+1=g2xi,t,forγ<x<1−μ,0<t≤T,−εσ3ui−1j+1−2uij+1+ui+1j+1h2+p3xiui+1j+1−ui−1j+12h+q3xiuij+1=g3xi,t,for1−μ≤x≤1,0<t≤T.To obtain the value of the fitting factor, multiply both sides of the equations in (23) by h and evaluating limit as h⟶0.
(24)limh⟶0−σ1ρui−1j+1−2uij+1+ui+1j+1+limh⟶0p1xiui+1j+1−ui−1j+12=0,where ρ=h/ε.The solution of Equation (18) comes from the theory of singular perturbation [16] given as
(25)ux,t=u0x,t+ae−pv11−x/ε,v=1,2,3,where u0x,t is outer solution.This implies the approximate solution atxi,tj is given as:
(26)limh⟶0uij=uoj0+ae−pv1/εepv1iρ.Using Equation (26), we have
(27a)limh⟶0ui+1j+1−ui−1j+1=ae−pv1/εepv1iρepv1ρ−e−pv1ρ,and
(27b)limh⟶0ui−1j+1−2uij+1+ui+1j+1=ae−pv1/εepv1iρepv1ρ−2+e−pv1ρ.Thus, using Equation (24), (27a), and (27b), we get
(28)σv=pviρ2cothpv1ρ2,v=1,2,3.Substituting Equation (28) into Equation (23), we get the following tri-diagonal system of equations that can be solved using Thomas algorithm.
(29)Ei,vj+1ui−1j+1+Fi,vj+1uij+1+Gi,vj+1ui+1j+1=Hi,vj+1,for i=1,2,⋯,m−1,j=0,1,⋯,n−1v=1,2,3, where
(30)Ei,vj+1=−εσvh2−pvxi2hFi,vj+1=2εσvh2+qvxi,Gi,vj+1=−εσvh2+pvxi2hHi,vj+1=gvxi,t.Sinceqvxi is nonnegative, then the system of Equation (29) becomes diagonally dominant.
(31)i.e.Fi,vj+1>Ei,vj+1+Gi,vj+1.Thus, the present method have convergent solution.


## 3.3. Convergence Analysis

Lemma 5 (Discrete Maximum Principle).
Suppose that the mesh functionwj+1xi satisfies wj+1x0≥0 and wj+1xm≥0. If Lεwj+1xi≥0 for 1≤i≤m−1, then wj+1xi>0, for all i,0≤i≤m.Proof.
Letwj+1xy=min1≤i≤m−1wj+1xi and suppose that wj+1xy<0, then
(32)Lεwj+1xy=−εσvwj+1xy−1−2wj+1xy+wj+1xy+1h2+pvxywj+1xy+1−wj+1xy−12h+qvwj+1xy<0,which contradicts our assumption.☐Thus,wj+1xi>0, for all i,0≤i≤m.Lemma 6 (Stability Estimate).
The solutionUij+1 of the discrete method satisfies the bound
(33)Uij+1≤θ−1maxLεUij+1+maxφl0,tj+1,φr1,tj+1,where qxi≥θ>0.Proof.
Define the barrier functionψi,j+1± as
(34)ψi,j+1±=θ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±Uij+1.
In the barrier function at the boundary condition, we obtain(35)ψj+1±0=θ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±φl0,tj+1≥0,ψj+1±1=θ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±φr1,tj+1≥0.
The barrier function at spatial domainxi, 1≤i≤m−1; we obtain
(36)Lεψi,j+1±=Lεθ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±LεUij+1=qvxiθ−1max∣LεUij+1∣+max∣φl0,tj+1,φr1,tj+1∣±LεUij+1≥0.
Using Lemma5, we have that ψi,j+1±≥0 for all xi,tj+1∈Ω¯. Thus, the required bound in Equation (33) is satisfied.Consistency of the method: Consider the Taylor series expansion of
(37)ui+1j+1=uij+1+hux,ij+1+h22uxx,ij+1+h33!uxxx,ij+1+⋯ui−1j+1=uij+1−hux,ij+1+h22uxx,ij+1−h33!uxxx,ij+1+⋯uij=uij+1−kut,ij+1+k22utt,ij+1−k33!uttt,ij+1+⋯.Using the above expansion of Equation (22) becomes
(38)1kuij+1−uij+1−kut,ij+1+k22utt,ij+1−⋯−εh2uij+1+hux,ij+1+h22uxx,ij+1+h33!uxxx,ij+1+⋯−2uij+1+uij+1−hux,ij+1+h22uxx,ij+1−h33!uxxx,ij+1+⋯+Axi2huij+1+hux,ij+1+h22uxx,ij+1+h33!uxxx,ij+1+⋯−uij+1−hux,ij+1+h22uxx,ij+1−h33!uxxx,ij+1+⋯+Bxiuij+1+Cxiuxi−γ,tj+1+Dxiuxi+μ,tj+1−Fxi,tj+1.After simplifying the terms, the truncation error becomes(39)T.E.=ui,tj+1−εuxx,ij+1+Axiux,ij+1+Bxiuij+1+Cxiuxi−γ,tj+1+Dxiuxi+μ,tj+1uij+1−Fxi,tj+1−k2utt,ij+1−εh212uxxxx,ij+1+Axih26uxxx,ij+1,which gives
(40)T.E.=−k2utt,ij+1−εh212uxxxx,ij+1+Axih26uxxx,ij+1.Thus,limh,k⟶0,0T.E.=0, which shows that the present method is consistent. Since the method is consistent and stable, then using Lax equivalence theorem, the present method is convergent.Theorem 7.
Let the exact solution and numerical solution of Equation (1), respectively, are u and U. Then,
(41)sup0<ε≤1maxi,j∣uxi,tj−Uij∣≤Mh2+k.Proof.
For the proof, one can refer [17].☐


## 4. Richardson Extrapolation Approach

The Richardson extrapolation approach has been described, and it is designed to improve the accuracy of the computed solutions in the basic scheme.LetD2n⊆D22n, where D22n is the mesh obtained from bisecting the step size k. Denote the numerical solution obtained from D22n by U¯jx, we have
(42a)ujx−Ujx≤Ck+Rjnx,x,tj∈D1×D2n,(42b)ujx−U¯jx≤Ck2+Rj2nx,x,tj∈D1×D22n,where Rjnx and Rj2nx are the remainder terms of the error. Now, subtracting the inequality (42b) from (42a) to obtain the extrapolation formula.
(43)ujx−Ujx−2ujx−U¯jx=Rjnx−Rj2n,which gives that
(44)Ujextx=2U¯jx−Ujx,is an approximate solution [12].Theorem 8.
Letuxi,tj+1 and Ui,j+1ext be the solution of problems in (1) and (44), respectively, then the proposed scheme satisfies the following error estimate
(45)sup0<ε<1maxxi,tj+1∣uxi,tj+1−Ui,j+1ext∣≤Ch2+k2.Proof.
Using the error for the temporal and spatial discretization gives the required bound.☐


## 5. Numerical Examples

To determine the efficacy of the current scheme, we looked at model problems that had been addressed in the literature and had approximate solutions that could be compared.We used the double-mesh principle to estimate the absolute maximum error of the current approach when the exact solution for the given problem was unknown. We use the following formula to approximate the absolute maximum error at the selected mesh points:Case 1.
If the exact solution is known,(46)EεM,N=maxxi,tj∈Ωuxi,tj−uiext,j.Case 2.
If the exact solution is unknown,(47)EεM,N=maxxi,tj∈Ωuiext,jM,N−uiext,j2M,2N.We also evaluate the corresponding rate of convergence.(48)RM,N=logEεM,N−logEε2M,2Nlog2.Example 1.
LetAx=2−x2,Bx=x−3,Cx=−2,Dx=−1,Fx,t=10t2e−tx1−x, where x,t∈0,1×0,1, and with initial boundary condition,
(49)ux,t=0,∀x,t∈Ωl=x,t:−γ≤x≤0,and0≤t≤1,ux,t=0,∀x,t∈Ωr=x,t:1≤x≤1+μ,and0≤t≤1,ux,t=0,∀x,t∈Ωb=x,t:0≤x≤1,and0≤t≤1.Example 2.
LetAx=1+x+x2,Bx=1+x2,Cx=−0.25+0.5x2,Dx=−0.25,Fx,t=sinπx1−x, where x,t∈0,1×0,1, with initial and boundary condition,
(50)ux,t=0,∀x,t∈Ωl=x,t:−γ≤x≤0,and0≤t≤1,ux,t=0,∀x,t∈Ωr=x,t:1≤x≤1+μ,and0≤t≤1,ux,t=0,∀x,t∈Ωb=x,t:0≤x≤1,and0≤t≤1.Example 3.
LetAx=1−x2/2,Bx=x+6,Cx=−4,Dx=−1,Fx,t=x1−x, where x,t∈0,1×0,3, with initial and boundary condition,
(51)ux,t=0,∀x,t∈Ωl=x,t:−γ≤x≤0,and0≤t≤3,ux,t=0,∀x,t∈Ωr=x,t:1≤x≤1+μ,and0≤t≤3,ux,t=0,∀x,t∈Ωb=x,t:0≤x≤1,and0≤t≤3.


## 6. Discussions and Results

We have presented the method for solving spatial delayed singularly perturbed parabolic partial differential equation. The basic mathematical procedures are defining the model problem, decomposing into three equations, approximating time variable using implicit Euler’s method, approximating the delay term using Taylor series expansion of order two, approximating the spatial variable using the central difference method, and finding fitting factor. Finally, apply Richardson extrapolation method to accelerate the accuracy of the method.Three model examples are used to exemplify the performance of the proposed method. The maximum error and rate of convergence are shown in Tables1–3 with different values of ε, delay parameters, and mesh length. The physical behavior of the solution are shown in Figures 1–4. We examined the suggested numerical scheme for stability, consistency, and ε uniform convergence. As shown in the result, the current method is second-order convergent with respect to time and spatial variables, the rate of convergence is two and more accurate than some of the methods that appear in the literature.Table 1
Maximum absolute point-wise error and rate of convergence for Example1 before and after Richardson extrapolation method applied, where γ=0.5ε, μ=0.6ε, and L≥8.
ε↓ M,N ⟶32,3264,64128,128256,256512,512Before10−25.5909e-032.7792e-031.3847e-036.9081e-043.4501e-0410−46.2590e-033.3063e-031.5449e-037.7016e-043.8448e-0410−66.2656e-033.1096e-031.5465e-037.7096e-043.8488e-0410−86.2657e-033.1096e-031.5465e-037.7097e-043.8488e-0410−106.2657e-033.1096e-031.5465e-037.7097e-043.8488e-04⋮⋮⋮⋮⋮⋮10−L6.2657e-033.1096e-031.5465e-037.7097e-043.8488e-04Rate of convergence10−L1.01071.00771.00431.0022After10−21.8542e-045.5542e-051.5119e-053.9402e-061.0055e-0610−42.2699e-046.6325e-051.7858e-054.6297e-061.1786e-0610−62.2742e-046.6437e-051.7886e-054.6369e-061.1804e-0610−82.2743e-046.6438e-051.7887e-054.6370e-061.1804e-0610−102.2743e-046.6438e-051.7887e-054.6370e-061.1804e-06⋮⋮⋮⋮⋮⋮10−L2.2743e-046.6438e-051.7887e-054.6370e-061.1804e-06Rate of convergence10−L1.77531.89311.94761.9739Table 2
Maximum absolute point-wise error rate of convergence for Example2 before and after Richardson extrapolation method applied, where γ=0.5ε, μ=0.6ε, and L≥8.
ε↓ M,N ⟶16,1632,3264,64128,128256,256Before10−28.9343e-034.4707e-032.2286e-031.1098e-035.5384e-0410−49.3800e-034.7241e-032.3512e-031.1687e-035.8202e-0410−69.3844e-034.7269e-032.3526e-031.1693e-035.8233e-0410−89.3845e-034.7269e-032.3526e-031.1694e-035.8233e-0410−109.3845e-034.7269e-032.3526e-031.1694e-035.8233e-04⋮⋮⋮⋮⋮⋮10−L9.3845e-034.7269e-032.3526e-031.1694e-035.8233e-04Rate of convergence10−L0.98941.00661.00851.0058After10−21.8103e-035.1433e-041.5297e-043.9939e-051.0095e-0510−42.0571e-036.2361e-041.9584e-045.1831e-051.3148e-0510−62.0599e-036.2482e-041.9636e-045.1977e-051.3185e-0510−82.0599e-036.2483e-041.9637e-045.1979e-051.3186e-0510−102.0599e-036.2483e-041.9637e-045.1979e-051.3186e-05⋮⋮⋮⋮⋮⋮10−L2.0599e-036.2483e-041.9637e-045.1979e-051.3186e-05Rate of convergence10−L1.72101.66991.91761.9789Table 3
Maximum absolute point-wise error and rate of convergence for Example3 before and after Richardson extrapolation method applied, where γ=0.1ε, μ=0.2ε, and L≥6.
ε↓ M,N ⟶8,816,1632,3264,64128,128Before10−26.1105e-034.0192e-032.4412e-031.4312e-037.9252e-0410−46.5144e-034.2051e-032.5534e-031.4901e-038.2602e-0410−66.5190e-034.2076e-032.5547e-031.4908e-038.2641e-0410−86.5190e-034.2076e-032.5547e-031.4908e-038.2641e-0410−106.5190e-034.2076e-032.5547e-031.4908e-038.2641e-04⋮⋮⋮⋮⋮⋮10−L6.5190e-034.2076e-032.5547e-031.4908e-038.2641e-04After10−21.9128e-039.9156e-044.3385e-042.1411e-044.8046e-0510−42.0383e-031.0236e-034.5811e-042.1921e-045.4179e-0510−62.0407e-031.0240e-034.5837e-042.1924e-045.4255e-0510−82.0407e-031.0240e-034.5837e-042.1924e-045.4256e-0510−102.0407e-031.0240e-034.5837e-042.1924e-045.4256e-05⋮⋮⋮⋮⋮⋮10−L2.0407e-031.0240e-034.5837e-042.1924e-045.4256e-05Figure 1
The physical behavior of the solutions for Example1 at m=n=32, ε=10−2, γ=0.5ε, and μ=0.6ε.Figure 2
The physical behavior of the solutions for Example1 at m=n=32, ε=10−2, γ=0.5ε, and μ=0.6ε at different time level.Figure 3
The physical behavior of the solutions for Example2 at m=n=64, ε=10−2, γ=0.5ε, and μ=0.6ε.Figure 4
The physical behavior of the solutions for Example3 at m=n=32, ε=10−2, γ=0.1ε, and μ=0.2ε.



---

*Source: 1661661-2021-10-22.xml*