Abstract
This study deals with the numerical solution of parabolic convection-diffusion problems involving two small positive parameters and arising in modeling of hydrodynamics. To approximate the solution, the backward Euler method for time stepping and fitted trigonometric-spline scheme for spatial discretization are considered on uniform meshes. The resulting scheme is shown to be uniformly convergent and its rate of convergence is one in the time variable and two in the space variable. The accuracy and rate of convergence are enhanced by using the Richardson extrapolation. To support the theoretically shown convergence analysis, we have taken some numerical examples and compared the absolute maximum error of the current method with some methods existing in the literature.
1. Introduction
Consider the following time-dependent singularly perturbed one-dimensional parabolic convection-diffusion problem with initial and boundary conditions of Dirichlet type on the domain .constrained to boundary conditionswhere are two small positive parameters. We assume that the coefficients , , and are sufficiently regular functions with the conditions for all . We also assume that the given data satisfy sufficient smoothness and compatibility conditions on the corner of the domains (0, 0) and (1, 0), i.e., . These conditions confirm the existence of a unique solution to the problem.
Physical occurrences, where energy, heat, or fluid is transformed inside a physical system due to the movement of molecules within the system (convection) and the spread of particles through the random motion from regions of higher concentration to regions of lower concentration (diffusion) are described by parabolic partial differential equations called convection-diffusion problems. In the case when the spatial derivatives are multiplied by small positive parameters (such as the partial differential equations of hydrodynamics), the problems are said to be singularly perturbed problems (SPPs) [1, 2]. The solution to these problems can have a steep gradient where boundary layers occur [3]. Using an analytic method or asymptotic expansions may be impossible to construct a solution or may fail to simplify the given problem. As a result, numerical approximations are the only option.
Classical numerical methods are inefficient to approximate their solution as the perturbation parameter goes to zero. In recent years, more due attention has been paid by different scholars on the special purpose numerical methods which are not affected by this parameter. For instance, Clavero et al. [4] defined a higher-order finite difference scheme on a piecewise-uniform Shishkin mesh, Jha and Kadalbajoo [5] constructed a robust layer adapted difference method with Shishkin-Bakhvalov mesh, and Zahra et al. [6] have developed an exponential-spline scheme on spatial piecewise-uniform mesh of Shishkin type. Das and Natesan [7] and Das and Mehrmann [8] obtained a numerical solution via adaptively generated mesh with the concept of equidistribution, Beckett and Mackenzie [9] solved singularly perturbed reaction-diffusion problem using grid equidistribution, Kadalbajoo and Jha [10] have used exponentially fitted cubic spline to solve singularly perturbed convection-diffusion problem, and Mittal and Jain [11] defined cubic B-spline collocation method for solving convection-diffusion equations. Kadalbajoo et al. [12] have solved the convection-diffusion problem with dominating convection term using a fitted collocation method. More recently, fitted operators on uniform meshes have been used to generate parameter-uniform methods like fitting the two small parameters through a finite difference [13], fitting the parameters through nonstandard method of Micken’s [14, 15], B-spline technique [16], central difference approximation [17], and piece-wise uniform mesh [18] for singularly perturbed parabolic differential equation of convection-diffusion type with two small parameters.
The numerical solution of two-point boundary-value problems which depend only on the spatial variable has been considered by Aziz and Khan [19], Bawa [20], Phaneendra et al. [21], Sirisha et al. [22], and Kumara Swamy et al. [23]. Time dependent parabolic problems subjected to initial and boundary conditions have been treated by Mohanty et al. [24] using nonpolynomial splines involving trigonometric functions. All the problems considered by these authors are with one perturbation parameter and each of the proposed methods is nonfitted. With this enthusiasm, in this work, we discuss trigonometric spline with fitted operator finite difference method on a uniform mesh to use it as a tool to solve time-dependent singularly perturbed one-dimensional initial-boundary value problem (IBVP) (1), involving two small positive perturbation parameters. This approach has the advantage over finite difference methods in which it provides continuous approximations not only for the solution but also for its considered derivatives at every point of the range of integration. Also, the continuous differentiability of the trigonometric part of the splines considered here reduces the loss of smoothness inherited by polynomial splines [21].
The manuscript is organized as follows: in Section 2, the properties of the continuous problem are given. Error estimates and detailed procedure for the construction of the trigonometric spline scheme are presented in Section 3. In Section 4, the parameter-uniform convergence analysis has been carried out. Numerical examples and results are illustrated by tables and graphs in Section 5. At last, discussion of the results obtained and conclusion are given in Section 6.
2. Properties of the Continuous Problem
In this section, we consider some properties of the continuous problem through maximum principle and bounds on the solution and on its derivatives that enable us to see a priori estimates for the solution and its derivatives.
Lemma 1 (Maximum principle). Let . If , for all on the boundary of the domain , and , for all then , for all .
Proof. The proof is given in [16].
Lemma 2 (Bounds of the solution). Let be the solution of equation (1) and , then we have the estimatewhere is the maximum norm.
Proof. The proof is given in [16].
Lemma 3 (Bounds of the derivatives). The derivatives of the solution satisfy the following bound for all non-negative integers such that .where the constant C is independent of and and depends only on the bounded derivatives of the coefficients and the source term.
Proof. The proof is given in [25].
3. Construction of the Method
3.1. Temporal Discretization
For the semidiscretization of the IBVP (1), we use the backward Euler method with uniform time-steps in the time-interval . This gives a time-mesh and
In the th time step, the semidiscrete solution of equation (1) satisfies the following differential equation of the space variable .where .
Lemma 4 (Semidiscrete maximum principle). At time level, let . If , and , then .
Proof. Let such thatsuppose that . Then, we have and from derivative tests , and . Then, , which contradicts the hypothesis. Hence, it follows that .
Now, it is observed that the operator satisfies the semidiscrete maximum principle; hence, by [4] the semidiscrete scheme (6) is stable. Furthermore, in order to analyze the convergence of the semidiscrete scheme, let us see the following error estimate.
Lemma 5 (Error estimate). Suppose and be solutions of equations (1) and (6), respectively. If , then the local error in the temporal direction satisfiesand the global error satisfieswith a constant .
Proof. From Taylor’s series expansion we haveSince is also expected to satisfy equation (6), we haveand satisfies equation (6), then we can write it asSimplifying further, we can arrive atfrom which we can deduceMoreover, the global error is given byTherefore, the time semidiscretization process is uniformly convergent of order one.
Lemma 6. Let , then the solution of the boundary-value problem equation (6) satisfies the estimatefor all .
Proof. For the proof, one can see [26].
3.2. Spatial Discretization
We consider a uniform space mesh with nodal point on [0, 1] such that with .
Let us introduce a fitting factor which controls the effect of the perturbation parameter . Then, at , equation (6) can be written as
Let be an approximate solution to of equation (17) obtained by the segment passing through the points and . The segment is supposed to satisfy the interpolation condition, and its first-order derivative is continuous at the interior nodes. Assume each segment has the formwhere and are constants and is a free parameter used to maintain consistency of the method. Let us denote
Substituting equations (18) into (19) and assigning , we can have
Solving equation (20) for the coefficients we get
From the continuity of first derivative, we have , whereThen,from which by simple arithmetic manipulation, we getwhere . Applying L’Hospital’s rule as , i.e., yields a scheme matching that of the ordinary cubic spline in , since
The consistency relation for equation (24) leads us to choose and in such a way that it satisfies the relation [19, 20]. Taking the Taylor expansion of and its first derivative about only up to second-order gives
Similarly,
From equation (26), we get
Substituting equations (28) into (27) givesand from equation (17) using spline’s second derivatives, we have
Substituting equations (31)–(33) into (24), we arrive at the following difference scheme:where
Therefore, the required scheme developed in equation (34) is termed as a fitted operator finite difference method obtained through trigonometric-spline and used to solve the problem in equation (1). It can be solved by the matrix inversion method after determining the value of the fitting factor in the next section.
3.3. Determination of the Fitting Factor
By fixing the time variable and considering the homogeneous part of equation (6) with constant coefficients that are the minimum values of the corresponding variable coefficients, we havewhere . Its characteristic equation on the th time level is , whose solutions are given by the relation
Assume there are two real solutions and that describes the boundary layers at and , respectively, based on the following two cases [27]: Case 1: if , as , Equation (1) has two boundary layers which behaves like the reaction-diffusion case with each of width at and . Then, the complementary function of equation (6) is where and are real constant numbers. Case 2: if , as ,
In this case, equation (1) has two boundary layers in the neighborhood of and with different width and , respectively. Then, the complementary function of equation (6) iswhere and are real constant numbers. The current problem is the generalization of problems with one perturbation parameter which is studied intensively as convection-diffusion or reaction-diffusion.
From the theory of singular perturbations given in [2] and using Taylor’s series expansion at layer regions and restriction to their first terms, we get the asymptotic solution aswhere is the solution of the reduced problem (when ) of equation (17) that is given by
By considering small enough, the discretized form of equation (42) is
Now, use equations (44) into (34) and restrict in each Taylor expansion to the first terms. Then, taking limits both sides as and solving for , we getwhere .
To this end, using equations (34) and (45), it is observable that
This shows the coefficient matrix associated to the difference operator is irreducibly diagonally dominant.
4. Parameter-Uniform Convergence Analysis
So far, we have shown that the continuous solution and its derivatives are bounded and errors due to the introduction of the discrete approximation in the time variable can be estimated (controlled). To realize the stability and consistency of the developed scheme, we further see the error estimate in the spatial variable and the total discrete scheme.
Lemma 7 (Discrete maximum principle). Assume the discrete function satisfies , and on . Then, at each point of .
Proof. To follow the proof by contradiction, let there exists a point where such thatand suppose that , then we have . But by using equation (25), the assumptions and , and the series representation into equation (34) we arrive atThis contradicts the assumption on , and then it completes the proof.
From Lemma 7, it follows that the discrete operator satisfies the discrete maximum principle and then the coefficient matrix associated to it is monotone. Moreover, as it is also irreducibly diagonally dominant, then it is an M-matrix. This guarantees for the existence of unique discrete solution. In the next lemma, we discuss the uniform stability of the discrete solution.
Lemma 8. The solution of the discrete scheme equation (13) satisfies the following boundwhere .
Proof. Let and define barrier functions asThe values of these barrier functions on the boundary points areand on the discretized domain isThen, by applying the discrete maximum principle given in Lemma 7, the proof is immediate.
As a result, the method is uniformly stable in the maximum norm. To establish a parameter-uniform convergence of the discrete scheme equation (34), let the truncation error be given by
Expanding each term by Taylor series method about up to third order and simplifying yields
After Taylor expansion of each coefficient up to first term and arithmetic manipulations, we reach at
By Lemma 6, is bounded, then the order of the method in the spatial direction is .
Theorem 9. Let be the solution of the problem in equation (1) at each grid point and be its numerical solution obtained by the proposed scheme given in equation (13). Then, the error estimate for the fully discrete scheme is given by
Proof. Applying the triangular inequality, the proof is immediate from Lemmas 5 and 6. Therefore, the proposed method is convergent independent of the perturbation parameters and its rate of convergence is , one in the time variable and two in the space variable.
It is well known that the more refined the step sizes, the better accurate results obtained. But this is computationally tough. So, using the Richardson extrapolation which depends on two already calculated approximations with different embedded meshes, we can minimize it and enhance the accuracy of our result with less refined mesh.
4.1. Richardson Extrapolation
Let be the approximate solutions obtained on the tensor products , respectively. The rate of convergence of the proposed method is shown to be by Theorem 9. On different meshes, the approximate can be given as
The results obtained through the Richardson extrapolation denoted and defined by better approximate the exact solution than . From equations (60) and (61), we get the error estimate
Furthermore, we have the following error estimate.
Theorem 10. Let be the solution of equation (1) and be the numerical solution obtained through the Richardson extrapolation as defined in equation (18). Then, the error after extrapolation satisfies
Proof. The proof is given in [28].
5. Numerical Examples and Results
Under this section, we apply the current method to the following numerical examples with and verify the theoretical results experimentally.
Example 1. For the problem in [29],
Example 2. For the problem in [6],whereand the exact solution is given by
It is difficult to get the exact solution of Example 1. Therefore, to illustrate the performance of the proposed scheme through the estimation of the maximum point-wise error, using the double mesh principle is indispensable. We define the absolute maximum errors before the Richardson extrapolation asand after the Richardson extrapolation aswhere and are the numbers of mesh points in and directions with and step sizes, respectively. is the approximate solution obtained using and number of meshes and is the approximate solution obtained using a double number of meshes and with half step sizes. and are given in equations (60) and (61), respectively. But for Example 2, since we have an exact solution, the absolute maximum errors before and after the Richardson extrapolation are defined, respectively, aswhere is the exact solution at the grid points. The rate of convergence for both examples before and after extrapolation is, respectively, obtained by
6. Discussion and Conclusion
We have described a trigonometric-spline method for solving time-dependent singularly perturbed convection-diffusion problems using the fitted operator technique. Applying the backward Euler method, we have first transformed the continuous time-dependent problem into an ordinary differential equation (ODE) of space variable . And then, in turn, using a trigonometric-spline discrete scheme in the space variable, we have obtained a fully discrete scheme. In our numerical experimentation, the calculated maximum point-wise errors before and after extrapolation are presented in Table 1, and the corresponding rates of convergence are given in Table 2. The overall current results for Examples 1 and 2 are shown in Tables 3 and 4, respectively.
The outcomes displayed in Tables 1, 3, and 4 clearly show that, as the perturbation parameter goes to zero, the absolute maximum point-wise error remains constant. This indicates the developed scheme given in. Equation (34) and the extrapolated technique are insensitive to the perturbation parameters and gives better results than the results in the literature. The rate of convergence of the current method is increased from one to two as a consequence of the Richardson extrapolation supporting the theoretically asserted hypothesis. This is illustrated in Tables 2 and 5. As the number of mesh points increases, the absolute maximum point-wise errors decrease, and the rate of convergence increases. These show us that there is no oscillation or unexpected change in the solution, that is, the effect of the perturbation parameter is controlled by the fitted numerical scheme obtained. From the graphs presented in Figures 1 and 2, one can observe that the problem given in Example 1 has a solution which exhibits a parabolic boundary layer. Besides this, Figure 2(a) indicates that the two boundary layers have equal width, as discussed in Case 1, whereas Figure 2(b) shows the boundary layers have different width ratifying in Case 2. From Figure 3, we can observe the profile of the numerical and exact solution of Example 2. Furthermore, since in Example 2, the boundary layer appears at the right end of the domain depending on the sign of the convection coefficient which is illustrated in Figure 4. We have plotted the maximum point-wise errors in the log-log scale for the sake of revealing the numerical order of convergence before and after the Richardson extrapolation technique in Figure 5, and again it endorses the theoretical order of convergence and the error being constant. Generally, the numerical results obtained by our method confirm the theoretical results. In our future works, we consider problems of the type treated in [32, 33].
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
Data Availability
No data were used to support the findings of this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.