Hybrid Fitted Numerical Scheme for Singularly Perturbed Convection-Diffusion Problem with a Small Time Lag

. In this article, a singularly perturbed convection-di ﬀ usion problem with a small time lag is examined. Because of the appearance of a small perturbation parameter, a boundary layer is observed in the solution of the problem. A hybrid scheme has been constructed, which is a combination of the cubic spline method in the boundary layer region and the midpoint upwind scheme in the outer layer region on a piecewise Shishkin mesh in the spatial direction. For the discretization of the time derivative, the Crank-Nicolson method is used. Error analysis of the proposed method has been performed. We found that the proposed scheme is second-order convergent. Numerical examples are given, and the numerical results are in agreement with the theoretical results. Comparisons are made, and the results of the proposed scheme give more accurate solutions and a higher rate of convergence as compared to some previous ﬁ ndings available in the literature.


Introduction
The delay differential equations are versatile in mathematical modeling of processes where they provide a realistic simulation of the real-world phenomena.The real-world operations/interactions that take time to complete can be utilized to simulate the time lag experience such as gestation time, incubation period, and transportation delays.In the model, if a small parameter multiplies the highest order derivative term, involving at least one shift term in the temporal variable, we call it singularly perturbed time delay differential equations (SPTDDE).These problems arise in the varied area of science and engineering models, for instance, in population dynamics, in epidemiology, in respiratory system, and in tumor growth [1][2][3][4][5][6].For more additional models, one can refer [7,8].Due to the appearance of the boundary layer in the solution of a singularly perturbed differential equation, classical numerical methods on equidistant grids are inadequate and fail to provide a reliable approximation, when the perturbation parameter tends to zero, unless otherwise, if one uses an unacceptably large number of grid points.Several articles have been written on the solution method for singularly perturbed delay differential equations, to cite a few [9][10][11][12][13].Among the recently conducted studies on SPTDDE of the convection-diffusion type, having a right end boundary layer, in [14], the authors used an implicit-trapezoidal scheme on uniform mesh for temporal discretization, and for spatial discretization, a hybrid scheme, which is a combination of the midpoint upwind scheme and the central difference scheme on Shishkin type meshes, is applied.In [8], the scheme is constructed using the Crank-Nicolson method for temporal discretization, and a midpoint upwind finite difference scheme on a fitted piecewise-uniform mesh in spatial discretization is applied.In [15], the scheme is devised using backward Euler's scheme on uniform mesh for temporal discretization and a new stable finite difference scheme on Shishkin mesh for spatial discretization.In [16], the problem is solved using the Crank-Nicolson method in temporal discretization, and in the spatial discretization, an exponentially fitted operator finite difference method on uniform mesh is used.All these developed schemes can work for both small and large time delays.However, except [14], the results of the above studies are first-order convergent, and there are some exceptional properties which work only for a small time delay ðδ < εÞ that cannot be assessed by these authors.
So, when we come to a singularly perturbed convectiondiffusion problem of small time lag, the following researchers address the case which works only for small time lag.In [17], the authors used the backward Euler scheme for temporal discretization and a central difference scheme with an adaptive mesh selection strategy for spatial discretization.In [18], the authors used the Euler method to discretize the time derivative and a B-spline collocation scheme for the spatial discretization on a uniform mesh.In this paper, both small and large time delays are considered, and the developed scheme is first-order convergent.In [19], the scheme is developed using the backward Euler method in the discretization of the time derivative, and a higher-order finite difference method is employed for the approximation of the spatial derivative.In this scheme, an exponential fitting factor is introduced, and the resulting scheme is first-order convergent.
From these, we are motivated to construct and analyze a higher order ε-uniform numerical scheme, for the problem considered in [17].The proposed hybrid scheme is a combination of the cubic spline method and the midpoint upwind scheme on piecewise Shishkin mesh in the spatial direction and the Crank-Nicolson method in the temporal direction.The advantage of the present method over the other methods in [17][18][19] is that it is a second-order accurate in both time and space variables, as well as its accuracy.In addition, in the paper [17][18][19], Taylor's series expansion is applied at the beginning without any restriction on the domain.In such a case, the advantage of the interval condition for the delay term in the approximation is meaningless.So we incorporate the condition where Taylor's series expansion is applied without losing the interval condition.
Notations: all through this paper, C and its subscripts denote generic positive constants independent of the perturbation parameter ε and mesh sizes.Also, k•k denotes the standard supremum norm, defined as for a function f defined on some domain D.
Lemma 3. The solution to the problem (1) satisfies the following bound: Proof.Please refer [22].
where r and y are nonnegative integers.

Numerical Scheme
The problem defined in (1) can be rewritten as where with initial and boundary conditions, D 1 = ð0, 1Þ × ð0, δÞ and D 2 = ð0, 1Þ × ðδ, T. Since the delay parameter (δ) is smaller than the singular perturbation parameter (ε), it is possible to apply Taylor's series expansion as follows: Substituting (11) into (9), we obtain where 4.1.Temporal Discretization.This section is devoted to the discretization of the temporal variable.So on the time domain ½0, T, we use uniform mesh Ω M t = ft j : t j = jΔt, for j = 0, 1, ⋯, M, Δt = T/Mg, where M is the number of mesh intervals and Δt is the time step.On Ω M t , problem ( 12) is discretized by using the Crank-Nicolson method as follows: After some rearrangement, ( 14) is written in the form where 3 Abstract and Applied Analysis and U j+1 ðxÞ = Uðx, t j+1 Þ is the semidiscrete approximation to the exact solution uðx, tÞ of (1) at the ðj + 1Þ th time level.
which is a contradiction.Hence, the minimum of ξ j+1 ðxÞ is nonnegative.
To prove that the approach is ε-uniform, more specific information on the behavior of the exact solution is required.This is done by decomposing the solution uðx, tÞ into a smooth component v ε and a singular component w ε as follows: u = v ε + w ε ; for more detail on the decomposition, one can refer [25].Lemma 9.The regular and the singular components of the solution of the discrete scheme (15) satisfy the following bounds: Proof.Follow the approach given in [21,23,24].

Discretization in the Spatial Direction
4.2.1.Shishkin Mesh.The piecewise-uniform mesh having N ≥ 4 mesh intervals on ½0, 1 is generated by dividing the interval ½0, 1 into two subintervals as ½0, τ and ½1 − τ, 1, where the transition parameter τ separates the coarse and fine region and is given by where τ 0 is a constant satisfying τ 0 ≥ 1/jα 0 j.Hence, the piecewise-uniform mesh is given by Abstract and Applied Analysis where 4.3.Cubic Spline Difference Scheme.In this subsection, we derive the cubic spline scheme on the meshes, ðxÞ at the nodal points x 0 , x 1 , ⋯, x N , satisfying the following properties: (i) Q j+1 ðxÞ coincides with a polynomial of degree three on each subintervals Then, the cubic spline function can be written as where From the properties of spline [26], For the approximation of dU j+1ðxÞ /dx and d 2 U j+1ðxÞ /dx 2 , we use the Taylor series approximations as follows: Multiplying (31) by h 2 i+1 /h 2 i and then subtracting the resulting equation from (30), we get In the same fashion multiplying (30) by h 2 i+1 /h 2 i and then adding the resulting equation into (31), we get Using (32) and (33), in we obtain the following approximation: Then, express ( 14) at x = x k in the form 5 Abstract and Applied Analysis Insert (32), (35), and (36) at ðj + 1Þ th and j th time level into (37), and substitute the resulting equation into (29).After some simplification, we obtain where ĥi = h i + h i+1 , 4.4.Midpoint Upwind Scheme.For any mesh function

and central difference D +
x D − x operators as follows: Using the midpoint upwind scheme, ( 15) can be discretized as where After some simplification, we obtain the following threeterm recurrence relation: Abstract and Applied Analysis where 4.4.1.Hybrid Scheme.The upwind scheme is stable and produces a first-order accurate solution [27][28][29].The cubic spline method satisfies the discrete maximum norm only in the fine mesh region, whereas in the coarse mesh region, the solution is unstable [28,29].The benefit of constructing a hybrid scheme is that it gives higher-order parameter uni-form convergent schemes that are oscillation-free across the region.So, using a piecewise-uniform Shishkin mesh, we develop a second-order parameter uniform convergent scheme that combines the cubic spline approach in the fine mesh region and the midpoint upwind scheme in the coarse mesh region.Thus, the discrete problem of (1) is given by Therefore, we have the following totally discrete scheme: where the coefficients s − i,j , s c i,j , and s + i,j are described in (39); also, r − i,j , r c i,j , and r + i,j are described in (44).

Stability and Error Analysis
In this section, we study the stability and ε-uniform convergence of the proposed scheme.For the analysis, we follow 7 Abstract and Applied Analysis the approach given in [23,25].The solution U N,M on D is decomposed as U = V + W, where V and W are the numerical solutions of the regular and singular components, respectively; for detail on the decomposition, please refer [25].
Lemma 10.Assume that N ≥ N 0 , where N 0 is the smallest positive integer satisfying Then, we have Proof.Divide the interval into two cases.
Case 1.For 1 ≤ i ≤ N/2, since the mesh is piecewise uniform, so that it is possible to fix the mesh size in the inner layer region by h i = h = 2τ 0 εlnN/lnN and substitute (47) into (39) and simply, we get jα 0 j < kαk, h < 1/N and (48).Similarly, s + < 0, and from (39), it is clear that Case 2. For N/2 < i < N, the mesh size in the outer layer region is h i = ð2ð1 − τ 0 ε log NÞ/NÞ < 2/N.Next, using (47) and ( 48) into (44), it is easy to see that s − < 0 and s + < 0. Simple mathematical calculation gives Hence, the proof is completed.
Note: Lemma 10 confirms that under the hypotheses given in ( 47) and (48), the tridiagonal matrix corresponding to the operator L N,M hy defined in ( 46) is an M-matrix.This implies that the operator L N,M hy satisfies the next discrete maximum principle.Assume that Proof.Assume that there exists a mesh point ði * , j + 1Þ for , and assume that The uniform stability estimate of the discrete solution is provided by the subsequent lemma.Lemma 12. Let U i,j+1 be the solution of the discrete problem (48).Then, Proof.Let us construct the barrier function as follows: and using Lemma 11, we can obtain the required bound.
Case 2. In the outer layer region, applying a similar procedure like the inner layer region, we obtain Theorem 13.Let uðx i , t j+1 Þ and U i,j+1 be the solution of problems ( 1) and (45), respectively; then, the scheme satisfies the following error estimate in the spatial direction: Proof.The proof is divided into two cases as a singular part and a regular part.
Case 1.The bound for the regular component of the solution U i,j of (1).
(i) On ½0, τ In this region, ℏ = 2τ/N = 2τ 0 εlnN/N, and Lemma 9 gives It is easy to see As a result, by applying Lemma 11, we obtain the estimate In this region, ℏ ≤ 2/N, and Lemma 9 gives and following the procedure like ½0, τ, we get Therefore, combining (66) and (73), we can obtain the desired result.Theorem 14.Let uðx i , t j Þ and U i,j be the solution of problems (1) and (45), respectively; then, the proposed scheme satisfies the following bound: Proof.Combining Theorems 7 and 13 gives the required result.

Numerical Illustration
In this section, the performance of the proposed scheme is tested through numerical examples.The exact solution of the following examples is not known, so to compute the maximum point-wise errors, we use the double mesh principle given by the following formula: where U N,M ðx i , t j Þ denote the numerical solution.The corresponding rate of convergence is computed by using the following formula: In addition, the ε-uniform maximum point-wise error E N,M is computed as and the corresponding ε-uniform rate of convergence R N,M is given by Example 1.Consider the following problem [19]:

Conclusion
A singularly perturbed convection-diffusion problem of small time lag is treated, via a hybrid fitted mesh scheme for the space discretization and the Crank-Nicolson method on a uniform mesh for time derivative.Due to the presence of a small perturbation parameter, the problem exhibits the left side boundary layer at x = 0.The error analysis of the proposed scheme is proved, and the proposed scheme is second order ε-uniform convergent.The numerical results in Tables 1 and 2 of Example 1 and Example 2, respectively, confirm that the tabulated numerical results are in agreement with the theoretical error estimates.In addition, we observed that the proposed scheme is more accurate and gives a higher rate of convergence as compared to some studies available in the literature.Since the approach is mesh dependent, when the number of mesh points increases, the effectiveness of the proposed scheme also increases.Due to the time limitations, we only discuss small time delays; however, with some slight adjustments, the scheme can also work for large time delays.

Lemma 11 .
Let ζ N,M be any mesh function defined on D N,M .

Figure 1 :
Figure 1: Surface plot of numerical solution of Example 1 at N = M = 64.
Numerical solution at t = 2, when N = 64

Figure 2 :
Figure 2: One-dimensional plot and log-log plot for Example 1.

Figure 3 :
Figure 3: Surface plot of numerical solution of Example 2 at N = M = 64.
Figures 1, 2(b), 3, and 4(b) of the surface plot and one-dimensional plot for the numerical solution of the problems in Example 1 and Exam-ple 2 clearly demonstrate the behavior of the boundary layer as ε ⟶ 0. The maximum point-wise error of Example 1 and Example 2 is plotted in Figures 2(a) and 4(a), respectively, in the log-log plot.