High-order blended compact difference schemes for the 3D elliptic partial differential equation with mixed derivatives and variable coefficients

Blended compact difference (BCD) schemes with fourth- and sixth-order accuracy are proposed for approximating the three-dimensional (3D) variable coefficients elliptic partial differential equation (PDE) with mixed derivatives. With truncation error analyses, the proposed BCD schemes can reach their theoretical accuracy, respectively, for the interior gird points and require 19 points compact stencil. They fully blend the implicit compact difference (CD) scheme and the explicit CD scheme together to make the derivation method and programming easier. The BCD schemes are also decoupled, which means the unknown function and its derivatives are separately resolved by different finite difference equations. Moreover, the sixth-order schemes are developed to solve the first-order derivatives, the second-order derivatives and the second-order mixed derivatives on boundaries. Several test problems are applied to show that the present BCD schemes are more accurate than those in the literature.


Introduction
In the paper, we study the steady 3D elliptic PDE au xx + bu yy + cu zz + pu x + qu y + ru z + d 1 u xy + d 2 u yz + d 3 where a, b, c, p, q, r, d 1 , d 2 , d 3 , s and f are sufficiently smooth functions and have the required partial derivatives on . The computational domain is a union of rectangular solid. ∂ is the boundary of . The elliptic PDE like (1) is of primary importance in various fields of engineering and science [1,2]. The numerical solution to the elliptic PDE has important interest in numerical analysis. In the past few decades, a lot of numerical methods, which include meshless methods [3], spline collocation methods [4], finite element methods [5], fast domain de-composition methods [6], Sinc-Galerkin methods [7] and finite difference methods [2,, have been proposed by many authors. Among the methods above, the finite difference method has been widely used in scientific research and engineering practice because of its simple structure, it being easy to understand and needing only a small amount of calculation.
The traditional difference schemes generally have low accuracy and calculation instability. And the non-compact high-order finite difference schemes are computational not efficient for a conventional type problem. However, high-order CD schemes have the advantages of having small discrete stencil, high-order accuracy, smaller element sensitivity and good numerical stability, which make it attractive in the fields of partial differential equations and computational fluid dynamics.
According to the difference of discrete objects, the compact difference scheme is divided into an explicit CD scheme and an implicit CD scheme. The former mainly discretize the differential equation (or model equation). The latter only discretizes derivatives of the unknown function involved in the model equation (or differential equation). Various special techniques have been developed rapidly based on implicit compact difference schemes [8][9][10][11][12][13][14][15][16][17][18]. In 1975, Adam [8], Kreiss [9] and Hirsh [10] proposed fourth-order CD schemes for different types of partial differential equations based on the Hermite formula, respectively. In Ref. [10], Hirsh solved the Burgers equation, the boundary layer problem and the driven cavity problem. In 1992, Lele [11] independently proposed a new type of CD method with quasi-spectral resolution and obtained a variety of symmetric compact difference schemes. Ma and Fu [12] proposed the upwind CD scheme with arbitrary-order accuracy in order to suppress the non-physical high frequency oscillation of the numerical solution near the shock wave. Deng and Zhang [13] developed the nonlinear CD scheme and overcame the non-physical oscillation of the symmetric compact difference scheme, which directly improved the quality of the stability simulation. Chu and Fan [14] derived a sixth-order combined CD scheme which is known as (CCD) scheme. In 1998, Mahesh [15] also proposed a CCD scheme of first-order and second-order derivatives with quasi-spectral resolution. Subsequently, the combined compact difference scheme is applied to resolve the wave propagation problem and the solution of the Navier-Stokes equation by Sengupta et al. [16,17]. On the basis of Ref. [14], Lee et al. [18] constructed a new CCD scheme (CCD2) to solve two-dimensional (2D) elliptic problems with mixed derivative. Meanwhile, many explicit CD schemes have been also developed. Using Taylor series expansion and the method of undetermined coefficients, Gupta et al. [19] developed a fourth-order polynomial CD scheme for the 2D convection-diffusion equation with variable coefficients. Dennis and Hudson [20] developed a fourth-order CD scheme for solving the 2D convection-diffusion equation by using Numerov approximation and applied it to solve the Navier-Stokes equations. Spotz and Carey [21], Li et al. [22] and Gupta [23] proposed a nine-point fourth-order CD method for the Navier-Stokes equations in the form of incompressible vorticity-stream function, respectively. Ananthakrishnaiah et al. [24] also developed a fourth-order CD scheme for the 3D elliptic equations by using a Taylor series expansion and the method of undetermined coefficients. Chen et al. [25] developed a perturbation fourth-order exponential compact difference scheme for convection-diffusion equations with constant and variable coefficients. Tian and Dai [26] and Radhakrishna-Pillai [27] proposed fourth-order exponential CD schemes for the 1D and 2D steady convection-diffusion equations, respectively. Mohanty [28] constructed a fourth-order CD scheme for the 3D nonlinear elliptic PDE. Gupta and Kouatchou [29] derived a 19-point fourth-order and a 27-point sixth-order CD schemes for the 3D Poisson equation by using the symbol software named Mathematica. Ge and Zhang [30] developed a 19-point fourth-order CD scheme for resolving the 3D linear elliptic equations by using Maple software. For the 3D steady convection-diffusion equation, Tian and Cui [31] and Zhang [32] proposed an explicit fourth-order CD scheme, respectively; Ma and Ge [33] developed a sixth-order CD scheme by using extrapolation technology; Mohamed et al. [34] proposed a fourth-order exponential CD scheme. Recently, Ma and Ge [35] developed a new class of BCD schemes to solve the 2D elliptic PDE with mixed derivatives by combining the advantages of the explicit CD and implicit CD. The BCD schemes cannot only achieve higher-order accuracy, but they also reduce the complexity of algebraic equations, so that they can be solved iteratively by the decoupling method.
As far as we know, there are few reports about the BCD schemes for the 3D elliptic equations with mixed derivatives, especially those with sixth-order accuracy. The main aim of the present paper is to extend our research work for 2D elliptic equations [35] to the 3D cases with variable coefficients and mixed derivatives to derive fourth-and sixth-order BCD schemes. The outline of this paper is organized as follows. Section 2 presents the two kinds of BCD schemes for the 3D elliptic PDE. In Sect. 3, we give truncation error analyses of the BCD schemes, respectively. Next, in Sect. 4, we compare our schemes with other schemes in the literature when numerical tests are conducted. Finally, some concluding remarks are given in Sect. 5.

The blended compact difference (BCD) schemes
In this section, the development of BCD formulations for Eq. (1) is briefly discussed. For convenience, the derivatives are symbolized as u, u x , u y , u z , u xx , u yy , u zz , respectively. The general Dirichlet boundary condition is considered. Here we use the notations

Inner grid points
Assume the problem domain to be cubical and construct on it a uniform Cartesian mesh of steps h x , h y and h z in the x, y and z directions, respectively. For convenience, we use a local coordinate system [30]. The approximate values of functions u, u x , u y , u z , u xx , u yy , u zz , The approximate values of its immediate 18 neighboring points are denoted by u l , u xl , u yl , u zl , u xxl , u yyl , u zzl , ∂ x ∂ y u l , ∂ y ∂ z u l , ∂ z ∂ x u l , l = 1, 2, . . . , 18, as in Fig. 1. We use the Taylor series expansions at point (x i , y j , z k ) for u x , u y , u z . We obtain here δ x , δ y and δ z are central difference operators for u x , u y , u z . Substituting Eqs.
(2)-(4) into Eq. (1), we obtain Differentiating Eq. (1) with respect to x, y and z, we get substituting Eqs. (6)-(8) into Eq. (5), we get where In order to get a fourth-order formulation for Eq. (9), all the derivatives {u xx , u yy , u zz , are approximated as follows: Here Eqs. (28)- (30) have been studied in Refs. [43,44]. Substituting Eqs. (28)-(40) into Eq. (9), we have Substituting δ x , δ y , δ z , δ 2 x , δ 2 y , δ 2 z , δ x δ y , δ y δ z , δ z δ x given in Appendix A at grid point 0 into Eq. (41). Then eliminating the sixth-order truncation error terms, we obtain the required fourth-order BCD scheme with 19 grid points Notice that there are four unknowns {u, u x , u y , u z } to be determined for fourth-order scheme. In order to match the system, in the inner field, u x , u y and u z are computed, respectively, by and to get a sixth-order compact formulation for Eq. (9), consider the following approximations for all the derivatives: Substituting Eqs. (46)-(61) into Eq. (9) and rearranging it, we have Substituting at grid point 0 into (62), and neglecting the truncation error terms, the 19-point sixth-order BCD scheme can be derived as follows: Notice that all the derivatives are separately computed and are demanded to achieve sixth-order accuracy. In the study, we directly use the sixth-order schemes of u x , u y , u z , u xx , u yy and u zz in Ref. [14] as well as the sixth-order CCD2 schemes of ∂ x ∂ y u, ∂ y ∂ z u and ∂ z ∂ x u in Ref. [18]. The sixth-order schemes of u x , u y , u z , u xx , u yy and u zz in x-, y-and z-directions are as follows: The mixed derivatives {∂ x ∂ y u, ∂ y ∂ z u, ∂ z ∂ x u} are computed with the nine-point sixth-order schemes in Ref. [18], which are expressed as follows: In the explicit fourth-order CD scheme [32], the third-and fourth-order derivatives of the truncation errors are represented by the original differential equation. But in the BCD schemes, the fifth-and sixth-order derivatives are represented by an unknown function, its first-and second-order derivatives, while the third-order derivatives are represented by means of the method in Ref. [32].

Boundary formulas for derivatives
It is of great significance to construct higher-order boundaries schemes for the BCD schemes at inner points. If the boundaries schemes are not well constructed, the accuracy and stability of the numerical solution will be affected. For the fourth-order BCD scheme, we adopt the consistent fourth-order boundary schemes [37], which have good stability and accuracy (see Appendix B).
For the sixth-order BCD scheme, because the first-order, second-order and secondorder mixed derivatives are unknown at the boundaries, so we need to construct sixthorder boundaries schemes for them. Firstly, we consider the left boundary. In order to obtain the sixth-order scheme of the left boundary (i = 0 and i = N x in x-direction) of u x , we assume that the unknown function and its first-order derivative u x have the following relationship: The left boundary format discrete template is shown in Fig. 2. Next, the unknown function u and its first-order derivative u x in Eq. (73) are expanded by a Taylor series at point i = 0. By matching the coefficients of the unknown function and its derivatives, we can obtain the following constrained linear equations: With the help the Matlab software, we obtain the sixth-order scheme of the left boundary: similarly, the remaining sixth-order schemes of the boundaries (right, down, up, rear and front) for {u x , u y , u z } are given as follows: u N x ,j,k -1877 300 u N x -1,j,k + 8u N x -2,j,k + 5 3 u N x -3,j,k -7u N x -4,j,k + 47 12 u N x -5,j,k -5 4 u N x -5,j,k (j = 0, 1, . . . , N y ; k = 0, 1, . . . , N z ), With a similar method, we can get the boundary schemes for {u xx , u yy , u zz } on boundaries. All sixth-order boundaries schemes (left, right, down, up, rear and front) for {u xx , u yy , u zz } are given as follows: (u xx ) 0,j,k -6(u xx ) 1,j,k (u yy ) i,0,k -6(u yy ) i,1,k (u yy ) i,N y ,k -6(u yy ) i,N y -1,k Finally, all sixth-order boundaries schemes for {∂ x ∂ y u, ∂ y ∂ z u, ∂ z ∂ x u} are derived, similar to deriving the boundaries schemes of {u x , u y , u z } (see Appendix C). Taking ∂ x ∂ y u as an example, the ∂ x ∂ y u can be regarded as (u x ) y or (u y ) x . All sixth-order boundaries schemes (left, right, down, up) for ∂ x ∂ y u are given as follows: With a similar method, we are able to obtain the boundaries schemes for {∂ y ∂ z u, ∂ z ∂ x u}.
All sixth-order boundaries schemes (rear, front, down, up,) are given for ∂ y ∂ z u as follows: All sixth-order boundaries schemes are given for ∂ x ∂ z u as follows: In the section, we obtain all sixth-order boundaries schemes which require more than 3 grid points. Although they are not compact in the traditional sense, the major compact structure is from the interior difference equations.

Computing process
Next, we will give computing process of the sixth-order BCD scheme. Computation of the fourth-order BCD is close to it, so we will not talk about it more.
Step 5: Computing u (m) by the present sixth-order BCD scheme (63) to compute the transport variable u.
If the condition F -Lu (m) < ε is satisfied, then stop. Otherwise, let m = m + 1 and repeat steps 2-5. Here m -the iteration number; L -the linear difference operator of the scheme (63); · -L ∞ norm; ε -the convergence tolerance.

Truncation error analysis
For simplicity, it is supposed that h x , h y and h z is equal to h. Truncation error analysis is given as follows. Using the Taylor series expansions at point (x i , y j , z k ) Differentiating Eq. (1) with respect to x, y and z, we get Substituting Eqs. (102)-(104) into Eq. (1) and rearranging it: Here, 16 , A 17 and F are defined in (9). In order to get a sixth-order compact scheme for Eq. (42), consider the following approximations for all the derivatives. These derivatives are approximated as follows: Notice that all the derivatives {u x , u y , u z , u xx , u yy , u zz , ∂ x ∂ y u, ∂ y ∂ z u, ∂ z ∂ x u} are calculated independently. Detailed derivation of the above discretization can be found in [20]. By substituting the truncation errors of the derivatives, we can get truncation error of the sixthorder BCD scheme as follows: Similarly, we can obtain a truncation error of the fourth-order BCD scheme as follows:

Numerical experiments
In order to verify the accuracy and reliability of the present BCD schemes, numerical experiments with five test problems which have analytical solutions are carried out in this section. All the test problems, where the right-hand function and the Dirichlet boundary conditions can be given using the analytical solutions, are defined on the unit cube domain = (0, 1)×(0, 1)×(0, 1). The hybrid biconjugate gradient stabilized method (BiCGstab (2)) is selected to resolve the resulting linear systems in the problems. The errors and the con-vergence order of the method are obtained according to the following definitions: Here u(x i , y j , z k ) is the exact solution. Error1 and Error2 are the maximum absolute errors estimated for two different grid step sizes h 1 and h 2 . For comparison, the proposed BCD schemes are used to compute the numerical solutions of all problems and the results with those computed by the explicit fourth-order CD scheme [24]. The maximum absolute errors and convergence orders with different values of h are listed in Tables 1-5, respectively. These results show clearly the BCD4 scheme, the explicit fourth-order CD scheme [24] and the BCD6 scheme can reach their fourth-and sixth-order accuracy and the BCD4 scheme gets a slightly better accurate solution than the explicit fourth-order CD scheme. However, the BCD6 scheme produces a much better     accurate solution than both the BCD4 scheme and the explicit fourth-order CD scheme [24]. (1) as

Conclusions
In this paper, we have constructed the fourth-and sixth-order BCD schemes for the 3D variable coefficients elliptic PDE with mixed derivatives. Firstly, based on Taylor series expansion and truncation error remainder, combined with the fourth-order Padé schemes of the first-order derivatives, a new fourth-order BCD scheme is constructed. In this new scheme, the unknown function and its first-order derivatives are regarded as the unknown variables in the calculation. Then, on the basis of the fourth-order BCD scheme proposed above, a new sixth-order BCD scheme is proposed by replacing the fifth-and sixth-order derivatives of the truncation errors with the linear combination of the unknown function, the first-and second-order as well as the second-order mixed derivatives. In other words, in the sixth-order BCD scheme, the unknown function, its first-and second-order derivatives as well as the second-order mixed derivatives are regarded as the unknown variables. At the same time, the sixth-order boundaries schemes of the first-order derivatives, the second-order derivatives and the second-order mixed derivatives are proposed. Finally, numerical results indicate that the present BCD schemes exhibit a very good resolution and high accuracy for all test problems.  The sixth-order approximation of left boundary of mixed derivative may be obtained from a relation of the form ∂ x ∂ y u 0,j,k + α∂ x ∂ y u 1,j,k = (u y ) x 0,j,k + α (u y ) x 1,j,k = a 0 (u y ) 0,j,k + a 1 (u y ) 1,j,k + a 2 (u y ) 2,j,k + a 3 (u y ) 3,j,k + a 4 (u y ) 4,j,k + a 5 (u y ) 5,j,k + a 6 (u y ) 6,j,k /h x . ( C . 1 ) Here j = 0, 1, . . . , N y ; k = 0, 1, . . . , N k the coefficients a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 (for the subscript see Fig. 3) and α are derived by matching the Taylor series coefficients of various orders.
The detailed derivation process is given below.