# A Parameter-Uniform Numerical Scheme for Solving Singularly Perturbed Parabolic Reaction-Diffusion Problems with Delay in the Spatial Variable

Authors: Ababi Hailu Ejere; Gemechis File Duressa; Mesfin Mekuria Woldaregay; Tekle Gemechu Dinka
Journal: International Journal of Mathematics and Mathematical Sciences
 (2023)

---

## Abstract

The objective of this research work is to develop and analyse a numerical scheme for solving singularly perturbed parabolic reaction-diffusion problems with large spatial delay. The presence of the small positive parameter on the term with the highest order of derivative exhibits two strong boundary layers in the solution of the problem, and the large delay term gives rise to a strong interior layer. The layers’ behavior makes it difficult to solve the problem analytically. To treat such a problem, we developed a numerical scheme using the Crank–Nicolson method in the time direction and the central difference method in the spatial direction via nonstandard finite difference methods on uniform meshes. Stability and convergence analyses for the obtained scheme have been established, which show that the developed numerical scheme is uniformly convergent. To confirm the theoretical analysis, model numerical examples are considered and demonstrated.


---

## Body

## 1. Introduction

A delay differential equation is an equation in which the evolution of the dependent variable depends on a certain time retarded to a previous state [1]. Different model problems in science and engineering are described using delay differential equations. For instance, in the modeling of biological mathematics [2], problems in control theory [3], in modeling finite-time event-triggered synchronization for stochastic neural networks along with the passivity and passification conditions [4], human pupil light reflex [5], HIV infections [6], in modeling the upper bound for bidirectional associative memory neural networks under global asymptotic robust stability condition [7], COVID-19 pandemic [8], in modeling sampled-data state estimation for delayed Markovian jamp neural networks based on passive theory [9] and others.A class of differential equation in which a small positive parameterε (often 0<ε≪1) is multiplied by the term with the highest order derivative and consists of one or more negatively shifted terms is said to be a singularly perturbed delay differential equation. The solution involves two boundary layers because of the influence of ε and interior layer, which is caused by the large delay in the spatial variable. The layer behavior brings a challenging factor to obtain an analytical solution to the problem. Also, classical numerical methods do not give satisfactory numerical solutions as they do not consider the boundary or interior layers into account unless extremely large mesh points are considered, which leads to a massive computational cost, or it may fail to give an accurate result.For the response to the previously mentioned challenges in solving singularly perturbed differential equations with delay parameters, different research works have been reported over the last few decades. For instance, Cimen and Amiraliyev [10] solved singularly perturbed delay differential equation with layer behavior by constructing a finite difference scheme using the fitted mesh method and obtained an essentially first-order parameter-uniformly convergent result. Duressa [11] treated a singularly perturbed problem by constructing a numerical scheme using a fitting parameter introduced in a finite difference approximation. Kumar et al. [12] solved singularly perturbed problems involving a term with a large negative shift by developing an exponentially fitted numerical scheme using Numerov’s method. Debela and Garoma [13] solved a singularly perturbed problems involving large and small delays by introducing a nonpolynomial spline approach. In reference [14], the singularly perturbed problem with small shift was treated using the average weight method in the time direction and fitted operator method in the spatial direction. In reference [15], a singularly perturbed delay differential equation is treated by decomposing the domain into subdomains and obtaining a convergent numerical result. In reference [16], a singularly perturbed Fredholm integrodifferential equation is considered. The authors constructed a new difference scheme by the method of integral identities using interpolating quadrature rules and obtained a uniformly convergent numerical method in the maximum norm.Singularly perturbed problems involving large delays have been also treated by few researchers. For instance, in reference [17], singularly perturbed delay differential equation subjected to integral boundary conditions have been treated using a fitted mesh method. In reference [18], a singularly perturbed parabolic problem with a large delay was treated by proposing a scheme applying the Crank–Nicolson method on a uniform time mesh and the central difference method on a piecewise uniform spatial mesh. Bansal and Sharma [19] and Swaminathan et al. [20] solved singularly perturbed parabolic problem involving significant delay in the spatial variable by constructing a scheme applying implicit Euler techniques on uniform time grids and central difference method on piecewise uniform spatial grids. In [21], a uniformly convergent numerical scheme has been developed using the fitted mesh method for a time-dependent singularly perturbed delay differential equation.In this study, we are motivated to treat a time-dependent singularly perturbed parabolic reaction-diffusion problem involving a unit delay in the spatial variable. The influence of the perturbation parameter and the delay in the reaction term cause difficulties in solving the problem analytically. To overcome such challenges, we developed a more accurate numerical scheme on uniformly discretized time and spatial domains using nonstandard finite difference methods. The convergence analysis is also established for the proposed numerical scheme, from which we demonstrate the uniform convergence of the method.The other sections of the paper are organized as follows: Section2 deals with the description of the analytical results. The numerical method of the study is briefly discussed in Section 3. To demonstrate the theoretical findings, we present numerical examples and discussion in Section 4. The conclusion of the study is presented in Section 5.Notation: we usedC as a generic constant that does not depends on the perturbation parameter and mesh elements. N and M represent the number of mesh sizes in the spatial and temporal directions, respectively. For a given continuous function vx, we used vℜ¯ to denote the supreme norm of the function.


## 2. Continuous Problem: Preliminary and Priori Estimates

We considered a singularly perturbed differential equation onℜ¯=0,2×0,T as follows:(1)Lεy≡∂y∂t−ε∂2y∂x2+axyx,t+bxyx−1,t=gx,t,subjected to the initial and interval-boundary conditions(2a)yx,0=y0x,x∈Γ¯,(2b)yx,t=ϕx,t,x,t∈ℜ−,(2c)y2,t=γt,2,t∈ℜ+,where 0<ε≪1, Γ¯=0,2, ℜ−=x,t:x∈−1,0;t∈0,T, and ℜ+=2,t:t∈0,T for a finite time T and Lε is a differential operator independent of ε. We assume the functions involved in the continuous problem to be sufficiently smooth. Moreover, for an arbitrary positive constant λ, the coefficient functions are considered to satisfy(3)ax+bx≥2λ>0andbx<0,x∈Γ¯.Depending on the interval condition, equation (1) can be equivalently written as follows:(4)Lεy≡∂y∂t−ε∂2y∂x2+axyx,t=gx,t−bxϕx−1.t,x,t∈0,1×0,T,∂y∂t−ε∂2y∂x2+axyx,t+bxyx−1,t=gx,t,x,t∈1,2×0,1,with the conditions given in equations (2a)–(2c). Setting ε=0, (4) is reduced to(5)Lεy0≡∂y0∂t+axy0x,t=gx,t−bxϕx−1.t,x,t∈0,1×0,T,∂y0∂t+axy0x,t+bxy0x−1,t=gx,t,x,t∈1,2×0,T,in the reduced problem (5), since y0x,t need not necessarily satisfy the conditions y00,t=ϕ0,t, y02,t=γt, y01+,t=y01−,t and ∂y01+,t/∂x=∂y01−,t/∂x, in general the solution yx,t exhibits two strong boundary layers at the two ends of the domain and the strong interior layer at x=1 [22]. Also, we impose compatibility conditions as follows:(6)y00,0=ϕ0,0,y02,0=γ2,0,∂ϕ0,0∂t−ε∂2y00∂x2+a0y00=g0,0−b0ϕ−1.0,∂γ2,0∂t−ε∂2y02∂x2+a2y22+b2y01,0=g2,0,which show that the considered data held at the corner points (0, 0) and (2, 0). From the previous assumptions, the existence of C independent ε is guaranteed and the problem has a unique solution [23]. Thus, ∀x,t∈ℜ¯ by [24], we can obtain(7)yx,t−yx,0≤Ct.Ifε small enough, the solution yx,t of (1)–(2c) closes to y0x,t. Also, without loss of generality assuming that all the data given in (2a)–(2c) are identically zero [25], so that the following properties of yx,t are true.Lemma 1 (Bound of the continuous solution).
The solutionyx,t of (1)–(2c) can be bounded as yx,t≤C for all x,t∈ℜ¯.Proof.
From (7), we can obtain yx,t−y0x≤yx,t−y0x≤Ct, which gives yx,t≤Ct+y0x, ∀x,t∈ℜ¯. Thus, as y0x is bounded and t∈0,T, we obtain yx,t≤C for all x,t∈ℜ¯.Lemma 2 (Continuous maximum principle).
Assume thatωx,t is any continuous function in ℜ¯. If ωx,t≥0,x,t∈∂ℜ and Lεωx,t≥0,x,t∈ℜ, then ωx,t≥0,x,t∈ℜ¯.Proof.
Letξ,ι∈ℜ¯ such that ωξ,ι=minℜ¯ωx,t and let ωξ,ι<0. By the given hypothesis, we have ξ,ι∉∂ℜ and ωxξ,ι=0, ωtξ,ι=0, ωxxξ,ι≥0.Case 1: on (0, 1], we obtainLε,1ωξ,ι=∂ω/∂tξ,ι−ε∂2ω/∂x2ι,ι+aξωξ,ι=−ε∂2ω/∂x2ξ,ι+aξωξ,ι≤0.Case 2: on (1, 2), we obtainLε,2ωξ,ι=∂ω/∂tξ,ι−ε∂2ω/∂x2ξ,ι+aξωξ,ι+bξωξ−1,ι≤0.
From the two cases, we obtained a contradiction to the hypothesis. Thus, our assumption fails andωξ,ι≥0, so that ωx,t≥0,x,t∈ℜ¯.Lemma 3 (Continuous stability estimate).
Letyx,t be the solution of (1). Then, it can be bounded as follows:(8)yx,t≤maxy∂ℜ,1λLεy.Proof.
For the proof, we define barrier functions asω±x,t=maxy∂ℜ,1/λLεy±yx,t. Then, we haveForx=0,ω±0,t=maxy∂ℜ,1/λLεy±y0,t≥1/λLεy±y0,t≥0.Forx=2,ω±2,t=maxy∂ℜ,1/λLεy±y2,t≥1/λLεy±y2,t≥0.Forx∈ (0.1] ,Lεω±x,t=ωt±−εωxx±+axω±x,t+bxϕ±x−1,t≥axmaxy∂ℜ,1/λg±gx,t≥0,x,t∈ℜ.Forx∈1,2,Lεω±x,t=ωt±−εωxx±+axω±x,t+bxy±x−1,t≥2λmaxy∂ℜ,1/λg±gx,t≥0,x,t∈ℜ.
Thus, applying Lemma2, we obtain ω±x,t≥0 for all x,t∈ℜ¯, which implies the required stability estimate.Lemma 4 (Bound on the derivatives of the solution).
The derivatives of the solutionyx,t of (1) can be bounded as follows:(9)∂i+jyx,t∂xi∂tj≤C1+ε−i/2exp−λεx+exp−λε1−xforx,t∈ℜ¯−,C1+ε−i/2exp−λεx−1+exp−λε2−xforx,t∈ℜ¯+,for all non-negative integer i,j such that i+j∈0,4.Proof.
To determine the bound on the derivatives of the solution for the continuous problem, we follow the procedures in references [26, 27]. For x,t∈ℜ¯, consider the case when i=0, which implies differentiating yx,t with respect to time. For j=0, it implies Lemma 1. For j=1, rearranging (1) gives(10)∂y∂t=ε∂2y∂x2−axyx,t−bxyx−1,t+gx,t.
Without loss of generality assuming thatyx,0=y0=0, yx,t=ϕx,t=0, y2,t=γt=0 (10), becomes(11)∂yx,0∂t=ε∂2yx,0∂x2−axyx,0−bxyx−1,0+gx,0.
Thus, alongx=0 and x=2, we can obtain yx,t=0 which gives ∂y/∂t=0 and along t=0, we obtain yx,0=0, so that ∂2y/∂x2=0,t=0. On the interval [0, 1], we have yx−1,0=ϕx−1,0=0 and on the interval [1, 2], we have yx−1,0=y0x−1,0=0. Combining these results give yx−1,0=0, so that equation (11) becomes ∂y/∂tx,0=gx,0 and since g is a smooth function, we have yx,0≤C on ∂ℜ, which implies that(12)∂y∂t≤Conℜ¯.
In a similar way forj=2, assuming that all the data given in (2a)–(2c) are identically zero and differentiating with respect to t gives ∂2y/∂t2−ε∂3y/∂x2∂t+ax∂y/∂t+bx∂yx−1,t/∂t=gtx,t. and along t=0, it becomes(13)∂2yx,0∂t2−ε∂3yx,0∂x2∂t+ax∂yx,0∂t+bx∂yx−1,0∂t=gtx,0.
However, from∂yx,0/∂t=gx,0, it follows that ∂3yx,0/∂2x∂t=gxxx,0 and ∂yx−1,0/∂t=0, so that (13) becomes ∂2yx,0/∂t2=εgxxx,0−axgx,0+gtx,0. Therefore, ∂2yx,0/∂t2≤C on ∂ℜ and since g is a smooth function, we have ∂2yx,t/∂t2≤C. Thus, for large value of C, we obtain ∂2yx,t/∂t2≤C, x,t∈ℜ¯.
On the other hand, to obtain the bound of derivatives of the solution with respect to the spatial variable, we consider fori=0,1,2,3,4. For i=0, it is the immediate result of Lemma 1. For i=1, let x∈Γ and consider a neighborhood I=ι,ι+ε, ∀x∈I. Then, for some ι in I¯=ι,ι+ε and t∈0,T, by the mean value theorem, we obtain(14)∂yι,t∂x=yι+ε,t−yι,tε=1εyι+ε,t−yι,t≤Cε−1/2y.
Forx∈I¯, we have(15)∂y∂x=∂yι,t∂x+ε−1∫ιx∂ys,t∂t+asys,t+bsys−1,t−gs,tds,which implies that(16)∂y∂x≤∂yι,t∂x+Cε−1∫ιxy+gds=∂yι,t∂x+Cε−1y+gε1/2.
Using equation (14) in the previous result yields ∂y/∂xx,t≤Cε−1/2, ∀x,t∈ℜ¯. Rearranging the terms in equation (1) gives(17)∂2y∂x2=ε−1∂y∂t+axyx,t+bxyx−1,t−gx,t.
For fixedt in 0,T and for x∈0,2, since y≤C, ∂y/∂t≤C and g is a smooth function, we can obtain ∂2y/∂x2x,t≤Cε−1. If we differentiate both sides of (15) with respect to t, then it becomes(18)∂3y∂x2∂t=ε−1∂2y∂t2+ax∂y∂t+bx∂yx−1,t∂t−gtx,t.
Since∂y/∂t≤C, ∂2y/∂t2≤C and gtx,t≤C, we have ∂3yx,t/∂x2∂t≤Cε−1. By a similar procedure for the remaining values of i and j satisfying 0≤i+2j≤4, the derivatives of the solution can be bounded. For the derivative of the solution with Shishkin decomposition, we refer [28].


## 3. Numerical Method

### 3.1. Temporal Semidiscretization

Let0,T be a given time domain for some finite time T, which can be subdivided into M equal intervals forming uniform temporal mesh as ℜtM=tj=jΔt,j+1=1,2,…,M−1,T=MΔt. From the rules upon which nonstandard schemes are based [29], the time derivative in equation (1) takes the form dy/dt⟶yj+1-yj/µΔt, where μΔt is the denominator function satisfying μΔt=Δt+OΔt2. Thus, for the subequation dy/dt=−λy, where ax≥λ>0, the standard forward Euler discretization is yj+1-yj/Δt=−λyj, which can be written as follows:(19)yj+1=1−λΔtyj.However, it is well known that equation (19) leads to instability if λΔt>1 [30]. Instead, from the fact that 1−λΔt=exp−λΔt+OλΔt2 for sufficiently small λΔt (i.e., 0<λΔt≪1), it is possible to make the replacement 1−λΔt≈exp−λΔt, so that Δt≈1-exp−λΔt/λ. Adopting this for the variable coefficient, we have a denominator function(20)μΔt=1−exp−axΔtax.Now, applying the Crank–Nicolson method in the temporal direction incorporating the denominator function (20), we obtain a semidiscrete scheme as follows:(21)LεMyj+1x≡yj+1x−yjxμΔt−ε2∂2yj+1∂x2+∂2yj∂x2+ax2yj+1x+yjx+bx2yj+1x−1+yjx−1=12gj+1x+gjx,and rearranging terms, we can rewrite as follows:(22)LεMyj+1x≡−ε∂2yj+1∂x2+2μΔt+axyj+1x=ε∂2yj∂x2+2μΔt−axyjx−2bxϕj+1/2x−1+2gj+1/2x,forx∈0,1,−ε∂2yj+1∂x2+2μΔt+axyj+1x+2bxyj+1x−1=ε∂2yj∂x2+2μΔt−axyjx+2gj+1/2x,forx∈1,2,subjected to yj+1x,0=y0j+1x,x∈Γ, yj+1x=ϕx,tj+1 and yj+12=γ2,tj+1 for j+1=1,2,…,M−1, where yj+1x is the approximate of yx,tj+1 at j+1th level. We can rewrite equation (22) as follows:(23)2μΔt+LεMyj+1x=ϑj+1x,x∈Γ,where(24)LεMyj+1x=−ε∂2yj+1∂x2+axyj+1x,x∈0,1,−ε∂2yj+1∂x2+axyj+1x+bxyj+1x−1,x∈1,2,ϑj+1x=ε∂2yj∂x2−axyjx−2bxϕj+1/2x−1+2gj+1/2x+2μΔt,x∈0,1,ε∂2yj∂x2−axyjx−bxyjx−1+2gj+1/2x+2μΔt,x∈1,2,where ϕj+1/2x−1=1/2ϕj+1x−1+ϕjx−1 and gj+1/2x=1/2gj+1x+gjx.Lemma 5 (Semidiscrete maximum principle).
For a sufficiently smooth functionωj+1 in ℜ¯, if ωj+10≥0, ωj+12≥0 and Lεωj+1x≥0 for all x∈Γ, then ωj+1x≥0, ∀x∈Γ¯.Proof.
Letr∈Γ¯ such that ωj+1r=minΓ¯ωj+1x and suppose ωj+1r<0. From the given conditions, we have r∉∂Γ and ωxj+1r=0, ωxxj+1r≥0.Case 1: forr∈0,1, we have(25)Lε,1Mωj+1r=−ε∂2ωj+1r∂x2+2μΔt+arωj+1r≤0.Case 2: forr∈1,2, we have(26)Lε,2Mωj+1r=−ε∂2ωj+1r∂x2+2μΔt+arωj+1r+brωj+1r−1≤−ε∂2ωj+1∂x2r+2μΔt+ar+brωj+1r≤0.
From the two cases, we observe that the given hypothesis is contradicted and hence our assumption fails, so thatωj+1x≥0, ∀x∈Γ¯. Thus, the maximum principle is satisfied by the operator LεM. Also, for the semidiscrete operator, we have(27)Lε,xM−1≤11+2λΔt,which is used to estimate truncation errors in the semidiscrete scheme.Lemma 6 (Stability estimate of the semidiscrete problem).
Let the solution of equation (22) be yj+1x. Then. it can be estimated as yj+1x≤ϑj+1/1+2λΔt+maxyj+10,yj+12,∀x∈Γ¯.Proof.
The proof can be easily shown following similar procedure as Lemma3.Definition 1.
At thej+1th time level, we can define local truncation error as the difference between the solution yx,tj+1 and the numerical solution yj+1x of the semidiscrete problem (23) and given ej+1x=yx,tj+1−yj+1x and the global error Ej+1 as the contribution of local truncation error up to the j+1th time level.Lemma 7 (Estimate of the local truncation error).
Suppose it is given that∂kyx,t/∂tk≤C for all x,t∈ℜ¯ and k=0,1,2. Then, at j+1th time level, we have ej+1≤CΔt3.Proof.
Let us choose the center pointx,tj+1/2 of the temporal mesh as a convenient base point where the Taylor series is expanded. Then, yx,tj+1 and yx,tj can be expanded as follows:(28a)yj+1x=yj+1/2x+Δt2∂yj+1/2∂t+Δt28∂yj+1/2∂t2+OΔt3,(28b)yjx=yj+1/2x−Δt2∂yj+1/2∂t+Δt28∂yj+1/2∂t2+OΔt3.
Now, subtracting (28b) from (28a) gives(29)yj+1x−yjxμΔt=∂yj+1/2∂t+OΔt3.
Using (1) into (29) at x,tj+1/2, we obtain(30)yj+1x−yjxμΔt=ε∂2yj+1/2∂x2−axyj+1/2x−bxyj+1/2x−1+gx,tj+1/2+OΔt3.
Putting (30) into the time semidiscrete problem (23) gives(31)ε∂2yj+1/2∂x2−axyj+1/2x−bxyj+1/2x−1+gx,tj+1/2+OΔt3−ε2∂2yj+1∂x2+∂2yj∂x2+ax2yj+1x+yjx+bx2yj+1x−1+yjx−1=12gj+1x+gjx,which can be written equivalently as follows:(32)LεMyx,tj+1−yj+1x=OΔt3.
Using (27) into (32), we obtain a boundary value problem of the type LMej+1x=OΔt3 with ej+10,tj+1=0 and ej+12=0. Applying Lemma 5, we can obtain an estimate for the local truncation error as ej+1≤CΔt3.Lemma 8 (Estimation of the global error).
Assuming that Lemma7 holds, the global error in the time direction can be estimated as Ej+1≤CΔt2, ∀j such that jΔt≤T.Proof.
Using Lemma7 up to the j+1th level, we have(33)Ej+1=∑k=1jek,jΔt≤T≤e1+⋯+ej≤CΔt2.
Thus,Ej+1≤CΔt2, where C is independent of ε and the mesh numbers.Lemma 9 (Bound of the derivatives of the solution).
Let the solution of (23) be yj+1x. Then its derivatives can be bounded as follows:(34)dnyj+1dxn≤C1+ε−n/2exp−λεx+exp−λε1−x,x∈Γ¯,n=014,C1+ε−n/2exp−λεx−1+exp−λε2−x,x∈Γ¯,n=014.Proof.
The proof can be carried out following the procedures given in Theorem 3 of [26].


### 3.2. Fully Discretized Scheme

Let [0, 2] be subdivided intoN equal steps of length h and form a mesh ℜxN=0=x0,x1,…,xN/2=1,xN/2+1,…,xN=2,xi=ih,i=0,1,…,N,h=2/N. Applying the techniques and rules of Mickens [31], we develop a nonstandard finite difference method for (18). Let β1 and β2 be denominator functions for the discrete difference schemes corresponding to the constant coefficient subequations(35a)−εd2yj+1dx2+λyj+1x=0,(35b)εd2yjdx2−λyjx=0,respectively, where ax≥λ>0.The characteristic equations of each of (35a) and (35b) has two distinct roots r1=expλ/εx and r2=exp−λ/εx. Using the theory of finite differences methods for (35a) and (35b) respectively, we have(36a)−εyi−1j+1−2yij+1+yi+1j+1β12+λyij+1=0,(36b)εyi−1j−2yij+yi+1jβ22−λyij=0.Now, if we denoteyij+1 to be the approximation of (35a) and (35b) at the grid point xi, then it is given as yij+1=c1eλ/εxi+c2e−fracλεxi. Following the approaches used in [30], we can obtain(37)detyi−1j+1eλ/εxi−1e−λ/εxi−1yi−1j+1eλ/εxie−λ/εxiyi+1j+1eλ/εxi+1e−λ/εxi+1=0.Simplifying this equation gives us(38)yi−1j+1−2coshλεhyij+1+yi+1j+1=0.By the hyperbolic identitysinhd/2=±coshd−1/2, we have(39)coshλεh=2sinh2λεh2+1.Putting (39) into (38) gives(40)yi−1j+1+−4sinh2λεh2−2yij+1+yi+1j+1=0.Comparing the coefficients in (36a) and (40), we obtain the denominator function as β12=4ε/λsinh2λ/εh/2. By a similar procedure, from (35b) we can obtain β22=4ε/λ∗sinh2λ∗/εh/2, where ax≥λ∗>0. If we adopt these for the variable coefficient problem, then we have β12=4ε/aj+1xisinh2aj+1xi/εh/2 and β22=4ε/ajxisinh2ajxi/εh/2. Using these denominator functions into the semidiscrete scheme (23), we obtain the discrete scheme as follows:(41)2μΔt+LεN,Myij+1=ϑij+1,j+1=1,2,…,M,where(42)LεN,Myij+1=−εyi−1j+1−2yij+1+yi+1j+1β12+aiyij+1,i=1,2,…,N2,−εyi−1j+1−2yij+1+yi+1j+1β12+aiyij+1+biyi−N/2j+1,i=N2+1,…,N−1,ϑij+1=εyi−1j−2yij+yi+1jβ22−aiyij−2biϕi−N/2j+1/2+2μΔt+2gij+1/2,i=1,2,…,N2,εyi−1j−2yij+yi+1jβ22−aiyij−biyi−N/2j+2μΔt+2gij+1/2,i=N2+1,…,N−1,with the discrete initial and boundary conditions yi0=y0xi,i=0,1,…,N, yij+1=ϕij+1,i=−N/2,−N/2+1,…,0,j+1=1,2,…,M and yNj+1=γj+1,j+1=1,2,…,M. From (41), we can obtain a system(43)Aiyi−1j+1+Biyij+1+Aiyi+1j+1=Diyi−1j+Eiyij+Diyi+1j+Gij+1,where(44)Ai=−ai4sinh2aij+1/εh/2,Bi=ai2sinh2ai/εh/2+ai+2ai1−exp−aiΔt,Di=ai4sinh2ai/εh/2,Ei=−ai2sinh2ai/εh/2−ai+2ai1−exp−aiΔt,Gij+1=2gij+1/2−2biϕi−N/2j+1/2,i=1,2,…,N2,j+1=1,2,…,M,2gij+1/2−2biyi−N/2j+1/2,i=N2+1,…,N−1,j+1=1,2,…,M.Lemma 10 (Discrete maximum principle).
Supposeωij+1 be any mesh function. If ω0j+1≥0, ωNj+1≥0 and 2/μΔt+LεN,Mωij+1≥0, i=1,2,…,N−1, then wij+1≥0, i=0,1,…,N,j+1=1,2,…,M.Proof.
Assume that there existι∈1,…,N−1 such that ωιj+1=mini∈0,Nωij+1<0. From the hypothesis, we see that ι∉0,N. Assuming that ωι+1j+1−ωιj+1>0 and ωιj+1−ωι−1j+1<0,Case 1: forι=1,2,…,N/2, we have(45)2μΔt+LεN,Mωιj+1=2ωιj+1μΔt−εωι−1j+1−2ωιj+1+ωj+1j+1β12+aιωιj+1<0.Case 2: forι=N/2+1,N/2+2,…,N−1, we have(46)2μΔt+LεN,Mωιj+1=2ωιj+1μΔt−εωι−1j+1−2ωιj+1+ωj+1j+1β12+aιωιj+1+bιωι−N/2j+1≤2ωιj+1μΔt−εωι−1j+1−2ωιj+1+ωj+1j+1β12+aιωιj+1+bιωιj+1<0.
By the previously mentioned cases forι=1,2,…,N−1, we obtained 2/μΔt+LεN,Mωιj+1≤0, which is a contradiction to the hypothesis. Thus, the assumption ωιj+1<0 is wrong, which implies that ωij+1≥0,∀i=0,1,…,N.Lemma 11 (Uniform stability estimate).
The numerical solutionyij+1 of (41), i=0,1,…,N, j+1=1,2,…,M is bounded as yij+1≤ϑij+1/1+2λΔt+maxy0j+1,yNj+1.Proof.
Define a barrier functionωi,j+1±=ϑij+1/1+2λΔt+maxy0j+1,yNj+1±yij+1.Fori=0, ω0,j+1±=ϑ0j+1/1+2λΔt+maxy0j+1,yNj+1±y0j+1≥ϑ0j+1/1+2λΔt≥0.Fori=N, ωN,j+1±=ϑNj+1/1+2λΔt+maxy0j+1,yNj+1±yNj+1≥ϑNj+1/1+2λΔt≥0.Fori=1,2,…,N/2, we have(47)2μΔt+LεN,Mωi,j+1±=−εωi−1,j+1±+2εωi,j+1±−εωi+1,j+1±β12+aiωi,j+1±≥aimaxy0j+1,yNj+1≥0.Fori=N/2+1,…,N−1, we have(48)2μΔt+LεN,Mωi,j+1±=−εωi−1,j+1±+2εωi,j+1±−εωi+1,j+1±β12+aiωi,j+1±+biωi−N/2,j+1±≥ai+bimaxy0j+1,yNj+1≥0.
Applying the discrete maximum principle givesωi,j+1±≥0,i=0,1,…,N, which implies the discrete stability estimate.Theorem 1.
Letyj+1xi be the solution of (23) and yij+1 be the solution of (41). Then, the truncation error is bounded as follows:(49)2μΔt+LεN,Myj+1xi−yij+1≤Ch21+supxi∈0,1d1xi,1+supxi∈1,2d2xi,where d1xi=exp−axi/εxi+exp−axi/ε1−xi/ε3/2 and d2xi=exp−axi/εxi−1+exp−axi/ε2−xi/ε3/2.Proof.
From the theory of finite differences, we define the finite difference operators asD+yij+1=yi+1j+1−yij+1/h, D−yij+1=yij+1−yi−1j+1/h and D+D−yij+1=D+−D−yij+1/h2. Then, using the differential and difference operators, we have(50)2μΔt+LεN,Myj+1xi−yij+1=−εd2yj+1xidx2−h2D+D−yij+1βι2≤εd2yj+1xidx2−D+D−yij+1+εh2βι2−1D+D−yij+1≤Cεh2d4yj+1xidx4+εh2βι2−1D+D−yij+1.
Here, claim thatεh2/βι2−1≤Ch2. For the denominator function βι2=4ε/axisinh2axi/εh/2 from the nonstandard finite difference, let us define ξ∈0,∞, such that ξ=axih2/ε. Then, we have(51)βι2=4h2ξsinh2ξ2=h2ξeξ−2+e−ξ.
Using this, we obtain(52)εh2βι2−1=εξeξ−2+e−ξ−1=axih2ξ+2−eξ−e−ξξeξ−2+e−ξ.
LetZξ=ξ+2−eξ−e−ξ/ξeξ−2+e−ξ. Using L’Hospital’s rule repeatedly, we can find that(53)limξ⟶0Zξ=−112andlimξ⟶∞Zξ=0.
This implies thatZξ≤C, so that εh2/βι2−1≤Caxih2≤Ch2, which is the claim. Thus, we have(54)2μΔt+LεN,Myj+1xi−yij+1≤Cεh2d4yj+1xidx4+Ch2d2yj+1xidx2.
Thus, applying the bound of derivatives of the solution, we obtain the required estimate.Lemma 12.
For a given meshN, assuming ε⟶0 gives(55)limε⟶0maxxi∈0,1exp−axi/εxi+exp−axi/ε1−xiεk/2=0,maxxi∈1,2exp−axi/εxi−1+exp−axi/ε2−xiεk/2=0,for all i=1,2,…,N−1 and k=1,2,3,….Proof.
For the proof, we refer to Lemma 3.3 of [32].Theorem 2.
Ifyxi,tj+1 is the solution of (1) and yij+1 is the solution of (41), then the uniform error is given by sup0<ε≪1limi=0,1,…,Nyxi,tj+1−yij+1≤CΔt2+N−2,j+1=1,2,…,M.Proof.
Using Lemma12 into Theorem 1 yields 2/μΔt+LεN,Myj+1xi−yij+1≤Ch2. Then by Lemma 10, we obtain(56)yj+1xi−yij+1≤Ch2.
Sinceyxi,tj+1−yij+1=yxi,tj+1−yj+1xi+yj+1xi−yij+1≤yxi,tj+1−yj+1xi+yj+1xi−yij+1, by Lemma 8 and equation (56), we can conclude the result.


## 3.1. Temporal Semidiscretization

Let0,T be a given time domain for some finite time T, which can be subdivided into M equal intervals forming uniform temporal mesh as ℜtM=tj=jΔt,j+1=1,2,…,M−1,T=MΔt. From the rules upon which nonstandard schemes are based [29], the time derivative in equation (1) takes the form dy/dt⟶yj+1-yj/µΔt, where μΔt is the denominator function satisfying μΔt=Δt+OΔt2. Thus, for the subequation dy/dt=−λy, where ax≥λ>0, the standard forward Euler discretization is yj+1-yj/Δt=−λyj, which can be written as follows:(19)yj+1=1−λΔtyj.However, it is well known that equation (19) leads to instability if λΔt>1 [30]. Instead, from the fact that 1−λΔt=exp−λΔt+OλΔt2 for sufficiently small λΔt (i.e., 0<λΔt≪1), it is possible to make the replacement 1−λΔt≈exp−λΔt, so that Δt≈1-exp−λΔt/λ. Adopting this for the variable coefficient, we have a denominator function(20)μΔt=1−exp−axΔtax.Now, applying the Crank–Nicolson method in the temporal direction incorporating the denominator function (20), we obtain a semidiscrete scheme as follows:(21)LεMyj+1x≡yj+1x−yjxμΔt−ε2∂2yj+1∂x2+∂2yj∂x2+ax2yj+1x+yjx+bx2yj+1x−1+yjx−1=12gj+1x+gjx,and rearranging terms, we can rewrite as follows:(22)LεMyj+1x≡−ε∂2yj+1∂x2+2μΔt+axyj+1x=ε∂2yj∂x2+2μΔt−axyjx−2bxϕj+1/2x−1+2gj+1/2x,forx∈0,1,−ε∂2yj+1∂x2+2μΔt+axyj+1x+2bxyj+1x−1=ε∂2yj∂x2+2μΔt−axyjx+2gj+1/2x,forx∈1,2,subjected to yj+1x,0=y0j+1x,x∈Γ, yj+1x=ϕx,tj+1 and yj+12=γ2,tj+1 for j+1=1,2,…,M−1, where yj+1x is the approximate of yx,tj+1 at j+1th level. We can rewrite equation (22) as follows:(23)2μΔt+LεMyj+1x=ϑj+1x,x∈Γ,where(24)LεMyj+1x=−ε∂2yj+1∂x2+axyj+1x,x∈0,1,−ε∂2yj+1∂x2+axyj+1x+bxyj+1x−1,x∈1,2,ϑj+1x=ε∂2yj∂x2−axyjx−2bxϕj+1/2x−1+2gj+1/2x+2μΔt,x∈0,1,ε∂2yj∂x2−axyjx−bxyjx−1+2gj+1/2x+2μΔt,x∈1,2,where ϕj+1/2x−1=1/2ϕj+1x−1+ϕjx−1 and gj+1/2x=1/2gj+1x+gjx.Lemma 5 (Semidiscrete maximum principle).
For a sufficiently smooth functionωj+1 in ℜ¯, if ωj+10≥0, ωj+12≥0 and Lεωj+1x≥0 for all x∈Γ, then ωj+1x≥0, ∀x∈Γ¯.Proof.
Letr∈Γ¯ such that ωj+1r=minΓ¯ωj+1x and suppose ωj+1r<0. From the given conditions, we have r∉∂Γ and ωxj+1r=0, ωxxj+1r≥0.Case 1: forr∈0,1, we have(25)Lε,1Mωj+1r=−ε∂2ωj+1r∂x2+2μΔt+arωj+1r≤0.Case 2: forr∈1,2, we have(26)Lε,2Mωj+1r=−ε∂2ωj+1r∂x2+2μΔt+arωj+1r+brωj+1r−1≤−ε∂2ωj+1∂x2r+2μΔt+ar+brωj+1r≤0.
From the two cases, we observe that the given hypothesis is contradicted and hence our assumption fails, so thatωj+1x≥0, ∀x∈Γ¯. Thus, the maximum principle is satisfied by the operator LεM. Also, for the semidiscrete operator, we have(27)Lε,xM−1≤11+2λΔt,which is used to estimate truncation errors in the semidiscrete scheme.Lemma 6 (Stability estimate of the semidiscrete problem).
Let the solution of equation (22) be yj+1x. Then. it can be estimated as yj+1x≤ϑj+1/1+2λΔt+maxyj+10,yj+12,∀x∈Γ¯.Proof.
The proof can be easily shown following similar procedure as Lemma3.Definition 1.
At thej+1th time level, we can define local truncation error as the difference between the solution yx,tj+1 and the numerical solution yj+1x of the semidiscrete problem (23) and given ej+1x=yx,tj+1−yj+1x and the global error Ej+1 as the contribution of local truncation error up to the j+1th time level.Lemma 7 (Estimate of the local truncation error).
Suppose it is given that∂kyx,t/∂tk≤C for all x,t∈ℜ¯ and k=0,1,2. Then, at j+1th time level, we have ej+1≤CΔt3.Proof.
Let us choose the center pointx,tj+1/2 of the temporal mesh as a convenient base point where the Taylor series is expanded. Then, yx,tj+1 and yx,tj can be expanded as follows:(28a)yj+1x=yj+1/2x+Δt2∂yj+1/2∂t+Δt28∂yj+1/2∂t2+OΔt3,(28b)yjx=yj+1/2x−Δt2∂yj+1/2∂t+Δt28∂yj+1/2∂t2+OΔt3.
Now, subtracting (28b) from (28a) gives(29)yj+1x−yjxμΔt=∂yj+1/2∂t+OΔt3.
Using (1) into (29) at x,tj+1/2, we obtain(30)yj+1x−yjxμΔt=ε∂2yj+1/2∂x2−axyj+1/2x−bxyj+1/2x−1+gx,tj+1/2+OΔt3.
Putting (30) into the time semidiscrete problem (23) gives(31)ε∂2yj+1/2∂x2−axyj+1/2x−bxyj+1/2x−1+gx,tj+1/2+OΔt3−ε2∂2yj+1∂x2+∂2yj∂x2+ax2yj+1x+yjx+bx2yj+1x−1+yjx−1=12gj+1x+gjx,which can be written equivalently as follows:(32)LεMyx,tj+1−yj+1x=OΔt3.
Using (27) into (32), we obtain a boundary value problem of the type LMej+1x=OΔt3 with ej+10,tj+1=0 and ej+12=0. Applying Lemma 5, we can obtain an estimate for the local truncation error as ej+1≤CΔt3.Lemma 8 (Estimation of the global error).
Assuming that Lemma7 holds, the global error in the time direction can be estimated as Ej+1≤CΔt2, ∀j such that jΔt≤T.Proof.
Using Lemma7 up to the j+1th level, we have(33)Ej+1=∑k=1jek,jΔt≤T≤e1+⋯+ej≤CΔt2.
Thus,Ej+1≤CΔt2, where C is independent of ε and the mesh numbers.Lemma 9 (Bound of the derivatives of the solution).
Let the solution of (23) be yj+1x. Then its derivatives can be bounded as follows:(34)dnyj+1dxn≤C1+ε−n/2exp−λεx+exp−λε1−x,x∈Γ¯,n=014,C1+ε−n/2exp−λεx−1+exp−λε2−x,x∈Γ¯,n=014.Proof.
The proof can be carried out following the procedures given in Theorem 3 of [26].


## 3.2. Fully Discretized Scheme

Let [0, 2] be subdivided intoN equal steps of length h and form a mesh ℜxN=0=x0,x1,…,xN/2=1,xN/2+1,…,xN=2,xi=ih,i=0,1,…,N,h=2/N. Applying the techniques and rules of Mickens [31], we develop a nonstandard finite difference method for (18). Let β1 and β2 be denominator functions for the discrete difference schemes corresponding to the constant coefficient subequations(35a)−εd2yj+1dx2+λyj+1x=0,(35b)εd2yjdx2−λyjx=0,respectively, where ax≥λ>0.The characteristic equations of each of (35a) and (35b) has two distinct roots r1=expλ/εx and r2=exp−λ/εx. Using the theory of finite differences methods for (35a) and (35b) respectively, we have(36a)−εyi−1j+1−2yij+1+yi+1j+1β12+λyij+1=0,(36b)εyi−1j−2yij+yi+1jβ22−λyij=0.Now, if we denoteyij+1 to be the approximation of (35a) and (35b) at the grid point xi, then it is given as yij+1=c1eλ/εxi+c2e−fracλεxi. Following the approaches used in [30], we can obtain(37)detyi−1j+1eλ/εxi−1e−λ/εxi−1yi−1j+1eλ/εxie−λ/εxiyi+1j+1eλ/εxi+1e−λ/εxi+1=0.Simplifying this equation gives us(38)yi−1j+1−2coshλεhyij+1+yi+1j+1=0.By the hyperbolic identitysinhd/2=±coshd−1/2, we have(39)coshλεh=2sinh2λεh2+1.Putting (39) into (38) gives(40)yi−1j+1+−4sinh2λεh2−2yij+1+yi+1j+1=0.Comparing the coefficients in (36a) and (40), we obtain the denominator function as β12=4ε/λsinh2λ/εh/2. By a similar procedure, from (35b) we can obtain β22=4ε/λ∗sinh2λ∗/εh/2, where ax≥λ∗>0. If we adopt these for the variable coefficient problem, then we have β12=4ε/aj+1xisinh2aj+1xi/εh/2 and β22=4ε/ajxisinh2ajxi/εh/2. Using these denominator functions into the semidiscrete scheme (23), we obtain the discrete scheme as follows:(41)2μΔt+LεN,Myij+1=ϑij+1,j+1=1,2,…,M,where(42)LεN,Myij+1=−εyi−1j+1−2yij+1+yi+1j+1β12+aiyij+1,i=1,2,…,N2,−εyi−1j+1−2yij+1+yi+1j+1β12+aiyij+1+biyi−N/2j+1,i=N2+1,…,N−1,ϑij+1=εyi−1j−2yij+yi+1jβ22−aiyij−2biϕi−N/2j+1/2+2μΔt+2gij+1/2,i=1,2,…,N2,εyi−1j−2yij+yi+1jβ22−aiyij−biyi−N/2j+2μΔt+2gij+1/2,i=N2+1,…,N−1,with the discrete initial and boundary conditions yi0=y0xi,i=0,1,…,N, yij+1=ϕij+1,i=−N/2,−N/2+1,…,0,j+1=1,2,…,M and yNj+1=γj+1,j+1=1,2,…,M. From (41), we can obtain a system(43)Aiyi−1j+1+Biyij+1+Aiyi+1j+1=Diyi−1j+Eiyij+Diyi+1j+Gij+1,where(44)Ai=−ai4sinh2aij+1/εh/2,Bi=ai2sinh2ai/εh/2+ai+2ai1−exp−aiΔt,Di=ai4sinh2ai/εh/2,Ei=−ai2sinh2ai/εh/2−ai+2ai1−exp−aiΔt,Gij+1=2gij+1/2−2biϕi−N/2j+1/2,i=1,2,…,N2,j+1=1,2,…,M,2gij+1/2−2biyi−N/2j+1/2,i=N2+1,…,N−1,j+1=1,2,…,M.Lemma 10 (Discrete maximum principle).
Supposeωij+1 be any mesh function. If ω0j+1≥0, ωNj+1≥0 and 2/μΔt+LεN,Mωij+1≥0, i=1,2,…,N−1, then wij+1≥0, i=0,1,…,N,j+1=1,2,…,M.Proof.
Assume that there existι∈1,…,N−1 such that ωιj+1=mini∈0,Nωij+1<0. From the hypothesis, we see that ι∉0,N. Assuming that ωι+1j+1−ωιj+1>0 and ωιj+1−ωι−1j+1<0,Case 1: forι=1,2,…,N/2, we have(45)2μΔt+LεN,Mωιj+1=2ωιj+1μΔt−εωι−1j+1−2ωιj+1+ωj+1j+1β12+aιωιj+1<0.Case 2: forι=N/2+1,N/2+2,…,N−1, we have(46)2μΔt+LεN,Mωιj+1=2ωιj+1μΔt−εωι−1j+1−2ωιj+1+ωj+1j+1β12+aιωιj+1+bιωι−N/2j+1≤2ωιj+1μΔt−εωι−1j+1−2ωιj+1+ωj+1j+1β12+aιωιj+1+bιωιj+1<0.
By the previously mentioned cases forι=1,2,…,N−1, we obtained 2/μΔt+LεN,Mωιj+1≤0, which is a contradiction to the hypothesis. Thus, the assumption ωιj+1<0 is wrong, which implies that ωij+1≥0,∀i=0,1,…,N.Lemma 11 (Uniform stability estimate).
The numerical solutionyij+1 of (41), i=0,1,…,N, j+1=1,2,…,M is bounded as yij+1≤ϑij+1/1+2λΔt+maxy0j+1,yNj+1.Proof.
Define a barrier functionωi,j+1±=ϑij+1/1+2λΔt+maxy0j+1,yNj+1±yij+1.Fori=0, ω0,j+1±=ϑ0j+1/1+2λΔt+maxy0j+1,yNj+1±y0j+1≥ϑ0j+1/1+2λΔt≥0.Fori=N, ωN,j+1±=ϑNj+1/1+2λΔt+maxy0j+1,yNj+1±yNj+1≥ϑNj+1/1+2λΔt≥0.Fori=1,2,…,N/2, we have(47)2μΔt+LεN,Mωi,j+1±=−εωi−1,j+1±+2εωi,j+1±−εωi+1,j+1±β12+aiωi,j+1±≥aimaxy0j+1,yNj+1≥0.Fori=N/2+1,…,N−1, we have(48)2μΔt+LεN,Mωi,j+1±=−εωi−1,j+1±+2εωi,j+1±−εωi+1,j+1±β12+aiωi,j+1±+biωi−N/2,j+1±≥ai+bimaxy0j+1,yNj+1≥0.
Applying the discrete maximum principle givesωi,j+1±≥0,i=0,1,…,N, which implies the discrete stability estimate.Theorem 1.
Letyj+1xi be the solution of (23) and yij+1 be the solution of (41). Then, the truncation error is bounded as follows:(49)2μΔt+LεN,Myj+1xi−yij+1≤Ch21+supxi∈0,1d1xi,1+supxi∈1,2d2xi,where d1xi=exp−axi/εxi+exp−axi/ε1−xi/ε3/2 and d2xi=exp−axi/εxi−1+exp−axi/ε2−xi/ε3/2.Proof.
From the theory of finite differences, we define the finite difference operators asD+yij+1=yi+1j+1−yij+1/h, D−yij+1=yij+1−yi−1j+1/h and D+D−yij+1=D+−D−yij+1/h2. Then, using the differential and difference operators, we have(50)2μΔt+LεN,Myj+1xi−yij+1=−εd2yj+1xidx2−h2D+D−yij+1βι2≤εd2yj+1xidx2−D+D−yij+1+εh2βι2−1D+D−yij+1≤Cεh2d4yj+1xidx4+εh2βι2−1D+D−yij+1.
Here, claim thatεh2/βι2−1≤Ch2. For the denominator function βι2=4ε/axisinh2axi/εh/2 from the nonstandard finite difference, let us define ξ∈0,∞, such that ξ=axih2/ε. Then, we have(51)βι2=4h2ξsinh2ξ2=h2ξeξ−2+e−ξ.
Using this, we obtain(52)εh2βι2−1=εξeξ−2+e−ξ−1=axih2ξ+2−eξ−e−ξξeξ−2+e−ξ.
LetZξ=ξ+2−eξ−e−ξ/ξeξ−2+e−ξ. Using L’Hospital’s rule repeatedly, we can find that(53)limξ⟶0Zξ=−112andlimξ⟶∞Zξ=0.
This implies thatZξ≤C, so that εh2/βι2−1≤Caxih2≤Ch2, which is the claim. Thus, we have(54)2μΔt+LεN,Myj+1xi−yij+1≤Cεh2d4yj+1xidx4+Ch2d2yj+1xidx2.
Thus, applying the bound of derivatives of the solution, we obtain the required estimate.Lemma 12.
For a given meshN, assuming ε⟶0 gives(55)limε⟶0maxxi∈0,1exp−axi/εxi+exp−axi/ε1−xiεk/2=0,maxxi∈1,2exp−axi/εxi−1+exp−axi/ε2−xiεk/2=0,for all i=1,2,…,N−1 and k=1,2,3,….Proof.
For the proof, we refer to Lemma 3.3 of [32].Theorem 2.
Ifyxi,tj+1 is the solution of (1) and yij+1 is the solution of (41), then the uniform error is given by sup0<ε≪1limi=0,1,…,Nyxi,tj+1−yij+1≤CΔt2+N−2,j+1=1,2,…,M.Proof.
Using Lemma12 into Theorem 1 yields 2/μΔt+LεN,Myj+1xi−yij+1≤Ch2. Then by Lemma 10, we obtain(56)yj+1xi−yij+1≤Ch2.
Sinceyxi,tj+1−yij+1=yxi,tj+1−yj+1xi+yj+1xi−yij+1≤yxi,tj+1−yj+1xi+yj+1xi−yij+1, by Lemma 8 and equation (56), we can conclude the result.


## 4. Numerical Experiments and Discussion

To demonstrate the application of the developed numerical method, we treat model problems in this section. For the considered examples, since the exact solutions are not given, we use [33] as eεN,M=lim1≤i,j≤N,Myi,j+1N,M−yi,j+12N,4M. The absolute maximum error is determined as eN,M=maxεeεN,M. Also, we can determine the maximum pointwise convergence rate as rεN,M=logeεN,M−logeε2N,4M/log2. And the uniform convergence rate is obtained as rN,M=maxεrεN,M.Example 1.
Consider−∂y/dt+ε∂2y/∂x2−5yx,t+2yx−1,t=−2, yx,0=sinπx for x∈0,2, yx,t=0 for x,t∈x,t:x∈−1,0,0≤t≤2 and y2,t=0 for 2,t∈2,t:0≤t≤2 [19].Example 2.
Consider−∂y/dt+ε∂2y/∂x2−x+6yx,t+x2+1yx−1,t=−3, yx,0=0 for x∈0,2, yx,t=0 for x,t∈x,t:x∈−1,0,0≤t≤2 and y2,t=0 for 2,t∈2,t:0≤t≤2 [18].Example 3.
Consider−∂y/dt+ε∂2y/∂x2−3yx,t+yx−1,t=−1, yx,0=0 for x∈0,2, yx,t=0 for x,t∈x,t:x∈−1,0,0≤t≤2 and y2,t=0 for 2,t∈2,t:0≤t≤2 [18, 19].For each model example, maximum pointwise errors, rate of convergences, and numerical solutions are computed using MATLAB R2019a software package. We have computed the maximum pointwise error and convergence rate at a point using a variant of the double mesh principle, as shown in Tables1–3. From Tables 1–3, we see that increasing the mesh numbers decreases the truncation error, and by decreasing ε, we obtain stabled pointwise errors. This shows that our scheme is ε-uniformly convergent. The accuracy of the results of the proposed scheme is compared with other published works for each example, as given in Table 4.Table 1
Maximum absolute error and rate of convergence for Example1.
εN⟶3264128256512↓M⟶83212851220482−045.9335e − 021.0399e − 021.4848e − 031.9482e − 044.1549e − 052−067.9072e − 021.4990e − 022.3170e − 033.1664e − 044.1042e − 052−088.1385e − 022.0199e − 024.1768e − 036.9064e − 041.2562e − 042−108.1361e − 022.0405e − 026.2087e − 031.8411e − 034.6589e − 042−128.1240e − 022.0054e − 025.2325e − 032.9563e − 038.3560e − 042−148.1230e − 022.0042e − 024.7207e − 031.7769e − 035.1954e − 042−168.1229e − 022.0041e − 024.7136e − 031.1282e − 033.0203e − 042−188.1229e − 022.0041e − 024.7134e − 031.1209e − 032.7989e − 042−208.1228e − 022.0041e − 024.7134e − 031.1209e − 032.7197e − 042−228.1229e − 022.0041e − 024.7134e − 031.1209e − 032.7196e − 042−248.1229e − 022.0041e − 024.7134e − 031.1209e − 032.7196e − 04eN,M8.1385e − 022.0405e − 026.2087e − 032.9563e − 038.3560e − 04rN,M1.99581.71661.07051.8229—Table 2
Maximum absolute error and rate of convergence for Example2.
εN⟶183672144288↓M⟶1872288115246082−042.8392e − 026.7928e − 031.1932e − 031.8358e − 043.2692e − 052−063.0557e − 029.9272e − 032.7617e − 036.2088e − 041.4208e − 042−083.0410e − 029.3899e − 033.0633e − 038.6458e − 042.0542e − 042−103.0405e − 029.1485e − 032.8443e − 037.2955e − 041.8578e − 042−123.0404e − 029.1460e − 032.5253e − 036.0965e − 041.8431e − 042−143.0404e − 029.1460e − 032.5229e − 036.9753e − 041.9546e − 042−163.0404e − 029.1460e − 032.5229e − 036.9489e − 042.0669e − 042−183.0404e − 029.1460e − 032.5229e − 036.9489e − 042.0335e − 042−203.0404e − 029.1460e − 032.5229e − 036.9489e − 042.0335e − 042−223.0404e − 029.1460e − 032.5229e − 036.9489e − 042.0335e − 042−243.0404e − 029.1460e − 032.5229e − 036.9489e − 042.0335e − 04eN,M3.0557e − 029.9272e − 033.0633e − 038.6458e − 042.0669e − 04rN,M1.62201.69631.82502.0645—Table 3
Maximum absolute error and rate of convergence for Example3.
εN⟶3264128256512↓M⟶83212851220482−041.2490e − 022.1437e − 032.8863e − 047.9803e − 052.1163e − 052−062.1674e − 024.0560e − 036.0966e − 048.2225e − 051.2670e − 052−082.7470e − 027.3034e − 031.4872e − 032.7381e − 045.9638e − 052−102.7310e − 028.3566e − 031.9371e − 034.8577e − 041.1996e − 042−122.7287e − 027.6132e − 031.9513e − 034.8654e−041.1612e − 042−142.7287e − 027.5769e − 032.0051e − 034.5717e − 041.1799e − 042−162.7287e − 027.5768e − 031.9603e − 035.4465e − 041.2315e − 042−182.7287e − 027.5768e − 031.9602e − 034.9447e − 041.2319e − 042−202.7287e − 027.5768e − 031.9602e − 034.9442e − 041.2393e − 042−222.7287e − 027.5768e − 031.9602e − 034.9442e − 041.2389e − 042−242.7287e − 027.5768e − 031.9602e − 034.9442e − 041.2389e − 04eN,M2.7470e − 028.3566e − 032.0051e − 035.4465e − 041.2393e − 04rN,M1.71692.05921.88032.1358—Table 4
Comparison of the proposed method with other published results in the literature.
Rates of convergence for Example1N:64128256512M:321285122048Proposed schemeε−162.08812.06281.90131.8972ε−182.08812.07212.00172.0036ε−202.08812.07212.04312.0032Results in reference [19]ε−161.79081.83141.51211.6261ε−181.79081.83411.51011.6258ε−201.79081.83541.50911.6257Rate of convergence forExample 3Proposed schemeε−161.95051.84772.14491.8972ε−201.95051.98702.00501.8970Results in reference [19]ε−161.33691.54021.63331.6722ε−201.33691.54021.63331.6722Maximum error and rate of convergence forExample 2N=M:183672144Proposed schemeeN1.3081e-037.8288e-044.7719e-042.3539e-04rN0.74060.71301.02070.9999Results in reference [18]eN,M1.3600e-021.0900e-025.5600e-032.2800e-03rN,M0.31930.97121.28610.8787Maximum error and rate of convergence forExample 3Proposed schemeeN9.7731e-034.5542e-031.9746e-036.7654e-04rN1.10161.20561.54531.3814Results in reference [18]eN,M1.1200e-027.0100e-032.9700e-031.1400e-03rN,M0.67601.23901.38141.4895The numerical solutions for each model example are simulated in Figures1–6. Figures 1, 3, and 5 show line plots of the solutions at different time levels, and to show the physical behaviors of the solutions, surface plots are simulated in Figures 2, 4, and 6. From all figures, the effect of ε can be easily observed. That is, if we decrease the value of ε, the width of the two boundary layers and the interior layer decreases, indicating that the errors in the layers have been resolved. To visualize the shape of the solution profile for each example, contour lines are also drawn with the surface plots. Figure 7 depicts log-log plots of the maximum pointwise errors corresponding to N for each example, demonstrating the ε-uniform convergence of the developed scheme.Figure 1
Line plots for the numerical solution of Example1 for N=128 and M=64 at different time levels. (a) [ε = 20], (b) [ε = 2−16].
(a)(b)Figure 2
Surface plots for the numerical solution of Example1 for N=128 and M=64. (a) [ε = 20], (b) [ε = 2−16].
(a)(b)Figure 3
Line plots for the numerical solution of Example2 for N=128 and M=64 at different time levels. (a) [ε = 20], (b) [ε = 2−14].
(a)(b)Figure 4
Surface plots for the numerical solution of Example2 for N=128 and M=64. (a) [ε = 20], (b) [ε = 2−14].
(a)(b)Figure 5
Line plots for the numerical solution of Example3 for N=128 and M=64. (a) [ε = 20], (b) [ε = 2−14].
(a)(b)Figure 6
Surface plots for the numerical solution of Example3 for N=128 and M=64. (a) [ε = 20], (b) [ε = 2−14].
(a)(b)Figure 7
Log-log plots ofN versus maximum absolute error for (a) Example 1, (b) Example 2, and (c) Example 3.
(a)(b)(c)


## 5. Conclusion

In this research work, we focused on developing and analysing a numerical method for singularly perturbed parabolic reaction-diffusion problems with large delays. Due to the presence ofε, the solution involves two boundary layers, and the large spatial delay causes the interior layer at x=1. The original differential equation is transformed into a discrete finite difference scheme by using the Crank–Nicolson method and the central difference method in the temporal and spatial direction, respectively, via nonstandard finite difference methods. Stability and uniform convergence analysis are investigated and proved. To validate the theoretical findings, model numerical examples are considered. For each example, the numerical solution is computed and given in the graphical form. Maximum pointwise errors and rate of convergences of the proposed method are displayed in tabular forms, for example. From the theoretical and numerical results, we can deduce that the method is ε-uniformly convergent. For future work, our target is to extend the present method for meaningful mathematical models of higher dimensional singularly perturbed problems involving large delays.



---

*Source: 7215106-2023-03-02.xml*