A Family of Sixth-Order Compact Finite-Difference Schemes for the Three-Dimensional Poisson Equation

We derive a family of sixth-order compact finite-difference schemes for the three-dimensional Poisson’s equation. As opposed to other research regarding higher-order compact difference schemes, our approach includes consideration of the discretization of the source function on a compact finite-difference stencil. The schemes derived approximate the solution to Poisson’s equation on a compact stencil, and thus the schemes can be easily implemented and resulting linear systems are solved in a high-performance computing environment. The resulting discretization is a one-parameter family of finite-difference schemes which may be further optimized for accuracy and stability. Computational experiments are implemented which illustrate the theoretically demonstrated truncation errors.


Introduction
In this article, we derive a family of sixth-order compact finite-difference schemes for Poisson's equation. Let Ω be an open, bounded regular hexahedron in R 3 , and consider −Δu f, in Ω, u g, on ∂Ω. 1.1 In order to obtain the high computational efficiency and the performance of higherorder methods, a complete characterization of the truncation error in the respective multivariable form must be formulated and minimized. The schemes are designed based on local multivariate Taylor expansion of the solution 1 constrained by higher-order derivatives of the equation about the center of the three-dimensional compact difference stencil. Higherorder compact finite-difference schemes for approximating elliptic equations have been well studied, 2-5 , since they achieve high-order accuracy without significant increase in the 2 Advances in Numerical Analysis resolution of the grid points. The technique of minimizing the truncation error has been extended to develop higher-order compact schemes 6 and for other application problems 7, 8 . In other approaches 1, 4 , the univariate Taylor series expansion is used to derive the finite-difference approximations of the individual derivative in terms of the differential equation and then coupled to obtain the numerical schemes for multiple spatial dimensions. Subsequently, the truncation errors are formulated to assess the accuracy of the schemes. Our approach in deriving finite difference schemes utilizes the fully multivariate Taylor series expansion rather than univariate expansions in each coordinate direction. First, the multivariate approximations to the unknown and the source function are substituted into the partial differential equation about the center of the local compact stencil. Then, the formal error for discretizing the equation is formulated using the discrete approximations of the unknowns, the sources, and weight parameters to mimic the derivatives in the equation. By determining the parameters to annihilate the leading coefficients of the error, the parameterbased fourth-order compact schemes are derived. By setting the parameter to zero, the traditional fourth-order compact scheme 1 is recovered. Numerical experiments show that the resulting schemes are much more stable and robust when the remaining free parameter is chosen in an effective manner.
In order to discuss compact stencil schemes in three spatial dimensions, we must first number the points on the compact stencil. For this article, the stencil will be labeled according to the diagram in Figure 1. In order to describe schemes for Poisson's equation in three spatial dimensions, we notice that the numerical approximation of 1.1 results in the following matrix equation: where the entries of H are determined by the stencil weights associated with the Laplacian operator Δ : ∂/∂x 2 ∂/∂y 2 ∂/∂z 2 , and the entries of Q are determined by the stencil weight associated with the right-hand side function f. We denote the collection of weights associated with a particular stencil as w 0 , w 1 , . . . , w 26 , where the weights are labeled according to Figure 1. Also, notice that in our derivation, we use stencil weights which are symmetric with respect to the coordinate axes and equal in each spatial direction. Thus it makes sense to say that any compact stencil in three spatial dimensions is determined by one of four values, the stencil weight at the center of the compact stencil denoted w 0 , the stencil weights in the directions of the coordinate axes i.e., steps of size h are taken in only one direction denoted w 1 · · · w 6 , the stencil weights where steps of size h are taken in two coordinate directions denoted w 7 · · · w 18 , and the stencil weights where steps of size h are taken in all three coordinate directions denoted w 19 · · · w 26 . By this logic, three-dimensional compact stencil for Poisson's equation may be described by eight values, w 0 , w 1 , w 7 , and w 19 for each of the matrices H and Q. For simplicity, we label these values α 0 , α 1 , α 7 , α 19 weights for the matrix H, and values β 0 , β 1 , β 7 , β 19 weights for the matrix Q. This nomenclature is illustrated in Figures 2 and 3.
In this article, we derive several schemes for 1.1 , these schemes we will label HOC4, HOC6 1 , HOC6 2 , and HOC6. HOC4 is defined as the traditional fourth-order scheme for Poisson's equation, that is,  We then derive the scheme HOC6 1 , which is defined as the sixth-order scheme in which the weights of H are selected in the usual way, and the weights of Q are chosen in a special way in order to remove the sixth-order error term. This scheme is given by where w 7 is determined by the formula, Advances in Numerical Analysis where ξ 1 : ξ 3 : ∂u 6 ∂x 2 ∂y 2 ∂z 2 .

1.6
This scheme shows that the stencil can be selected in a special way in order to annihilate the fourth-order error term and therefore produce a sixth-order scheme. Following this, we derive the scheme HOC6 2 in which the weights are rearranged from HOC6 1 in order so the fourth order error term may be eliminated a priori. In order to do this, we must alter the coefficients of the matrix H. This scheme is defined as follows: where w 7 is given by the formula Finally, we give the most general sixth-order scheme, denoted HOC6, which is a one parameter family of sixth-order schemes. This scheme is derived in the same way as HOC6 2 Advances in Numerical Analysis 5 and utilizes all 27 stencil weights for each of the matrices H and Q. The scheme HOC6 is given as follows: HOC6 : where w 7 and w 19 satisfy the relationship: Along with the definition of HOC6 in 1.9 -1.10 , we show how the weights can be effectively approximated so that the values of the fourth order partial derivatives of the right-hand side function f do not have to be calculated analytically in order to obtain a sixth-order scheme. The paper is outlined as follows. In Section 2, a derivation for the fourth-order compact scheme HOC4 is given which begins with the intuitive second-order scheme and eliminates the second-order error term. This illustrates the manner in which we will derive the sixthorder schemes later in the paper. In Section 3, a one-parameter family of fourth-order schemes is given and information about the truncation error is used in order to derive the schemes HOC6 1 and HOC6 2 . In Section 4, the complete one-parameter family of sixth-order schemes is presented; HOC6 is derived from a three-parameter family of fourth-order schemes. In Section 5, the computational implementation is discussed and numerical experiments are included which illustrate the convergence of the complete family of schemes, HOC4, HOC6 1 , HOC6 2 , and HOC6. In Section 6, stability of the schemes is discussed in terms of the conditioning of the mass matrix and in terms of approximating the solution of a generalized eigenvalue problem. Finally, in Section 7, conclusions are given.

The Fourth Order Scheme, HOC4
In order to illustrate the ideas in the sequel, we provide a derivation of the fourth-order HOC scheme for Poisson's equation by approximating the second-order error term on the local stencil and utilizing this approximation in the scheme.
From this figure, we immediately see that a second-order accurate scheme for Poisson's equation is defined by approximating each of the second-order derivatives in 1.1 , We will define the truncation error for this scheme and all subsequent schemes in this article by subtracting terms approximating u from the terms approximating f. For the second-order difference scheme 2.1 , the truncation error is given by In order to derive the fourth order scheme, we simply approximate the second-order terms in the truncation error on the local stencil. We, however, notice that the fourth order partial derivatives, u xxxx , u yyyy , u zzzz , cannot be approximated on the local stencil. However, we can utilize the differential equation 1.1 in order to rewrite the second-order error in a way which can be approximated on the local stencil.
Using the relationships, we see that the second-order error term may be rewritten as Next, using the facts that multiplying by their coefficients in 2.4 and moving to the other side of 2.2 , we have that a fourth order scheme is given by We will refer to this scheme as HOC4. This is the scheme found most prevalently

A One-Parameter Family of Fourth Order Schemes and the Sixth-Order Scheme HOC6 1
In this section, we consider a one parameter family of fourth order schemes and show that the schemes can be selected in a grid-dependent fashion which eliminates the fourth order error term and therefore determines a sixth-order scheme. An alternate stencil for the mass matrix and the one-parameter family of stencils for the mass matrix for the fourth order schemes appear in Figures 6 and 7 respectively. We may obtain the following expression for the error:

Advances in Numerical Analysis
where ξ 1 , ξ 2 , and ξ 3 are defined as in 1.6 . Eliminating the fourth order error term in the above equation, we have

A Priori Error Minimization
Obviously, although the scheme presented in the previous section certainly is valid, it will not be particularly useful unless we can derive a scheme which can be determined utilizing only information about the source function f. In order to derive a more suitable sixth-order scheme, we perform manipulations which will allow us to rewrite almost all of the terms in the error equation in terms of f, and the remaining terms will be approximated on the compact stencil. First, recall that the fourth order term for the error was given by First, making the substitutions, we obtain the following equivalent form for 3.3 : Next, making the substitutions we obtain the following equivalent form for 3.3 : Advances in Numerical Analysis Noticing that we can approximate the final term of the expression on the compact stencil, we obtain a sixth-order approximation scheme by setting f xxxx f yyyy f zzzz f xxyy f xxzz f yyzz 3.8 and using the following approximation: Adding this approximation to the stencil for u, we obtain the stencil which appears in Figure 8.

A Two-Parameter Family of Sixth-Order Schemes, HOC6
In the previous section, we illustrated how one can derive a sixth-order scheme in which stencil weights are selected in a grid-dependent fashion based only upon values of the source function f.

4.1
Notice that this is a three-parameter family of schemes, which may be selected in any particular fashion. It is our intention, however, to select the coefficients α 19 , β 7 , β 19 in such a way as to eliminate the fourth order term from the truncation error.
For this particular family of schemes, the error expression is given by Using similar analysis as in the previous section, we immediately obtain that values for the parameters for which the fourth order error term is eliminated are given by α 19 1/30 i.e., the stencil given in Figure 8 and

Computational Implementation
Notice that the term in 4.3 is prohibitive; since in order to define weights according to this formula, various fourth order partial derivatives of the source function f must be known exactly. However, for our computational implementation, we approximate each of the fourth order partial derivatives of f using difference approximations. The fourth order difference approximations of the mixed derivative terms f xxyy , f xxzz , f yyzz may be approximated in the usual way; see, for example, 2.6 , in order to obtain the formula To approximate the other fourth order partial derivative terms f xxxx , f yyyy , f zzzz , we notice that these quantities cannot be approximated on the compact stencil. However, using the fact that, for example, and utilizing the value h h/2, we obtain the following approximation formula: where f 0 : f x, y, z ,

5.4
The implementation was carried out using the C programming language, and matrices were assembled in compressed row sparse matrix form and solved using the Preconditioned BiConjugate Gradient Stabilized Method implemented in the Sparselib sparse matrix library available from NIST. Computations were performed on a Dell desktop computer with a 3.19 GHz Intel Processor and 1.99 GB of RAM. Notice that in order to solve the three-dimensional Poisson equation on a unit box with h 1/64, a 250, 047 × 250, 047 sparse matrix must be solved. However, using the hardware and software described, we were able to obtain solutions when h 1/64 in about 13 minutes.

Computational Experiments
In this section, we illustrate the experimental convergence rates for each of our schemes HOC4, HOC6 1 , HOC6 2 , and HOC6. It is interesting to note that the scheme HOC6 allows for a free parameter, so that the schemes may be optimized for additional accuracy or computational stability. Also, notice that in each of the cases, HOC6 1 represents the minimum error obtained. However, recall from our derivation that weights for HOC6 1 are determined according to the solution u, whereas the weights for HOC6 2 and HOC6 are minimized a priori.

5.7
Results are given for six cases, HOC4, HOC6 1 , HOC6 2 β 19 0 , HOC6 β 7 0 , HOC6 β 7 β 19 , and HOC6 β 7 2β 19 . Table 2 summarizes our results, illustrating the experimental convergence rates for each of the schemes. Notice that in both of the experiments the choice β 7 2β 19 yields the best a priori error for the scheme HOC6. This choice was selected in such a way that the error contributions from β 7 and β 19 are identical see 4.3 .

Stability of Schemes
In this section, we illustrate that the families of fourth-order schemes and thus the sixth-order scheme will exhibit increased computational stability when compared to the traditional fourth order scheme, a property which makes the alternative schemes desirable for implementation in applied problems. In order to do this, we test the fourth order methods by solving the corresponding eigenvalue problem −Δu λu in Ω, u 0, on ∂Ω.
6.1 Proceeding as before, we approximate the eigenvalue problem 6.1 by Hu λQu, 6.2 where H and Q are the stiffness and mass matrices formed using the schemes described above. We utilize the weights given in 4.1 .
An initial numerical concern is the diagonal dominance of the H matrix. We note that the restriction on α 19 for row diagonal dominance is The restrictions on β 7  Next, we investigate the computational accuracy and the conditioning of solving the approximate matrix equation 6.2 . Our preliminary results show that α 19 0 i.e., 19 points for H can adequately describe families of fourth-order schemes for solving the generalized eigenvalue problem 6.1 . In Figure 9, we illustrate the conditioning of the mass matrix Q depending on the choices of β 7 and β 19 as described in 6.6 . We observe that Q is best conditioned when β 7 β 19 1/48 and is the worst for the case β 7 0, β 19 0. Moreover, any selections of β 7 , β 19 away from this case yield an improvement in the conditioning of matrix Q. Now, we analyze the stability of the proposed schemes by solving 6.1 with Ω : 0, 1 × 0, 1 × 0, 1 . We compare the traditional fourth-order scheme HOC4 with 7 grid points for the mass matrix Q β 7 0, β 19 0 and the case with 19 grid points for Q β 7 1/48, β 19 0 in Figure 10.
The analysis shows that the alternative fourth-and sixth-order schemes provide credible stable alternative schemes for large-scale application problems 9, 10 to the traditional fourth order discretization with β 7 0, β 19 0 . The mass matrix for the traditional discretization becomes more and more poorly conditioned as the size of Q increases as illustrated in Figure 9. On the other hand, for instance, the condition number of Q with β 7 1/48, β 19 0 stays approximately close to 1.8 almost independent of the size of Q and provides comparable accuracies as also illustrated in Figure 10 to that of the traditional discretization.

Conclusion
In this article, we have derived and demonstrated a family of sixth-order compact finitedifference schemes for Poisson's equation in three spatial dimensions. Such schemes can be extended in order to derive sixth-order schemes for stationary and transient propagation problems in two or three spatial dimensions. It will be the subject of future work to reframe the concept of compact finite differencing in terms of a finite volume scheme for conservative form flow equations, and demonstrate the utility of these schemes in applied problems.