# Fitted Tension Spline Method for Singularly Perturbed Time Delay Reaction Diffusion Problems

Authors: Ermias Argago Megiso; Mesfin Mekuria Woldaregay; Tekle Gemechu Dinka
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

A uniformly convergent numerical method is presented for solving singularly perturbed time delay reaction-diffusion problems. Properties of the continuous solution are discussed. The Crank–Nicolson method is used for discretizing the temporal derivative, and an exponentially fitted tension spline method is applied for the spatial derivative. Using the comparison principle and solution bound, the stability of the method is analyzed. The proposed numerical method is second-order uniformly convergent. The theoretical analysis is supported by numerical test examples for various values of perturbation parameters and mesh size.


---

## Body

## 1. Introduction

Delay differential equations in which its highest order derivative term is multiplied by a small parameterε are known as singularly perturbed delay differential equations (SPDDEs). A good example of SPDDEs is the mathematical model for the control system of the furnace used to produce metal sheets which is given by(1)∂u∂t−ε∂2u∂x2=vgux,t−τ+cfux,t−τ−ux,t,where u is the temperature distribution in a metal sheet moving at an instantaneous material strip velocity v and heated by a distributed temperature source given by the function f; both v and f are dynamically changing by a controller monitoring the current temperature distribution. The controller’s finite speed causes a fixed delay of length τ [1].When the perturbation parameterε tends to zero, the smoothness of the solution of SPDDEs descends, and it forms a boundary layer [2]. Since most of the problems have no known analytical/exact solution, formulating numerical methods becomes mandatory [3]. It is of theoretical and practical interest to consider numerical methods for solving SPDDEs. Developing a numerical method whose convergence does not depend on the perturbation parameter has great importance [4]. Owning this, different authors have developed numerical schemes for the solution of SPDDEs. Chakravarthy et al. [5–7] used an adaptive mesh method and fitted operator methods for solving singularly perturbed differential-difference equations. Kumar and Kadhalbajoo [8] computed the numerical solution of SPDDEs using a B-spline collocation method on Shishkin mesh. Kanth and Kumar in [9–11] used the tension spline method to solve singularly perturbed convection-dominated differential equations. Woldaregay et al. [3, 12–16] developed an exponentially fitted numerical method for solving different singularly perturbed differential-difference equations. In [14, 17], they proposed the nonstandard finite difference techniques and the fitted mesh methods, respectively. In [18, 19], an initial value technique with an exponential fitting factor is applied to solve convection-dominated delay differential equations. The difference method on Shishkin for solving a singularly perturbed linear second-order delay differential equation is developed in [20].Govindarao et al. [2] solved a singularly perturbed delay parabolic reaction-diffusion problem using the implicit Euler scheme for time derivative on the uniform mesh, and they used the central difference scheme on the Shishkin mesh for spatial discretization. Kumar and Chandra Sekhara Rao [21] developed a numerical scheme using a suitable combination of fourth order compact difference scheme and central difference scheme on generalized Shishkin mesh. Gowrisankar and Natesan [22] proposed a numerical scheme using the backward-Euler method for temporal discretization and the upwind finite difference scheme for spatial derivative. Erdogan and Amiraliyev [23] treated singularly perturbed second-order delay differential equations using an exponentially fitted difference scheme on a uniform mesh.In some cases, the cubic spline does not provide an accurate picture of the shape of the graph. In such cases, the tension spline method is preferable. A tension spline is a cubic spline that has a tension factor applied to it to stretch the graph closer to the given points. A tension factorδ is a number that is used to generate different conditions for the spline. When the tension factor is set to zero, it gives the usual cubic spline. The objective of this study is to formulate a uniform numerical method for solving singularly perturbed time delay reaction-diffusion problems. Furthermore, we need to establish the stability and uniform convergence of the method. To the best of our knowledge, the exponentially fitted tension spline method has not been used for treating SPDDEs. In this paper, we employ the Crank–Nicolson method to temporal discretization and the exponentially fitted tension spline method to spatial discretization.


## 2. Continuous Problem

On the domainD=Ωx×Ωt=0,1×0,T, we consider singularly perturbed time delay reaction-diffusion problems of the form(2)Lεux,t≡∂u∂t−ε∂2u∂x2+axux,t,=−bx,tux,t−τ+fx,t,with the interval condition(3)ux,t=ψbx,t,x,t∈Ω¯x×−τ,0,=0,1×−τ,0,and the boundary conditions(4)u0,t=ψlt,u1,t,=ψrt,t∈Ω¯t=0,T,where ε∈0,1 is the perturbation parameter, and τ is the delay parameter. The functions a,b,ψl,ψr,ψb, and f are assumed to be sufficiently smooth and bounded that satisfy ax≥α>0,bx,t≥β>0x,t∈D¯=Ω¯x×Ω¯t.


### 2.1. Properties of the Analytical Solution

To guarantee the existence of unique solution for (2)–(4), it is presumed that given data ax,bx,t, and fx,t are Holder’s continuous, and the compatibility conditions at the corner points 0,0,1,0,0,−τ, and 1,−τ are fulfilled. That is(5)ψb0,0=ψl0,ψb1,0=ψr0,∂ψl0∂t−ε∂2ψb0,0∂x2+a0ψb0,0,=−b0,0ψb0,−τ+f0,0,∂ψr0∂t−ε∂2ψb1,0∂x2+a1ψb1,0,=−b1,0ψb1,−τ+f1,0.Solution of the problem in (2)–(4) satisfies the following maximum principle.Lemma 1 (see [24]).
(The maximum principle) Assume thata,b∈C0D¯ and let u∈C2D∩C0D¯. Suppose that u≥0 on ∂D=D¯−D. Then, Lεu≥0 in D implies that u≥0 in D¯.
The stability of the solution of (2)–(4) in the maximum norm is established by the following lemma:Lemma 2 (see [25]).
Letu be the solution of the problem in (2)–(4). Then, it satisfies the ε-uniform upper bound(6)u≤1+αTmaxLεu,u∂D,where the constant α=maxD¯0,1−α≤1 and the norm . is denoted for the maximum norm which is defined as u=maxD¯ux,t.Lemma 3 (see [26]).
Derivatives of the solution of (2)–(4) satisfy the following bound:(7)∂k+lu∂xk∂tl≤C1+ε−k/2exp−αεx+exp−αε1−x,for k=0,1,2,3,4 and l=0,1,2.


## 2.1. Properties of the Analytical Solution

To guarantee the existence of unique solution for (2)–(4), it is presumed that given data ax,bx,t, and fx,t are Holder’s continuous, and the compatibility conditions at the corner points 0,0,1,0,0,−τ, and 1,−τ are fulfilled. That is(5)ψb0,0=ψl0,ψb1,0=ψr0,∂ψl0∂t−ε∂2ψb0,0∂x2+a0ψb0,0,=−b0,0ψb0,−τ+f0,0,∂ψr0∂t−ε∂2ψb1,0∂x2+a1ψb1,0,=−b1,0ψb1,−τ+f1,0.Solution of the problem in (2)–(4) satisfies the following maximum principle.Lemma 1 (see [24]).
(The maximum principle) Assume thata,b∈C0D¯ and let u∈C2D∩C0D¯. Suppose that u≥0 on ∂D=D¯−D. Then, Lεu≥0 in D implies that u≥0 in D¯.
The stability of the solution of (2)–(4) in the maximum norm is established by the following lemma:Lemma 2 (see [25]).
Letu be the solution of the problem in (2)–(4). Then, it satisfies the ε-uniform upper bound(6)u≤1+αTmaxLεu,u∂D,where the constant α=maxD¯0,1−α≤1 and the norm . is denoted for the maximum norm which is defined as u=maxD¯ux,t.Lemma 3 (see [26]).
Derivatives of the solution of (2)–(4) satisfy the following bound:(7)∂k+lu∂xk∂tl≤C1+ε−k/2exp−αεx+exp−αε1−x,for k=0,1,2,3,4 and l=0,1,2.


## 3. Numerical Scheme

### 3.1. Temporal Semidiscretization

LetM be the number of mesh points in the discretization of Ω¯t=0,T, and m is the number of mesh points in the discretization of the interval −τ,0. Note that T=kτ for some positive integer k. The time domain 0,T is discretized using a uniform mesh with time step Δt as ΩtM=tj=jΔt,j=0,1,2,…M,tM=T,Δt=T/M and ΩtM=tj=jΔt,j=0,1,2,…M,tM=T,Δt=T/M. We approximated the temporal derivative term of (2)–(4) using the averaged Crank–Nicolson method, which gives a system of BVPs(8)1+Δt2LεΔtUj+1x=1−Δt2LεΔtUjx−Δtbj+1/2xψbx+Δtfj+1/2x,forj=0,1,…,m,1−Δt2LεΔtUjx−Δtbj+1/2xUj−mx+Δtfj+1/2x,forj=m+1,…,M−1,with the boundary conditions(9)Uj+10=ψltj+1,Uj+11=ψrtj+1,where LεΔtUj+1x=−εUxxj+1x+axUj+1x. Here, Uj+1x is denoted for the approximation of ux,t at the j+1th time level and bj+1/2x=bx,tj+1+bx,tj/2, similarly for fj+1/2x also.Lemma 4.
For eachj=0,1,2,…M−1, the local truncation error in the temporal direction is given by(10)LTEj+1≤C1Δt3,where C1 is a constant independent of ε and mesh size Δt.Proof.
Using Taylor’s series approximation forux,tj and ux,tj+1 centering at tj+1/2, we obtain(11)ux,tj=ux,tj+1/2−Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3,ux,tj+1=ux,tj+1/2+Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3.
From (11), we obtain(12)ux,tj+1−ux,tjΔt=utx,tj+1/2+OΔt2.
Substituting (12) into (2), we obtain(13)ux,tj+1−ux,tjΔt=εuxxx,tj+1/2−axux,tj+1/2+fj+1/2x+OΔt2−bj+1/2xψbx,forj=0,1,…,m,bj+1/2xux,tj+1/2−m,forj=m+1,…,M−1,
Since the errorLTEj+1x:=ux,tj+1−Uj+1x satisfies the semidiscrete difference scheme(14)1+Δt2LεΔtLTEj+1x=OΔt3,LTEj+10=0=LTEj+11.
Hence, by applying the maximum principle, we obtain(15)LTEj+1≤C1Δt3.
Next, we need to show the bound for the global error of temporal discretization. Let us denoteGTEj+1 as the global error up to the j+1th time step.Lemma 5.
The global error in temporal discretization attj+1 time step is given by(16)GTEj+1≤C1Δt2,j=1,2,…,M−1.Proof.
Using the local truncation error up to thej+1th time step given in the above lemma, we obtain the global error up to the j+1th time step as(17)GTEj+1=∑l=1MLTEl+1≤LTE1+LTE2+⋯+LTEj+1≤C1TΔt2sincej+1Δt≤T=CΔt2,C1T=C,where C is a constant independent of ε and Δt.Lemma 6.
LetUj+1x be a solution of the semidiscrete problem in (7), (8) at the j+1 time level. Then, its derivatives satisfy the bound(18)dkUj+1xdxk≤C1+ε−k/2expαεx+expαε1−x,where ax≥α>0 for k=0,1,...,4.Proof.
It follows from Lemma3.


### 3.2. Spatial Discretization

A functionSx,δ=Sx is a class of C2Ω¯x, which interpolates ux at the mesh points xi depending on the parameter δ, reduces to cubic spline in Ω¯x as δ⟶0, and it is termed as a parametric cubic spline function [27, 28]. The spatial domain is discretized into N equal number of subintervals, each of length h=1/N. Let 0=x0,xN=1 and xi=ih,i=0,1,2,…,N be the mesh points. For spatial discretization, we apply an exponentially fitted tension spline method which helps us control the influence of the singular perturbation parameter. In xi,xi+1, the spline function Sx satisfies the differential equation(19)S′′x−δSx=S′′xi−δSxixi+1−xh+S″xi+1−δSxi+1x−xih,where Sxi=uxi and δ>0 is termed as a cubic spline in compression. Solving the linear second-order differential equation in (19) and determining the arbitrary constants from the interpolation conditions Sxi+1=uxi+1,Sxi=uxi, we get(20)Sx=h2λ2sinhλMi+1sinhλx−xih+Misinhλxi+1−xh,−h2λ2Mi+1−λ2h2uxi+1x−xih+Mi−λ2h2uxixi+1−xh,where λ=hδ1/2 and Mk=u″xk for k=i,i±1. Now, differentiating (19) and letting x⟶xi, on the interval xi,xi+1, we obtain(21)S′xi+=uxi+1−uxih−hλ2Mi+11−λsinhλ+Miλcothλ−1.Similarly, on the intervalxi−1,xi , we obtain(22)S′xi−=uxi−uxi−1h+hλ2Miλcothλ−1+Mi−11−λsinhλ.Equating the left and right hand derivatives atxi gives(23)uxi−uxi−1h+hλ2Miλcothλ−1+Mi−11−λsinhλ,=uxi+1−uxih−hλ2Mi+11−λsinhλ+Miλcothλ−1.Rearranging, we obtain the tridiagonal system(24)λ1Mi−1+2λ2Mi+λ1Mi+1=uxi−1−2uxi+uxi+1h2,for i=1,2,3,…,N−1 where λ1=1/λ2λ/sinhλ−1 and λ2=1/λ11−λcothλ.The condition of continuity in (24) ensures the continuity of the first derivatives of the spline Sx at interior nodes. Using into (24), (8) gives(25)−εMk+akUj+1xk=εMk−akUjxk+Δt−bkj+1/2ψbxk+fkj+1/2,forj=0,1,2,…,m,εMk−akUjxk+Δt−bkj+1/2Uj−mxk+fkj+1/2,forj=m+1,…,M−1,for k=i,i±1, where Mk=u″xk,ak=axk and bkj+1/2=bj+1/2xk.To handle the effect of the perturbation parameterε, we use the exponential fitting factor(26)σ=h2/4εaisinhh/2ai/ε2,that stabilizes the discrete scheme. Substituting the fitting factor into (26), (25) we obtain the difference equation as(27)1+Δt2LεΔt,hUij+1:≡vi−Ui−1j+1+vicUij+1+vi+Ui+1j+1,=zi−Ui−1j+zicUij+zi+Ui+1j+Kj,where(28)vi−=−Δt2εσh2−λ1ai−1+λ1,zi−=Δt2εσh2−λ1ai−1+λ1,vic=Δtεσh2+λ2ai+2λ2,zic=−Δtεσh2+λ2ai+2λ2,vi+=−Δt2εσh2−λ1ai+1+λ1,zi+=Δt2εσh2−λ1ai+1+λ1,Kj=λ1Δt−bi−1j+1/2ψbxi−1+fi−1j+1/2+2λ2Δt−bkj+1/2ψbxi+fij+1/2,+λ1Δt−bi+1j+1/2ψbxi+1+fi+1j+1/2,forj=0,1,2,…,m,λ1Δt−bi−1j+1/2Ui−1j−m+fi−1j+1/2+2λ2Δt−bij+1/2Uij−m+fij+1/2,+λ1Δt−bi+1j+1/2Ui+1j−m+fi+1j+1/2,forj=m+1,…,M−1.Lemma 7.
(Discrete comparison principle) Let there exist a comparison functionQij+1 such that if U0j+1≤Q0j+1 and UNj+1≤QNj+1, then 1+Δt/2LεΔt,hUij+1≤1+Δt/2LεΔt,hQij+1, i=1,2,…,N−1, implies that Uij+1≤Qij+1,i=0,1,2,…,N.Proof.
The matrix1+Δt/2LεΔt,hUij+1 is a size of N+1×N+1 with its entries for i=1,2,…,N−1 which are vi−<0,vic>0, and vi+<0. As a result, the coefficient matrix satisfies the M matrix property. So, there is an inverse matrix that is non-negative. This proves that the discrete solution exists and is unique.Lemma 8.
(Uniform stability estimate) The solutionUij+1 of the discrete scheme (21) satisfies the bound(29)Uij+1≤1+Δt/2LεΔt,hUij+11+Δtα/2+maxψltj+1,ψrtt+1,where α>0 is lower bound of ai=axi.Proof.
Let us consider a barrier functionμi,j+1±  =  1+Δt/2LεΔt,hUij+1/1+Δtα/2+maxψltj+1,ψrtt+1±Uij+1. We can easily show that μ0,j+1±≥0 and μN,j+1±≥0.(30)1+Δt2LεΔt,hμi,j+1±=μi,j+1±+Δt2−εσD+D−μi,j+1±+aiμi,j+1±,=1+Δt2LεΔt,hai1+Δt/2LεΔt,hUij+11+Δtα/2+maxψltj+1,ψrtt+1,=1+Δt2ai1+Δt/2LεΔt,hUij+11+Δtα/2+maxψltj+1,ψrtt+1,±Uij+1±Δt2−εσρD+D−Uij+1+aiUij+1,±1+Δt2aiUij+1≥0.
Using the discrete comparison principle, we obtainμi,j+1±≥0,i=0,1,2…,NLemma 9 (see [29, 30]).
For a fixed number of meshN and for positive integer k as ε⟶0(31)limε⟶0max1≤i≤N−1exp−Cxi/εεk/20,limε⟶0max1≤i≤N−1exp−C1−xi/εεk/2=0.Lemma 10.
IfQij+1 be any mesh function such that Q0j+1=QNj+1=0. Then,(32)Qij+1≤1αmax1≤k≤N−1LεΔt,hQkj+1.Proof.
Consider two barrier functions of the formμi,j+1±=S±Qkj+1, where S=1/αmax1≤k≤N−1LεΔt,hQkj+1. It is easily shown that μ0,j+1±≥0,μN,j+1±≥0. Next, we show that(33)LεΔt,hμi,j+1±=LεΔt,hS±Qij+1,=±LεΔt,hQij+1+aiαmax1≤k≤N−1LεΔt,hQkj+1≥0.
Hence, using the discrete comparison principles givesQij+1≤1/amax1≤k≤N−1LεΔt,hQkj+1Theorem 1.
The discrete solution satisfies the following error bound:(34)Uj+1x−Uij+1≤CN−2.Proof.
The truncation error of the scheme (27) is given by(35)LεΔtUj+1x−LεΔt,hUij+1=−εUxxj+1x+axUj+1x−εσh2Ui+1j+1−2Uij+1+Ui−1j+1+θUi−1j+1+θUi+1j+1+2λ2+λ2aiUij+1,where θ=λ11+Δt/2ai and 2λ2=1−2λ1. Using Taylor’s series expansion of the terms Ui+1j+1 and Ui−1j+1, we have(36)Ui±1j+1=Uj+1xi±hUxj+1xi+h22Uxxj+1xi±h36Uxxxj+1xi+h424Uxxxxj+1xi±h5120Uxxxxxj+1xi+⋯
Using the truncated Taylor’s series expansions of the termsUi+1j+1 and Ui−1j+1 yields(37)LεΔtUj+1x−LεΔt,hUij+1=Δt2ai−εΔt2Uxxj+1xi−−εΔtσ2h2h2Uxxj+1x+h212Uxxxxj+1x,+θh2Uxxj+1xi+θh412Uxxxxj+1xi.
Next, we use a truncated Taylor’s series expansion of the exponential fitting factorσ of order five(38)σ=1−ζi2h212+ζi4h4240.
Now, substituting (38) into (37), we obtain(39)LεΔtUj+1xi−LεΔt,hUij+1=ε12Uxxj+1xiζi2+Uxxxxj+1xi+θUxxj+1xih2+,εζi4240Uxxj+1xi−ζi4144Uxxj+1xi+θ12Uxxxxj+1xih4+εζi42880Uxxxxj+1xih6.
Using Lemma6 together with Lemma 9 and the relation N−2>N−4>N−6…, the discrete scheme satisfies the bound(40)LεΔtUj+1xi−LεΔt,hUij+1≤CN−2.
Using the bound in Lemma10, we obtain Uj+1x−Uij+1≤CN−2.Theorem 2.
Letu and U be the solution of (2)–(4) and (27), respectively. Then, the following uniform error bound holds(41)u−U≤CN−2+Δt2.Proof.
The combination of temporal and spatial error bounds gives the required result.


## 3.1. Temporal Semidiscretization

LetM be the number of mesh points in the discretization of Ω¯t=0,T, and m is the number of mesh points in the discretization of the interval −τ,0. Note that T=kτ for some positive integer k. The time domain 0,T is discretized using a uniform mesh with time step Δt as ΩtM=tj=jΔt,j=0,1,2,…M,tM=T,Δt=T/M and ΩtM=tj=jΔt,j=0,1,2,…M,tM=T,Δt=T/M. We approximated the temporal derivative term of (2)–(4) using the averaged Crank–Nicolson method, which gives a system of BVPs(8)1+Δt2LεΔtUj+1x=1−Δt2LεΔtUjx−Δtbj+1/2xψbx+Δtfj+1/2x,forj=0,1,…,m,1−Δt2LεΔtUjx−Δtbj+1/2xUj−mx+Δtfj+1/2x,forj=m+1,…,M−1,with the boundary conditions(9)Uj+10=ψltj+1,Uj+11=ψrtj+1,where LεΔtUj+1x=−εUxxj+1x+axUj+1x. Here, Uj+1x is denoted for the approximation of ux,t at the j+1th time level and bj+1/2x=bx,tj+1+bx,tj/2, similarly for fj+1/2x also.Lemma 4.
For eachj=0,1,2,…M−1, the local truncation error in the temporal direction is given by(10)LTEj+1≤C1Δt3,where C1 is a constant independent of ε and mesh size Δt.Proof.
Using Taylor’s series approximation forux,tj and ux,tj+1 centering at tj+1/2, we obtain(11)ux,tj=ux,tj+1/2−Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3,ux,tj+1=ux,tj+1/2+Δt2utx,tj+1/2+Δt28uttx,tj+1/2+OΔt3.
From (11), we obtain(12)ux,tj+1−ux,tjΔt=utx,tj+1/2+OΔt2.
Substituting (12) into (2), we obtain(13)ux,tj+1−ux,tjΔt=εuxxx,tj+1/2−axux,tj+1/2+fj+1/2x+OΔt2−bj+1/2xψbx,forj=0,1,…,m,bj+1/2xux,tj+1/2−m,forj=m+1,…,M−1,
Since the errorLTEj+1x:=ux,tj+1−Uj+1x satisfies the semidiscrete difference scheme(14)1+Δt2LεΔtLTEj+1x=OΔt3,LTEj+10=0=LTEj+11.
Hence, by applying the maximum principle, we obtain(15)LTEj+1≤C1Δt3.
Next, we need to show the bound for the global error of temporal discretization. Let us denoteGTEj+1 as the global error up to the j+1th time step.Lemma 5.
The global error in temporal discretization attj+1 time step is given by(16)GTEj+1≤C1Δt2,j=1,2,…,M−1.Proof.
Using the local truncation error up to thej+1th time step given in the above lemma, we obtain the global error up to the j+1th time step as(17)GTEj+1=∑l=1MLTEl+1≤LTE1+LTE2+⋯+LTEj+1≤C1TΔt2sincej+1Δt≤T=CΔt2,C1T=C,where C is a constant independent of ε and Δt.Lemma 6.
LetUj+1x be a solution of the semidiscrete problem in (7), (8) at the j+1 time level. Then, its derivatives satisfy the bound(18)dkUj+1xdxk≤C1+ε−k/2expαεx+expαε1−x,where ax≥α>0 for k=0,1,...,4.Proof.
It follows from Lemma3.


## 3.2. Spatial Discretization

A functionSx,δ=Sx is a class of C2Ω¯x, which interpolates ux at the mesh points xi depending on the parameter δ, reduces to cubic spline in Ω¯x as δ⟶0, and it is termed as a parametric cubic spline function [27, 28]. The spatial domain is discretized into N equal number of subintervals, each of length h=1/N. Let 0=x0,xN=1 and xi=ih,i=0,1,2,…,N be the mesh points. For spatial discretization, we apply an exponentially fitted tension spline method which helps us control the influence of the singular perturbation parameter. In xi,xi+1, the spline function Sx satisfies the differential equation(19)S′′x−δSx=S′′xi−δSxixi+1−xh+S″xi+1−δSxi+1x−xih,where Sxi=uxi and δ>0 is termed as a cubic spline in compression. Solving the linear second-order differential equation in (19) and determining the arbitrary constants from the interpolation conditions Sxi+1=uxi+1,Sxi=uxi, we get(20)Sx=h2λ2sinhλMi+1sinhλx−xih+Misinhλxi+1−xh,−h2λ2Mi+1−λ2h2uxi+1x−xih+Mi−λ2h2uxixi+1−xh,where λ=hδ1/2 and Mk=u″xk for k=i,i±1. Now, differentiating (19) and letting x⟶xi, on the interval xi,xi+1, we obtain(21)S′xi+=uxi+1−uxih−hλ2Mi+11−λsinhλ+Miλcothλ−1.Similarly, on the intervalxi−1,xi , we obtain(22)S′xi−=uxi−uxi−1h+hλ2Miλcothλ−1+Mi−11−λsinhλ.Equating the left and right hand derivatives atxi gives(23)uxi−uxi−1h+hλ2Miλcothλ−1+Mi−11−λsinhλ,=uxi+1−uxih−hλ2Mi+11−λsinhλ+Miλcothλ−1.Rearranging, we obtain the tridiagonal system(24)λ1Mi−1+2λ2Mi+λ1Mi+1=uxi−1−2uxi+uxi+1h2,for i=1,2,3,…,N−1 where λ1=1/λ2λ/sinhλ−1 and λ2=1/λ11−λcothλ.The condition of continuity in (24) ensures the continuity of the first derivatives of the spline Sx at interior nodes. Using into (24), (8) gives(25)−εMk+akUj+1xk=εMk−akUjxk+Δt−bkj+1/2ψbxk+fkj+1/2,forj=0,1,2,…,m,εMk−akUjxk+Δt−bkj+1/2Uj−mxk+fkj+1/2,forj=m+1,…,M−1,for k=i,i±1, where Mk=u″xk,ak=axk and bkj+1/2=bj+1/2xk.To handle the effect of the perturbation parameterε, we use the exponential fitting factor(26)σ=h2/4εaisinhh/2ai/ε2,that stabilizes the discrete scheme. Substituting the fitting factor into (26), (25) we obtain the difference equation as(27)1+Δt2LεΔt,hUij+1:≡vi−Ui−1j+1+vicUij+1+vi+Ui+1j+1,=zi−Ui−1j+zicUij+zi+Ui+1j+Kj,where(28)vi−=−Δt2εσh2−λ1ai−1+λ1,zi−=Δt2εσh2−λ1ai−1+λ1,vic=Δtεσh2+λ2ai+2λ2,zic=−Δtεσh2+λ2ai+2λ2,vi+=−Δt2εσh2−λ1ai+1+λ1,zi+=Δt2εσh2−λ1ai+1+λ1,Kj=λ1Δt−bi−1j+1/2ψbxi−1+fi−1j+1/2+2λ2Δt−bkj+1/2ψbxi+fij+1/2,+λ1Δt−bi+1j+1/2ψbxi+1+fi+1j+1/2,forj=0,1,2,…,m,λ1Δt−bi−1j+1/2Ui−1j−m+fi−1j+1/2+2λ2Δt−bij+1/2Uij−m+fij+1/2,+λ1Δt−bi+1j+1/2Ui+1j−m+fi+1j+1/2,forj=m+1,…,M−1.Lemma 7.
(Discrete comparison principle) Let there exist a comparison functionQij+1 such that if U0j+1≤Q0j+1 and UNj+1≤QNj+1, then 1+Δt/2LεΔt,hUij+1≤1+Δt/2LεΔt,hQij+1, i=1,2,…,N−1, implies that Uij+1≤Qij+1,i=0,1,2,…,N.Proof.
The matrix1+Δt/2LεΔt,hUij+1 is a size of N+1×N+1 with its entries for i=1,2,…,N−1 which are vi−<0,vic>0, and vi+<0. As a result, the coefficient matrix satisfies the M matrix property. So, there is an inverse matrix that is non-negative. This proves that the discrete solution exists and is unique.Lemma 8.
(Uniform stability estimate) The solutionUij+1 of the discrete scheme (21) satisfies the bound(29)Uij+1≤1+Δt/2LεΔt,hUij+11+Δtα/2+maxψltj+1,ψrtt+1,where α>0 is lower bound of ai=axi.Proof.
Let us consider a barrier functionμi,j+1±  =  1+Δt/2LεΔt,hUij+1/1+Δtα/2+maxψltj+1,ψrtt+1±Uij+1. We can easily show that μ0,j+1±≥0 and μN,j+1±≥0.(30)1+Δt2LεΔt,hμi,j+1±=μi,j+1±+Δt2−εσD+D−μi,j+1±+aiμi,j+1±,=1+Δt2LεΔt,hai1+Δt/2LεΔt,hUij+11+Δtα/2+maxψltj+1,ψrtt+1,=1+Δt2ai1+Δt/2LεΔt,hUij+11+Δtα/2+maxψltj+1,ψrtt+1,±Uij+1±Δt2−εσρD+D−Uij+1+aiUij+1,±1+Δt2aiUij+1≥0.
Using the discrete comparison principle, we obtainμi,j+1±≥0,i=0,1,2…,NLemma 9 (see [29, 30]).
For a fixed number of meshN and for positive integer k as ε⟶0(31)limε⟶0max1≤i≤N−1exp−Cxi/εεk/20,limε⟶0max1≤i≤N−1exp−C1−xi/εεk/2=0.Lemma 10.
IfQij+1 be any mesh function such that Q0j+1=QNj+1=0. Then,(32)Qij+1≤1αmax1≤k≤N−1LεΔt,hQkj+1.Proof.
Consider two barrier functions of the formμi,j+1±=S±Qkj+1, where S=1/αmax1≤k≤N−1LεΔt,hQkj+1. It is easily shown that μ0,j+1±≥0,μN,j+1±≥0. Next, we show that(33)LεΔt,hμi,j+1±=LεΔt,hS±Qij+1,=±LεΔt,hQij+1+aiαmax1≤k≤N−1LεΔt,hQkj+1≥0.
Hence, using the discrete comparison principles givesQij+1≤1/amax1≤k≤N−1LεΔt,hQkj+1Theorem 1.
The discrete solution satisfies the following error bound:(34)Uj+1x−Uij+1≤CN−2.Proof.
The truncation error of the scheme (27) is given by(35)LεΔtUj+1x−LεΔt,hUij+1=−εUxxj+1x+axUj+1x−εσh2Ui+1j+1−2Uij+1+Ui−1j+1+θUi−1j+1+θUi+1j+1+2λ2+λ2aiUij+1,where θ=λ11+Δt/2ai and 2λ2=1−2λ1. Using Taylor’s series expansion of the terms Ui+1j+1 and Ui−1j+1, we have(36)Ui±1j+1=Uj+1xi±hUxj+1xi+h22Uxxj+1xi±h36Uxxxj+1xi+h424Uxxxxj+1xi±h5120Uxxxxxj+1xi+⋯
Using the truncated Taylor’s series expansions of the termsUi+1j+1 and Ui−1j+1 yields(37)LεΔtUj+1x−LεΔt,hUij+1=Δt2ai−εΔt2Uxxj+1xi−−εΔtσ2h2h2Uxxj+1x+h212Uxxxxj+1x,+θh2Uxxj+1xi+θh412Uxxxxj+1xi.
Next, we use a truncated Taylor’s series expansion of the exponential fitting factorσ of order five(38)σ=1−ζi2h212+ζi4h4240.
Now, substituting (38) into (37), we obtain(39)LεΔtUj+1xi−LεΔt,hUij+1=ε12Uxxj+1xiζi2+Uxxxxj+1xi+θUxxj+1xih2+,εζi4240Uxxj+1xi−ζi4144Uxxj+1xi+θ12Uxxxxj+1xih4+εζi42880Uxxxxj+1xih6.
Using Lemma6 together with Lemma 9 and the relation N−2>N−4>N−6…, the discrete scheme satisfies the bound(40)LεΔtUj+1xi−LεΔt,hUij+1≤CN−2.
Using the bound in Lemma10, we obtain Uj+1x−Uij+1≤CN−2.Theorem 2.
Letu and U be the solution of (2)–(4) and (27), respectively. Then, the following uniform error bound holds(41)u−U≤CN−2+Δt2.Proof.
The combination of temporal and spatial error bounds gives the required result.


## 4. Numerical Results and Discussions

We considered three examples for validating the theoretical analysis. We computed the maximum absolute error and the rate of convergence of the scheme for different mesh lengths and the values ofε. For the test examples, the exact solution is not given. The double mesh technique is applied to compute the maximum absolute error. Let Ui,Nj,M be the computed solution on N and M number of mesh points in space and time discretization, respectively, and Ui,2Nj,2M be the computed solution on the double number of mesh points 2N and 2M by including the midpoints xi+1/2=xi+1+xi/2 and tj+1/2=tj+1+tj/2 into the mesh points. The maximum absolute error is calculated by using the formula(42)EεN,M=maxi,jUi,Nj,M−Ui,2Nj,2M,and the uniform error is calculated by using the formula EN,M=maxε∈0,1EεN,M The rate of convergence is calculated as(43)rεN,M=logEεN,M−logEε2N,2Mlog2,and the uniform rate of convergence is calculated as(44)rN,M=logEN,M−logE2N,2Mlog2.Example 1.
[24] Consider the following problem with the retarded argument(45)∂ux,t∂t−ε∂2ux,t∂x2+1+x22ux,t+ux,t−τ=t3,ux,t=0,x,t∈0,1×−1,0,u0,t=u1,t=0,t∈0,2.Example 2.
[25] Consider the following problem with the retarded argument(46)∂ux,t∂t−ε∂2ux,t∂x2+x2ux,t=t3−ux,t−τ,ux,t=0,x,t∈0,1×−1,0,u0,t=u1,t=0,t∈0,2.Example 3.
[31] Consider the following problem with the retarded argument(47)∂ux,t∂t−ε∂2ux,t∂x2+1.1+x2ux,t=t3−ux,t−τ,ux,t=0,x,t∈0,1×−1,0,u0,t=u1,t=0,t∈0,2.
In Tables1, 2, and 3, respectively, we show the maximum point-wise absolute error and the rate of convergence of Examples 1, 2, and 3. In each column of the tables of 1–3 (or for each N and M as ε⟶0, the maximum absolute error becomes uniform; this shows that the convergence of the method is independent of the perturbation parameter. One can observe in these tables that the method converges with an order of convergence two for each values of the perturbation parameter. The comparison in Table 4 shows that the proposed method is more accurate than the scheme in [24]. In Figures 1 and 2, the boundary layer formation on the solution is given as ε goes small.Table 1
Example1, the maximum absolute error and the rate of convergence of the scheme.
ε↓M = 641282565121024N = 2040801603202−64.2439e-038.8457e-043.5017e-041.8423e-049.4623e-052.26231.33690.926550.961250.97892−82.2925e-024.6993e-031.0750e-032.2823e-049.6683e-052.28642.12812.23581.23920.982282−108.9288e-022.4218e-024.9511e-031.1946e-032.7467e-041.88242.29032051221208223352−128.9288e-022.4218e-024.9511e-031.1946e-032.7467e-041.88242.29032.05122.12082.2335EN,M8.9288e-022.4218e-024.9511e-031.1946e-032.7467e-04rN,M1.88242.29032.05122.12082.2335Table 2
Example2, the maximum absolute error and the rate of convergence of the scheme.
ε↓M = 641282565121024N = 2040801603202−63.9903e-038.1168e-044.1671e-042.2002e-041.1291e-042.29750.961870.92140.962460.981542−82.2316e-024.56683e-031.0299e-032.2688e-041.1671e-042.2882.14872.18250.959000.98032−108.8696e-022.3903e-024.8842e-031.1725e-032.6934e-041.89172.29101.42692.12212.16992−122.1312e-019.2178e-022.4753e-025.0491e-031.2449e-031.20921.89682.29352.02002.0749EN,M2.1312e-019.2178e-022.4753e-025.0491e-031.2449e-03rN,M1.20921.89682.29352.02002.0749Table 3
Example3, the maximum absolute error and the rate of convergence of the scheme.
ε↓M = 641282565121024N = 2040801603202−67.2293e-031.6486e-033.5961e-041.4320e-047.3755e-052.13262.19671.32840.957220.973662−83.3115e-027.7596e-031.8118e-034.2553e-049.3382e-052.09342.09862.0902.18801.26732−101.0114e-013.4597e-028.0418e-031.8959e-034.6360e-041.54762.10512.08462.03192.0927e2−122.0856e-011.0410e-013.5372e-028.1868e-031.9433e-031.00251.55732.11122.07482.0097EN,M2.0856e-011.0410e-013.5372e-028.1868e-031.9433e-03rN,M1.00251.55732.11122.07482.0097Table 4
Example1 comparison of ε-uniform error EN,M and ε-uniform rate of convergence rN,M .
Methods↓M = 641282565121024N = 204080160320Proposed methodEN,M8.9288e-022.4218e-024.9511e-031.1946e-032.7467e-04rN,M1.88242.29032.05122.12082.2335Result in [24]EN,M5.3429e-022.7108e-021.3653e-026.8512e-033.4318e-03rN,M0.97890.98950.99470.99740.9987Figure 1
Solution of Example2, with boundary layer formation, (a)ε=2−6 (b)ε=2−12.
(a)(b)Figure 2
Solution of Example3, with boundary layer formation (a)ε=2−4, (b)ε=2−12.
(a)(b)


## 5. Conclusion

In this paper, the exponentially fitted tension spline method is developed for solving the singularly perturbed time delay reaction-diffusion problems. The solution of considered problems exhibits exponential boundary layers on the left and right side of the spatial domain. We employed the Crank–Nicolson method in temporal discretization and an exponentially fitted tension spline method in spatial discretization to set up the numerical method. The bounds and properties of the analytical solution are explored. An exponential fitting parameter is induced to stabilize the influence of the perturbation parameter on the discrete solution. The stability and convergence analysis of the method are discussed and proved. The performance of the method is illustrated through the numerical examples. It is proved that the method gives second-order convergence in space and in time.



---

*Source: 8669718-2022-10-12.xml*