Abstract
A numerical treatment via a difference scheme constructed by the Crank–Nicolson scheme for the time derivative and cubic spline in tension for the spatial derivatives on a layer resolving nonuniform Bakhvalov-type mesh for a singularly perturbed unsteady-state initial-boundary-value problem with two small parameters is presented. Error analysis of the constructed scheme is discussed and shown to be parameter-uniformly convergent with second-order convergence. Numerical experimentation is taken to confirm the theoretical findings.
1. Introduction
The following class of unsteady-state singularly perturbed one-dimensional parabolic convection-diffusion problemson the domain , together with prescribed initial and boundary conditions of Dirichlet-type,where and are two small positive parameters considered. The functions , and are assumed to be sufficiently smooth satisfying . Furthermore, we assume the compatibility condition at the two corner points, i.e., . Within this frame, a unique solution exists for the problem and exhibits a boundary layer on both the left and right lateral boundaries of the spatial domain with significantly different width depending on the relation between and [1–6].
Singularly perturbed problems (SPPs) widely appear in modeling physical problems such as fluid dynamics, heat and mass transfer in chemical engineering, quantum mechanics [7], elasticity, the theory of plates and shells, oil and gas reservoir simulation, and magnetic-hydrodynamic flow in general [2] and two-parameter singularly perturbed problems (TPSPPs) occur in chemical flow reactor theory [5] as well as in the case of boundary layers controlled by suction (or blowing) of some fluid [8] in particular.
Several articles dealt with the solution of an SPP and showed that it is revealed to have a sudden change in the boundary layer region and behaves normally in the outer region. Due to this multiscale character, the classical numerical methods cease to give satisfactory numerical results and are not competitive computationally. Therefore, to overcome the shortcomings, researchers are forced to construct numerical methods that are flexible enough to provide a solution whose accuracy is independent of the perturbation parameter(s), and its singular nature can easily be captured. The sounding numerical methods on hand are the specially designed (fitted) ones. Just to list some of these, the authors of [9–15] used the fitted operator method, the authors of [16–18] implemented the fitted mesh method, Gowrisankar and Natesan [19] constructed a layer-adapted scheme through the equidistribution of a positive integrable monitor function, and Khan et al. [20] applied variable mesh in which more mesh points are in the layer region. On the other part, the authors of [21, 22] used unfitted methods and the results yielded are not uniform.
In this paper, we construct a numerical scheme for the parabolic convection-diffusion problem (1-2) through the fitted mesh method in which Crank–Nicolson discretization is used for the time variable and cubic spline in tension on layer-adapted meshes of Bakhvalov-type for the spatial variable discretization. The scheme constructed from this composition is in its first use for the problem under consideration.
2. Properties of the Continuous Problem
In this section, we consider some properties of the continuous problem through the 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. Let . If , for all on the boundary of the domain , and , for all , then , for all .
Proof. Assume to the contrary that there exists a point such thatIt is obvious that and implies . Applying the differential operator in equation (1) on at the critical point yieldsBy the partial derivative test, we have , and . Incorporating our assumptions, with , we can arrive at . This contradicts the given statement. Therefore, we can conclude that the minimum of is nonnegative.
Lemma 2. A solution of equation (1) satisfies the estimatewhere is the maximum norm.
Proof. For the proof, one can refer [23].
Lemma 3. The derivatives of the solution satisfy the following bound for all nonnegative 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 [2].
3. Discretization and Mesh Generation
Consider a uniform mesh on the time domain with mesh intervals spaced by and given by . Applying the Crank–Nicolson method for the time variable in equations (1) and (2), we obtain the following two-level semidiscrete scheme:where and
A characteristic equation corresponding to equation (7) that is helpful to describe the boundary layers occurring at the two end points of is
This defines two continuous functions
Let
Lemma 4. For a fixed number , a certain order , and the solution of equation (3), the following holds:where is a constant independent of .
Proof. For the proof, readers can refer [24].
Lemma 5. The error estimate in the temporal direction satisfies
Proof. The proof is given in [11].
3.1. Graded Mesh
Let be the number of graded mesh points, , of the interval that are arbitrarily spaced by yielding a space mesh . Now, divide the domain into three subintervals, namely, , and , such that the second subinterval is with equidistant mesh subintervals and the rest two are with gradually graded mesh subintervals. The two transition points and are chosen such that in and in and given by
Since we are using nonuniform meshes in the spatial direction, we do not assume and to be . Furthermore, the relation is given in [24].
For the two boundary layers, following the technique in [3], we have two generating functionswhich satisfy . Analogous to [24], the mesh points in the spatial direction can be defined by
This gives condensed and refined meshes in the layer regions and points where the mesh switches from fine to coarse and vice versa are chosen to be and . The step sizes are related by
Lemma 6. The mesh sizes satisfy
Proof. For ,However, from the mean value theorem, there exists such thatThis gives .
For , from equation (15), we have . The proof in the third subinterval can be shown in a similar way as in the first subinterval.
3.2. Total Difference Scheme
Let interpolate the solution and pass through the points and such thatwhere and is a free parameter used to maintain consistency. The first derivatives of from the right and left of are given, respectively, by
Then, from the continuity condition, we have and obtainwhere . Application of the L’Hopital’s rule gives
At using the notation , from equation (7), we can havewhere, from [25],
Now, adopting the cubic spline in tension used in [25] yields
Substituting equations (24) and (25) into equation (26) and writing it in compact form, we arrive atwhere ,
The system (27) together with the given conditions is uniquely solvable. Furthermore, the matrix associated with this system of equations satisfies the property of an M-matrix that is given in [26].
4. Convergence and Stability Analysis
Obtaining stability conditions for the Crank–Nicolson by using either the Fourier analysis or the matrix method is subjected to severe restrictions. This is worse for problems with variable coefficients like in our case. Thus, the maximum principle analysis can be taken as an alternative means [27]. So, to establish the parameter uniform convergence of the proposed scheme, let us see the following lemmas.
Lemma 7. Assume that
Then, under these assumptions, the matrix associated with equation (27) at each time level is an M-matrix. Hence, the operator satisfies the following discrete maximum principle.
Lemma 8. (discrete maximum principle). Let the discrete function satisfy and . Then, implies .
Proof. To follow the proof by contradiction, let there exist a mesh point for such thatand suppose . Then, this indicates , and for , we obtainwhich contradicts the assumption . This follows that .
Lemma 9. For any discrete function defined on the tensor product such that for all on the boundary of , thenfor .
Proof. The proof follows from the application of the discrete maximum principle to a mesh function given by
Lemma 10. The truncation error in the spatial direction, denoted by , satisfies
Proof. Let the truncation error in the spatial variable with fixed time be given byThe Taylor series expansion of each term in space up to the third order derivative and collecting , and givewhere upon simplification and restricting the expansion of the coefficients to their first term yieldUsing the relation , Lemma 4 and equations (17) and (23), it can be further simplified toAs the operator satisfies the discrete maximum principle on the tensor product and the error is estimated, the current method is uniformly stable in the maximum norm. Applying the triangular inequality from equations (12) and (38), we can have the following theorem.
Theorem 1. Let be the solution of equations (1) and (2) and , be the solution of equation (13). Then, the error estimate for the fully discrete scheme is given by
Therefore, the presented method is convergent independent of the perturbation parameters and its rate of convergence is two.
5. Computational Results and Discussion
To verify the theoretical results experimentally, we solve the following examples.
Example 1. For the problem in [28],
Example 2. Consider the following singular perturbation initial-boundary-value problem in [28],
Example 3. Consider the following singular perturbation initial-boundary-value problem in [2],As the considered examples have no exact solution, we calculate the absolute maximum errors using the double mesh principle as follows:where is an approximate solution obtained using and step sizes in the and directions, respectively, and is an approximate solution obtained by bisecting the step sizes. As well, the corresponding rate of convergence is defined byThe maximum absolute error and rate of convergence for Examples 1–3 are given in Tables 1–3. From these results, one can observe that the current method converges independently of the perturbation parameters and gives more accurate numerical results than the existing results available in the literature. These results are also in excellent agreement with the theoretical findings. The solution profiles illustrated in Figures 1–3 indicate the meshes in the layer regions are more condensed and refined than the outer layer. From the log-log plots in Figures 4–6, we can conclude that our method is parameter uniform.
6. Conclusions
The present study provides the numerical solutions for a singularly perturbed unsteady-state initial-boundary-value problem with two small parameters. The method comprises of developing a scheme through discretization of the time variable by the Crank–Nicolson method on a uniform mesh and the space variable by cubic spline in tension with Bakhvalov-type mesh. The method is flexible and successful in capturing the sudden change of the solution behavior of the problem in the boundary layer region.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.