Abstract
This paper deals with numerical treatment of nonstationary singularly perturbed delay convection-diffusion problems. The 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. The resulting time-dependent singularly perturbed convection-diffusion problems are solved using Crank-Nicolson method for temporal discretization and hybrid method for spatial discretization. The 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. The present method gives accurate and nonoscillatory solutions in regular and boundary layer regions of the solution domain. The stability and the uniform convergence of the scheme are proved. The scheme converges uniformly with almost second-order rate of convergence.
1. 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–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:with suitable initial and boundary conditions. In (1), is the pressure and and are the velocity components in the and directions, respectively. Parameter 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 highest-order 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., ). But for the case where (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. The 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. The 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. They first treated the delays using Taylor’s series approximation and used upwind FDM on Shishkin mesh for spatial derivative approximation. The 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. They 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. The 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. The symbolsandare denoted for the number of mesh intervals in spatial and temporal discretization, respectively. Symbol(in some cases indexed) denotes positive constant independent ofand. The normrepresents the maximum norm.
2. 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 formwhere for some fixed number , with smooth boundary , and is singular perturbation parameter and is delay parameter satisfying . The functions , and are assumed to be sufficiently smooth, bounded, and independent of parameter to guarantee the existence of unique solution. The coefficient functions and are assumed to satisfywhere and are the lower bounds for and , respectively.
2.1. 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
For small values of , (2) and (5) are asymptotically equivalent, since the difference between the two is . We have since is lower bound of .
The problem in (5) exhibits regular boundary layer and the position of the boundary layer depends on the following condition: If , left boundary layer exists and, for , right boundary layer exists. In case changes sign, inside layer will exist [15, 16]. For the case where , the problem reduces to singularly perturbed parabolic PDEs, in which for small it exhibits boundary layer. The layer is maintained for but it is sufficiently small. We assume also that implies occurrence of right boundary layer.
Lemma 1. Let; there exists a constantindependent ofsuch that the solutionsatisfies
Remark 1. There does not exist a constantindependent ofsuch that, since the boundary layer is on the right side of the domain.
The problem obtained by setting in (5) is called reduced problem and is given as
For small values of , the solution of (5) is very close to the solution of (5).
To show the bounds of the solution of (5), we assume without loss of generality that [17]. Since is sufficiently smooth and using the property of norm, we prove the following lemma.
Lemma 2. The solutionof (5) is bounded as
Proof. From the inequality , we have which implies that ; since and is bounded, this implies that there exists a constant such that .
The differential operator is denoted for the differential equations in (4) as
Lemma 3 (the maximum principle). Let be a sufficiently smooth function defined on which satisfies . Then implies that .
Proof. Suppose that there exists such that . It is clear that ; that is, . Since which implies that , and . Giving that which is contradiction to the assumption made . Therefore, .
Lemma 4. Suppose that Lemmas 2 and 3 hold. Then, the bound on the derivatives of with respect to is given by
Proof. See [18].
Lemma 5 (stability estimate). The solution of the continuous problem in (5) is bounded as
Proof. Define two barrier functions as . At the initial stage, we haveOn the boundaries, we haveOn domain , we haveSince and , by the hypothesis of the maximum principle, we obtain , which gives the required bound.
Lemma 6. The derivative of the solution of the problem in (5) satisfies the bound
Proof. See [18].
3. Numerical Scheme Formulation
3.1. Temporal Semidiscretization
Subdivide the time domain into M intervals as , where . Let denote the approximation of at the time level discretization. Using Crank-Nicolson method for semidiscretizing the problem in (4), we obtain a system of boundary value problems:wherefor and .
Next, we state the discrete maximum principle for semidiscretized problem as follows.
Lemma 7. (semidiscrete maximum principle). For each , let be sufficiently smooth function on domain . If and , then .
Proof. Suppose that is such that . From assumption, it is clear that implies that . Applying the property of extrema values in calculus gives , and , given that which is a contradiction to . Therefore, we conclude that . Hence, operator satisfies the semidiscrete maximum principle.
The local error for each time step is denoted by and given by .
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 , we obtainSubstituting (19) into (5), we obtainwhere . By simplifying, (20) can be rewritten asand, from (8), we haveFrom (21) and (22), we obtainwith the boundary conditions and . Hence, applying the maximum principle gives .
Next, we need to show the bound for the global error of the temporal discretization. Let denote the global error estimate up to the time step.
Lemma 9 (Global error estimate). The global error estimate up to time step is given by
Proof. Using the local error estimate up to the time step given in Lemma 8, we obtain the global error estimate at the time step aswhere is a constant independent of and .
Lemma 10. For , let be a solution of (16); then the solution satisfies the following bound:where .
Proof. Construct barrier functions . We need to show that if and , this implies that . We easily obtain that and , andUsing the semidiscrete maximum principle, we conclude that . Hence, the required bound is satisfied.
3.2. Decomposition of Solution
Since the bounds on the derivatives of the solution are not sharp enough for the proof of uniform convergence, we need to derive more stronger error bounds. This can be achieved by decomposing the solution into regular and singular component and doing error bound separately for each. Let
The regular component is the solution to the nonhomogeneous problem:and the singular component is the solution to the homogeneous problem:
Lemma 11. For , the regular component of the solution and its derivatives satisfies the boundIn general, the derivatives are bounded as
Proof. See [18].
Lemma 12. For , the singular component of the solution and its derivatives satisfies the bound:
Proof. See [18].
The next lemma gives a bound for the derivatives of solution.
Lemma 13. The solutions of the boundary value problems in (16) satisfy the bound
Proof. The results from Lemmas 11 and 12 give the required bound.
3.3. Spatial Discretization and the Hybrid FDM
Let be the discretized domain of the spatial domain, where be the number of grid points and must be even positive integer. For each , we define . Since the considered problem exhibits right boundary layer, we set a mesh transition parameter which divides the domain into outer layer region and boundary layer region , wherewhere 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 in with mesh spacing and similarly a uniform mesh is placed in with the mesh spacing . So, the mesh point is given as
In a similar fashion, one can generate a Shishkin mesh for the left boundary layer problem.
For the mesh functions at the grid points , the backward and forward difference operators are given asrespectively. The central difference operator is
Similarly, the second-order derivative is given by
Remark 2. The approximationis second-order consistent on any mesh, just as on equidistant meshes, but this is not the case for: a term, 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]. Whenis small, the central difference approximationtoon 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 aswhereand the notation stands for and similarly for others.
3.4. Properties of the Discrete Scheme and Uniform Convergence
Assume that . We assign each of the subintervals and with number of equidistant grid points. Let be the mesh width in the subinterval and let be the mesh width in . These mesh widths satisfy
To attain stability condition while using central difference scheme on the singular component, we need the following condition [19]:
Lemma 14 (discrete maximum principle). The discrete schemeand for givenandfor fixedhas a solution. Ifforand ifand, then,.
Proof. The matrix associated with operator is of size and satisfies the property of -matrix with the given condition in (43). See the detailed proof in [20].
Lemma 15. Let . Then, there exists a positive constant such that .
Proof. Using the discrete operator on the mesh function , for , we have
For the case where , we have
Lemma 16. Forand fixed, define a mesh functionwith the convention that if , then . Then, for , we havefor some constant .
Proof. Since , we have . Now, putting these on the discrete operator, we haveFor , we haveHence, the required result holds. For the case where ,
Lemma 17. Letbe any smooth function defined on. For, define truncation errortoatwhich is to be; then there exists a constantsuch thatfor .
Proof. See Lemma 4.1 in [16].
Lemma 18. For each , we have
Proof. For each , we haveNow multiplying for , we obtain
Lemma 19. For Shishkin mesh defined above, there exists a constant such that
Proof. Let by Lemma 1 in [20].Therefore, we haveSince is bounded for , the desired result is obtained.
3.5. Decomposition of Numerical Solution
Here, we decompose the numerical solution 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 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). The error in the regular componentsatisfies the following estimate at thetime level:
Proof. For the error in the regular component,Let ; then, . Using the barrier function and applying Lemmas 14 and 15, we have
Theorem 2 (error in singular component). The error in the singular component satisfies the following bound:
Proof. Let us analyse the error in singular component. We know thatand and for . Now, let us use Lemma 16 and considering the barrier function asfor by discrete maximum principle, we getTherefore, using Lemmas 18 and 12, we geton the coarse mesh part, we haveand, 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 rather than . Now, taking in (71), we getNow, modifying Lemma 17 for and using Lemma 12, we haveBy Lemma 18 and using for , now, consider a barrier functionand, using the discrete maximum principle, we obtainHence, the proof is complete.
Theorem 3. The 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 Theorems 1 and 2 gives the following result.
Theorem 4. Let , , and 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 Theorem 3 gives the required bound.
4. Numerical Results and Discussion
In this section, we consider numerical examples to illustrate the theoretical findings of the developed hybrid numerical scheme.
Example 1. Consider the problem is subject to the initial condition and the interval-boundary conditions .
Example 2. Consider the problemwith subject to the initial condition and the interval-boundary conditions .
In the considered examples, the exact solution of the problems is not known. We use the double mesh technique to calculate maximum absolute error. Let denote the computed solution of the problem on number of mesh points in and directions, respectively, and denotes the computed solution on double number of mesh points by including the mid-points and into the mesh points. The maximum point-wise absolute error is given by
The -uniform error estimate is calculated using
The rate of convergence of the scheme is calculated usingand the -uniform rate of convergence is calculated using
In Figure 1, solution profiles of Examples 1 and 2 are plotted for different values of time level by taking and . The solutions of the problems in the considered examples exhibit right boundary layer of thickness . 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 at . In Figure 5, the absolute errors of the proposed scheme for , and are given. As is seen on the figures, the absolute error is dominant in the boundary layer region.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
In Tables 1 and 2, the maximum point-wise absolute error, -uniform error, and the -uniform rate of convergence of the scheme are given. The results are given for the boundary layer and regular regions separately. In each column (or for each and ), as , the maximum absolute error becomes uniform. This 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 . The 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. The 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.
5. Conclusion
This paper deals with the numerical treatment of singularly perturbed parabolic delay convection-diffusion problems. The solution of the considered problem exhibits boundary layer behaviour. The delay term is approximated using Taylor’s series and resulted in asymptotically equivalent singularly perturbed parabolic partial differential equations. The developed numerical scheme is constituted of Crank-Nicolson method in temporal discretization with hybrid finite difference on Shishkin mesh for the spatial discretization. The 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. The -uniform convergence of the scheme is proved. The applicability of the scheme is investigated by considering test examples exhibiting boundary layer. The 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. The 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
The authors declare no conflicts of interest.