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

In this article, a numerical solution is proposed for singularly perturbed delay parabolic reaction-diffusion problem with mixedtype 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−2 ln2NÞ, 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.


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 = f0g × ð0, T and Γ r =f1g × ð0, T. In this work, we consider the following singularly perturbed delay parabolic reaction-diffusion problem subject to the initial condition and boundary conditions 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 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] where Ω ⊂ ℝ is a bounded domain and ðx, tÞ ∈ Ω × ð0,∞Þ. The state variable pðx, 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][4][5][6][7][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][15][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 differentialdifference equation arising in computational neuroscience in [22]. Recently, scholars in [23][24][25][26] studied timedependent singularly perturbed parabolic partial differential equations.
This article is structured as follows. The properties of continuous problem are given in Section 2. 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.
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.
Let vðx, tÞ be any function in the domain of problem. Then, we have the bound 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 H€ older continuous and also satisfy an appropriate compatibility conditions at the corner points ð0, 0, 1, 0, 0,−τÞ and ð1,−τÞ. The boundary functions ϕ l , ϕ r ∈ C k ð½0, TÞ, ϕ b ∈ C 1,k ðΓ b Þ are said to satisfy the k th order compatibility condition at the initial function if 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 function a, b, f ∈ C ð2+α,1+α/2Þ ð DÞ, ϕ l , ϕ r ∈ C ð3+αÞ/2 ½ð0, TÞ, ϕ b ∈ C ð4+α,2+α/2Þ ðΓ b Þ, α ∈ ð0, 1Þ. Assume that the compatibility conditions for k = 0, 1, 2 are fulfilled. Then, the problem has a unique solution and the derivatives of the solution u satisfy the bound 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 Theorem 3 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 function a, b, f ∈ C ð4+α,2+α/2Þ ð DÞ, ϕ l , ϕ r ∈ C ð5+αÞ/2 ð½0, TÞ, ϕ b ∈ C ð6+α,3+α/2Þ ðΓ b Þ, α ∈ ð0, 1Þ. Under the smoothness and compatibility conditions, we have the bounds for i, Proof. The details of the proof are found in [11,27].

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.

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 where M denotes the number of mesh elements in temporal direction. Uniform meshes with step size Δt, Ω M t and Ω s t 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, t j = jΔt, j ≥ −s. We obtain the following system of ordinary differential equations subject to the boundary conditions where U j+1 ðxÞ = Uðx, t j+1 Þ and U j+1 is the numerical solution at the ðj + 1Þth time level. For each time step, equations (12)-(13) can be rewritten as subject to the initial and boundary conditions, respectively, Here, p j+1 = a j+1 + ð1/ΔtÞ, p j+1 ≥ α > 0 and Abstract and Applied Analysis By using the initial condition, we can evaluate the righthand side as The local truncation error of an implicit Euler scheme for the temporal semidiscretization is given by e j+1 = uðx, t j+1 Þ − U j+1 ðxÞ. This error measures the contribution of each time step to the global error of the time semidiscretization.
then the local error bound in the temporal direction is given by 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 by e j = uð x, t j Þ − U j ðxÞ.

Lemma 6.
Under the hypothesis of Lemma5, the global error estimate at t j is given by We conclude that time semidiscretization is first-order uniformly convergent.

Spatial Discretization.
We first construct nonequidistant (layer adapted) Shishkin mesh as follows. We divide the three nonoverlapping subintervals ½0, σ, ðσ, 1 − σÞ, and ½1 − σ, 1 into N/4, N/2, and N/4 equidistant subintervals. We define the transition parameter σ as σ = min f1/4, 2 be the set of mesh points. Now, we define piecewise uniform mesh points as with piecewise uniform mesh spacing h i = 4σ/N, if i = 1ð1Þð N/4Þ, i = ðð3N/4Þ + 1Þð1ÞN, and h i = ð2ð1 − 2σÞÞ/N, if i = ðð N/4Þ + 1, 1Þð3N/4Þ. Now, we apply the extended cubic Bspline collocation method to find the approximate solution for problem (14)- (16). Let Δ : 0 = x 0 < ⋯ < x N = 1 be the spatial domain ½0, 1 with a piecewise uniform mesh spacing where I i = ½x i−2 , x i−1 and −mðm − 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 Bspline 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 C 2 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 E N+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 Bspline on Ω. Since each E i ðxÞ is also a piecewise cubic with knots at Ω, each E i ðxÞ ∈ φ 3 ðΩÞ. Suppose that E 3 ðΩÞ =span Since the functions E i ′ s are linearly independent on ½0, 1, E 3 ðΩÞ is an ðN + 3Þ-dimensional. Let Sðx, λÞ be the B-spline interpolating function for uðx, tÞ at the nodal points and Sðx, λÞ ∈ E 3 ðΩÞ. Therefore, we seek an approximate solution Sðx, λÞ of the problem (14)-(16) which is given by where γ i are unknown real coefficients to be determined by requiring that Sðx, λÞ satisfies (14)- (16) at N + 1 collocation points and boundary conditions. The values of extended Bsplines E i ðx, λÞ and its derivatives at the nodal points can be calculated from (23) and depicted in Table 1.

Abstract and Applied Analysis
An approximate solution over typical subinterval ½x i , x i+1 can be defined as Now, substituting the values of E i for S i ðλÞ and E i ′ ′ for S i ′ ′ðλÞ as stated in Table 1 in equation (14), we get ðN + 1Þ linear equations in ðN + 3Þ unknowns as where the coefficients are given by Boundary conditions in (16) at x 0 and x N 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: where the coefficients are given by a 1 The entries of the column vector G are given as Since p j+1 i > 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 Sðx, λÞ.

Error Analysis
This section proves the ε-uniformly convergence of the proposed method in the spatial direction. For this, we use the following lemma.
Proof. We know that Extended cubic B-spline E i ðx, λÞ is nonzero at only three nodal points. Thus, at any nodal value x i , from Table 1 we obtain Table 1: Values of E i ðx, λÞ and its derivatives at nodal points.
Let YðxÞ be the unique cubic spline interpolate from an approximate solution Sðx, λÞ of the problem (14)- (16) to the solution uðx, tÞ which is given by Lemma 8. Let YðxÞ ∈ C 2 ð½0, 1Þ be the cubic spline interpolant associated with a solutionũðxÞ. IfũðxÞ ∈ C 4 ð½0, 1Þ, it follows from the estimate of Hall [30] that the standard cubic spline interpolation error estimate holds, for where λ n are constants independent of h i and N. Lemma 9. [31]. If matrix A is strictly diagonally dominant by rows and constant ζ = min Then, we have the bound Theorem 10. Let Sðx, λÞ be an extended cubic B-spline collocation approximation from the space of extended B-splines φ 3 ðΩÞ to the solution of (14)- (16) andûðx i Þ is the analytical solution to the problem. IfgðxÞ ∈ C 2 ð½0, 1Þ , then the parameter-uniform error estimate satisfies the bound where C is a constant independent of ε and N.
Proof. To estimate the error |ûðx i Þ − Sðx i , λÞ | , it follows immediately from the estimates in Lemma 8 and equation (37) that from (37).
Case 1. When σ = 1/4, the mesh is uniform with spacing 1/N, that is, h i = 1/N and 2 ffiffi ε p ln N ≥ 1/4 gives ε −1/2 ≤ C ln N. From this, we get ε −1 ≤ ðC ln NÞ 2 . 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 Since CN −2 ðln NÞ 4 ≤ Cðln NÞ 2 , we obtain the following estimate: Case 2. If Ω i lies in the boundary layer regions, then mesh spacing h i ≤ Cε 1/2 N −1 ln N. Using the bound in the layer regions together with the estimates in equation (42), we have Since N −2 ðln NÞ 4 ≤ ðln NÞ 2 , we obtain the following estimate: Abstract and Applied Analysis On the other hand, for the subinterval ½σ, 1 − σ, that is, for the outer region, the mesh spacing is Using this in equation (42) together with the bounds gives us Since N −2 ðln NÞ 4 is very small number, we obtain the following estimate: Combining the above estimates for both the cases, we have We where Since It can be seen that for reasonable large N, the matrix M is strictly diagonally dominant and thus nonsingular. From the estimate in Lemma 9, we get Combining the bounds in equations (51)-(54), we obtain Let e = ðe 0 ,⋯,e N Þ T , where e i = γ i − γ i . Now, from (51), we have Using (53) and (54) in (56), we have the following estimate: We have where a 1 , b 1 , and c 1 are defined in (27). From this, it is a simple task to obtain |γ −1 − γ −1 | ≤CN −2 ðln NÞ 2 and |γ N+1 − γ N+1 | ≤CN −2 ðln NÞ 2 . Therefore, we have the following estimation from boundary conditions as Therefore, the above inequality enables us to estimate |S ðx, λÞ − YðxÞ | as  Abstract and Applied Analysis Using (56) and Lemma 7, we obtain Using triangle inequality and the results in (49) and (61) gives Hence, this completes the proof.
Theorem 11. Let Sðx i , λÞ be the extended B-spline collocation approximation to the solution uðx, tÞ of the problem (14)- (16). Then, the parameter-uniform error estimate of fully discrete scheme is given by where C is a constant independent of mesh parameters and ε.
Proof. The result of Lemma 6 and Theorem 10 proves this theorem.

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 where U N,Δt ðx i , t j Þ denote the numerical solution obtained at ðN, ΔtÞ mesh points whereas U 2N,Δt/2 ðx 2i , t 2j Þ denote the numerical solution at ð2N, 2MÞ mesh points. The uniform rate of convergence is calculated by the formula Example 1. Consider the following singularly perturbed delay problem [11].
The computed maximum point-wise errors E N,Δt ε are given in Table 2. From this result, it is clear that the proposed method gives an ε-uniform convergence.
Computational results in Table 3 confirm that the present method has improved the finite difference method in the literature. Table 4 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. Figure 1 depicts the numerical simulation of solution profile at N = 2 7 , Δt = 0:1/2 4 , ε = 2 −12 , and λ = −0:55, which indicates parabolic boundary layers at x = 0 and x = 1.
The maximum point-wise error values are plotted using log-log scale as can be seen in Figure 2.

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 reactiondiffusion problem subject to mixed boundary conditions. The proposed method is shown to be accurate of OðΔt + N −2 ln 2 NÞ 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.

Data Availability
No data were used to support the findings of this study.

Conflicts of Interest
The authors declare that there is no conflict of interest regarding the publication of this manuscript.

10
Abstract and Applied Analysis