# A Nonstandard Fitted Operator Method for Singularly Perturbed Parabolic Reaction-Diffusion Problems with a Large Time Delay

Authors: Awoke Andargie Tiruneh; Getachew Adamu Derese; Dagnachew Mengstie Tefera
Journal: International Journal of Mathematics and Mathematical Sciences
 (2022)

---

## Abstract

In this paper, we design and investigate a higher orderε-uniformly convergent method to solve singularly perturbed parabolic reaction-diffusion problems with a large time delay. We use the Crank–Nicolson method for the time derivative, while the spatial derivative is discretized using a nonstandard finite difference approach on a uniform mesh. Furthermore, to improve the order of convergence, we used the Richardson extrapolation technique. The designed scheme converges independent of the perturbation parameter (ε-uniformly convergent) and also achieves fourth-order convergent in both time and spatial variables. Two model examples are considered to demonstrate the applicability of the suggested method. The proposed method produces better accuracy and a higher rate of convergence than some methods that appear in the literature.


---

## Body

## 1. Introduction

Many practical problems which are modelled by parameter dependent differential equations and delay singularly perturbed differential equations have many application in engineering and applied mathematics such as the effect of various physical parameter on heat and mass transfer [1, 2], heat transport in the presence of magnetic field and modulation [3, 4], ground water flow, extraction of oil from reservoirs, population dynamics, and fluid flow. The solution of time delayed singularly perturbed partial differential equation is not only determined by its present state but also by a certain past state. One of the difficulties in solving singularly perturbed problems is developing ε− uniformly convergent methods. As far as finite difference methods are concerned, two approaches have ordinarily been used by numerical analysis to obtain ε− uniformly convergent methods for solving singularly perturbed problems, and these are namely fitted operator and fitted mesh finite difference methods [5].The majority of researchers focused on the numerical solution of singularly perturbed partial differential equations without a delay parameter [6–10]. Numerical approaches for solving singularly perturbed parabolic time delay reaction-diffusion equations that are ε− uniformly convergent are not well established. Gowrisankar and Natesan [11] presented a numerical scheme for singularly perturbed parabolic delay convection diffusion equation using piecewise fitted mesh approach, and the method have first-order in time and first-order up to a logarithmic factor in space. Singh et al. [12] analyzed a domain decomposition method for solving singularly perturbed parabolic reaction-diffusion problems with time delay, and the method have rate of convergence one. Yulan et al. [13] examined time delayed singularly perturbed reaction-diffusion problem using the barycentric interpolation collocation method but they did not show the convergence analysis. Kumar and Kumari [14] introduced a cubic B-spline method overfitted mesh for a class of singularly perturbed delay parabolic partial differential equations first-order in t and second-order accurate in x. Dagnachew et al. [15] designed a fitted operator method for singularly perturbed parabolic convection diffusion equation with small space delays, and the method have second order in both spatial and temporary variables. Kumar [16] developed a collocation method consisting of cubic B-spline basis functions over piecewise-uniform mesh for singularly perturbed parabolic convection diffusion problems with a time delay. Govindarao et al. [17] presented for the solution of singularly perturbed delay parabolic reaction-diffusion initial-boundary-value problem using Shishkin mesh strategy. Mesfin and Gemechis [18] introduced numerical solution for singularly perturbed parabolic PDEs with shift parameter using the Crank–Nicolson method in the temporal discretization and nonstandard finite difference method in the spatial discretization, and the method have order of convergence OΔx+Δt2. Fasica and Gemechis [19] proposed a nonstandard finite difference method for singularly perturbed parabolic reaction-diffusion subject to Robin boundary conditions, and the method have the rate of convergence one.In this study, we developed a nonstandard finite difference method (NSFDM) for solving singularly perturbed parabolic reaction-diffusion problems with a large time delay. We applied a nonstandard finite difference method (Mickens rule) over the Crank–Nicolson finite difference method. To accelerate the order of convergence, the Richardson extrapolation technique is applied. The current method is more accurate and has faster rate of convergence than those methods that appear in the literature. The present method has a fourth-order of convergence in both the spatial and temporal variables after Richardson extrapolation technique applied.Notation. Throughout this work,M and C are positive constants that are independent of the perturbation parameter ε and m and n. The norm ∥.∥ denotes the maximum norm.


## 2. Description of the Problem

We consider the singularly perturbed parabolic reaction-diffusion initial boundary value problem with time delay overD=Dx×Dt with Dx=0,1, Dt=0,T, and Ω=Ωl∪Ωb∪Ωr of the form:(1)ut−ϵuxx+px,tux,t+qx,tux,t−τ=gx,t,x,t∈D,ux,t=ψbx,t,onΩb=x,t:0≤x≤1,−τ≤t≤0,u0,t=ψlt,onΩl=0,t:0≤t≤T,u1,t=ψrt,onΩr=1,t:0≤t≤T,where 0<ϵ≪1 and τ>0 are the given constants. The functions px,t,qx,t,gx,t over D¯ and ψbx,t,ψlt,ψrtoverΩ¯ are sufficiently smooth and bounded functions that satisfy(2)px,t≥α>0,x,t∈D¯.For some positive integerb, the terminal time T is considered to satisfy the condition T=bτ. The existence of a unique solution of equation (1) is guaranteed with the assumption that the data are smooth and satisfy appropriate compatibility conditions at the corner points 0,0,1,0,0,−τ,and1,−τ. The delay terms and the required compatibility condition at the corner points are ψb0,0=ψl0,ψb1,0=ψr0, and(3)dψl0dt−ϵ∂2ψb0,0∂x2+p0,0ψb0,0+q0,0ψb0,−τ=g0,0,dψr0dt−ϵ∂2ψb1,0∂x2+p1,0ψb1,0+q1,0ψb1,−τ=g1,0.Equation (1) admits a unique solution under the above assumptions and compatibility conditions, and the solution exhibits boundary layers along x=0,x=1 in [20].(4)ux,t−ψbx,0≤Ct,x,t∈D,ux,t≤C,x,t∈D¯,where C is the positive real constant. For the proof of equation (4), see Govindarao and Mohapatra [21].Lemma 1.
The solution of the equation (1) satisfies the following estimate:Lemma 2.
Continuous maximum principle: we assume thatθx,t is any sufficiently smooth function satisfying θx,t≥0, ∀x,t∈∂D. Then, Lϵθx,t≥0, ∀x,t∈D implies that θx,t≥0, ∀x,t∈D¯, where Lϵθx,t=θtx,t−ϵθxxx,t+px,tθx,t.Proof.
We assumex∗,t∗∈D¯−∂D such that θx∗,t∗=minx,t∈D¯θx,t and suppose θx∗,t∗<0. This gives that θtx∗,t∗=0, θxx∗,t∗=0, and θxxx∗,t∗≥0, and this implies(5)Lϵθx∗,t∗=θtx∗,t∗−ϵθxxx∗,t∗+px∗,t∗x∗,t∗=−ϵθxxx∗,t∗+px∗,t∗θx∗,t∗≤0.
We haveLϵθx,t≤0 which contradicts our assumption.
Thus,θx∗,t∗≥0 which leads to θx,t≥0, ∀x,t∈D¯.Lemma 3.
Stability estimate: letux,t be the solution of the continuous problem for equation (1), and we have(6)ux,t≤α−1Lϵu+uΩ,where Lϵu=ut−ϵuxx+px,tux,t.Proof.
Let the two barrier functionπ± be defined as
π±=α−1Lϵu+uΩ±ux,t, x,t∈D¯, we have(7)π±0,t=α−1Lϵu+uΩ±u0,t≥uΩ±u0,t≥0,π±1,t=α−1Lϵu+uΩ±u1,t≥uΩ±u1,t≥0.
Also, forx,t∈Ωb,(8)π±x,t=α−1Lϵu+uΩ±ux,t≥uΩ±ux,t≥0.
Furthermore,∀x,t∈D(9)Lϵπ±x,t=πt±x,t−ϵπxx±x,t+px,tπ±x,t=px,tα−1Lϵu+uΩ±Lϵux,t≥Lϵu+αuΩ±Lϵux,t≥Lϵu±Lϵux,t≥0.
Therefore, using Lemma2, we obtained the desired outcome.


## 3. Derivation of the Numerical Method

In this section, we develop a nonstandard finite difference method over Crank–Nicolson scheme to solve equation (1).


### 3.1. Time Discretization

We first perform a semidiscretization of the problem in time direction by using the average of implicit and explicit Euler method with uniform time-stepsk=T/n in the time-interval 0,T. This results a time-mesh Dtn≡tj=jk,j=0,1,2,…,n, and the method results in a system of boundary value problems.(10)uj+1x−ujxk−ϵ2uxxj+1x+uxxjx+12pj+1xuj+1x+pjxujx=−12qj+1xψbx,tj+1−1+qjxψbx,tj−1+12gj+1x+gjx,for0<x<1and0≤j≤nb,uj+1x−ujxk−ϵ2uxxj+1x+uxxjx+12pj+1xuj+1x+pjxujx=−12qj+1xuj+1−n/bx+qjxuj−n/bx+12gj+1x+gjx,for0<x<1andnb<j≤n,where ux,tj+1=uj+1x, px,tj+1=pj+1x, qx,tj+1=qj+1x, gx,tj+1=gj+1x, and b=T/τ.Rearranging equation (10), we get(11)−ϵuxxj+1x+p¯j+1xuj+1x=ϵuxxjx−p¯¯jxujx+g¯j+1x+g¯jx,for0<x<1and0≤j≤nb,−ϵuxxj+1x+p¯j+1xuj+1x=ϵuxxjx−p¯¯jxujx+qj+1xuj+1−n/bx+qjxuj−n/bx+gj+1x+gjx,for0<x<1andnp<j≤n,where p¯j+1x=pj+1x+2/k, p¯¯jx=pjx−2/k, and g¯j+1x=gj+1x−qj+1xψbx,tj+1−1.In the operator form, equation (11) can be written as(12)L∗uj+1x=g1x,tj+1,0<x<1and0≤j≤nb,L∗uj+1x=g2x,tj+1,0<x<1andnb<j≤n,with initial and boundary condition(13)ux,0=ψbx,0,u0,t=ψlt,u1,t=ψrt,where(14)L∗uj+1x=−ϵuxxj+1x+p¯j+1xuj+1x,g1x,tj+1=ϵuxxjx−p¯¯jxujx+g¯j+1x+g¯jx,g2x,tj+1=ϵuxxjx−p¯¯jxujx−qj+1xuj+1−n/bx−qjxuj−n/bx+gj+1x+gjx.Lemma 4 (Maximum principle for continuous semidiscretization).
We assumeθj+1x∈C2D¯x. If θj+10≥0, θj+11≥0, and L∗θj+1x≥0∀x∈Dx, then θj+1x≥0∀x∈D¯x.Proof.
Suppose thatx∗∈D¯x such that(15)θj+1x∗=minx∈Dx¯θj+1x<0.
Sinceθj+10≥0, θj+11≥0, then x∗∉0,1 and x∗∈Dx.
The differential operatorL∗θj+1x at the point x∗ value becomes(16)L∗θj+1x∗=−ϵθxxj+1x∗+p¯j+1x∗θj+1x∗,<0,sinceθxxj+1x∗≥0,which contradicts our assumption.
Therefore,θj+1x≥0,∀x∈D¯x.The local truncation error of the semidiscretization method of equation (11) is given by ej+1≡uj+1−u¯j+1, where u¯j+1 is the computed solution of the boundary value problem.Lemma 5 (Estimation of local error).
Suppose that(17)∂v∂tvux,t≤C,x,t∈D¯,0≤v≤2.The local error estimate in the temporal direction is given by(18)ej+1∞=C1k3.Proof.
Using Taylor series expansion, we have(19)ux,tj+1=ux,tj+1/2+k2utx,tj+1/2+k28utt+Ok3,(20)ux,tj=ux,tj+1/2−k2utx,tj+1/2+k28utt+Ok3.
Subtracting equations (19) and (20), we have(21)ux,tj+1−ux,tjk=utx,tj+1/2+Ok2,ux,tj+1−ux,tjk=ϵuxxx,tj+k2−px,tj+k2ux,tj+k2+g1x,tj+k2,0<x<1,0≤j≤np,ux,tj+1−ux,tjk=ϵuxxx,tj+k2−px,tj+k2ux,tj+k2+g2x,tj+k2,0<x<1,np<j≤n.where(22)ux,tj+k2=ux,tj+ux,tj+12+Ok2,gx,tj+k2=gx,tj+gx,tj+12+Ok2.
The solution to the local error is(23)L∗ej+1=Ok3,ej+10=0=ej+11.
Thus, using the maximum principle, the required result is satisfied.Lemma5 at each time step plays significant role to the global error in the temporal discretization which is defined in the following lemma.Lemma 6 (Global error estimate).
The global error estimate atjth time level is given by(24)Ej∞≤C3k2.Proof.
Using Lemma5, the global error estimate at jth time level is given by(25)Ej∞=∑l=1jel∞≤e1∞+e2∞+⋯+ej∞≤Cjkk2,using Lemma5≤CTk2,sincejk≤T=C3k2,where C3 is the positive real constant independent of ϵ and k.Lemma 7 (Stability estimate for semidiscretization).
Letyx be the solution of the continuous problem for equation (11). Then, we have(26)yx≥α−1L∗yx+yxΩ.Proof.
This lemma can be proved using the same approach as Lemma3.


### 3.2. Spatial Discretization with NSFDM

The development of the exact finite difference approach led to the discovery of a nonstandard discrete modeling method. LetDxm be the partition of the spatial domain 0,1 into m subinterval such that(27)xi=ih,fori=0,1,2,…,m,where h=1/m is the step length in spatial direction such that h>ϵ. The resulted boundary value problems (11) are treated using the nonstandard finite difference scheme rules of Mickens in [22] to find the solution as follows:(28)−ϵui−1j+1−2uij+1+ui+1j+1ϕ2i,j+1+p¯ij+1uij+1=ϵui−1j−2uij+ui+1jϕ2i,j−p¯¯ijuij+g¯ij+1+g¯ij,for0<i<mand0≤j≤nb,−ϵui−1j+1−2uij+1+ui+1j+1ϕ2i,j+1+p¯ij+1uij+1=ϵui−1j−2uij+ui+1jϕ2i,j−p¯¯ijuij+qij+1uij+1−nb+qijuij−nb+gij+1+gij,for0<i<mandnb<j≤n,with initial and boundary conditions(29)uxi,0=ψbxi,0,u0,tj=ψltj,u1,tj=ψrtj,where fori=1,2,…,m−1andj=1,2,…,n.(30)ϕ2i,j=4ρ2i,jsinh2ρi,jh2,withρi,j=p¯ijϵ.The scheme (28) results in a tridiagonal system of linear equations(31)Eij+1ui−1j+1+Fij+1uij+1+Gij+1ui+1j+1=Hij+1,where(32)Eij+1=−ϵϕ2i,j+1Fij+1=2ϵϕ2i,j+1+p¯ij+1Gij+1=−ϵϕ2i,j+1∀i,jandHij+1=ϵui−1j−2uij+ui+1jϕ2i,j−p¯¯ijuij+g¯ij+1+g¯ij,∀i,0≤j≤nb,ϵui−1j−2uij+ui+1jϕ2i,j−p¯¯ijuij−qij+1uij+1−nb−qijuij−nb+gij+1+gij,∀i,nb<j≤n.The discrete operator of the scheme (31) satisfies the following maximum principle.Lemma 8.
(Discrete maximum principle). Suppose thatLm,n is the discrete operator of the scheme (31) and ωij be any mesh function satisfying, ωi0≥0,ω0j≥0,ωmj≥0,∀i,j. If Lm,nωij≥0 in Dmn, then ωij≥0 in D¯mn.Proof.
Letωsl=mini,jωij, for all ωij∈D¯mn and ωsl<0. It follows that ωs+1l−ωsl>0 and ωs−1l−ωsl>0. Thus, Lm,nωsl<0, which contradicts our assumptions.
Therefore,ωij≥0∈D¯mn.


## 3.1. Time Discretization

We first perform a semidiscretization of the problem in time direction by using the average of implicit and explicit Euler method with uniform time-stepsk=T/n in the time-interval 0,T. This results a time-mesh Dtn≡tj=jk,j=0,1,2,…,n, and the method results in a system of boundary value problems.(10)uj+1x−ujxk−ϵ2uxxj+1x+uxxjx+12pj+1xuj+1x+pjxujx=−12qj+1xψbx,tj+1−1+qjxψbx,tj−1+12gj+1x+gjx,for0<x<1and0≤j≤nb,uj+1x−ujxk−ϵ2uxxj+1x+uxxjx+12pj+1xuj+1x+pjxujx=−12qj+1xuj+1−n/bx+qjxuj−n/bx+12gj+1x+gjx,for0<x<1andnb<j≤n,where ux,tj+1=uj+1x, px,tj+1=pj+1x, qx,tj+1=qj+1x, gx,tj+1=gj+1x, and b=T/τ.Rearranging equation (10), we get(11)−ϵuxxj+1x+p¯j+1xuj+1x=ϵuxxjx−p¯¯jxujx+g¯j+1x+g¯jx,for0<x<1and0≤j≤nb,−ϵuxxj+1x+p¯j+1xuj+1x=ϵuxxjx−p¯¯jxujx+qj+1xuj+1−n/bx+qjxuj−n/bx+gj+1x+gjx,for0<x<1andnp<j≤n,where p¯j+1x=pj+1x+2/k, p¯¯jx=pjx−2/k, and g¯j+1x=gj+1x−qj+1xψbx,tj+1−1.In the operator form, equation (11) can be written as(12)L∗uj+1x=g1x,tj+1,0<x<1and0≤j≤nb,L∗uj+1x=g2x,tj+1,0<x<1andnb<j≤n,with initial and boundary condition(13)ux,0=ψbx,0,u0,t=ψlt,u1,t=ψrt,where(14)L∗uj+1x=−ϵuxxj+1x+p¯j+1xuj+1x,g1x,tj+1=ϵuxxjx−p¯¯jxujx+g¯j+1x+g¯jx,g2x,tj+1=ϵuxxjx−p¯¯jxujx−qj+1xuj+1−n/bx−qjxuj−n/bx+gj+1x+gjx.Lemma 4 (Maximum principle for continuous semidiscretization).
We assumeθj+1x∈C2D¯x. If θj+10≥0, θj+11≥0, and L∗θj+1x≥0∀x∈Dx, then θj+1x≥0∀x∈D¯x.Proof.
Suppose thatx∗∈D¯x such that(15)θj+1x∗=minx∈Dx¯θj+1x<0.
Sinceθj+10≥0, θj+11≥0, then x∗∉0,1 and x∗∈Dx.
The differential operatorL∗θj+1x at the point x∗ value becomes(16)L∗θj+1x∗=−ϵθxxj+1x∗+p¯j+1x∗θj+1x∗,<0,sinceθxxj+1x∗≥0,which contradicts our assumption.
Therefore,θj+1x≥0,∀x∈D¯x.The local truncation error of the semidiscretization method of equation (11) is given by ej+1≡uj+1−u¯j+1, where u¯j+1 is the computed solution of the boundary value problem.Lemma 5 (Estimation of local error).
Suppose that(17)∂v∂tvux,t≤C,x,t∈D¯,0≤v≤2.The local error estimate in the temporal direction is given by(18)ej+1∞=C1k3.Proof.
Using Taylor series expansion, we have(19)ux,tj+1=ux,tj+1/2+k2utx,tj+1/2+k28utt+Ok3,(20)ux,tj=ux,tj+1/2−k2utx,tj+1/2+k28utt+Ok3.
Subtracting equations (19) and (20), we have(21)ux,tj+1−ux,tjk=utx,tj+1/2+Ok2,ux,tj+1−ux,tjk=ϵuxxx,tj+k2−px,tj+k2ux,tj+k2+g1x,tj+k2,0<x<1,0≤j≤np,ux,tj+1−ux,tjk=ϵuxxx,tj+k2−px,tj+k2ux,tj+k2+g2x,tj+k2,0<x<1,np<j≤n.where(22)ux,tj+k2=ux,tj+ux,tj+12+Ok2,gx,tj+k2=gx,tj+gx,tj+12+Ok2.
The solution to the local error is(23)L∗ej+1=Ok3,ej+10=0=ej+11.
Thus, using the maximum principle, the required result is satisfied.Lemma5 at each time step plays significant role to the global error in the temporal discretization which is defined in the following lemma.Lemma 6 (Global error estimate).
The global error estimate atjth time level is given by(24)Ej∞≤C3k2.Proof.
Using Lemma5, the global error estimate at jth time level is given by(25)Ej∞=∑l=1jel∞≤e1∞+e2∞+⋯+ej∞≤Cjkk2,using Lemma5≤CTk2,sincejk≤T=C3k2,where C3 is the positive real constant independent of ϵ and k.Lemma 7 (Stability estimate for semidiscretization).
Letyx be the solution of the continuous problem for equation (11). Then, we have(26)yx≥α−1L∗yx+yxΩ.Proof.
This lemma can be proved using the same approach as Lemma3.


## 3.2. Spatial Discretization with NSFDM

The development of the exact finite difference approach led to the discovery of a nonstandard discrete modeling method. LetDxm be the partition of the spatial domain 0,1 into m subinterval such that(27)xi=ih,fori=0,1,2,…,m,where h=1/m is the step length in spatial direction such that h>ϵ. The resulted boundary value problems (11) are treated using the nonstandard finite difference scheme rules of Mickens in [22] to find the solution as follows:(28)−ϵui−1j+1−2uij+1+ui+1j+1ϕ2i,j+1+p¯ij+1uij+1=ϵui−1j−2uij+ui+1jϕ2i,j−p¯¯ijuij+g¯ij+1+g¯ij,for0<i<mand0≤j≤nb,−ϵui−1j+1−2uij+1+ui+1j+1ϕ2i,j+1+p¯ij+1uij+1=ϵui−1j−2uij+ui+1jϕ2i,j−p¯¯ijuij+qij+1uij+1−nb+qijuij−nb+gij+1+gij,for0<i<mandnb<j≤n,with initial and boundary conditions(29)uxi,0=ψbxi,0,u0,tj=ψltj,u1,tj=ψrtj,where fori=1,2,…,m−1andj=1,2,…,n.(30)ϕ2i,j=4ρ2i,jsinh2ρi,jh2,withρi,j=p¯ijϵ.The scheme (28) results in a tridiagonal system of linear equations(31)Eij+1ui−1j+1+Fij+1uij+1+Gij+1ui+1j+1=Hij+1,where(32)Eij+1=−ϵϕ2i,j+1Fij+1=2ϵϕ2i,j+1+p¯ij+1Gij+1=−ϵϕ2i,j+1∀i,jandHij+1=ϵui−1j−2uij+ui+1jϕ2i,j−p¯¯ijuij+g¯ij+1+g¯ij,∀i,0≤j≤nb,ϵui−1j−2uij+ui+1jϕ2i,j−p¯¯ijuij−qij+1uij+1−nb−qijuij−nb+gij+1+gij,∀i,nb<j≤n.The discrete operator of the scheme (31) satisfies the following maximum principle.Lemma 8.
(Discrete maximum principle). Suppose thatLm,n is the discrete operator of the scheme (31) and ωij be any mesh function satisfying, ωi0≥0,ω0j≥0,ωmj≥0,∀i,j. If Lm,nωij≥0 in Dmn, then ωij≥0 in D¯mn.Proof.
Letωsl=mini,jωij, for all ωij∈D¯mn and ωsl<0. It follows that ωs+1l−ωsl>0 and ωs−1l−ωsl>0. Thus, Lm,nωsl<0, which contradicts our assumptions.
Therefore,ωij≥0∈D¯mn.


## 4. Convergence Analysis

We shall omit the time level index for the sake of simplicity. The following is how we get the local truncation error of the present nonstandard finite difference method from equation (28):(33)Lϵm,nUi−ui=Lx,ϵ−Lx,ϵmui,using nonstandard finite difference method(34)Lϵm,nUi−ui=−ϵui″+ϵui−1−2ui+ui+1ϕi2.Using Taylor series expansion ofui−1 and ui+1 in (34), we have that(35)Lϵm,nUi−ui=−ϵui″+ϵϕi2h2ui″+h412u4ξi,ξi∈xi−1,xi+1.The truncated Taylor series expansion of1/ϕi2 is given by(36)1ϕi2=ρi244ρi2h2−13+ρi2h260.Using equation (36), equation (34) becomes(37)Lϵm,nUi−ui=ϵ12u4ξi−ρi2ui″h2+ϵρi2ρi2240ui″−1144u4ξih4+ϵρi42880u4ξih6.Thus,(38)Lϵm,nUi−ui≤Mh2,where M is independent to ϵ.Theorem 1.
LetUij be the numerical solution of (31) and ux,tj be the solution of (11) both at time level j. Then,(39)max0≤i≤mUij−uij≤Mh2.Therefore, the nonstandard finite difference method is a second order uniformly convergent scheme in space variable.Theorem 2.
LetUij be the numerical solution of (31) and ux,t be the solution of (1). Then,(40)max0≤i≤m,0≤j≤nUij−uij≤Mh2+k2.Therefore, the nonstandard finite difference method is a second order uniformly convergent scheme in both space and time variables.


## 5. Richardson Extrapolation Method

The Richardson extrapolation method is intended to increase the accuracy of the basic scheme’s computed solutions. From Theorem2, we have(41)uxi,tj+1−Uij+1≤Mh2+k2,where uxi,tj+1 and Uij+1 are the exact and approximate solution for Dmn mesh interval, respectively, and M is a constant independent of ϵ, k, and h.LetD2m2n be the mesh found by doubling the mesh interval of Dmn and approximate solution using D2m2n mesh interval is U¯ij+1. Equation (41) with mesh size h,k≠0 becomes(42)uxi,tj+1−Uij+1=Mh2+k2+Rmn,xi,tj+1∈Dmn,and equation (41) with mesh size h/2,k/2≠0 becomes(43)uxi,tj+1−U¯ij+1=Mh24+k24+R2m2n,xi,tj+1∈D2m2n,where the remainders Rmn and R2m2n are Oh4+k4. Subtracting equation (42) from equation (43) and obtained extrapolation formula,(44)Uij+1ext=4U¯ij+1−Uij+13,where Uij+1ext is the approximate solution of uxi,tj+1 using the Richardson extrapolation method with truncation error.(45)uxi,tj+1−Uij+1ext≤Mh4+k4.Theorem 3.
LetUij+1ext be the numerical solution of (31) after the Richardson extrapolation method applied and ux,t be the solution of (1). Then,(46)max0≤i≤m,0≤j<nuij−Uij+1ext≤Mh4+k4.Therefore, after the Richardson extrapolation method applied, the proposed method becomes fourth order uniformly convergent scheme in both space and time variables.


## 6. Numerical Examples and Results

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 is unknown [23]. We use the following formula to approximate the absolute maximum error at the selected mesh points:Case 1.
If the exact solution is known,(47)Eϵm,n=maxxi,tj∈Ωuxi,tj−uijext.Case 2.
If the exact solution is unknown,(48)Eϵm,n=maxxi,tj∈Ωuij−u2i−12j−1ext.where uij is the approximation solution obtained by using m and n number of subintervals of space and time variables, respectively, and u2i−12j−1ext is the approximation solution obtained by using 2m and 2n number of subintervals of space and time variables, respectively, after applying the Richardson extrapolation method.
We also evaluate the corresponding rate of convergence:(49)Rϵm,n=logEϵm,n−logEϵ2m,2nlog2.Example 1.
We consider the following singularly perturbed parabolic delay initial-boundary-value problem:(50)utx,t−ϵuxxx,t+1.1+x2−ux,t−1=t3,x,t∈0,1×0,2,ux,t=0,x,t∈0,1×−1,0,u0,t=0=u1,t,t∈0,2,whose exact solution is unknown.Example 2.
We consider the following singularly perturbed parabolic delay initial-boundary-value problem in [24]:(51)utx,t−ϵuxxx,t+1+x22=−ux,t−1+t3,x,t∈0,1×0,2,ux,t=0,x,t∈0,1×−1,0,u0,t=0=u1,t,t∈0,2,whose exact solution is unknown.


## 7. Discussion and Conclusion

We proposed a nonstandard finite difference method for solving singularly perturbed parabolic delay partial differential equations of the reaction-diffusion type. The considered problems’ solution exhibits twin boundary layers on the right and left sides of the spatial domain. The basic mathematical procedures are defining the model problem, approximating the time variable using the Crank–Nicolson method, approximating the spatial variable using nonstandard finite difference scheme using Mickens rule, and reducing it into a tridiagonal systems of equation; this can be solved by using the Thomas algorithm. Lastly, to enhance the accuracy of the method, we use the Richardson extrapolation technique.To validate the applicability of the proposed method, we used two model examples. The maximum absolute error and rate of convergence are shown in Tables1 and 2 with different values of ϵ,m,andn. Comparison of maximum absolute error and rate of convergence with other methods that appear in the literature are shown in Tables 3 and 4. The behavior of maximum absolute error is shown in Figure 1. The physical behavior of the solution is shown in Figures 2 and 3.Table 1
Maximum absolute pointwise error and rate of convergence for Example1 before and after the Richardson extrapolation method applied, where L≥6.
ε↓m, n⟶10, 1020, 2040, 4080, 80160, 160After10−65.7318e − 063.3273e − 072.0558e − 081.2776e − 097.9666e − 1110−85.7318e − 063.3273e − 072.0558e − 081.2776e − 097.9666e − 1110−105.7318e − 063.3273e − 072.0558e − 081.2776e − 097.9666e − 1110−125.7318e − 063.3273e − 072.0558e − 081.2776e − 097.9666e − 1110−145.7318e − 063.3273e − 072.0558e − 081.2776e − 097.9666e − 11Rate of convergence10−L4.10664.01664.00814.0034Before10−68.3948e − 032.3134e − 036.0360e − 041.5393e − 043.8852e − 0510−88.3948e − 032.3134e − 036.0360e − 041.5393e − 043.8852e − 0510−108.3948e − 032.3134e − 036.0360e − 041.5393e − 043.8852e − 0510−128.3948e − 032.3134e − 036.0360e − 041.5393e − 043.8852e − 0510−148.3948e − 032.3134e − 036.0360e − 041.5393e − 043.8852e − 05Rate of convergence10−L1.85951.93831.97131.9862Table 2
Maximum absolute pointwise error and rate of convergence for Example2 before and after the Richardson extrapolation method applied, where L≥6.
ε↓m, n⟶16, 1632, 3264, 64128, 128256, 256After10−61.5252e − 065.6891e − 082.8561e − 091.5286e − 109.4897e − 1210−81.5252e − 065.6891e − 082.8561e − 091.5808e − 109.4897e − 1210−101.5252e − 065.6891e − 082.8561e − 091.5808e − 109.4897e − 1210−121.5252e − 065.6891e − 082.8561e − 091.5808e − 109.4897e − 1210−141.5252e − 065.6891e − 082.8561e − 091.5808e − 109.4897e − 12Rate of convergence10−L4.74464.31614.22374.0581Before10−64.7003e − 031.1960e − 033.0108e − 047.5494e − 051.8899e − 0510−84.7003e − 031.1960e − 033.0108e − 047.5494e − 051.8899e − 0510−104.7003e − 031.1960e − 033.0108e − 047.5494e − 051.8899e − 0510−124.7003e − 031.1960e − 033.0108e − 047.5494e − 051.8899e − 0510−144.7003e − 031.1960e − 033.0108e − 047.5494e − 051.8899e − 05Rate of convergence10−L1.97461.93831.99571.9980Table 3
Comparison of absolute maximum errors and rate of convergenceRm,n for Example 1.
ε↓m, n⟶32, 864, 16128, 32256, 64Present method10−61.3999e − 058.6123e − 075.0742e − 083.1249e − 0910−71.3999e − 058.6123e − 075.0742e − 083.1249e − 0910−81.3999e − 058.6123e − 075.0742e − 083.1249e − 09Rm,n4.02284.08514.0213Singh et al. [12]10−67.95e − 012.03e − 015.12e − 021.28e − 0210−77.95e − 012.03e − 015.12e − 021.28e − 0210−87.95e − 012.03e − 015.12e − 021.28e − 02Rm,n1.971.992.00Table 4
Comparison of absolute maximum errors and rate of convergenceRm,n for Example 2.
ε↓m, n⟶32, 3264, 64128, 128256, 256Present method10−65.6891e − 082.8561e − 091.5286e − 109.4897e − 1210−75.6891e − 082.8561e − 091.5286e − 109.4897e − 1210−85.6891e − 082.8561e − 091.5286e − 109.4897e − 12Rm,n4.31614.22374.0581Govindarao and Mohapatra [21]10−62.6347e − 032.4293e − 041.6525e − 051.0538e − 0610−72.6347e − 032.4293e − 041.6525e − 051.0538e − 0610−82.6347e − 032.4293e − 041.6525e − 051.0538e − 06Rm,n3.43903.87713.9710Figure 1
Log-log plot of maximum absolute errors before and after the Richardson extrapolation method applied for Example1.Figure 2
The physical behavior of the solutions for Example2 for m = 20, n = 20, and ε=10−10.Figure 3
The behavior of spatial vs solution graph with different time level for Example2 for m = 20, n = 20, and ε=10−10.The convergence of proposed scheme is insensitive to the perturbation parameter. After applying the Richardson extrapolation technique, the present method has order of four rate of convergence. The performance of the proposed scheme is assessed by comparing the results, and it is demonstrated that the numerical results are more accurate than that of existing difference schemes in the literature.As future directions of this research, we will apply the scheme for nonlinear singularly perturbed parabolic reaction and convection diffusion problems.



---

*Source: 5625049-2022-07-30.xml*