A higher-order blended compact difference (BCD) method for solving the general 2D linear second-order partial differential equation

A higher-order blended compact difference (BCD) scheme is proposed to solve the general two-dimensional (2D) linear second-order partial differential equation. The distinguishing feature of the present method is that methodologies of explicit compact difference and implicit compact difference are blended together. Sixth-order accuracy approximations for the firstand second-order derivatives are employed, and the original equation is also discretized based on a 9-point stencil, which is different from the work of Lee et al. (J. Comput. Appl. Math. 264:23–37, 2014). A truncation error analysis is performed to show that the scheme is of sixth-order accuracy for the interior grid points. Simultaneously, sixth-order accuracy schemes are proposed to compute the grid points on the boundaries for the firstand second-order derivatives. Numerical experiments are conducted to demonstrate the accuracy and efficiency of the present method.


Introduction
In this paper, we consider the general two-dimensional (2D) linear partial differential equation in the form a(x, y)u xx + b(x, y)u yy + c(x, y)u xy + p(x, y)u x + q(x, y)u y + r(x, y)u = f (x, y). (1) Here the unknown function u, the variable coefficient functions a, b, c, p, q, r, and the forcing function f are assumed to be continuously differentiable and have the required partial derivatives on Ω. Ω is a continuous rectangular domain. A suitable Dirichlet condition is prescribed on the boundary (∂Ω) in this paper.
Equation (1) is widely used in the fields of porous media flow [2,3] and when coordinate transformations are applied to a convection-diffusion equation on non-rectangular domains [4], which also generate Eq. (1). So it is both theoretically and practically important to investigate their accurate, stable and efficient numerical methods. In the past decades, high-order compact (HOC) difference methods have attracted increasing interests for solving such equations [1,[5][6][7][8]. The main feature of HOC methods is that they have higher-order accuracy (usually fourth-order, even higher) and higher spectral resolution with relatively few grid points (compact grid stencil) than the conventional difference methods. There are two kinds of compact difference methods. One is called the explicit compact difference method and the other is called the implicit compact difference method. For the explicit compact difference method, all the derivatives in partial differential equations are explicitly approximated by the nodal values of the target function. Usually, 9 grid points are used for solving the 2D problems and no more than 27 grid points (usually 19 points) for solving the three-dimensional (3D) problems. For instance, in Ref. [9], Gupta et al. derived a fourth-order compact finite-difference scheme for the 2D convectiondiffusion equation. Then they extended the method to the 2D elliptic equation without the mixed derivative term [5]. Karaa extended the work of Gupta et al. to the 2D elliptic and parabolic problems with mixed derivatives [8]. The basic idea of their methods lies in utilizing the Taylor expansions by 2D power series of all the functions involved in the differential equation at the reference grid point (i, j), substituting these expansions into the original equation and comparing the coefficients of x i y j to get the linear constraints on the unknown coefficients. In Ref. [8], the author considered the model equations with the mixed derivative terms, but the coefficients of the second derivatives u xx and u yy were equivalent to constants. In other words, if they are functions and a(x, y) = b(x, y), and c(x, y) = 0 in Eq. (1), there do not exist any explicit fourth-order compact difference schemes with 9-point grid stencil as was declared in Ref. [6]. And the authors of [6] also pointed out that explicit sixth-order compact difference schemes exist only if a(x, y) = b(x, y) = 1, c(x, y) = 0, and p(x, y) = q(x, y). Special cases are the Poisson equation and the Helmholtz equation. Some fourth-and sixth-order compact difference schemes are reported in Refs. [10][11][12][13][14][15][16]. Additionally, some fourth-order compact difference schemes for the convection-diffusion equations and the steady incompressible Navier-Stokes equations are reported in Refs. [17][18][19][20][21][22][23][24][25].
Different from the explicit compact difference method, the implicit compact difference method is to treat the derivatives and the target function as unknowns simultaneously, with all derivatives involved in computation. The advantage of the method is that they are easier to attain higher-order accuracy by using additional information of the derivatives that are possibly unnecessary for reality use than the explicit compact difference method. In the 1970s, Kreiss [26], Hirsh [27], and Adam [28] have proposed Hermitian compact techniques using fewer grid nodes (three instead of five) at each coordinate direction to solve partial differential equations (PDEs). Later on, Lele [29] proposed a class of compact finite-difference schemes with a range of spatial scales (spectral-like resolution). Chu and Fan [30] proposed a combined compact difference (CCD) scheme, which can be regarded as an extension of the standard Padé schemes discussed by Lele [29]. In [30], the CCD scheme is derived by using local Hermitian polynomials and is sixth-order accuracy with 3 points for one dimension (see Appendix 1). Fourier analysis shows that the CCD scheme has better spectral resolution than many other existing compact or non-compact highorder schemes. Fu et al. [31] developed an upwind fifth-order compact scheme. Deng and Zhang [32] developed a high-order accurate weighted compact nonlinear scheme. More recently, implicit compact difference schemes were used to solve hyperbolic equations and various kinds of fluid dynamics and engineering problems [33][34][35][36][37][38][39][40][41][42][43].
As we discussed above, for the general 2D linear second-order equations with variable coefficients and the mixed derivative term like (1), it is impossible to get an explicit fourthorder compact difference scheme with 9 grid points except under certain conditions mentioned above. But it is possible to get an implicit fourth-order, even sixth-order compact difference scheme for such equations with no more than 9 grid points. Very recently, Lee et al. [1] extended the work in [30] and proposed a combined compact difference (CCD2) scheme to solve Eq. (1). The highlight of this work is a sixth-order accurate difference scheme for the mixed derivative with a 9-point compact stencil at the interior region. It is a pity that just the fourth-order accuracy schemes on the boundaries are proposed in [1], which results in the convergence order of the CCD2 scheme is less than six for some problems. So, in this paper, we are aiming at developing a new sixth-order compact difference scheme for solving the general 2D linear partial differential equation (1). Our method makes use of the explicit compact difference scheme in which the derivatives of the original differential equation are discretized and the implicit compact difference scheme in which all derivatives with their own schemes are involved in computation independently, which is named the blended compact difference (BCD) scheme. Namely, the original governing equation is discretized with the 9-point compact stencil (like the explicit compact difference scheme) and a sixth-order compact difference scheme is proposed, but it still includes the representations of the first-and second-order derivatives (like the implicit compact difference scheme), which are computed by the sixth-order CCD schemes [30]. This scheme differs from the explicit compact difference method because the computations of the first-and second-order derivatives are also needed in the scheme. And it also differs from the implicit CCD2 scheme [1] because the governing equation is required to be discretized simultaneously in the present method. Although the computational cost increases, it is more accurate and more convenient to conduct computation than the method described in [1,30]. We will illustrate this conclusion in the following sections.
This paper is organized as follows. In Sect. 2, we present the method for deriving a sixthorder compact difference scheme for solving Eq. (1). Simultaneously, the sixth-order accuracy schemes are proposed to compute the grid points on the boundaries. In Sect. 3, truncation error analysis is done to show that the scheme is sixth-order accuracy for the interior grid points. Next, numerical tests are conducted in Sect. 4, and we compare the present scheme with other existing numerical methods in the literature. Finally, we give a brief conclusion in Sect. 5.

Methodology of sixth-order BCD scheme
In this section, we briefly discuss the development of BCD formulation for Eq. (1). Assume the problem domain to be rectangular Ω = [0, L x ] × [0, L y ], and we discrete the domain with uniform grids and set h x and h y to be the grid step-length in the x-and y-directions, respectively. The central mesh point (x i , y j ) is denoted by 0 and the other 8 mesh points at (x i ± h x , y j ), (x i , y j ± h y ), (x i ± h x , y j ± h y ), are denoted by numbers 1 -8 (see Fig. 1). For convenience, we denote the derivatives by {u, u x , u xx , u y , u yy , u xy }, respectively. We consider the general Dirichlet boundary condition. In the following sections, we will discuss how to implement the BCD scheme on a rectangular region.

Inner grid points
For each interior grid point (x i , y j ), x i = ih x , y j = jh y , i = 1, 2, . . . , N x -1, j = 1, 2, . . . , N y -1. There are 6 unknowns {u, u x , u y , u xx , u yy , u xy } to be determined and consequently we should provide 6 independent difference equations. Note that the original 3-point sixthorder CCD schemes for the first-and second-order derivatives as well as the 9-point sixthorder scheme of the mixed derivative u xy are proposed in [30] and [1], respectively. We simply borrow the sixth-order schemes for the five derivatives {u x , u y , u xx , u yy , u xy } that are listed in Appendix 1.
In order to get the BCD scheme, we use the Taylor series expansions at point (x, y).
Here δ x and δ y (see Appendix 2) are standard central difference operators for the firstorder derivatives. Substituting (2) and (3) into (1),we obtain the following modified differential equation corresponding to Eq. (1): We now replace all the third-order derivatives in Eq. (4). Differentiating Eq. (1) with respect to x and y, respectively: Substituting Eqs. (5)-(6) into Eq. (4),and rearranging it, we obtain Au xx +Bu yy +Cu xy +Du x +Ēu y +Ḡu xyy +Hu xxy + (r + pδ x + qδ y )u wherē To get a sixth-order compact scheme for Eq. (7), we use the following approximations for all the derivatives {u xx , u yy , u xxy , u xyy , u x (5) , u x (6) , u y (5) , u y (6) }: Substituting Eqs. (16)-(23) into Eq. (7) and rearranging it, we have Substituting the standard finite-difference operators (see Appendix 2) into (24) and omitting the truncation error terms, we obtain the following BCD scheme: All the coefficients are explicitly given as follows: Equation (25) is the new BCD scheme we developed. It involves the values of the unknown function on 9 grid points (like explicit compact difference scheme) and their first-and second-order derivatives on 3 grid points on each coordinate direction (like implicit compact difference scheme). All the derivatives are separately required to be approximated up to sixth-order accuracy. Namely, the first-order derivatives u x and u y are approximated by the sixth-order CCD interior schemes (59) and (61) in Appendix 1, respectively. The second-order derivatives u xx and u yy are approximated by the sixth-order CCD interior schemes (60) and (62) in Appendix 1, respectively. The However, the BCD scheme is decoupled by means of solving the unknown function and its various derivatives separately with an iteration process, so the derivation of the schemes, the algorithm design and programming are simple and convenient to operate.

Boundary grid points
In the above section, we have developed difference equations for each interior grid point.
Many applications involve computations in domains with non-periodic boundaries. In this section, we will introduce approximations for the first-and second-order derivatives for the boundary nodes. For non-periodic problems, these approximations are, of necessity, non-central or one-sided. Additional sixth-order expressions are needed to compute nodes at the boundaries nodes to close the system. Consider left boundary i = 0, the sixth-order schemes of the first-order derivative may be obtained from a relation of the form The 5 (for the subscript see Fig. 2) and α are derived by matching the Taylor series coefficients of various orders. By some tedious calculations, we are able to obtain the linear equations as shown below: And resolving it by Matlab software, we can get the results as follows: So, we can get the left boundary condition of u x for j = 0, 1, . . . , N y : (u x ) 0,j + 5(u x ) 1,j = -197 60 u 0,j -5 12 u 1,j + 5u 2,j -5 3 u 3,j + 5 12 u 4,j -1 20 The derivation of the following boundary approximations for the first-order and secondorder derivatives is exactly analogous to the above process. The boundary schemes are summarized now.
Right boundary condition of u x for j = 0, 1, . . . , N y : Bottom boundary condition of u y for i = 0, 1, . . . , N x : Top boundary condition of u y for i = 0, 1, . . . , N x : (u y ) i,N y + 5(u y ) i,N y -1 = 197 60 u i,N y + 5 12 Left boundary condition of u xx for j = 0, 1, . . . , N y : (u xx ) 0,j -6(u xx ) 1,j = -403 18 u 0,j + 33u 1,j - Right boundary condition of u xx for j = 0, 1, . . . , N y : Bottom boundary condition of u yy for i = 0, 1, . . . , N x : Top boundary condition of u yy for i = 0, 1, . . . , N x : (u yy ) i,N y -6(u yy ) i,N y -1 = -403 18 In the end, we will derive boundary conditions of the mixed derivative. We can regard the second mixed derivative u xy as the first partial derivative of u x (or u y ) with respect to y (or x). Similar to the first-order derivative boundary scheme, by using the 2D Taylor expansions and some tedious calculations (see Appendix 3), we are able to deduce some difference schemes for approximating u xy for grid points on the boundaries. A group of sixth-order accurate equations on the four boundaries can be given here.
Left boundary condition of u xy for j = 0, 1, . . . , N y : Right boundary condition of u xy for j = 0, 1, . . . , N y : Top boundary condition of u xy for i = 0, 1, . . . , N x : Bottom boundary condition of u xy for i = 0, 1, . . . , N x : With more than 4 grid points along one direction, the above boundary difference equations are not strictly compact in the traditional sense, but they are necessary to retain the high-order accuracy of the BCD scheme. Nevertheless, we remark that these equations are only used for boundary grid points. The interior difference equations contribute to the major compact structure.
Note that if the Dirichlet boundary condition is replaced by the Neumann boundary condition, the equation (1) can also be solved because the function u is unknown at the boundaries while the values of the first-order derivatives u x and u y at the boundaries are known. According to Eqs. (28)-(31), the function u at the boundaries can be calculated by the sixth-order schemes as follows: u 0,j = 60 -h x (u x ) 0,j -5h x (u x ) 1,j -5 12 u 1,j + 5u 2,j -5 3 u 3,j + 5 12 u 4,j -1 20 u 5,j 197, (40)

Iterative procedure
The BCD scheme is decoupled by means of solving the unknown function u and its various derivatives u x , u y , u xx , u yy , u xy separately with an iteration process. A pseudo code of the BCD scheme is listed as follows: Give any initial guess u (0) , u (0) x , u (0) y , u (0) xx , u (0) yy and u 0 xy For k = 1, 2, . . . do Compute u (k) x by Scheme (59) in Appendix 1 and Schemes (28) and (29). Compute u (k) y by Scheme (61) in Appendix 1 and Schemes (30) and (31). Compute u (k) xx by Scheme (60) in Appendix 1 and Schemes (32) and (33).
Here k is iteration number. L is the linear difference operator for representation of the scheme (25). · is a sort of norm. η is the convergence tolerance.

Truncation error analysis
A brief truncation error analysis is given here. For the sake of simplicity, we suppose that h x and h y are equal to h, using the Taylor series expansions at point (x, y).

Numerical experiments
In this section, we conduct numerical experiments with three test problems chosen from the literature to test the high-order accuracy of the BCD scheme. All problems are defined on the unit square domain Ω = (0, 1) × (0, 1). For all three problems, the right-hand functions f and the Dirichlet boundary conditions are prescribed to satisfy the given analytic solution. We select the hybrid biconjugate gradient stabilized method (BiCGStab(2)) [44] to solve the resulting linear systems in all test problems. All iterative procedures are started with zero initial guesses and are terminated when the Euclidean norm of the residual vector is reduced by 10 10 . The code is written in Fortran 77 programming language with double precision arithmetic. All computations are run on a personal computer with an Intel (R) core (TM) i3-5005U double 2.0 GHz CPU and 4 GB memory. The error and the convergence rate of the method are computed with the following definition: where u(x i , y j ) is the exact solution. Error1 and Error2 are the maximum absolute errors estimated for the two different grid step sizes h 1 and h 2 .

Problem 1: With variable coefficients [1, 45]
Consider the following differential equation in the presence of a source term: The analytic solution is u = x 3 y 2 + x sin(x) cos(xy).
We notice that Problem 1 has a non-periodic boundary condition. For comparison, the results of the CCD2 scheme [1] are also given. Firstly, we use different mesh sizes from 1/8 to 1/64 to evaluate the computed accuracy order. Table 1 gives the maximum absolute errors and the convergence rate for Problem 1. We can clearly see that the BCD scheme obtains sixth-order accuracy and gets a more accurate solution than the CCD2 scheme. The CCD2 scheme cannot reach sixth-order accuracy because the fourth-order scheme is used for the boundaries computation, which influences the whole computed accuracy.

Problem 2: With high anisotropy ratio [1]
Consider the following differential equation in the presence of a source term: The analytic solution is u = sin(πx) sin(πy).
The second model problem is an anisotropy problem. ε is chosen to reflect the magnitude of the anisotropy. Note that ε = 1 gives an isotropic equation, and the degree of anisotropy increases when ε becomes smaller. The maximum absolute errors, the convergence rate and CPU time for ε = 10 -1 and 10 -3 , using CCD2 scheme and the BCD scheme, are given in Table 2 and Table 3, respectively. Numerical results show that the maximum  absolute errors of the BCD scheme and the CCD2 scheme are almost identical. They are seventh-order accurate and their excellent performances are obviously independent of the anisotropy ratio ε = 10 -1 and ε = 10 -3 . When ε = 10 -1 , the BCD scheme yields a slightly better solution than the CCD2 scheme. However, when ε = 10 -3 , the CCD2 scheme gets a slightly better solution than the BCD scheme. We notice that values on the boundaries are zero (periodic boundary conditions), for the CCD2 scheme, the fourth-order scheme for periodic boundary conditions does not influence the results of inner grid points. In Table 2 and Table 3 it can be shown that the CPU time is related to the number of grid points and anisotropy ratio. We find the same anisotropy ratio in which the number of grid points becomes higher when the CPU time is higher and the same number of grid points with which anisotropy ratio is higher when the CPU time is higher.
For comparison, we use the BCD, the FOC [8], and the CCD2 schemes [1] to compute the numerical solutions of Problem 3. Tables 4-6 give the maximum absolute errors, the convergence rate and CPU time, when Re = 10 2 , 10 4 , 10 6 , respectively. It is seen that the BCD and CCD2 schemes are almost not influenced by the increase of Reynolds number, while the FOC scheme gradually loses its fourth-order accuracy. Numerical results also show that the maximum absolute errors of the BCD scheme and the CCD2 scheme are almost identical. With seventh-order accuracy their excellent performances are obviously independent of the Re numbers. The BCD scheme yields slightly more accurate solution than the CCD2 scheme does. We can also see that the CPU time is related to the Re numbers, which means that Re becomes bigger when the CPU time is higher with the same number of grid points.    Finally, the maximum absolute error and the convergence rate of Problem 1, Problem 2 with ε = 10 -1 and Problem 3 with Re = 10 are given in Table 7 when h x is not equal to h y . Numerical results also show that the BCD scheme can reach its theoretical sixth-order accuracy.

Conclusions
In this paper, the BCD scheme is proposed to solve the general two-dimensional linear second-order partial differential equation. A truncation error analysis is done to show that the BCD scheme is sixth-order accuracy for the interior grid points. Besides, the sixthorder accuracy difference schemes are also proposed to compute the first-and secondorder derivatives for grid points on the boundaries. The superiority of the present method is that it fully exploits the merits of the explicit compact difference and implicit compact difference methods. At least three important conclusions are obtained. (i) The present method reaches sixth-order accuracy for smooth solution problems, even for those with large first derivative terms or anisotropy problems. (ii) We perform a sixth-order computation for the grid points on the boundaries, while Ref. [1] uses the fourth-order scheme, which decreases the accuracy order of solution for those problems with non-periodic boundary conditions. (iii) Our method is decoupled, i.e., we can solve the unknown function u and its various derivatives separately with an iteration process. Its advantage is that the derivation of the scheme, the algorithm design and programming are simple and easy to operate in the extension to high-dimensional problems. Especially for the problems of complex flow and heat transfer, it is relatively easy to construct a discrete scheme with high precision and for their boundary conditions. Numerical experiments are conducted to demonstrate the accuracy of the present scheme. It is shown that the present method is more accurate than those in the literature.

Appendix 1: The sixth-order schemes of the first-and second-order derivatives as well as the mixed derivative
For the first-and second-order derivatives along the x-and y-directions [30] For the mixed derivative u xy [1] (u xy ) 0 + 1 16 , δ x δ y u 0 = u 5u 6 + u 7u 8 4h x h y .

Appendix 3: The sixth-order schemes of the mixed derivative on the boundaries
The sixth-order scheme of left boundary of the mixed derivative may be obtained from a relation of the form (u xy ) 0,j + α(u xy ) 1,j = ∂u y ∂x 0,j + α ∂u y ∂x 1,j = a 0 (u y ) 0,j + a 1 (u y ) 1,j + a 2 (u y ) 2,j + a 3 (u y ) 3,j + a 4 (u y ) 4,j + a 5 (u y ) 5,j + a 6 (u y ) 6,j /h x , where j = 0, 1, . . . , N y , the coefficients a 0 , a 1 , a 2 , a 3 , a 4 , a 5 (for the subscript see Fig. 2) and α are derived by matching the Taylor series coefficients of various orders. The detailed derivation process is given here.