# A Uniformly Convergent Collocation Method for Singularly Perturbed Delay Parabolic Reaction-Diffusion Problem

Authors: Fasika Wondimu Gelu; Gemechis File Duressa
Journal: Abstract and Applied Analysis
 (2021)

---

## Abstract

In this article, a numerical solution is proposed for singularly perturbed delay parabolic reaction-diffusion problem with mixed-type boundary conditions. The problem is discretized by the implicit Euler method on uniform mesh in time and extended cubic B-spline collocation method on a Shishkin mesh in space. The parameter-uniform convergence of the method is given, and it is shown to beε-uniformly convergent of OΔt+N−2ln2N, where Δt and N denote the step size in time and number of mesh intervals in space, respectively. The proposed method gives accurate results by choosing suitable value of the free parameter λ. Some numerical results are carried out to support the theory.


---

## Body

## 1. Introduction

LetΩ=0,1, D=Ω×0,T, and Γ=Γl∪Γb∪Γr, where Γl and Γr are the left and right sides of the rectangular domain D corresponding to x=0 and x=1, respectively, and Γb is the base of the domain and given by Γb=0,1×−τ,0. Note that Γl=0×0,T and Γr=1×0,T. In this work, we consider the following singularly perturbed delay parabolic reaction-diffusion problem
(1)Lu≡ut−εuxx+ax,tux,t=−bx,tux,t−τ+fx,t,x,t∈D,subject to the initial condition
(2)ux,t=ϕbx,t,x,t∈Γb,and boundary conditions
(3)Blux,t≡ux,t−εuxx,t=ϕlt,x,t∈Γl,Brux,t≡ux,t+εuxx,t=ϕrt,x,t∈Γr,whereε0<ε<<1 is a diffusion parameter whose presence makes the problem singularly perturbed and τ>0 is a delay parameter. For the uniqueness of the solution, we assume that the functions involved in problem (1)–(3) are sufficiently smooth and bounded functions satisfying the following conditions
(4)ax,t≥α>0,bx,t≥β>0,x,t∈D¯.This type of problem arises in various areas of science and engineering. The problem with delay occurs in the field of many biological models like epidemiology and population ecology. A typical example of delay parabolic partial differential equation is the temporal Wazewska-Czyzewska and Lasota equation which describes the survival of red blood cells in animals. This equation may be extended by incorporating a spatial component. The spatiotemporal delay reaction-diffusion equation of the following form is discussed in [1, 2]
(5)∂px,t∂t=d∂2px,t∂x2−δpx,t+qe−apx,t−τ,where Ω⊂ℝ is a bounded domain and x,t∈Ω×0,∞. The state variable px,t denotes the number of red blood cells located at x at time t>0. The constant time delay τ denotes the time needed to produce blood cells. The parameter δ is the death rate of red blood cells, where the parameters q and a are related to the generation of red blood cells.The numerical solution of delay partial differential equation depends not only on the solution at a present stage but also at some past stages. Many researchers have proposed different numerical methods to solve singularly perturbed time delay parabolic reaction-diffusion equations; for instance, see [3–8]. Singularly perturbed two-parameter time delay parabolic problems are studied in [9] using hybrid method based on a uniform mesh in time and a layer-adapted Shishkin mesh in space. In [10], the authors studied singularly perturbed differential-difference convection-diffusion equations using a higher order numerical method. But all the mentioned authors considered singularly perturbed time delay parabolic reaction-diffusion equations subject to Dirichlet boundary conditions. The authors in [11] studied singularly perturbed time delay parabolic reaction-diffusion equation subject to mixed-type boundary conditions, yet to the best of our knowledge, no further study has been done. In this work, we apply the extended form of cubic B-spline collocation method for singularly perturbed time delay parabolic reaction-diffusion problem subject to mixed boundary conditions.B-spline functions have emerged as powerful techniques in the numerical solution of linear and nonlinear partial differential equations. B-spline functions are piecewise polynomial or nonpolynomial functions. B-spline functions were first formulated in 1946. Extended cubic B-spline is an improved version of cubic B-spline, where its basis is constructed in such a way that one free parameter,λ, is included and the degree of the piecewise polynomial is increased but the continuity of the extended cubic B-splines remains in order three. In [12], the authors proposed an extension of cubic B-spline of degree four with one free parameter, λ, which is called a shape parameter. This parameter is introduced within the basis function so as to change the shape of the spline. The authors in [13] generalized the extension to degrees five and six. Different researchers used extended cubic B-spline basis function to solve linear and nonlinear ordinary and partial differential equations; for example, see [14–16]. The authors in [17, 18] studied singularly perturbed semilinear differential equation of reaction-diffusion and convection-diffusion type, respectively, using cubic B-spline collocation method on a piecewise Shishkin mesh. Recently, cubic B-spline collocation method has been developed for time-dependent singularly perturbed differential-difference equations; see [19, 20]. An extended cubic B-spline collocation method has been developed for time-dependent singularly perturbed partial differential equations with time lag in [21] and singularly perturbed parabolic differential-difference equation arising in computational neuroscience in [22]. Recently, scholars in [23–26] studied time-dependent singularly perturbed parabolic partial differential equations.This article is structured as follows. The properties of continuous problem are given in Section2. In Section 3, the description of numerical scheme and bounds of error for time semidiscretization followed by spatial discretization using extended B-spline collocation method are discussed. Error analysis in spatial direction and the overall error bounds are given in Section 4. Some numerical computations are given in Section 5. The conclusion is given in Section 6.


## 2. Properties of Continuous Problem

In the study of the numerical aspects of singularly perturbed problems, their analytical aspects play an important role. Setting the valueε=0, the reduced problem corresponding to (1)–(3) is
(6)∂u0x,t∂t+ax,tu0x,t=−bx,tu0x,t−τ+fx,t,x,t∈D,ux,t=ϕbx,t,x,t∈Γb.The reduced problem (6) is an initial value problem which will not make use of the two boundary conditions. As a result, the solution of problem in (1)–(3) will have left and right boundary layers. The characteristic curve of the reduced problem in equation (6) is the vertical lines x=constant, which implies that boundary layers arising in the solution are of parabolic type. The problem in (1)–(3) satisfies the following continuous maximum principle.Theorem 1.
Assume thata∈C0,0D¯and letψ∈U∗=C2,1D∩C1,0D∗∩C0,0D¯be a sufficiently smooth function defined onDsuch thatLψx,t≥0,x,t∈D, Blψx,t≥0,x,t∈Γl, Brψx,t≥0,x,t∈Γr, andψx,t≥0,x,t∈Γ0whereD∗=D∪Γl∪Γr,Γ0=x,0:x∈0,1. Then,ψx,t≥0, for allx,t∈D¯.Proof.
The details of the proof are found in [11].Stability andε-uniform bound for problem (1) is established in the following theorem in the sense of the maximum norm which follows from Theorem 1.Theorem 2.
Letvx,tbe any function in the domain of problem. Then, we have the bound(7)∥v∥D¯≤max∥Blvx,t∥,∥Brvx,t∥,∥vx,t∥,1α∥Lvx,t∥.Proof.
See the details of the proof in [11].The existence and uniqueness for a solution of (1)–(3) can be established under the assumption that the data are Ho¨lder continuous and also satisfy an appropriate compatibility conditions at the corner points 0,0,1,0,0,−τ and 1,−τ. The boundary functions ϕl,ϕr∈Ck0,T,ϕb∈C1,kΓb are said to satisfy the kth order compatibility condition at the initial function if
(8)∂k∂tkϕb−ε∂ϕb∂x0,0=dkϕl0∂tk,∂k∂tkϕb+ε∂ϕb∂x1,0=dkϕr1∂tk.Therefore, the problem in (1)–(3) will have a unique solution which exhibits parabolic boundary layers at x=0 and x=1; see the details in [3, 11]. Now, we establish the classical bounds on the solution and its derivatives.Theorem 3.
Let the coefficients and source functiona,b,f∈C2+α,1+α/2D¯,ϕl,ϕr∈C3+α/20,T,ϕb∈C4+α,2+α/2Γb,α∈0,1. Assume that the compatibility conditions fork=0,1,2are fulfilled. Then, the problem has a unique solution and the derivatives of the solutionusatisfy the bound(9)∥∂i+ju∂xi∂tj∥D¯≤Cε−i2,i,j≥0,0≤i+2j≤4,where the constant C is independent of ε.Proof.
The proof of the first part is given in [27]. The bounds on the solution and its derivatives are explained in [11].The classical bounds in Theorem3 are not adequate for the proof of ε-uniform error estimate. Thus, the nonclassical bounds in singular and regular components and its derivatives are established in the following theorem.Theorem 4.
Let the coefficients and source functiona,b,f∈C4+α,2+α/2D¯,ϕl,ϕr∈C5+α/20,T,ϕb∈C6+α,3+α/2Γb,α∈0,1. Under the smoothness and compatibility conditions, we have the bounds fori,j≥0,0≤i+2j≤4(10)∂i+jvx,t∂xi∂tj≤C1+ε1−i/2,∂i+jwl∂xi∂tj≤Cε−i2exp−xε,∂i+jwr∂xi∂tj≤Cε−i2exp−1−xε.Proof.
The details of the proof are found in [11, 27].


## 3. Description of the Numerical Method

In this section, we utilize the implicit Euler method to discretize time derivative and then we introduce a piecewise uniform Shishkin mesh to discretize the space derivative using the extended form of cubic B-spline collocation method for the linear differential equations resulted from the time semidiscretization.


### 3.1. Time Semidiscretization

We discretize time derivative in (1)–(3) by means of the implicit Euler scheme on a uniform mesh with step length of Δt defined by
(11)ΩtM=tj:tj=t0+jΔt,j=1,⋯,M,t0=0,Δt=TM,where M denotes the number of mesh elements in temporal direction. Uniform meshes with step size Δt,ΩtM and Ωts with M and s mesh elements are used on the interval 0,T and −τ,0, respectively. The mesh size Δt is chosen in such a way that the delay parameter τ=sΔt, where s is a positive integer, tj=jΔt, j≥−s. We obtain the following system of ordinary differential equations
(12)Uj+1x=ϕbx,tj+1,x∈Ω,−s+1≤j≤0,Uj+1−UjΔt−εUxxj+1+aUj+1=−bj+1Uj−s+1+fj+1,subject to the boundary conditions
(13)Uj+10−εUxj+10=ϕltj+1,0≤j≤M,Uj+11+εUxj+11=ϕrtj+1,0≤j≤M,where Uj+1x=Ux,tj+1 and Uj+1 is the numerical solution at the j+1th time level. For each time step, equations (12)–(13) can be rewritten as
(14)−εUxxj+1+pUj+1=Rj+1,subject to the initial and boundary conditions, respectively,
(15)Uj+1x=ϕbx,tj+1,x∈Ω,−s+1≤j≤0,(16)Uj+10−εUxj+10=ϕltj+1,0≤j≤M,Uj+11+εUxj+11=ϕrtj+1,0≤j≤M.Here,pj+1=aj+1+1/Δt,pj+1≥α>0 and Rj+1=Uj/Δt−bj+1Uj−s+1+fj+1.By using the initial condition, we can evaluate the right-hand side as(17)Rj+1x=UjΔt−bj+1ϕbx,tj−s+1+fj+1,j=0,1,⋯,s,UjΔt−bj+1Ux,tj−s+1+fj+1,j=s+1,⋯,M.The local truncation error of an implicit Euler scheme for the temporal semidiscretization is given byej+1=ux,tj+1−Uj+1x. This error measures the contribution of each time step to the global error of the time semidiscretization.Lemma 5.
If(18)∣∂iux,t∂xi∣≤C,x,t∈D¯,0≤i≤2,then the local error bound in the temporal direction is given by
(19)∥ej+1∥∞≤CΔt2.Proof.
For the proof of the lemma, the readers can refer to [28].
The global error is the measure of the contribution of the local error estimate at each time step and is given byej=ux,tj−Ujx.Lemma 6.
Under the hypothesis of Lemma5, the global error estimate attjis given by(20)∥Ej∥∞≤CΔt,j≤TΔt.We conclude that time semidiscretization is first-order uniformly convergent.


### 3.2. Spatial Discretization

We first construct nonequidistant (layer adapted) Shishkin mesh as follows. We divide the three nonoverlapping subintervals0,σ, σ,1−σ, and 1−σ,1 into N/4, N/2, and N/4 equidistant subintervals. We define the transition parameter σ as σ=min1/4,2εlnN. Let Ω¯xN=xii=0N be the set of mesh points. Now, we define piecewise uniform mesh points as
(21)xi=ihi,ifi=01N4,σ+i−N4hi,ifi=N4+113N4,1−σ+i−3N4hi,ifi=3N4+11N,with piecewise uniform mesh spacing hi=4σ/N, if i=11N/4,i=3N/4+11N, and hi=21−2σ/N, if i=N/4+1,13N/4. Now, we apply the extended cubic B-spline collocation method to find the approximate solution for problem (14)–(16). Let Δ:0=x0<⋯<xN=1 be the spatial domain 0,1 with a piecewise uniform mesh spacing hi=xi+1−xi. The extended form of cubic B-spline of degree 4, Eix,λ,i=−1,0,1,⋯,N+1, is defined by [21, 29]
(22)Eix,λ=124hi44hi1−λx−xi−23+3λx−xi−24,x∈Ii,4−λhi4+12hi3x−xi−1+6hi22+λx−xi−12−12hix−xi−13−3λx−xi−14,x∈Ii+1,4−λhi4+12hi3xi+1−x+6hi22+λxi+1−x2−12hixi+1−x3−3λxi+1−x4,x∈Ii+2,4hi1−λxi+2−x3+3λxi+2−x4,x∈Ii+3,0,otherwise,where Ii=xi−2,xi−1 and −mm−2≤λ≤1 is a free parameter which is used to change the shape of the B-spline curve and m is the degree of extended cubic B-spline. The variation in m gives different forms of extended cubic B-spline functions [29]. The extended cubic B-spline function has one free parameter λ, when the free parameter λ tends to zero the extended cubic B-spline reduced to convectional cubic B-spline functions. For λ∈−8,1, cubic B-spline and extended cubic B-spline share the same properties such as local support, nonnegativity, partition of unity, and C2 continuity; the parameter λ controls the tension of the solution curve. The shape of extended cubic B-spline functions forces us to add two fictitious points E−1 and EN+1 to satisfy the boundary conditions. Since B-splines of degree p are p−1 times continuously differentiable piecewise polynomials that form a basis of the space of splines, let φ3Ω be the space of twice continuously differentiable piecewise extended cubic B-spline on Ω. Since each Eix is also a piecewise cubic with knots at Ω, each Eix∈φ3Ω. Suppose that E3Ω=span E−1,E0,⋯,EN,EN+1. Since the functions Ei′s are linearly independent on 0,1, E3Ω is an N+3-dimensional. Let Sx,λ be the B-spline interpolating function for ux,t at the nodal points and Sx,λ∈E3Ω. Therefore, we seek an approximate solution Sx,λ of the problem (14)–(16) which is given by
(23)Sx,λ=∑i=−1N+1γiEix,λ,where γi are unknown real coefficients to be determined by requiring that Sx,λ satisfies (14)–(16) at N+1 collocation points and boundary conditions. The values of extended B-splines Eix,λ and its derivatives at the nodal points can be calculated from (23) and depicted in Table 1.Table 1
Values ofEix,λ and its derivatives at nodal points.
xi−2xi−1xixi+1xi+2Eix,λ04−λ248+λ124−λ240Ei′x,λ012hi0−12hi0Ei′′x,λ02+λ2hi2−2+λhi22+λ2hi20An approximate solution over typical subintervalxi,xi+1 can be defined as
(24)Sx,λ=∑j=i−1i+2γjEjx,λ.Now, substituting the values ofEi for Siλ and Ei′′ for Si′′λ as stated in Table 1 in equation (14), we get N+1 linear equations in N+3 unknowns as
(25)ri−γi−1+ricγi+ri+γi+1=R~i,fori=0,1,⋯,N,where the coefficients are given by
(26)ri−=−ε2+λ2hi2+4−λ24p~i,ric=ε2+λhi2+8+λ12p~i,ri+=−ε2+λ2hi2+4−λ24p~i,p~i=pij+1,R~i=Rij+1.Boundary conditions in (16) at x0 and xN must be imposed to the system of equations in (25) to obtain the unique solution. Thus, we obtain the approximate solution at two boundary points as follows:
(27)a1γ−1+b1γ0+c1γ1=ϕltj+1,c1γN−1+b1γN+a1γN+1=ϕrtj+1,where the coefficients are given by a1=4−λ/24−ε/2hi,b1=8+λ/12,c1=4−λ/24+ε/2hi. Equations (25) and (27) lead to N+3×N+3 linear systems with N+3 unknowns γ−1,γ0,⋯,γN+1. Excluding the unknowns γ−1 and γN+1 from (27) for i=0 and i=N, then (25) becomes solvable system of N+1×N+1 linear equations in N+1 unknowns γ0,⋯,γN in matrix form as
(28)Mγ=G,i=1,2,⋯,N−1,where the entries of the tridiagonal matrix M=mij are given by
(29)mij=r0c−b1a1r0−,i=j=0,r0+−c1a1r0−,i=0,j=1ri−,i=j+1,j=01N−2,ric,i=j=11N−1,ri+,i=j−1,j=21N,rN−−c1a1rN+,i=N,j=N−1,rNc−b1a1rN+,i=j=N,0,∀∣i−j∣>0.The entries of the column vectorG are given as
(30)G=R~0j+1−r0−a1ϕltj+1,j=0,R~ij+1,j=11N−1,R~Nj+1−rN+a1ϕrtj+1,j=N.Sincepij+1>0, it is easily seen that for λ>−2 the matrix M is strictly diagonally dominant and hence nonsingular. Since M is nonsingular, we can solve the system (28) for γ0,⋯,γN. Hence, the extended cubic B-spline collocation method applied to problem (14) has a unique solution Sx,λ.


## 3.1. Time Semidiscretization

We discretize time derivative in (1)–(3) by means of the implicit Euler scheme on a uniform mesh with step length of Δt defined by
(11)ΩtM=tj:tj=t0+jΔt,j=1,⋯,M,t0=0,Δt=TM,where M denotes the number of mesh elements in temporal direction. Uniform meshes with step size Δt,ΩtM and Ωts with M and s mesh elements are used on the interval 0,T and −τ,0, respectively. The mesh size Δt is chosen in such a way that the delay parameter τ=sΔt, where s is a positive integer, tj=jΔt, j≥−s. We obtain the following system of ordinary differential equations
(12)Uj+1x=ϕbx,tj+1,x∈Ω,−s+1≤j≤0,Uj+1−UjΔt−εUxxj+1+aUj+1=−bj+1Uj−s+1+fj+1,subject to the boundary conditions
(13)Uj+10−εUxj+10=ϕltj+1,0≤j≤M,Uj+11+εUxj+11=ϕrtj+1,0≤j≤M,where Uj+1x=Ux,tj+1 and Uj+1 is the numerical solution at the j+1th time level. For each time step, equations (12)–(13) can be rewritten as
(14)−εUxxj+1+pUj+1=Rj+1,subject to the initial and boundary conditions, respectively,
(15)Uj+1x=ϕbx,tj+1,x∈Ω,−s+1≤j≤0,(16)Uj+10−εUxj+10=ϕltj+1,0≤j≤M,Uj+11+εUxj+11=ϕrtj+1,0≤j≤M.Here,pj+1=aj+1+1/Δt,pj+1≥α>0 and Rj+1=Uj/Δt−bj+1Uj−s+1+fj+1.By using the initial condition, we can evaluate the right-hand side as(17)Rj+1x=UjΔt−bj+1ϕbx,tj−s+1+fj+1,j=0,1,⋯,s,UjΔt−bj+1Ux,tj−s+1+fj+1,j=s+1,⋯,M.The local truncation error of an implicit Euler scheme for the temporal semidiscretization is given byej+1=ux,tj+1−Uj+1x. This error measures the contribution of each time step to the global error of the time semidiscretization.Lemma 5.
If(18)∣∂iux,t∂xi∣≤C,x,t∈D¯,0≤i≤2,then the local error bound in the temporal direction is given by
(19)∥ej+1∥∞≤CΔt2.Proof.
For the proof of the lemma, the readers can refer to [28].
The global error is the measure of the contribution of the local error estimate at each time step and is given byej=ux,tj−Ujx.Lemma 6.
Under the hypothesis of Lemma5, the global error estimate attjis given by(20)∥Ej∥∞≤CΔt,j≤TΔt.We conclude that time semidiscretization is first-order uniformly convergent.


## 3.2. Spatial Discretization

We first construct nonequidistant (layer adapted) Shishkin mesh as follows. We divide the three nonoverlapping subintervals0,σ, σ,1−σ, and 1−σ,1 into N/4, N/2, and N/4 equidistant subintervals. We define the transition parameter σ as σ=min1/4,2εlnN. Let Ω¯xN=xii=0N be the set of mesh points. Now, we define piecewise uniform mesh points as
(21)xi=ihi,ifi=01N4,σ+i−N4hi,ifi=N4+113N4,1−σ+i−3N4hi,ifi=3N4+11N,with piecewise uniform mesh spacing hi=4σ/N, if i=11N/4,i=3N/4+11N, and hi=21−2σ/N, if i=N/4+1,13N/4. Now, we apply the extended cubic B-spline collocation method to find the approximate solution for problem (14)–(16). Let Δ:0=x0<⋯<xN=1 be the spatial domain 0,1 with a piecewise uniform mesh spacing hi=xi+1−xi. The extended form of cubic B-spline of degree 4, Eix,λ,i=−1,0,1,⋯,N+1, is defined by [21, 29]
(22)Eix,λ=124hi44hi1−λx−xi−23+3λx−xi−24,x∈Ii,4−λhi4+12hi3x−xi−1+6hi22+λx−xi−12−12hix−xi−13−3λx−xi−14,x∈Ii+1,4−λhi4+12hi3xi+1−x+6hi22+λxi+1−x2−12hixi+1−x3−3λxi+1−x4,x∈Ii+2,4hi1−λxi+2−x3+3λxi+2−x4,x∈Ii+3,0,otherwise,where Ii=xi−2,xi−1 and −mm−2≤λ≤1 is a free parameter which is used to change the shape of the B-spline curve and m is the degree of extended cubic B-spline. The variation in m gives different forms of extended cubic B-spline functions [29]. The extended cubic B-spline function has one free parameter λ, when the free parameter λ tends to zero the extended cubic B-spline reduced to convectional cubic B-spline functions. For λ∈−8,1, cubic B-spline and extended cubic B-spline share the same properties such as local support, nonnegativity, partition of unity, and C2 continuity; the parameter λ controls the tension of the solution curve. The shape of extended cubic B-spline functions forces us to add two fictitious points E−1 and EN+1 to satisfy the boundary conditions. Since B-splines of degree p are p−1 times continuously differentiable piecewise polynomials that form a basis of the space of splines, let φ3Ω be the space of twice continuously differentiable piecewise extended cubic B-spline on Ω. Since each Eix is also a piecewise cubic with knots at Ω, each Eix∈φ3Ω. Suppose that E3Ω=span E−1,E0,⋯,EN,EN+1. Since the functions Ei′s are linearly independent on 0,1, E3Ω is an N+3-dimensional. Let Sx,λ be the B-spline interpolating function for ux,t at the nodal points and Sx,λ∈E3Ω. Therefore, we seek an approximate solution Sx,λ of the problem (14)–(16) which is given by
(23)Sx,λ=∑i=−1N+1γiEix,λ,where γi are unknown real coefficients to be determined by requiring that Sx,λ satisfies (14)–(16) at N+1 collocation points and boundary conditions. The values of extended B-splines Eix,λ and its derivatives at the nodal points can be calculated from (23) and depicted in Table 1.Table 1
Values ofEix,λ and its derivatives at nodal points.
xi−2xi−1xixi+1xi+2Eix,λ04−λ248+λ124−λ240Ei′x,λ012hi0−12hi0Ei′′x,λ02+λ2hi2−2+λhi22+λ2hi20An approximate solution over typical subintervalxi,xi+1 can be defined as
(24)Sx,λ=∑j=i−1i+2γjEjx,λ.Now, substituting the values ofEi for Siλ and Ei′′ for Si′′λ as stated in Table 1 in equation (14), we get N+1 linear equations in N+3 unknowns as
(25)ri−γi−1+ricγi+ri+γi+1=R~i,fori=0,1,⋯,N,where the coefficients are given by
(26)ri−=−ε2+λ2hi2+4−λ24p~i,ric=ε2+λhi2+8+λ12p~i,ri+=−ε2+λ2hi2+4−λ24p~i,p~i=pij+1,R~i=Rij+1.Boundary conditions in (16) at x0 and xN must be imposed to the system of equations in (25) to obtain the unique solution. Thus, we obtain the approximate solution at two boundary points as follows:
(27)a1γ−1+b1γ0+c1γ1=ϕltj+1,c1γN−1+b1γN+a1γN+1=ϕrtj+1,where the coefficients are given by a1=4−λ/24−ε/2hi,b1=8+λ/12,c1=4−λ/24+ε/2hi. Equations (25) and (27) lead to N+3×N+3 linear systems with N+3 unknowns γ−1,γ0,⋯,γN+1. Excluding the unknowns γ−1 and γN+1 from (27) for i=0 and i=N, then (25) becomes solvable system of N+1×N+1 linear equations in N+1 unknowns γ0,⋯,γN in matrix form as
(28)Mγ=G,i=1,2,⋯,N−1,where the entries of the tridiagonal matrix M=mij are given by
(29)mij=r0c−b1a1r0−,i=j=0,r0+−c1a1r0−,i=0,j=1ri−,i=j+1,j=01N−2,ric,i=j=11N−1,ri+,i=j−1,j=21N,rN−−c1a1rN+,i=N,j=N−1,rNc−b1a1rN+,i=j=N,0,∀∣i−j∣>0.The entries of the column vectorG are given as
(30)G=R~0j+1−r0−a1ϕltj+1,j=0,R~ij+1,j=11N−1,R~Nj+1−rN+a1ϕrtj+1,j=N.Sincepij+1>0, it is easily seen that for λ>−2 the matrix M is strictly diagonally dominant and hence nonsingular. Since M is nonsingular, we can solve the system (28) for γ0,⋯,γN. Hence, the extended cubic B-spline collocation method applied to problem (14) has a unique solution Sx,λ.


## 4. Error Analysis

This section proves theε-uniformly convergence of the proposed method in the spatial direction. For this, we use the following lemma.Lemma 7.
The extended cubic B-splinesE−1x,λ,E0x,λ,⋯,ENx,λ,EN+1x,λdefined in (22) satisfy the following inequality(31)∑i=−1N+1∣Eix,λ∣≤74,x∈0,1.Proof.
We know that(32)∣∑i=−1N+1Eix,λ∣≤∑i=−1N+1Eix,λ.Extended cubic B-splineEix,λ is nonzero at only three nodal points. Thus, at any nodal value xi, from Table 1 we obtain
(33)∑i=−1N+1∣Eix,λ∣=∣Ei−1xi,λ∣+∣Eixi,λ∣+∣Ei+1xi,λ∣=4−λ24+8+λ12+4−λ24=1<74.From Table1, for x∈xi−1,xi, we have
(34)∣Eixi,λ∣≤8+λ12,∣Ei−1xi−1,λ∣≤8+λ12.Similarly, forx∈xi−1,xi, we have that
(35)∣Ei+1xi,λ∣≤4−λ24,∣Ei−2xi−1,λ∣≤4−λ24.Now, for any point forx∈xi−1,xi, we get
(36)∑i=−1N+1∣Eix,λ∣=∣Ei−1x,λ∣+∣Eix,λ∣+∣Ei+1x,λ∣+∣Ei−2x,λ∣=20+λ12.Since−8<λ<1, thus 20+λ/12≤7/4 and this completes the proof.LetYx be the unique cubic spline interpolate from an approximate solution Sx,λ of the problem (14)–(16) to the solution ux,t which is given by
(37)Yx=∑i=−1N+1γ¯iEix,λ.Lemma 8.
LetYx∈C20,1be the cubic spline interpolant associated with a solutionu~x. Ifu~x∈C40,1, it follows from the estimate of Hall [30] that the standard cubic spline interpolation error estimate holds, forx∈Ωi≔xi,xi+1∈Ω(38)∥Dnu^x−Yx∥≤λn∥u∧4x∥hi4−n,n=0,1,2,3,where λn are constants independent of hi and N.Lemma 9.
[31]. If matrixAis strictly diagonally dominant by rows and constantζ=mini∣ai,i∣−∑i≠j∣ai,j∣. Then, we have the bound(39)∥A−1∥∞<1ζ.Theorem 10.
LetSx,λbe an extended cubic B-spline collocation approximation from the space of extended B-splinesφ3Ωto the solution of (14)–(16) andu^xiis the analytical solution to the problem. Ifg~x∈C20,1, then the parameter-uniform error estimate satisfies the bound(40)sup0<ε≤1max0≤i≤N∥u^xi−Sxi,λ∥≤CN−2lnN2,where C is a constant independent of ε and N.Proof.
To estimate the error∣u^xi−Sxi,λ∣, it follows immediately from the estimates in Lemma 8 and equation (37) that
(41)∣L^u^xi−L^Yxi∣≤εu∧′′xi−Y′′xi+pxiu^xi−Yxi∣≤ελ2u∧4hi2+∥p∥λ0u∧4hi4,
from (37).Therefore, rearrangement of the above expression gives(42)L^u^xi−L^Yxi≤ελ2hi2+∥p∥∞λ0hi4u∧4.Since the argument depends on whetherσ=1/4 or σ=2εlnN<1/4, there arises two cases.Case 1.
Whenσ=1/4, the mesh is uniform with spacing 1/N, that is, hi=1/N and 2εlnN≥1/4 gives ε−1/2≤ClnN. From this, we get ε−1≤ClnN2. In this case, we use a classical analysis to prove convergence. Using the classical bound in Theorem 3, that is, ∣u∧4∣≤Cε−2 together with equation (42) yields
(43)∣L^u^xi−L^Yxi∣≤Cε−2ελ2N−2+∥p∥∞λ0N−4≤CN−2ε−1+ε−2N−2≤CN−2ClnN2+CN−2lnN4.SinceCN−2lnN4≤ClnN2, we obtain the following estimate:
(44)∣L^u^xi−L^Yxi∣≤CN−2lnN2.Case 2.
IfΩi lies in the boundary layer regions, then mesh spacing hi≤Cε1/2N−1lnN. Using the bound in the layer regions together with the estimates in equation (42), we have
(45)∣L^u^xi−L^Yxi∣≤Cε−2ελ2C2εN−2ln2N+∥p∥∞λ0C4ε2N−4ln4N≤CN−2lnN2+N−2lnN4.SinceN−2lnN4≤lnN2, we obtain the following estimate:
(46)∣L^u^xi−L^Yxi∣≤CN−2lnN2.On the other hand, for the subintervalσ,1−σ, that is, for the outer region, the mesh spacing is hi=2N−11−2σ=2N−1−CεN−1lnN≤Cε1/2N−1lnN. Using this in equation (42) together with the bounds gives us
(47)∣L^u^xi−L^Yxi∣≤Cε−2ελ2C2εN−2ln2N+∥p∥∞λ0C4ε2N−4ln4N≤CN−2lnN2+N−2lnN4.SinceN−2lnN4 is very small number, we obtain the following estimate:
(48)∣L^u^xi−L^Yxi∣≤CN−2lnN2.Combining the above estimates for both the cases, we have(49)∣L^u^xi−L^Yxi∣≤CN−2lnN2.Therefore, we have that(50)∣L^Sxi−L^Yxi∣=∣R^xi−L^Yxi∣=∣L^u^xi−L^Yxi∣≤CN−2lnN2.We know thatL^Uxi=R^xi,0≤i≤N with the boundary conditions Ux0−εU′x0=ϕltj+1 and UxN+εU′xN=ϕrtj+1 leads to the linear system Mγ=G. Assume that L^Yxi=R¯xi,0≤i≤N with the boundary conditions Yx0−εY′x0=ϕ¯ltj+1 and YxN+εY′xN=ϕ¯rtj+1 leads to the linear system Mγ¯=G¯. It follows that
(51)Mγ−γ¯=G−G¯,where γ−γ¯=γ0−γ¯0,γ1−γ¯1,⋯,γN−γ¯NT and
(52)G−G¯=R~0−R¯0+r0−a1ϕ¯ltj+1−ϕltj+1R~1−R¯1⋮R~N−1−R¯N−1R~N−R¯N+rN+a1ϕ¯rtj+1−ϕrtj+1.SinceMγ−Mγ¯=G−G¯ implies that L^Sxi−L^Yxi, from (50) we have
(53)∥G−G¯∥≤CN−2lnN2.It can be seen that for reasonable largeN, the matrix M is strictly diagonally dominant and thus nonsingular. From the estimate in Lemma 9, we get
(54)∥M−1∥∞≤C.Combining the bounds in equations (51)–(54), we obtain
(55)∥γ−γ¯∥≤CN−2lnN2.Lete=e0,⋯,eNT, where ei=γi−γ¯i. Now, from (51), we have
(56)e=M−1G−G¯.Using (53) and (54) in (56), we have the following estimate:
(57)∥e∥≤CN−2lnN2.We have(58)a1γ−1−γ¯−1+b1γ0−γ¯0+c1γ1−γ¯1=ϕltj+1,c1γN−1−γ¯N−1+b1γN−γ¯N+a1γN+1−γ¯N+1=ϕrtj+1,where a1,b1, and c1 are defined in (27). From this, it is a simple task to obtain ∣γ−1−γ¯−1∣≤CN−2lnN2 and ∣γN+1−γ¯N+1∣≤CN−2lnN2. Therefore, we have the following estimation from boundary conditions as
(59)max−1≤i≤N+1∣γi−γ¯i∣≤CN−2lnN2.Therefore, the above inequality enables us to estimate∣Sx,λ−Yx∣ as
(60)∣Sx,λ−Yx∣=∑i=−1N+1γi−γ¯iEix,λ.Using (56) and Lemma 7, we obtain
(61)max0≤i≤N∣Sxi,λ−Yxi∣≤CN−2lnN2.Using triangle inequality and the results in (49) and (61) gives
(62)sup0<ε≤1max0≤i≤N∣u^xi−Sxi,λ∣≤CN−2lnN2.Hence, this completes the proof.Theorem 11.
LetSxi,λbe the extended B-spline collocation approximation to the solutionux,tof the problem (14)–(16). Then, the parameter-uniform error estimate of fully discrete scheme is given by(63)uxi−Sxi,λ≤CΔt+N−2ln2N,where C is a constant independent of mesh parameters and ε.Proof.
The result of Lemma6 and Theorem 10 proves this theorem.


## 5. Numerical Results and Discussions

To see the applicability and efficiency of the proposed method, an example is considered from the literature. Computations are done for reasonable value of the free parameterλ=0.99∈−8,1, which gives minimum error. Since the exact solution for the test example is unknown, we use the double mesh principle to calculate absolute errors. For each ε, we can determine the maximum point-wise errors using the formula as
(64)EεN,Δt=max0≤i≤N;t∈0,TUN,Δtxi,tj−U2N,2Mx2i,t2j,where UN,Δtxi,tj denote the numerical solution obtained at N,Δt mesh points whereas U2N,Δt/2x2i,t2j denote the numerical solution at 2N,2M mesh points. The uniform rate of convergence is calculated by the formula
(65)RN,Δt=log2EN,ΔtE2N.Δt/2.Example 1.
Consider the following singularly perturbed delay problem [11].
(66)ut−εuxx+1+x22ux,t=−ux,t−1+t3,x,t∈0,1×0,2,ux,t=0,x,t∈0,1×−1,0,u0,t−εux0,t=−128/35π−1/2t7/2,t∈0,2,u1,t+εux1,t=−128/35π−1/2t7/2,t∈0,2.The computed maximum point-wise errorsEεN,Δt are given in Table 2. From this result, it is clear that the proposed method gives an ε-uniform convergence.Table 2
The maximum point-wise errorsEεN,Δt for various ε and λ=0.99.
ε↓N=32, Δt=0.1/2264,0.1/23128,0.1/24256,0.1/25512,0.1/261024,0.1/272−69.7604e−24.9405e−22.9193e−21.7022e−24.4560e−32.1246e−32−79.7647e−24.9411e−22.4853e−21.4146e−24.4571e−32.1266e−32−89.7679e−24.9415e−22.4853e−21.2463e−25.1580e−32.1270e−32−91.1064e−14.9417e−22.4854e−21.2463e−25.1580e−32.1270e−32−101.1980e−15.5853e−22.4854e−21.2463e−25.1580e−32.1270e−32−111.1990e−15.6763e−22.6206e−21.2463e−25.1580e−32.1270e−32−121.1990e−15.6763e−22.6206e−21.2463e−25.1580e−32.1270e−32−131.1990e−15.6763e−22.6206e−21.2463e−25.1580e−32.1270e−32−141.1990e−15.6763e−22.6206e−21.2463e−25.1580e−32.1270e−32−151.1990e−15.6763e−22.6206e−21.2463e−25.1580e−32.1270e−32−161.1990e−15.6763e−22.6206e−21.2463e−25.1580e−32.1270e−3⋮⋮⋮⋮⋮⋮⋮2−281.1990e−15.6763e−22.6206e−21.2463e−25.1580e−32.1270e−32−291.1990e−15.6763e−22.6206e−21.2463e−25.1580e−32.1270e−32−301.1990e−15.6763e−22.6206e−21.2463e−25.1580e−32.1270e−3Computational results in Table3 confirm that the present method has improved the finite difference method in the literature.Table 3
Comparison ofε-uniform values for λ=0.99 using the proposed method and finite difference method.
∀ε↓N=32, Δt=0.1/2264,0.1/23128,0.1/24256,0.1/25512,0.1/26Our resultEN,Δt1.1990e−15.6763e−22.6206e−21.2463e−25.1580e−3RN,Δt1.07881.11511.17221.27281.2780Result in [11]EN,Δt2.3070e−11.2002e−16.1238e−23.0931e−21.5545e−2RN,Δt0.94270.97090.98530.99260.9963Table4 displays the comparison of computational results using classical cubic B-spline method for λ=0 and extended B-spline method for λ=0.99. Both the particular selections of the extension parameter give numerical results, but the later value of λ gives better result than the former one as well as improved results of finite difference method in the literature. It is clear from Table 4 that extended B-spline performs good result than classical cubic B-spline.Table 4
Comparison ofε-uniform values for λ=0 and λ=0.99 using the proposed method.
∀ε↓N=32, Δt=0.1/2264,0.1/23128,0.1/24256,0.1/25512,0.1/26Forλ=0EN,Δt1.5717e−17.3329e−23.3118e−21.4555e−26.9985e−3RN,Δt1.09991.14681.18611.22321.2216Forλ=0.99EN,Δt1.1990e−15.6763e−22.6206e−21.2463e−25.1580e−3RN,Δt1.07881.11511.17221.27281.2780Figure1 depicts the numerical simulation of solution profile at N=27,Δt=0.1/24,ε=2−12, and λ=−0.55, which indicates parabolic boundary layers at x=0 and x=1.Figure 1
Numerical simulation atN=27,Δt=0.1/24,ε=2−12,λ=−0.55.The maximum point-wise error values are plotted using log-log scale as can be seen in Figure2.Figure 2
Plot of maximum-uniform errors in log-log scale.To compare the numerical simulation for both classical B-spline and extended B-spline, solution profile is plotted atN=27,Δt=0.1/24,ε=2−12, and λ=0 in Figure 3, which shows the solution profile of classical B-spline.Figure 3
Numerical simulation atN=27,Δt=0.1/24,ε=2−12,λ=0.Figure4 is the numerical simulation using line graph at N=24,Δt=0.1/2,ε=2−12, and λ=−0.55, which indicates parabolic-type boundary layers at x=0 and x=1.Figure 4
Line graph atN=24,Δt=0.1/2,ε=2−12,λ=−0.55.


## 6. Conclusion

This article proposed an implicit Euler method for time derivative with uniform mesh and extended cubic B-spline collocation method for space derivative on Shishkin mesh to solve singularly perturbed delay parabolic reaction-diffusion problem subject to mixed boundary conditions. The proposed method is shown to be accurate ofOΔt+N−2ln2N by preserving an ε-uniform convergence. The suitable choice of the extension parameter λ minimizes the error. To see the effect of delay and extension parameter λ in the boundary regions, graphs are plotted on the solution profile for λ=−0.55∈−8,1.



---

*Source: 8835595-2021-03-13.xml*