A New Way to Generate an Exponential Finite Difference Scheme for 2 D Convection-Diffusion Equations

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 equationwith 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 onemay gainmore 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.


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][3][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].
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][10][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][14][15][16][17][18][19][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 secondorder 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][24][25][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][24][25][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, finitevolume method, and exponential Petrov-Galerkin finiteelement 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 secondorder 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 convectiondominated ones.Some concluding remarks are summarized in Section 5.

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 ℎ = 1/ and  = 1/.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  0 , and those values at its four immediate neighboring mesh points are denoted by   ,  = 1, 2, 3, 4. 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 (, ) > 0 and (, ) > 0 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   ( = 0, 1, 2, 3, 4) and  01 ,  10 are coefficients to be determined later and  0 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  0 = ,  0 = ,  0 = ,  0 = , and so forth, just for brevity.Set 1 ,  2 ,  3 ,  4 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, (, ) =  0 , (, ) =  0 , 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,  (0,3) =  2  (0,1) +  (1,0) +  (1,1)   −  (2,0) −  (2,1) −  −  (0,1) .( 9) 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). (0,4) ,  (0,5) . . .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 V =  00  +  01  (0,1) +  10  (1,0) +  11  (1,1) + ⋅ ⋅ ⋅ +  00  +  01  (0,1) + ⋅ ⋅ ⋅ and comparing it with (3), we know clearly that it is necessary to set  00 = 0,  01 = 0,  10 = 0, and  00 = 1,  01 =  01 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   ,  = 0, 1, 2, 3, 4, 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,  11 = 0 holds too.Substituting ( 25)-( 29) into ( 20) and ( 24), we get, respectively, We denote the numerical solution at stencil points by   ( = 0, 1, 2, 3, 4); the FD approximation can then be obtained by dropping the truncation error  0 from (3) with   ( = 0, 1, 2, 3, 4) defined in ( 25)-( 29) and  01 ,  10 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  01 = 0 and  10 = 0; we will call the scheme EDS0 as below.

Discussions of EDS and EDS0
(1) 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.
(2) The scheme EDS0 is equivalent to central differencing applied to a modified differential equation: where 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   ( = 0, 1, 2, 3, 4) with only two terms to be reserved, that is, ), and  0 ≈ (2/ℎ 2 ) + (2/ 2 ), together with  01 = 0 and  10 = 0, then EDS (32) becomes CDS.
Furthermore, if (, ) > 0 and (, ) ≡ 0 are satisfied in Ω, this problem typically has an exponential boundary layer at  = 1 and parabolic boundary layers at  = 0 and  = 1 (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.

Local Truncation Error.
We begin with analyzing the truncation error of the scheme about EDS0.For brevity we only consider the case that (, ) > 0 and (, ) > 0; similar analysis can be done for other cases.

Convergence Analysis.
For continuing our analysis we show a discrete maximum principle next.
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: 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  ∈  4 (Ω) and conditions for Lemmas 1 and 4 are satisfied, then where  =  1 /.
If (, ) ∈  1 (Ω), that is,  (,) (, ), (,  = 0, 1,  +  ≤ 1) are bounded by some constant, say by  3 .Since  01 ,  10 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 (, ) ≥ 0 and (, ) ≥ 0, but, in fact, in case of (, ) ≤ 0, (, ) ≤ 0 (see that it will not change the sign of   ( = 0, 1, 2, 3, 4) from ( 25)-( 29)), the method is still effective and similar theoretical analysis as we have shown above could be done.

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 Ω = (0, 1) × (0, 1).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  (, ) = Re, (,)=Re, (,)=0.(54) The exact solution of this problem is Figure 2 shows that there is no boundary layer when Re = 1, but when Re is large the boundary layers will emerge near  = 1 and  = 1.Notice that since the source term (, ) = 0, 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 3−16, which is almost the minimum positive machine-precision number (the minimum positive machine-precision number is 2.22045 * 10 −16 in doubleprecision 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 ℎ =  = 1/32.
We continue to implement our experiments; the convergence criterion for the Gauss Seidel iteration is chosen to be 10 −10 for the following problems.
We prefer to apply the Richardson extrapolation technique here to compute a fourth-order accurate solution since  , the difference solution at any mesh point (  ,   ) ∈ Ω ℎ, ; the Richardson extrapolation can be written as where  = 2 is the order of accuracy before the extrapolation and the order of accuracy will be increased to  + 2 after the extrapolation.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.
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  = 0.01, ℎ =  = 1/32.The errors of different schemes along  = 0.5 and  = 0.5 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  = 0 to  = 1 along the horizontal  axis, but EDS and EDS-RE get much better results (Table 8).

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

Figure 2 :
Figure 2: Exact solution of Example 7 for different values of Re.

Figure 4 :
Figure 4: Exact and EDS numerical solution for Example 8.

Figure 5 :
Figure 5: Exact solution of Example 9 for different values of .

Table 2 :
Absolute errors for EDS solution of Example 7 with different values of Re.

Table 3 :
Maximum absolute errors and rate of convergence, Example 8 with  = 100.

Table 4 :
Maximum absolute errors and rate of EDS0-RE and EDS-RE, Example 8 with  = 100.

Table 5 :
Maximum absolute errors in the computed solution, Example 8 with  = 1000.
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.

Table 6 :
Maximum absolute errors in the computed solution of Example 9.

Table 7 :
Maximum absolute errors and rate in the computed solution, Example 10 with  = 0.1.Note.The convergence criterion of the Gauss Seidel iteration for EDS-RE scheme is chosen to be 3 × 10 −16 to obtain highly accurate results.

Table 8 :
[31]absolute error in the computed solution, Example 10 with  = 0.001.Note.The last column of data comes from[31]of the MDC finite-element solution.