Uniformly Convergent Hybrid Numerical Method for Singularly Perturbed Delay Convection-Diffusion Problems

*is paper deals with numerical treatment of nonstationary singularly perturbed delay convection-diffusion problems. *e solution of the considered problem exhibits boundary layer on the right side of the spatial domain. To approximate the term with the delay, Taylor’s series approximation is used. *e resulting time-dependent singularly perturbed convection-diffusion problems are solved using Crank-Nicolson method for temporal discretization and hybrid method for spatial discretization. *e hybrid method is designed using mid-point upwind in regular region with central finite difference in boundary layer region on piecewise uniform Shishkin mesh. Numerical examples are used to validate the theoretical findings and analysis of the proposed scheme. *e present method gives accurate and nonoscillatory solutions in regular and boundary layer regions of the solution domain. *e stability and the uniform convergence of the scheme are proved. *e scheme converges uniformly with almost second-order rate of convergence.


Introduction
Singularly perturbed differential equations are differential equations in which the highest-order derivative term is multiplied by a small perturbation parameter ε. Singularly perturbed problems appear in several research areas of applied mathematics [1][2][3], for example, in water quality problem in river networks, in simulation of oil extraction from underground reservoirs, in convective heat transport problem with large Peclet numbers, in atmospheric pollution, and in fluid flow at high Reynolds number.
One of the common examples of singularly perturbed problems is the Navier-Stokes equation of fluid dynamics: z zx u 2 + p + z zy (uv) + 1 Re with suitable initial and boundary conditions. In (1), p is the pressure and u and v are the velocity components in the x and y directions, respectively. Parameter Re is known as Reynolds number which is proportional to the length scale, velocity scale, and density and inversely proportional to the kinematic viscosity of the fluid [4]. Singularly Perturbed Delay Differential Equations (SPDDEs) are differential equations in which the highestorder derivative term is multiplied by small perturbation parameter ε and involves at least one delay term. SPDDEs model process for which the evaluation does not only depend on the current state of the system but also includes the past history. So, it is called problem with memory.
In general, when the perturbation parameter tends to zero, the smoothness of the solution of SPDDEs deteriorates and it forms boundary layer [5]. As a result, standard numerical methods like FDM, FEM, and collocation method give accurate results only when mesh size is less than the perturbation parameter (i.e., h ≤ ε). But for the case where h > ε (or ε is very small), it leads to oscillations in the computed solution while using the standard numerical methods [5]. To overcome these oscillations, while using the standard numerical methods, a large number of mesh points are required, which is not practical due to round-off error, computer processing ability, and computer memory issue.
To handle this drawback, recently, different authors have developed numerical schemes for solving singularly perturbed parabolic problems having delay on the spatial variable. Chakravarthy and Kumar [6] used Taylor's approximation for the delay and discretize the domain using uniform mesh in the temporal direction and adaptive mesh is generated using the concept of entropy function for the spatial direction. e method is based on central finite difference scheme. Kumar and Kadalbajoo [7] first treated the delay using Taylor's series approximation and developed numerical scheme using upwind for time derivative with B-spline collocation for spatial derivatives approximation on Shishkin mesh.
In [8], Kumar used Taylor's series approximation to tackle the shift terms and semidiscretized the time direction using the Crank-Nicolson method. e resulting set of ODEs is discretized using a mid-point upwind finite difference scheme on Shishkin mesh. Rai and Sharma [9] developed numerical scheme for turning layer problem. ey first treated the delays using Taylor's series approximation and used upwind FDM on Shishkin mesh for spatial derivative approximation. e scheme converges with rate of convergence of almost one. Ramesh and Kadalbajo [5] first treated the delays using Taylor's series approximation and developed numerical scheme using upwind for time derivative discretization with upwind and mid-point upwind for the spatial derivatives discretization on Shishkin mesh. Rao and Chakravarthy [10] developed numerical scheme using upwind for time derivative with fitted operator FDM for spatial derivatives discretization. Kumar et al. [11] developed complete flux-finite volume method for treating the considered problem. Woldaregay and Duressa [12] developed uniformly convergent scheme using exponentially fitted operator finite difference method. ey proved the uniform convergence of the scheme with linear order of convergence.
Numerical treatment of singularly perturbed parabolic delay differential equations still needs improvement; developing higher-order uniformly convergent numerical scheme is an active research [13]. In this paper, we formulate almost second-order uniformly convergent numerical scheme for treating the considered problem. e proposed scheme consists of Crank-Nicolson method for temporal discretization with hybrid FDM (mid-point upwind in outer layer and central difference in boundary layer) on Shishkin mesh for spatial discretization.

Notation 1.
e symbols N and M are denoted for the number of mesh intervals in spatial and temporal discretization, respectively. Symbol C(in some cases indexed) denotes positive constant independent of ε and N. e norm ‖.‖ represents the maximum norm.

Statement of the Problem
We consider a class of singularly perturbed parabolic delay convection-diffusion problems containing delay on the reaction term of the spatial variable of the form where D � Ω × Λ � (0, 1) × (0, T] for some fixed number T > 0, with smooth boundary zD � D\D, and ε, 0 < ε ≪ 1 is singular perturbation parameter and δ is delay parameter satisfying δ < ε. e functions a(x), α(x), β(x), f (x, t), u 0 (x), ϕ(x, t), and ψ(1, t) are assumed to be sufficiently smooth, bounded, and independent of parameter ε to guarantee the existence of unique solution. e coefficient functions α(x) and β(x) are assumed to satisfy where α * and β * are the lower bounds for α(x) and β(x), respectively.

Properties of the Continuous Solution and Bounds.
In case of δ < ε, using Taylor's series approximation for the terms containing delay is appropriate [14]. So, we approximate 2 International Journal of Differential Equations Using (3) into (2) gives For small values of δ, (2) and (5) are asymptotically equivalent, since the difference between the two is O(δ 3 ). We e problem in (5) exhibits regular boundary layer and the position of the boundary layer depends on the following condition: If a(x) − δα(x) < 0, left boundary layer exists and, for a(x) − δα(x) < 0, right boundary layer exists. In case a(x) − δα(x) changes sign, inside layer will exist [15,16]. For the case where δ � 0, the problem reduces to singularly perturbed parabolic PDEs, in which for small ε it exhibits boundary layer. e layer is maintained for δ ≠ 0 but it is sufficiently small. We assume also that a(x) − δα(x) ≥ a * − δα * > 0 implies occurrence of right boundary layer.
ere does not exist a constant C independent of ε such that |u(x, t) − ψ(1, t)| ≤ Cx, since the boundary layer is on the right side of the domain. e problem obtained by setting ε 2 − (δ 2 /2)α * � 0 in (5) is called reduced problem and is given as For small values of ε 2 − (δ 2 /2)α * , the solution u(x, t) of (5) is very close to the solution u 0 (x, t) of (5).
To show the bounds of the solution u(x, t) of (5), we assume without loss of generality that u 0 (x) � 0 [17]. Since u 0 (x) is sufficiently smooth and using the property of norm, we prove the following lemma.

Lemma 2.
e solution u(x, t) of (5) is bounded as Proof. From the inequality |u( and u 0 (x) is bounded, this implies that there exists a constant C 1 such that |u(x, t)| ≤ C 1 .
□ e differential operator L is denoted for the differential equations in (4) as International Journal of Differential Equations □ Lemma 4 Suppose that Lemmas 2 and 3 hold. en, the bound on the derivatives of u(x, t) with respect to t is given by Proof. See [18].
International Journal of Differential Equations Lemma 6.
e derivative of the solution u(x, t) of the problem in (5) satisfies the bound Proof. See [18].

Numerical Scheme Formulation
denote the approximation of u(x, t j+1 ) at the (j + 1)th time level discretization. Using Crank-Nicolson method for semidiscretizing the problem in (4), we obtain a system of boundary value problems: Next, we state the discrete maximum principle for semidiscretized problem as follows.

Lemma 8 (local error estimate). Considering the result in Lemma 4, the local error in the temporal semidiscretization is given by
Proof. Using Taylor's series expansion about t j+1/2 , we obtain Substituting (19) into (5), we obtain International Journal of Differential Equations where . By simplifying, (20) can be rewritten as and, from (8), we have From (21) and (22), we obtain with the boundary conditions applying the maximum principle gives ‖e j+1 ‖ ≤ C 2 (Δt) 3 . □ Next, we need to show the bound for the global error of the temporal discretization. Let E j+1 denote the global error estimate up to the (j + 1)th time step.

Lemma 9 (Global error estimate).
e global error estimate up to t j+1 time step is given by Proof. Using the local error estimate up to the (j + 1)th time step given in Lemma 8, we obtain the global error estimate at the (j + 1)th time step as where C i is a constant independent of ε 2 − (δ 2 /2)α * and Δt.
(28) e regular component is the solution to the nonhomogeneous problem: and the singular component is the solution to the homogeneous problem: (1), Lemma 11. For j � 1, 2, . . . , M − 1, the regular component of the solution and its derivatives satisfies the bound In general, the derivatives are bounded as Proof. See [18].
□ e next lemma gives a bound for the derivatives of solution.

Lemma 13.
e solutions of the boundary value problems in (16) satisfy the bound Proof. e results from Lemmas 11 and 12 give the required bound. be the discretized domain of the spatial domain, where N be the number of grid points and must be even positive integer.
Since the considered problem exhibits right boundary layer, we set a mesh transition parameter τ which divides the domain [0, 1] into outer layer region Ω 1 � [0, 1 − τ] and boundary layer region where N denotes the number of mesh points in spatial discretization and σ denotes a constant that represents the order of the scheme. We set a uniform mesh N/2 in Ω 1 with mesh spacing H � 2(1 − τ)/N and similarly a uniform mesh N/2 is placed in Ω 2 with the mesh spacing h � 2τ/N. So, the mesh point x i is given as In a similar fashion, one can generate a Shishkin mesh for the left boundary layer problem.
For the mesh functions U j+1 (x i ) at the grid points x i , the backward and forward difference operators are given as respectively. e central difference operator D 0 x is International Journal of Differential Equations Similarly, the second-order derivative D + x D − x is given by is second-order consistent on any mesh, just as on equidistant meshes, but this is not the case for , which arises in the consistency error analysis, is only of first order on arbitrary mesh. For more details, the reader is referred to [3,19]. When ε is small, the central difference approximation D 0 x U j+1 (x i ) to (d/dx)U j+1 (x) on equidistant grid leads to nonphysical oscillations of the computed solution which is due to the loss in stability, unless the grid size is very small.
To overcome the challenge indicated in Remark 2, hybrid numerical scheme containing central difference in the boundary layer region and mid-point upwind scheme on outer region is developed in order to get better accuracy and rate of convergence. Using the operators above, the discrete hybrid scheme can be formulated as and the notation d(x i− 1/2 ) stands for d(x i ) + d(x i− 1 )/2 and similarly for others.
e matrix associated with operator L Δt,N is of size (N + 1) × (N + 1) and satisfies the property of M-matrix with the given condition in (43). See the detailed proof in [20].
Proof. Using the discrete operator on the mesh function Lemma 16. For i � 0, 1, . . . , N and fixed j, define a mesh function International Journal of Differential Equations with the convention that if i � 0, then S 0 � 1.

en, for
for some constant C.

Decomposition of Numerical Solution.
Here, we decompose the numerical solution U i,j+1 into regular and singular components in a similar fashion to the continuous case: where the regular part satisfies the nonhomogeneous equation: and the singular component W i,j+1 satisfies the homogeneous equation: and now the error in the numerical solution can also be decomposed as We estimate the error in the regular and singular solution separately.

Theorem 1 (error in the regular component).
e error in the regular component V i,j+1 satisfies the following estimate at the (j + 1)th time level: Proof. For the error in the regular component, Using the barrier function Z j+1 (x i ) � CH((ε 2 − (δ 2 /2)α * ) + H)(1 + x i ) and applying Lemmas 14 and 15, we have □ Theorem 2 (error in singular component). e error in the singular component satisfies the following bound: Proof. Let us analyse the error in singular component. We know that International Journal of Differential Equations 13 and |W 1,j+1 | � |W j+1 (1)| ≤ C and L Δt,N W i,j+1 � 0 for i � 1, 2, . . . , N − 1. Now, let us use Lemma 16 and considering the barrier function Z j+1 (x i ) as for W i,j+1 by discrete maximum principle, we get erefore, using Lemmas 18 and 12, we get on the coarse mesh part, we have and, on the finer part, we compute the error as we did for the smooth component, that is, by using the consistency and barrier function, but we restrict our domain to [1 − τ, 1] rather than [0, 1]. Now, taking i � N/2 in (71), we get Now, modifying Lemma 17 for N/2 < i < N and using Lemma 12, we have By Lemma 18 and using sinhx ≤ Cx for x ∈ [0, 1], now, consider a barrier function and, using the discrete maximum principle, we obtain Hence, the proof is complete. □ Theorem 3. e error due to the spatial discretization of the computed solution satisfies the following estimate: Proof. Combining the error estimates given above for the regular and singular components in eorems 1 and 2 gives the following result.
□ Theorem 4. Let u(x i , t j+1 ), U j+1 (x i ), and U i,j+1 be the solutions of the problems in (5), (16), and (40), respectively; then the proposed scheme satisfies the following error estimate: Proof. Using the error for the temporal and spatial discretization in Lemma 9 and eorem 3 gives the required bound.

Numerical Results and Discussion
In this section, we consider numerical examples to illustrate the theoretical findings of the developed hybrid numerical scheme.   T � 2 is subject to the initial condition u(x, 0) � 0, x ∈ [0, 1] and the interval-boundary conditions ϕ( i,j denotes the computed solution on double number of mesh points 2N, 2M by including the mid-points x i+1/2 � ((x i+1 + x i )/2) and t j+1/2 � ((t j+1 + t j )/2) into the mesh points. e maximum point-wise absolute error is given by

Example 2. Consider the problem
e ε-uniform error estimate is calculated using e rate of convergence of the scheme is calculated using and the ε-uniform rate of convergence is calculated using In Figure  One observes in Figures 2 and 3 that as the perturbation parameter ε goes small, the boundary layer formation becomes visible. In these figures, the effect of the singular perturbation parameter on the solution is shown. In addition, for small ε, there are a sufficient number of mesh points in the boundary layer region, which shows the layer resolving nature of the scheme. In Figure 4, the effect of delay parameter on the behaviour of the solution is shown by using different values for delay parameter for the test examples by keeping ε � 2 − 5 at T � 1. In Figure 5, the absolute   International Journal of Differential Equations errors of the proposed scheme for ε � 2 − 20 , N � 2 7 , and M � 2 7 are given. As is seen on the figures, the absolute error is dominant in the boundary layer region.
In Tables 1 and 2, the maximum point-wise absolute error, ε-uniform error, and the ε-uniform rate of convergence of the scheme are given. e results are given for the boundary layer and regular regions separately. In each column (or for each N and M), as ε ⟶ 0, the maximum absolute error becomes uniform. is shows that the convergence of the scheme is independent of the perturbation parameter. In the last two rows of Tables 1 and 2, we have the ε-uniform error and the ε-uniform rate of convergence of the scheme. One observes in these tables that the scheme converges with order of convergence of almost two. In Tables 3 and 4, the maximum absolute error for different values of the delay parameter is given for ε � 2 − 10 . e results in these tables show that the convergence of the scheme is also independent of the delay parameter. In Table 5, the comparison of the ε-uniform error and the ε-uniform rate of convergence is given. e result in the proposed scheme is more accurate than those in the schemes in [5,10,21]. Based on the results given in Tables 1-5, the proposed scheme is accurate, stable, and ε-uniformly convergent with rate of convergence of almost two.

Conclusion
is paper deals with the numerical treatment of singularly perturbed parabolic delay convection-diffusion problems. e solution of the considered problem exhibits boundary layer behaviour. e delay term is approximated using Taylor's series and resulted in asymptotically equivalent singularly perturbed parabolic partial differential equations.
e developed numerical scheme is constituted of Crank-Nicolson method in temporal discretization with hybrid finite difference on Shishkin mesh for the spatial discretization. e stability of the scheme is investigated using construction of barrier function for the solution bound and maximum principle for the existence of unique discrete solution. e ε-uniform convergence of the scheme is proved. e applicability of the scheme is investigated by considering test examples exhibiting boundary layer. e numerical results are tabulated in terms of maximum absolute error and ε-uniform error and are observed to be in agreement with the theoretical results. Effect of the delay on the solution profile is depicted in graphs. e computational result shows uniform convergence of the scheme with order of convergence of almost two.

Data Availability
No external data were used in this paper.

Conflicts of Interest
e authors declare no conflicts of interest.