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

Authors: Yaw Kyei; John Paul Roop; Guoqing Tang
Journal: Advances in Numerical Analysis
 (2010)

---

## 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.


---

## Body

## 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 ℝ3, and consider(1.1)-Δu=f,inΩ,u=g,on∂Ω.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 Figure1.    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:(1.2)Hu=Qf,
where the entries of H are determined by the stencil weights associated with the Laplacian operator (Δ:=∂/∂x2+∂/∂y2+∂/∂z2), 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 w0,w1,…,w26, 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 w0), the stencil weights in the directions of the coordinate axes (i.e., steps of size h are taken in only one direction)  (denoted w1=⋯=w6), the stencil weights where steps of size h are taken in two coordinate directions  (denoted w7=⋯=w18), and the stencil weights where steps of size h are taken in all three coordinate directions  (denoted w19=⋯=w26).Figure 1
Stencil numbering for three-dimensional Poisson problem.By this logic, three-dimensional compact stencil for Poisson's equation may be described by eight values,w0, w1, w7, and w19 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.Figure 2
Generic stencil weights for the left-hand side matrix,H.Figure 3
Generic stencil weights for the right-hand side matrixQ.In this article, we derive several schemes for (1.1), these schemes we will label HOC4, HOC61, HOC62, and HOC6. HOC4 is defined as the traditional fourth-order scheme for Poisson's equation, that is,(1.3)HOC4:{α0=4,α1=-13,α7=-16,α19=0,β0=12,β1=112,β7=0,β19=0.We then derive the scheme HOC61, 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(1.4)HOC61:{α0=4,α1=-13,α7=-16,α19=0,β0=1+24w72,β1=1-48w712,β7=w7,β19=0,
where w7 is determined by the formula,(1.5)w7(HOC61)=(1/144)ξ1-(1/240)ξ2ξ1+3ξ3,
where(1.6)ξ1:=∂6u∂x4∂y2+∂6u∂x4∂z2+∂6u∂x2∂y4+∂6u∂x2∂z4+∂6u∂y2∂z4+∂6u∂y4∂z2,ξ2:=∂6u∂x6+∂6u∂y6+∂6u∂z6,ξ3:=∂u6∂x2∂y2∂z2.
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 HOC62 in which the weights are rearranged from HOC61 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:(1.7)HOC62:{α0=6415,α1=-716,α7=-110,α19=-130,β0=1+24w72,β1=1-48w712,β7=w7,β19=0,
where w7 is given by the formula(1.8)w7(HOC62)=190-1240(fxxxx+fyyyy+fzzzzfxxyy+fxxzz+fyyzz).Finally, we give themost 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 HOC62 and utilizes all 27 stencil weights for each of the matrices H  and  Q. The scheme HOC6  is given as follows:(1.9)HOC6:{α0=6415,α1=-716,α7=-110,α19=-130,β0=1+24w7+32w192,β1=1-48w7-48w1912,β7=w7,β19=w19,
where w7  and  w19 satisfy the relationship:(1.10)w7+2w19=190-1240(fxxxx+fyyyy+fzzzzfxxyy+fxxzz+fyyzz).
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 Section2, 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 HOC61  and HOC62. 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, HOC61, HOC62, 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),(2.1)6u0-u1-u2-u3-u4-u5-u6h2≈f0.
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 (2.2)6u0-u1-u2-u3-u4-u5-u6h2-f0=h212[∂4u∂x4+∂4u∂y4+∂4u∂z4]+𝒪(h4).
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, uxxxx, uyyyy, uzzzz, 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,(2.3)-uxxxx-uxxyy-uxxzz=fxx,-uxxyy-uyyyy-uyyzz=fyy,-uxxzz-uyyzz-uzzzz=fzz,
we see that the second-order error term may be rewritten as(2.4)h26(uxxyy+uxxzz+uyyzz)-h212(fxx+fyy+fzz).
Next, using the facts that (2.5)fxx≈f3-2f0+f1h2,fyy≈f4-2f0+f2h2,fzz≈f5-2f0+f6h2,(2.6)uxxyy≈u7+u8+u9+u10-2u1-2u3+4u0h4,uxxzz≈u11+u13+u15+u17-2u5-2u6+4u0h4,uyyzz≈u12+u14+u16+u18-2u2-2u4+4u0h4,
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(2.7)24u0-2(u1+⋯+u6)-(u7+⋯+u18)6h2≈6f0+f1+⋯+f612.
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.Figure 4
Stencil weights for the left-hand side matrixH for the fourth-order scheme, HOC4.Figure 5
Stencil weights for the right-hand side matrixQ for the fourth-order scheme, HOC4.


## 3. A One-Parameter Family of Fourth Order Schemes  and the Sixth-Order Scheme HOC61

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 Figures6 and 7 respectively.Figure 6
Alternate stencil weights for the right-hand side matrixQ for a second fourth-order scheme.Figure 7
Stencil weights for the right-hand side matrixQ for the one-parameter family of fourth-order schemes.We may obtain the following expression for the error:(3.1)h4[(1144-w7)ξ1-1240ξ2-3w7ξ3]+O(h6),
where ξ1, ξ2, and  ξ3 are defined as in (1.6). Eliminating the fourth order error term in the above equation, we have (3.2)w7=(1/144)ξ1-(1/240)ξ2ξ1+3ξ3.


### 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 functionf. 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 (3.3)(1144-w7)ξ1-1240ξ2-3w7ξ3.First, making the substitutions,(3.4)-uxxxxxx-uxxxxyy-uxxxxzz=fxxxx,-uxxyyyy-uyyyyyy-uyyyyzz=fyyyy,-uxxzzzz-uyyzzzz-uzzzzzz=fzzzz,
we obtain the  following equivalent form for (3.3):(3.5)(190-w7)ξ1+1240(fxxxx+fyyyy+fzzzz)-3w7ξ3.
Next, making the substitutions (3.6)-uxxxxyy-uxxyyyy-uxxyyzz=fxxyy,-uxxxxzz-uxxyyzz-uxxzzzz=fxxzz,-uxxyyzz-uyyyyzz-uyyzzzz=fyyzz,
we obtain the following equivalent form for (3.3):(3.7)(w7-190)(fxxyy+fxxzz+fxxyy)+1240(fxxxx+fyyyy+fzzzz)-130ξ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(3.8)w7=190-1240(fxxxx+fyyyy+fzzzzfxxyy+fxxzz+fyyzz)
and using the following approximation:(3.9)ξ3≈-8u0+4(u1+⋯+u6)-2(u7+⋯+u18)+(u19+⋯+u26)h6.
Adding this approximation to the stencil for u, we obtain the stencil which appears in Figure 8.Figure 8
Stencil weights for the left-hand side matrixH for the sixth-order scheme, HOC61.


## 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 functionf. 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 (3.3)(1144-w7)ξ1-1240ξ2-3w7ξ3.First, making the substitutions,(3.4)-uxxxxxx-uxxxxyy-uxxxxzz=fxxxx,-uxxyyyy-uyyyyyy-uyyyyzz=fyyyy,-uxxzzzz-uyyzzzz-uzzzzzz=fzzzz,
we obtain the  following equivalent form for (3.3):(3.5)(190-w7)ξ1+1240(fxxxx+fyyyy+fzzzz)-3w7ξ3.
Next, making the substitutions (3.6)-uxxxxyy-uxxyyyy-uxxyyzz=fxxyy,-uxxxxzz-uxxyyzz-uxxzzzz=fxxzz,-uxxyyzz-uyyyyzz-uyyzzzz=fyyzz,
we obtain the following equivalent form for (3.3):(3.7)(w7-190)(fxxyy+fxxzz+fxxyy)+1240(fxxxx+fyyyy+fzzzz)-130ξ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(3.8)w7=190-1240(fxxxx+fyyyy+fzzzzfxxyy+fxxzz+fyyzz)
and using the following approximation:(3.9)ξ3≈-8u0+4(u1+⋯+u6)-2(u7+⋯+u18)+(u19+⋯+u26)h6.
Adding this approximation to the stencil for u, we obtain the stencil which appears in Figure 8.Figure 8
Stencil weights for the left-hand side matrixH for the sixth-order scheme, HOC61.


## 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 functionf. 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(4.1)H:α0=4+8α19,α1=-1-12α193,α7=-1+12α196,α19=free,Q:β0=1+24β7+32β192,β1=1-48β7-48β1912,β7=free,β19=free.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(4.2)h4[(1144-β7-2β19)ξ1-1240ξ2+(α19-3β7-6β19)ξ3]+O(h6).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(4.3)β7+2β19=190-1240(fxxxx+fyyyy+fzzzzfxxyy+fxxzz+fyyzz),
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 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 (fxxyy,fxxzz,fyyzz) may be approximated in the usual way; see, for example, (2.6), in order to obtain the formula(5.1)fxxyy+fxxzz+fyyzz≈12f0+f7-18-4f1-6h4.
To approximate the other fourth order partial derivative terms (fxxxx, fyyyy, fzzzz), we notice that these quantities cannot be approximated on the compact stencil. However, using the fact that, for example,(5.2)fxxxx≈f(x+2ĥ,y)-2f(x+ĥ,y)+6f(x,y)-2f(x-ĥ,y)+f(x-2ĥ,y)ĥ4,
and utilizing the value ĥ=h/2, we obtain the following approximation formula:(5.3)β7+2β19≈190-115(18f0-4f1/2+f1-612f0+f7-18-4f1-6),
where(5.4)f0:=f(x,y,z),f1-6:=f(x+h,y,z)+f(x-h,y,z)+f(x,y+h,z)+f(x,y-h,z)+f(x,y,z+h)+f(x,y,z-h),f7-18:=f(x+h,y+h,z)+f(x+h,y-h,z)+f(x-h,y+h,z)+f(x-h,y-h,z)+f(x+h,y,z+h)+f(x+h,y,z-h)+f(x-h,y,z+h)+f(x-h,y,z-h)+f(x,y+h,z+h)+f(x,y+h,z-h)+f(x,y-h,z+h)+f(x,y-h,z-h),f1/2:=f(x+h2,y,z)+f(x-h2,y,z)+f(x,y+h2,z)+f(x,y-h2,z)+f(x,y,z+h2)+f(x,y,z-h2).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 withh=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.


### 5.1. Computational Experiments

In this section, we illustrate the experimental convergence rates for each of our schemes HOC4, HOC61, HOC62, 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, HOC61 represents the minimum error obtained. However, recall from our derivation that weights for HOC61 are determined according to the solution u, whereas the weights for HOC62 and HOC6 are minimized a priori.Experiment 1.
LetΩ:=(0,1)×(0,1)×(0,1). For this example, we implement the schemes HOC4, HOC61, HOC62, and HOC6 for u(x,y,z)=x5y2z3 taken as the true solution of (1.1) and the schemes HOC4, HOC61, HOC62, and HOC6 defined as in (1.3), (1.4), (1.7), and (1.9), respectively. The weights w7 satisfy the approximate formula
(5.5)w7+2w19=1144(3y2+z23x2+3y2+z2).
Results are given for six cases, HOC4, HOC61, HOC62 (β19=0), HOC6 (β7=0), HOC6 (β7=β19), and HOC6 (β7=2β19). Table 1 summarizes our results, illustrating the experimental convergence rates for each of the schemes.Table 1
Experimental error results for Experiment1.
h∥u-uh∥L∞  HOC4cvge.  rate∥u-uh∥L∞  HOC61cvge.  rate1/41.785562·10-43.492412·10-61/81.104208·10-54.025.703779·10-85.941/167.016054·10-73.989.095885·10-105.971/324.439135·10-83.981.421643·10-116.001/642.775527·10-94.002.228478·10-136.00HOC62 (β19=0.0)rateHOC6 (β7=0.0)rate1/47.048372·10-61.884849·10-51/81.194951·10-75.883.043228·10-75.951/161.881611·10-95.994.914854·10-95.951/322.949327·10-116.007.711285·10-115.991/644.622188·10-136.001.203933·10-126.00HOC6 (β7=β19)rateHOC6 (β7=2β19)rate1/41.021620·10-55.900059·10-61/81.682916·10-75.921.002760·10-75.881/162.677366·10-95.971.572341·10-96.001/324.202498·10-115.992.475537·10-115.991/646.567039·10-136.003.858164·10-136.00Experiment 2.
LetΩ:=(0,1)×(0,1)×(0,1). For this example, we implement the schemes HOC4, HOC61, HOC62, and HOC6 for u(x,y,z)=z5sin(xy) taken as the true solution of (1.1) and the schemes HOC4, HOC61, HOC62, and HOC6 defined as in (1.3), (1.4), (1.7), and (1.9), respectively. The weights w7 satisfy the approximate formula
(5.6)w7+2w19=-χ1720χ2,
where
(5.7)χ1:=(40y3xz4+40x3yz4)cos(xy)+(3y6z4-5y4x2z4+60y2z4+100y4z2-5x4y2z4+60x2z4+3x6z4+100x4z2-600y2-600x2)sin(xy),χ2:=(240yxz2-8y3xz4-8x3yz4)cos(xy)+(-20y4z2-60y2x2z2+120z2+120y2+y4x2z4-12y2z4+x4y2z4-12x2z4-20x4z2+120x2)sin(xy).
Results are given for six cases, HOC4, HOC61, HOC62 (β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)).Table 2
Experimental error results for Experiment2.
∥u-uh∥L∞  HOC4cvge.  rate∥u-uh∥L∞   HOC61cvge. rate1/42.750218·10-57.190362·10-71/81.722579·10-64.001.129125·10-85.991/161.100952·10-73.971.807232·10-105.971/326.893032·10-94.002.830760·10-126.001/644.310131·10-104.004.354850·10-146.02HOC62 (β19=0.0)rateHOC6 (β7=0.0)rate1/41.299453·10-63.569410·10-61/82.354981·10-85.796.288367·10-85.831/163.679920·10-106.009.817293·10-106.001/325.771845·10-126.001.549958·10-115.991/648.967826·10-146.012.431874·10-135.99HOC6 (β7=β19)rateHOC6 (β7=2β19)rate1/41.999183·10-61.214070·10-61/83.407251·10-85.871.966693·10-85.951/165.371719·10-105.993.181549·10-105.951/328.469725·10-125.994.971176·10-126.001/641.334904·10-135.997.897849·10-145.98


## 5.1. Computational Experiments

In this section, we illustrate the experimental convergence rates for each of our schemes HOC4, HOC61, HOC62, 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, HOC61 represents the minimum error obtained. However, recall from our derivation that weights for HOC61 are determined according to the solution u, whereas the weights for HOC62 and HOC6 are minimized a priori.Experiment 1.
LetΩ:=(0,1)×(0,1)×(0,1). For this example, we implement the schemes HOC4, HOC61, HOC62, and HOC6 for u(x,y,z)=x5y2z3 taken as the true solution of (1.1) and the schemes HOC4, HOC61, HOC62, and HOC6 defined as in (1.3), (1.4), (1.7), and (1.9), respectively. The weights w7 satisfy the approximate formula
(5.5)w7+2w19=1144(3y2+z23x2+3y2+z2).
Results are given for six cases, HOC4, HOC61, HOC62 (β19=0), HOC6 (β7=0), HOC6 (β7=β19), and HOC6 (β7=2β19). Table 1 summarizes our results, illustrating the experimental convergence rates for each of the schemes.Table 1
Experimental error results for Experiment1.
h∥u-uh∥L∞  HOC4cvge.  rate∥u-uh∥L∞  HOC61cvge.  rate1/41.785562·10-43.492412·10-61/81.104208·10-54.025.703779·10-85.941/167.016054·10-73.989.095885·10-105.971/324.439135·10-83.981.421643·10-116.001/642.775527·10-94.002.228478·10-136.00HOC62 (β19=0.0)rateHOC6 (β7=0.0)rate1/47.048372·10-61.884849·10-51/81.194951·10-75.883.043228·10-75.951/161.881611·10-95.994.914854·10-95.951/322.949327·10-116.007.711285·10-115.991/644.622188·10-136.001.203933·10-126.00HOC6 (β7=β19)rateHOC6 (β7=2β19)rate1/41.021620·10-55.900059·10-61/81.682916·10-75.921.002760·10-75.881/162.677366·10-95.971.572341·10-96.001/324.202498·10-115.992.475537·10-115.991/646.567039·10-136.003.858164·10-136.00Experiment 2.
LetΩ:=(0,1)×(0,1)×(0,1). For this example, we implement the schemes HOC4, HOC61, HOC62, and HOC6 for u(x,y,z)=z5sin(xy) taken as the true solution of (1.1) and the schemes HOC4, HOC61, HOC62, and HOC6 defined as in (1.3), (1.4), (1.7), and (1.9), respectively. The weights w7 satisfy the approximate formula
(5.6)w7+2w19=-χ1720χ2,
where
(5.7)χ1:=(40y3xz4+40x3yz4)cos(xy)+(3y6z4-5y4x2z4+60y2z4+100y4z2-5x4y2z4+60x2z4+3x6z4+100x4z2-600y2-600x2)sin(xy),χ2:=(240yxz2-8y3xz4-8x3yz4)cos(xy)+(-20y4z2-60y2x2z2+120z2+120y2+y4x2z4-12y2z4+x4y2z4-12x2z4-20x4z2+120x2)sin(xy).
Results are given for six cases, HOC4, HOC61, HOC62 (β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)).Table 2
Experimental error results for Experiment2.
∥u-uh∥L∞  HOC4cvge.  rate∥u-uh∥L∞   HOC61cvge. rate1/42.750218·10-57.190362·10-71/81.722579·10-64.001.129125·10-85.991/161.100952·10-73.971.807232·10-105.971/326.893032·10-94.002.830760·10-126.001/644.310131·10-104.004.354850·10-146.02HOC62 (β19=0.0)rateHOC6 (β7=0.0)rate1/41.299453·10-63.569410·10-61/82.354981·10-85.796.288367·10-85.831/163.679920·10-106.009.817293·10-106.001/325.771845·10-126.001.549958·10-115.991/648.967826·10-146.012.431874·10-135.99HOC6 (β7=β19)rateHOC6 (β7=2β19)rate1/41.999183·10-61.214070·10-61/83.407251·10-85.871.966693·10-85.951/165.371719·10-105.993.181549·10-105.951/328.469725·10-125.994.971176·10-126.001/641.334904·10-135.997.897849·10-145.98


## 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(6.1)-Δu=λuinΩ,u=0,on∂Ω.Proceeding as before, we approximate the eigenvalue problem (6.1) by(6.2)Hu=λQu,
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 theH matrix. We note that the restriction on α19 for row diagonal dominance is(6.3)0≤α19≤112.
The restrictions on β7 and β19 for row diagonal dominance of the mass matrix Q are(6.4)|1+24β7+32β192|≥6|1-48β7-48β1912|+12|β7|+8|β19|
with feasible positivity solutions as(6.5)0≤β7≤124,0≤β19≤116,β7+β19≤148,6β7+4β19≤14
or(6.6)0≤β7≤148,0≤β19≤148,β7+β19≤148.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.Condition number ofQ depending on β7 and β19. (a) Grid spacing h=1/10; (b) Grid spacing h=1/20.
(a)(b)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.Computational accuracies for the first25 eigenvalues of (6.1) for different choices of β7  and   β19 for Q. (a) Grid spacing h=1/10; (b) Grid spacing h=1/20.
(a)(b)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.


## 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.



---

*Source: 352174-2010-05-26.xml*