Abstract
In this article, a singularly perturbed convection-diffusion 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 findings available in the literature.
1. 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–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–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 convection-diffusion 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–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–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, and its subscripts denote generic positive constants independent of the perturbation parameter and mesh sizes. Also, denotes the standard supremum norm, defined as for a function defined on some domain .
2. Problem Formulation
Consider the following singularly perturbed delay parabolic problem [17, 19]: subject to where is the perturbation parameter, is the delay parameter, , , , , and ; the functions , and on and on are assumed to be smooth and bounded functions that satisfy the conditions and on Under the above assumption, problem (1) exhibits a boundary layer along . The existence and uniqueness of the solution for problem (1) can be guaranteed by the sufficient smoothness of , and , along with the following compatibility conditions:
The problem (1) has a unique solution with the above conditions and assumptions [20]. Let us use the notation
3. Bounds on the Solution and Its Derivatives
Lemma 1. Suppose . Assume that , , and . Then, .
Proof. Suppose such that and It is visible that and . At the point , and Then, , which contradicts . Therefore,
Lemma 2. The solution of problem (1) is bounded and satisfies the following estimate:
Proof. Please refer [21].
Lemma 3. The solution to the problem (1) satisfies the following bound:
Proof. Please refer [22].
Lemma 4. The solution of (1) satisfies the bound
Proof. Let us take Using Lemma 3, and the fact that is bounded and . Hence, we get the required result.
Lemma 5. Let be the solution of problem (1), then where and are nonnegative integers.
Proof. The reader can follow the procedure in [21, 23, 24].
4. Numerical Scheme
The problem defined in (1) can be rewritten as where with initial and boundary conditions, and . 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 , we use uniform mesh , where is the number of mesh intervals and is the time step. On , problem (12) is discretized by using the Crank-Nicolson method as follows:
After some rearrangement, (14) is written in the form where and is the semidiscrete approximation to the exact solution of (1) at the time level.
Lemma 6. Suppose that . The local truncation error associated to scheme (15) satisfies
Proof. Using Taylor’s series expansion, expanding and centered at , we get Substituting (18) in to (1), we obtain where From (19), the local truncation error is the solution of the following BVP:
Theorem 7 (global error estimate). The global error estimate in the time direction at time level is given by
Proof. Using the local error estimate in Lemma 6, we obtain
Lemma 8 (discrete maximum principle). Suppose . Assume that , , and . Then, .
Proof. Suppose such that ; let . This gives , which implies that the point . Also, we have and . which is a contradiction. Hence, the minimum of 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 into a smooth component and a singular component as follows: ; 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].
4.2. Discretization in the Spatial Direction
4.2.1. Shishkin Mesh
The piecewise-uniform mesh having mesh intervals on is generated by dividing the interval into two subintervals as and , where the transition parameter separates the coarse and fine region and is given by where is a constant satisfying . Hence, the piecewise-uniform mesh is given by where
4.3. Cubic Spline Difference Scheme
In this subsection, we derive the cubic spline scheme on the meshes, and . Assume is an interpolating cubic spline function corresponding to the values of a function , , …, of a function at the nodal points , satisfying the following properties: (i) coincides with a polynomial of degree three on each subintervals (ii)(iii)
Then, the cubic spline function can be written as where .
From the properties of spline [26],
For the approximation of and , we use the Taylor series approximations as follows:
Multiplying (31) by and then subtracting the resulting equation from (30), we get
In the same fashion multiplying (30) by and then adding the resulting equation into (31), we get
Using (32) and (33), in we obtain the following approximation:
Then, express (14) at in the form
Insert (32), (35), and (36) at and time level into (37), and substitute the resulting equation into (29). After some simplification, we obtain where ,
4.4. Midpoint Upwind Scheme
For any mesh function , let us define the backward , forward , and central difference operators as follows:
Using the midpoint upwind scheme, (15) can be discretized as where
After some simplification, we obtain the following three-term recurrence relation: where
4.4.1. Hybrid Scheme
The upwind scheme is stable and produces a first-order accurate solution [27–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 uniform 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 , and are described in (39); also, , and are described in (44).
5. Stability and Error Analysis
In this section, we study the stability and -uniform convergence of the proposed scheme. For the analysis, we follow the approach given in [23, 25]. The numerical solution on is decomposed as , where and are the numerical solutions of the regular and singular components, respectively; for detail on the decomposition, please refer [25].
Lemma 10. Assume that , where is the smallest positive integer satisfying Then, we have
Proof. Divide the interval into two cases.
Case 1. For , since the mesh is piecewise uniform, so that it is possible to fix the mesh size in the inner layer region by and substitute (47) into (39) and simply, we get
and (48).
Similarly, , and from (39), it is clear that
Case 2. For , the mesh size in the outer layer region is . Next, using (47) and (48) into (44), it is easy to see that and . 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 defined in (46) is an -matrix. This implies that the operator satisfies the next discrete maximum principle.
Lemma 11. Let be any mesh function defined on . Assume that . Then, implies that .
Proof. Assume that there exists a mesh point for such that , and assume that , so it is clear that . Then, for ; this contradicts with the assumption . Hence, the result follows.
The uniform stability estimate of the discrete solution is provided by the subsequent lemma.
Lemma 12. Let 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.
5.1. Truncation Error
Case 1. Truncation error in the inner layer region using (15) at in (46) is given by where and . Substituting Taylor’s series expansion in the spatial variable of , , , , , , , , , and into (56) where It can be easily seen that . Thus, we get where .
Case 2. In the outer layer region, applying a similar procedure like the inner layer region, we obtain , which is
Theorem 13. Let and be the solution of problems (1) and (45), respectively; then, the proposed 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 of (1).
From the truncation error of the spatial direction (60), we have
Since , for , using Lemma 9, we obtain
Now, construct the barrier function for each interval; let us take ,
It is easy to see , and . As a result, by applying Lemma 11, it is possible to obtain the required estimate, which is
Apply a similar procedure for ; then, we can get
Case 2. The bound for the singular component of the solution of (1). The proof is done separately in and . (i)On In this region, , and Lemma 9 gives Now, construct the barrier function , It is easy to see , and , for . As a result, by applying Lemma 11, we obtain the estimate (ii)on In this region, , and Lemma 9 gives and following the procedure like , we get Therefore, combining (66) and (73), we can obtain the desired result.
Theorem 14. Let and 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.
6. 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 denote the numerical solution. The corresponding rate of convergence is computed by using the following formula:
In addition, the -uniform maximum point-wise error is computed as and the corresponding -uniform rate of convergence is given by
Example 1. Consider the following problem [19]:
Example 2. Consider the following problem [17]:
7. 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 . 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 Example 2 clearly demonstrate the behavior of the boundary layer as . 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. 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.
(a)
(b)
(a) Log-log plot
(b) Numerical solution at , when
(a)
(b)
(a) Log-log plot
(b) Numerical solution at , when 4
Data Availability
To support the finding of this study, no data were used.
Conflicts of Interest
The authors declare that they have no conflicts of interest.