Research Article  Open Access
A Family of SixthOrder Compact FiniteDifference Schemes for the ThreeDimensional Poisson Equation
Abstract
We derive a family of sixthorder compact finitedifference schemes for the threedimensional Poisson's equation. As opposed to other research regarding higherorder compact difference schemes, our approach includes consideration of the discretization of the source function on a compact finitedifference 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 highperformance computing environment. The resulting discretization is a oneparameter family of finitedifference schemes which may be further optimized for accuracy and stability. Computational experiments are implemented which illustrate the theoretically demonstrated truncation errors.
1. Introduction
In this article, we derive a family of sixthorder compact finitedifference schemes for Poisson's equation. Let be an open, bounded regular hexahedron in , and consider
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 higherorder derivatives of the equation about the center of the threedimensional compact difference stencil. Higherorder compact finitedifference schemes for approximating elliptic equations have been well studied, [2–5], since they achieve highorder accuracy without significant increase in the resolution of the grid points. The technique of minimizing the truncation error has been extended to develop higherorder 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 finitedifference 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 fourthorder compact schemes are derived. By setting the parameter to zero, the traditional fourthorder 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 are determined by the stencil weights associated with the Laplacian operator , and the entries of are determined by the stencil weight associated with the righthand side function . We denote the collection of weights associated with a particular stencil as , 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 ), the stencil weights in the directions of the coordinate axes (i.e., steps of size are taken in only one direction) (denoted ), the stencil weights where steps of size are taken in two coordinate directions (denoted ), and the stencil weights where steps of size are taken in all three coordinate directions (denoted ).
By this logic, threedimensional compact stencil for Poisson's equation may be described by eight values, , , , and for each of the matrices and . For simplicity, we label these values , , , weights for the matrix , and values , , , weights for the matrix . This nomenclature is illustrated in Figures 2 and 3.
In this article, we derive several schemes for (1.1), these schemes we will label HOC, HOC, HOC, and HOC6. HOC4 is defined as the traditional fourthorder scheme for Poisson's equation, that is,
We then derive the scheme HOC, which is defined as the sixthorder scheme in which the weights of are selected in the usual way, and the weights of are chosen in a special way in order to remove the sixthorder error term. This scheme is given by where is determined by the formula, where This scheme shows that the stencil can be selected in a special way in order to annihilate the fourthorder error term and therefore produce a sixthorder scheme. Following this, we derive the scheme HOC in which the weights are rearranged from HOC6 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 . This scheme is defined as follows: where is given by the formula
Finally, we give the most general sixthorder scheme, denoted HOC6, which is a one parameter family of sixthorder schemes. This scheme is derived in the same way as HOC6 and utilizes all 27 stencil weights for each of the matrices and . The scheme HOC6 is given as follows: where and 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 righthand side function do not have to be calculated analytically in order to obtain a sixthorder scheme.
The paper is outlined as follows. In Section 2, a derivation for the fourthorder compact scheme HOC4 is given which begins with the intuitive secondorder scheme and eliminates the secondorder error term. This illustrates the manner in which we will derive the sixthorder schemes later in the paper. In Section 3, a oneparameter family of fourthorder schemes is given and information about the truncation error is used in order to derive the schemes HOC and HOC. In Section 4, the complete oneparameter family of sixthorder schemes is presented; HOC6 is derived from a threeparameter family of fourthorder 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, HOC, HOC, 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.
2. The Fourth Order Scheme, HOC4
In order to illustrate the ideas in the sequel, we provide a derivation of the fourthorder HOC scheme for Poisson's equation by approximating the secondorder error term on the local stencil and utilizing this approximation in the scheme.
From this figure, we immediately see that a secondorder accurate scheme for Poisson's equation is defined by approximating each of the secondorder derivatives in (1.1), We will define the truncation error for this scheme (and all subsequent schemes in this article) by subtracting terms approximating from the terms approximating . For the secondorder difference scheme (2.1), the truncation error is given by In order to derive the fourth order scheme, we simply approximate the secondorder terms in the truncation error on the local stencil. We, however, notice that the fourth order partial derivatives, , , , cannot be approximated on the local stencil. However, we can utilize the differential equation (1.1) in order to rewrite the secondorder error in a way which can be approximated on the local stencil.
Using the relationships, we see that the secondorder 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 in the literature for Poisson's equation. The stencils for HOC4 are illustrated in Figures 4 and 5.
3. A OneParameter Family of Fourth Order Schemes and the SixthOrder 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 griddependent fashion which eliminates the fourth order error term and therefore determines a sixthorder scheme. An alternate stencil for the mass matrix and the oneparameter 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: where , , and are defined as in (1.6). Eliminating the fourth order error term in the above equation, we have
3.1. 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 . In order to derive a more suitable sixthorder scheme, we perform manipulations which will allow us to rewrite almost all of the terms in the error equation in terms of , 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): Noticing that we can approximate the final term of the expression on the compact stencil, we obtain a sixthorder approximation scheme by setting and using the following approximation: Adding this approximation to the stencil for , we obtain the stencil which appears in Figure 8.
4. A TwoParameter Family of SixthOrder Schemes, HOC6
In the previous section, we illustrated how one can derive a sixthorder scheme in which stencil weights are selected in a griddependent fashion based only upon values of the source function . For completeness, in this section we give the most general family of fourth order HOC schemes for Poisson's equation, and from that derive the complete family of sixthorder HOC schemes.
In order to derive the complete family of sixthorder schemes for Poisson's equation, we will first state the complete family of fourth order schemes for Poisson's equation, which are
Notice that this is a threeparameter family of schemes, which may be selected in any particular fashion. It is our intention, however, to select the coefficients 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 (i.e., the stencil given in Figure 8) and which yields a one parameter family of sixthorder schemes.
5. 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 must be known exactly. However, for our computational implementation, we approximate each of the fourth order partial derivatives of using difference approximations.
The fourth order difference approximations of the mixed derivative terms ( 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 (, , ), we notice that these quantities cannot be approximated on the compact stencil. However, using the fact that, for example, and utilizing the value we obtain the following approximation formula: where
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 threedimensional Poisson equation on a unit box with , a sparse matrix must be solved. However, using the hardware and software described, we were able to obtain solutions when in about 13 minutes.
5.1. Computational Experiments
In this section, we illustrate the experimental convergence rates for each of our schemes HOC4, HOC, HOC, 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, HOC represents the minimum error obtained. However, recall from our derivation that weights for HOC are determined according to the solution , whereas the weights for HOC and HOC6 are minimized a priori.
Experiment 1. Let . For this example, we implement the schemes HOC4, HOC, HOC, and HOC6 for taken as the true solution of (1.1) and the schemes HOC4, HOC, HOC, and HOC6 defined as in (1.3), (1.4), (1.7), and (1.9), respectively. The weights satisfy the approximate formula Results are given for six cases, HOC4, HOC, HOC (), HOC6 (), HOC6 (), and HOC6 (). Table 1 summarizes our results, illustrating the experimental convergence rates for each of the schemes.

Experiment 2. Let . For this example, we implement the schemes HOC4, HOC, HOC, and HOC6 for taken as the true solution of (1.1) and the schemes HOC4, HOC, HOC, and HOC6 defined as in (1.3), (1.4), (1.7), and (1.9), respectively. The weights satisfy the approximate formula where Results are given for six cases, HOC4, HOC, HOC (), HOC6 (), HOC6 (), and HOC6 (). Table 2 summarizes our results, illustrating the experimental convergence rates for each of the schemes. Notice that in both of the experiments the choice yields the best a priori error for the scheme HOC6. This choice was selected in such a way that the error contributions from and are identical (see (4.3)).

6. Stability of Schemes
In this section, we illustrate that the families of fourthorder schemes and thus the sixthorder 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
Proceeding as before, we approximate the eigenvalue problem (6.1) by where and 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 matrix. We note that the restriction on for row diagonal dominance is The restrictions on and for row diagonal dominance of the mass matrix are with feasible positivity solutions as or
Next, we investigate the computational accuracy and the conditioning of solving the approximate matrix equation (6.2). Our preliminary results show that (i.e., points for ) can adequately describe families of fourthorder schemes for solving the generalized eigenvalue problem (6.1). In Figure 9, we illustrate the conditioning of the mass matrix depending on the choices of and as described in (6.6). We observe that is best conditioned when and is the worst for the case , . Moreover, any selections of , away from this case yield an improvement in the conditioning of matrix .
(a)
(b)
Now, we analyze the stability of the proposed schemes by solving (6.1) with . We compare the traditional fourthorder scheme HOC4 with grid points for the mass matrix (, ) and the case with grid points for (, ) in Figure 10.
(a)
(b)
The analysis shows that the alternative fourth and sixthorder schemes provide credible stable alternative schemes for largescale application problems [9, 10] to the traditional fourth order discretization with . The mass matrix for the traditional discretization becomes more and more poorly conditioned as the size of increases as illustrated in Figure 9. On the other hand, for instance, the condition number of with stays approximately close to almost independent of the size of and provides comparable accuracies as also illustrated in Figure 10 to that of the traditional discretization.
7. Conclusion
In this article, we have derived and demonstrated a family of sixthorder compact finitedifference schemes for Poisson's equation in three spatial dimensions. Such schemes can be extended in order to derive sixthorder 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.
Acknowledgment
This Research was partially supported by NOAA Grant NA06OAR4810187.
References
 J. C. Strikwerda, Finite Difference Schemes and Partial Differential Equations, Wadsworth & BrooksCole Advanced Books & Software, Pacific Grove, Calif, USA, 1989. View at: MathSciNet
 L. Ge and J. Zhang, “Symbolic computation of high order compact difference schemes for three dimensional linear elliptic partial differential equations with variable coefficients,” Journal of Computational and Applied Mathematics, vol. 143, no. 1, pp. 9–27, 2002. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 M. Piller and E. Stalio, “Finitevolume compact schemes on staggered grids,” Journal of Computational Physics, vol. 197, no. 1, pp. 299–340, 2004. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 W. F. Spotz and G. F. Carey, “Highorder compact scheme for the steady streamfunction vorticity equations,” International Journal for Numerical Methods in Engineering, vol. 38, no. 20, pp. 3497–3512, 1995. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 H. Sun, N. Kang, J. Zhang, and E. S. Carlson, “A fourthorder compact difference scheme on face centered cubic grids with multigrid method for solving 2D convection diffusion equation,” Mathematics and Computers in Simulation, vol. 63, no. 6, pp. 651–661, 2003. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 S. K. Lele, “Compact finite difference schemes with spectrallike resolution,” Journal of Computational Physics, vol. 103, no. 1, pp. 16–42, 1992. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 K. Ito, Z. Li, and Y. Kyei, “Higherorder, Cartesian grid based finite difference schemes for elliptic equations on irregular domains,” SIAM Journal on Scientific Computing, vol. 27, no. 1, pp. 346–367, 2005. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 Z. Li and K. Ito, “Maximum principle preserving schemes for interface problems with discontinuous coefficients,” SIAM Journal on Scientific Computing, vol. 23, no. 1, pp. 339–361, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
 H. C. Elman and D. P. O'Leary, “Efficient iterative solution of the threedimensional Helmholtz equation,” Journal of Computational Physics, vol. 142, no. 1, pp. 163–181, 1998. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 Y. A. Erlangga, C. Vuik, and C. W. Oosterlee, “On a class of preconditioners for solving the Helmholtz equation,” Applied Numerical Mathematics, vol. 50, no. 34, pp. 409–425, 2004. View at: Publisher Site  Google Scholar  Zentralblatt MATH  MathSciNet
Copyright
Copyright © 2010 Yaw Kyei et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.