Advances in Numerical Analysis

Volume 2010, Article ID 352174, 17 pages

http://dx.doi.org/10.1155/2010/352174

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

Department of Mathematics, North Carolina A & T State University, Greensboro, NC 27411, USA

Received 24 October 2009; Accepted 17 March 2010

Academic Editor: Yin Nian He

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.

#### Abstract

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.

#### 1. 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 , and consider

In order to obtain the high computational efficiency and the performance of higher-order 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. Higher-order 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 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 parameter-based 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 are determined by the stencil weights associated with the Laplacian operator , and the entries of are determined by the stencil weight associated with the right-hand 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, three-dimensional 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 fourth-order scheme for Poisson's equation, that is,

We then derive the scheme HOC, which is defined as the sixth-order 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 sixth-order 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 fourth-order error term and therefore produce a sixth-order 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* 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 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 right-hand side function 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 sixth-order 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 HOC and HOC. 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, 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 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 from the terms approximating . 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, , , , 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 in the literature for Poisson's equation. The stencils for HOC4 are illustrated in Figures 4 and 5.

#### 3. 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: 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 sixth-order 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 sixth-order 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 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 . 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 sixth-order HOC schemes.

In order to derive the complete family of sixth-order 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 three-parameter 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 sixth-order 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 three-dimensional 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 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

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 fourth-order 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 .

Now, we analyze the stability of the proposed schemes by solving (6.1) with . We compare the traditional fourth-order scheme HOC4 with grid points for the mass matrix (, ) and the case with grid points for (, ) 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 . 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 sixth-order compact finite-difference 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.

#### Acknowledgment

This Research was partially supported by NOAA Grant NA06OAR4810187.

#### References

- J. C. Strikwerda,
*Finite Difference Schemes and Partial Differential Equations*, Wadsworth & Brooks-Cole 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 · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - M. Piller and E. Stalio, “Finite-volume compact schemes on staggered grids,”
*Journal of Computational Physics*, vol. 197, no. 1, pp. 299–340, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus - W. F. Spotz and G. F. Carey, “High-order compact scheme for the steady stream-function vorticity equations,”
*International Journal for Numerical Methods in Engineering*, vol. 38, no. 20, pp. 3497–3512, 1995. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - H. Sun, N. Kang, J. Zhang, and E. S. Carlson, “A fourth-order 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 · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - S. K. Lele, “Compact finite difference schemes with spectral-like resolution,”
*Journal of Computational Physics*, vol. 103, no. 1, pp. 16–42, 1992. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - K. Ito, Z. Li, and Y. Kyei, “Higher-order, 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 · View at Google Scholar · View at Zentralblatt MATH · View at 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 · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet - H. C. Elman and D. P. O'Leary, “Efficient iterative solution of the three-dimensional Helmholtz equation,”
*Journal of Computational Physics*, vol. 142, no. 1, pp. 163–181, 1998. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at Scopus - 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. 3-4, pp. 409–425, 2004. View at Publisher · View at Google Scholar · View at Zentralblatt MATH · View at MathSciNet