Abstract
The idea of direction changing and order reducing is proposed to generate an exponential difference scheme over a five-point stencil for solving two-dimensional (2D) convection-diffusion equation with source term. During the derivation process, the higher order derivatives along y-direction are removed to the derivatives along x-direction iteratively using information given by the original differential equation (similarly from x-direction to y-direction) and then instead of keeping finite terms in the Taylor series expansion, infinite terms which constitute convergent series are kept on deriving the exponential coefficients of the scheme. From the construction process one may gain more insight into the relations among the stencil coefficients. The scheme is of positive type so it is unconditionally stable and the convergence rate is proved to be of second-order. Fourth-order accuracy can be obtained by applying Richardson extrapolation algorithm. Numerical results show that the scheme is accurate, stable, and especially suitable for convection-dominated problems with different kinds of boundary layers including elliptic and parabolic ones. The idea of the method can be applied to a wide variety of differential equations.
1. Introduction
“Singular perturbation” means that a small perturbation may cause a large impact in mathematical or physical problems. This terminology was first used by Friedrichs and Wasow in their paper [1] in 1946. Despite this fairly long history, the subject of singular perturbations is not a settled one and there are still a lot of open problems to be investigated. Some surveys on the computational techniques for different kinds of singularly perturbed problems can be seen in [2–4].
Among abundant mathematical models for singular perturbation problems, the convection-diffusion problems play a very important role as they arise in fluid flows, groundwater flows, reactive flows, traffic flows, and so forth. Stynes [5] presented a brief history of the development of numerical methods for steady-state convection-diffusion problems, and he also extensively discussed some common numerical methods in [6].
In this paper, the 2D linear convection-diffusion model is considered: for transport variable defined in a bounded domain , with appropriate boundary condition: The second-order derivatives in (1) model diffusion while the first-order derivatives are associated with convective or transport process; , and are sufficiently smooth functions in . When the magnitudes of and are large (convection-dominated case), boundary and interior layers will normally emerge in the solutions of such kind of problems. The term “layers” means some thin regions in which the solution fluctuates rapidly.
It is well known that global unphysical oscillations may occur if standard discretization schemes on general meshes are used [7, 8]. Hence, stabilized methods and/or a priori adapted meshes are widely developed in order to get discrete solutions with satisfactory accuracy. In this paper, we mainly focus on special finite difference (FD) methods on uniform meshes.
The defect-correction method, as investigated in [9–11], made the early attempt to combine the accuracy of the central difference (CD) method and the stability of the upwind difference (UD) method. But according to Segal’s report [12], the method is not useful for convection-dominated problems because of its slow convergence.
In recent years, high order compact FD methods have aroused renewed interest and a variety of techniques has been developed. A compact difference scheme is one that is restricted to cells surrounding any given node and does not extend further, so it is convenient for use since no special techniques are needed for points near the boundary,.
A lot of compact FD methods may be classified as polynomial FD schemes for the influence coefficients are connected to polynomials functions of the coefficients of the differential operator; interested readers are referred to [13–20] and the references therein for more details. Numerical experiments have showed that some high order compact polynomial FD schemes can yield high accuracy results. But each scheme has its range of application, just as Stynes have stated in [5], some polynomial FD schemes are difficult to develop due to the need for extensive algebraic manipulation, and the local mesh size should be small enough to make sure that the coefficients matrix satisfies the maximum principle.
Another kind of compact FD methods is the exponential ones that the coefficients of the schemes are connected to exponential functions of the coefficients of the differential operator. The exponential FD scheme was first introduced by Allen and Southwell [21] to solve the second-order partial differential equation governing the transport of vorticity. Later Il’in [22] derived in principle the same scheme. But both [21] and [22] were investigated only for 1D case. MacKinnon and Johnson [23] derived a fourth-order exponential FD scheme for one-dimensional convection-diffusion equation and extended the formulation to two dimensions. Chen et al. [24] developed a perturbational fourth-order exponential FD scheme with diagonally dominant coefficient matrix for the convection-diffusion equations based on a second-order exponential FD scheme. The authors in [25, 26] also developed a class of fourth-order compact exponential FD schemes for solving 1D and 2D convection-diffusion problems, respectively. Here we should point out that the coefficients of the high-order exponential schemes in [23–26] involve both exponential and polynomial functions so they are still conditionally stable. Most references lack theoretical analysis of the stability and convergence of the scheme.
We notice that coefficients of most of those exponential type FD schemes [23–26] have close contact with a famous scheme Il’in [22] for 1D convection-diffusion problems. Roos [27] described ten different approaches to generate the Il’in and related schemes for the 1D problems, including the compact exponentially fitted method, collocation method, finite-volume method, and exponential Petrov-Galerkin finite-element method. But Roos also pointed out that “there are still a lot of open questions and technical difficulties” as to the generalization of most of those methods to the 2D case.
As is known to all, the remainder reapproximation technique is an efficient procedure to increase the accuracy of approximations for many problems in numerical analysis. More specifically, a basic scheme (e.g., CDS) is given first and the truncation error is analyzed by using Taylor series expansions and then the high-order derivatives that appear in the leading term of the truncation error are reapproximated by using information given by the original differential equation. The details of the idea can be seen in [13, 14, 19, 23, 25, 26].
Just inspired by the remainder reapproximation method, we aim to provide completely new way to generate an exponential FD scheme for 2D convection-diffusion problem directly without the aid of schemes for 1D case. The stencil coefficients of the FD scheme are to be determined on a five-point stencil. Different from the previous methods in the literature, the main idea is direction changing and order reducing. We keep as many high-order derivatives in the Taylor series expansions as we can at the very beginning of the construction of the scheme. During the derivation process the higher order derivatives along -direction are removed to the derivatives along -direction (or from -direction to -direction) iteratively using information given by the original differential equation (direction changing) and result in terms with the orders of -directional (or -directional) derivatives less than second-order (order reducing). Then infinite terms which constitute convergent series are kept to form a system of linear equations to determine the final stencil coefficients of the scheme.
Here I should emphasize that the meaning of the paper lies not merely on the scheme over a five-point stencil but also on the broad prospect of the idea in application: it can be used to construct the scheme of high-order convergence over a nine-point stencil, it can also be used to construct schemes on nonuniform meshes and to solve many other different kinds of problems.
The outline of the rest of the paper is as follows. In Section 2, a detailed derivation of a five-point exponential difference scheme is included for the convection-diffusion problem (1) together with a thorough discussion. Proofs of stability and convergence of the scheme are given in Section 3. Firstly, we show that our scheme is of positive type and then based on estimates of the remainder stability and convergence of the difference scheme are analyzed using barrier functions and the maximum principle for elliptic differential equation. The error is proved to be of second-order accuracy. In Section 4, different kinds of problems are included for numerical experiments, fourth-order accuracy can be obtained by Richardson extrapolation algorithm. Results show that the method can cope with problems with and without boundary layers and especially fit for convection-dominated ones. Some concluding remarks are summarized in Section 5.
2. A New Way to Generate Exponential FD Scheme and Some Discussions
2.1. Derivation of the FD Scheme
Let and be two positive integers. Partition the region by the uniform rectangular grid of points with a spacing in the -direction and in the -direction, where and . Denote all the interior mesh points in by and the mesh points on the boundary by .
The value of a function at a reference mesh point is denoted by , and those values at its four immediate neighboring mesh points are denoted by . The discretized values , which will be used in the following part of the paper, have their obvious meanings too. Before we go on, we assume and firstly just for brevity, though we may see at the end of the section that this is not the necessary condition for our method.
The compact five-grid point stencil is shown in Figure 1.
The method will be based on a five-point approximation using Taylor series expansion at the reference mesh point numbered 0 in Figure 1. For simplicity, we introduce denotation for any smooth function . Assume that the exact solution to differential equation (1) satisfies the following expression: where and are coefficients to be determined later and is the term connected with the local truncation error that also will be shown later. Sometimes we may also omit the subscript 0 and write , and so forth, just for brevity.
Set can be defined at the reference mesh point 0, respectively, using Taylor’s theorem: From the differential equation (1) one can have In a next step we freeze the data and assume that the convection terms are constant in the neighbouring meshes of grid 0; that is, , and then using (8) as an auxiliary relation, we may remove all the higher order derivatives (equal to or more than second) of function along -direction to the derivatives of along -direction; for instance, Equation (9) is really the key to learn the idea of the construction of the differencing scheme. We have removed the higher order (equal to 3) derivatives of function along -direction to the derivatives of along -direction (direction changing), and result in the right side of (9) with the orders of -directional derivatives of function less than 2 (order reducing). can be treated in the same way. Generally, substituting (10) into (6) iteratively until all the derivatives of direction- are less than second-order. The idea of direction changing and order reducing is just reflected in this procedure which can be carried out easily by means of Mathematica software, and then we obtain Noticing that the infinite terms in the brackets of (11) constitute convergent series, we rewrite it as with Similarly, we get with Substituting (5), (12), and (14) into (4) yields Denoting and comparing it with (3), we know clearly that it is necessary to set , and in order to get the numerical scheme of (1), and after some easy calculation, we get five equations:
Instead of (8), if we change differential equation (3) into the form similar procedure can be carried out as above, (17) will be obtained together with the following formulae:
These seven equations, (17)–(19) and (22)-(23) which contain five unknown values , are compatible; in fact (17)–(19) and (22) can be easily solved by means of Mathematica software to give and then (23) holds immediately. In fact, holds too. Substituting (25)–(29) into (20) and (24), we get, respectively, We denote the numerical solution at stencil points by ; the FD approximation can then be obtained by dropping the truncation error from (3) with defined in (25)–(29) and given in (30) and (31). And the boundary points values are directly given by the boundary condition (2) as
Numerical methods like these, whose coefficients involve exponential functions, are known collectively as compact exponentially difference schemes. For brevity we call the scheme EDS (exponential difference scheme) in the sequel. If we freeze the data of to be constant locally too, then in (32) we should let and ; we will call the scheme EDS0 as below.
2.2. Discussions of EDS and EDS0
It can be easily seen from (25)–(29) that the coefficients matrix of difference operator is five-diagonal and of positive type: the matrix is diagonally dominant with positive diagonal elements and negative off-diagonal elements. And since upwind effect is preserved in the scheme so hopefully it will serve well for convection-dominated problems.
The scheme EDS0 is equivalent to central differencing applied to a modified differential equation: where Since the scheme EDS0 can be said to have artificial diffusion. Likewise, the scheme EDS is equivalent to central differencing applied to the following modified differential equation:
(3) If we implement Taylor expansion for with only two terms to be reserved, that is, , and , together with and , then EDS (32) becomes CDS.
(4) If considering and on some reference point , we have no difficulty to follow the whole derivation process from the beginning of this section to arrive at the difference operator similar to defined in (32) but with , and with no change in (25), (27), and (31). Correspondingly, . Easily we can see that EDS and EDS0 are still of positive type.
Furthermore, if and are satisfied in , this problem typically has an exponential boundary layer at and parabolic boundary layers at and (see [5, 28]). The parabolic layers raise interesting numerical issues. They cause numerical instabilities that are far less serious than those engendered by exponential layers; yet it is difficult for traditional methods to approximate them accurately. But EDS and EDS0 can cope with such kind of problems without difficulty, and we will provide numerical Example 10 in Section 4 to validate it.
As for at some points (or ) and , similar formulae can be derived.
(5) If and hold in , EDS becomes CDS for Poisson’s equation.
(6) We assume and at the beginning of the section just for brevity, as for or , the difference scheme (32) still holds.
3. Analysis of Stability and Convergence of the Difference Scheme
3.1. Local Truncation Error
We begin with analyzing the truncation error of the scheme about EDS0. For brevity we only consider the case that and ; similar analysis can be done for other cases.
Estimating the truncation error of EDS0 from (3) yields Substituting Taylor expansions (5)–(7) together with values given in (25)–(29) into (37), and making some arrangement, we have Using Taylor expansion again, we arrive at Substituting the two formulae above into (38), we obtain Assuming that the solution to differential equations (1) and (2) satisfies , that is to say, are bounded by some constant, we use the norm and assuming and , then we can easily see from (40) that the truncation error of the difference approximation satisfies Here the constant , which depends on the coefficients of differential operator defined in (1) but not on mesh steps , or . We summarize the conclusion as below.
Lemma 1. Assuming that the solution of differential equations (1) and (2) satisfies , and both and are satisfied in , the truncation error of EDS0 is of second-order which is shown in (42).
Remark 2. In all cases, (42) will hold for EDS0. For instance, if in , from (40) we can see that (42) still holds with . In particular, if in , then (1) becomes Poisson’s equation, and (42) holds with as well known to us and so do other cases.
Now we complete the local truncation error estimation for EDS0. As for EDS (32), if , we can see that EDS is just a slight modification of EDS0 without changing the latter’s accuracy order since .
3.2. Convergence Analysis
For continuing our analysis we show a discrete maximum principle next.
Lemma 3 (maximum principle). If is a mesh function such that for any , then attains its maximum (minimum) value for some .
Proof. From discussions in Section 2, the difference operator defined in (32) is of positive type, so it satisfies the discrete maximum principle. Detailed proof can be seen in [29].
The maximum principle leads to a stability estimate in the discrete maximum norm, as we will now demonstrate it. In this paper, we write discrete maximum norm for mesh functions:
Lemma 4. With defined in (32), if (or ) in , we have, for any mesh function ,
Proof. We only consider the case . Set and define the mesh function ; then and in . Furthermore, setting , we get and since for , it follows from Lemma 3 that . Since in , so This proves Lemma 4.
Lemma 4 immediately shows the existence and uniqueness of the solution of the EDS0. We are now ready for an error estimate.
Theorem 5. Let and be the solutions of difference scheme (32)-(33) and differential equation (1)-(2), respectively, and assuming and conditions for Lemmas 1 and 4 are satisfied, then where .
Proof. From (3) and (32), the error satisfies, at the interior reference mesh points,
Since , from (42), the result therefore follows by application of Lemma 4 to , since for .
If , that is, ) are bounded by some constant, say by . Since can be considered as just relatively minor modifications of EDS0 (equal to second-order), so the truncation error of EDS is still of second-order accuracy. The maximum principle and convergent theorem holds too for EDS.
Remark 6. So far all the derivation process and analysis of EDS are based on the fact that and , but, in fact, in case of (see that it will not change the sign of from (25)–(29)), the method is still effective and similar theoretical analysis as we have shown above could be done.
4. Numerical Experiments
In this section, we test our methods with four examples. Example 7 is a convection-diffusion equation with constant coefficients and has exponential boundary layers. Example 8 is the one with constant coefficients and has no boundary layers, but the convection terms have singular lines on the boundary; that is, and vanish on some part of the boundary. Example 9 is convection-dominated equation with variable coefficients and has corner exponential boundary layer (Figure 5). At last, we provide Example 10 with both exponential and parabolic boundary layers.
We discuss the accuracy and stability of the EDS in comparison with CDS and UDS methods and some results given by other authors.
Still the computational domain for all the following problems is a unit square . For the sake of simplicity, we assume that the step sizes and are equal. However, this is not a restriction at all for the methods proposed in this paper. The mesh points and the boundary points can be simplified to be denoted by and .
The exact solutions for these problems are known, so the maximum pointwise error at all mesh points can be calculated using the formula where is the exact solution and is the numerical solution. The rate of convergence is defined by
Example 7.
Consider convection-diffusion equation with constant coefficients:
and with boundary condition:
In order to prevent overflow of the computer, we rewrite (52) as
Comparison of (51) with (1) shows that
The exact solution of this problem is
Figure 2 shows that there is no boundary layer when , but when Re is large the boundary layers will emerge near and . Notice that since the source term , so EDS and EDS0 are virtually identical.
The Gauss-Seidel iteration is used to solve the resulting systems of equations. The convergence criterion for the iteration is chosen to be , which is almost the minimum positive machine-precision number (the minimum positive machine-precision number is in double-precision computer system). The maximum absolute errors for different values of and are given in Table 1.
From Table 1 we can see that EDS gives more accurate results compared with CDS and UDS. In fact, the result of the EDS is almost as accurate as the machine-precision just because EDS is the exact scheme for Example 7. An exact difference scheme [30] is one for which the computed solution agrees exactly with the true solution at the mesh points. It is easy to verify this fact because if we substitute the exact solution (55) and the values of the coefficients (25)–(29) into the truncation error (37) we will get
at any reference point.
Table 2 shows clearly that with Re increases, the results become more and more accurate which validates the unconditional numerical stability in solving these linear system of equations of the EDS scheme. Figure 3 shows the absolute errors of EDS for Example 7, with mesh steps .
We continue to implement our experiments; the convergence criterion for the Gauss Seidel iteration is chosen to be for the following problems.
We prefer to apply the Richardson extrapolation technique here to compute a fourth-order accurate solution since the technique will not destroy the unconditional stability of EDS0 and EDS. Now we denote by the difference solution at any mesh point ; the Richardson extrapolation can be written as
where is the order of accuracy before the extrapolation and the order of accuracy will be increased to after the extrapolation.
| (a) |
| (b) |
| (c) |
| (d) |
| (a) |
| (b) |
Example 8.
Consider
where boundary condition:
holds in . The source term is determined so as being a solution of (58).
The figure of the exact solution is presented on the left in Figure 4. The solution has no relation to and has no boundary layer.
From Table 3 we can see that EDS0 and EDS get better results than UDS for Example 8, and as step decreases, the rate of EDS0 and EDS get closer and closer to second-order, whereas the rate of UDS gets closer to 1.
In order to get much more accurate results, we carry out Richardson’s extrapolation. In this paper, we denote the method of the EDS0 (EDS) with Richardson’s extrapolation by EDS0-RE (EDS-RE). The results are given in Table 4. We can see that Richardson’s extrapolation really increases the accuracy a lot.
Table 5 shows that EDS is a little bit more accurate than EDS0, but EDS-RE is more accurate than EDS0-RE.
(a) Exact solution
| (b) EDS numerical solution when |
| (a) |
| (b) |
Example 9. Consider convention-dominated problem with variable convection coefficient:
with boundary condition
The source term is determined so as being a solution of (60).
Corner boundary layer occurs at and . Table 6 lists maximum absolute errors in the computed solution of the problem. Table 6 again shows that EDS method gets much more accurate results for convection-dominated problems.
Example 10. Consider a constant coefficient convection-diffusion equation (1) with , and ; that is,
The boundary condition is prescribed as
The exact solution [31] is
where .
This problem represents a convection dominated flow and was used as one of the test problems by Gupta et al. [18] and many other authors [31, 32]. Figure 6 shows the exact solution with and . For most part of the domain, the exact solution is smooth, but it has a steep boundary layer (exponential layer) along the downstream edge at and has parabolic layers along the edges at and .
Results of Example 10 from different methods are shown in Table 7. We see that EDS and Richardson’s extrapolation for EDS get more accurate results.
Figure 7 draws errors of EDS for Example 10 with . The errors of different schemes along and are plotted in Figure 8. We can see that UDS cannot approximate the exact solution quite well near the exponential layer. With fixed coordinate , the errors of all the methods increase with which increases from to along the horizontal axis, but EDS and EDS-RE get much better results (Table 8).
| (a) |
| (b) |
| (a) Errors of three methods for |
| (b) Errors of three methods for |
| (c) Errors of EDS for |
| (d) Errors of EDS for |
5. Conclusions
Idea of direction changing and order reducing has been used to construct a compact exponential FD scheme over a five-point stencil for 2D convection-diffusion problem in the paper. The scheme is of concise exponential type and unconditionally stable; the convergent accurate order is second. Four numerical examples show that EDS0 and EDS can get much highly accurate results for different kinds of problems, including those with boundary layers and without layers. Richardson’s extrapolation can improve the accurate order to four.
We should mention that if one takes a scheme for its 1D analogue to form the “tensor product” of the scheme in 2D case, one can directly arrive at a scheme on five-point stencil but not on a nine-point stencil. But our method can be applied to construct much more accurate scheme over a nine-point stencil. Also, the idea can be used to form schemes on nonuniform grids. Theoretically, idea of direction changing and order reducing can also be used to solve varieties of problems, for example, parabolic equation , and so forth. Those will be investigated in future studies.
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
The author would like to thank the reviewers for their helpful suggestions.