Analysis on Sixth-Order Compact Approximations with Richardson Extrapolation for 2 D Poisson Equation

By using Richardson extrapolation and fourth-order compact finite difference scheme on different scale grids, a sixth-order solution is computed on the coarse grid. Other three techniques are applied to obtain a sixth-order solution on the fine grid, and thus give out three kinds of Richardson extrapolation-based sixth order compact computation methods. By carefully analyzing the truncation errors respectively on 2D Poisson equation, we compare the accuracy of these three sixth order methods theoretically. Numerical results for two test problems are discussed.


Introduction
High order and high efficiency numerical computation for partial differential equations is very important in many scientific and engineering modeling problems.Compared to low order (second order) methods, high order methods can achieve satisfactory errors on much coarser grids and thus greatly curtail the computational cost [1].Although the most extensively used methods to obtain high order accuracy are spectral methods [2] [3], these methods have some limitations.The most significant one is that the discretization of PDE by spectral methods leads to the solution of large systems of linear or non-linear equations involving full matrices [4].In contrast, finite difference (and finite element) methods lead to systems with sparse matrices that can be handled by efficient iterative methods to reduce the computation complexity and computation cost substantially [5].Therefore, taking into account the computational efficiency, we focus on developing numerical algorithms based on high order finite difference methods here.Among various high order finite difference methods, high order compact (HOC) finite difference schemes are extremely noticeable.Compared to using straightforward central differences to obtain high order accuracy on a larger stencil which results the increase of bandwidth of the coefficient matrix and rises to a problem at the points close to the boundaries, the HOC schemes only use the center and adjacent points (i.e., a 9-points stencil is used in HOC schemes in 2D) which avoid extra special treatments for those points close to the boundaries and further improve computational efficiency.In the past two decades, HOC schemes have been well studied.Various fourth order compact (FOC) finite difference schemes have been developed for Poisson equations, convection-diffusion equations, and Navier-Stokes equations [6] [7] [8] [9].
Recently, there has been growing interest in developing sixth order compact finite difference schemes.By using Taylor series expansion, Soptz and Carey [10] developed a compact scheme which can achieve sixth order accuracy only when the derivatives of the source term can be determined analytically.Sutmann used Padé approximation discussed by Lele [11] on the Taylor expansion for the discretized Laplace operator to develop sixth order compact schemes [12].Although the schemes need less grid points than the straightforward expansion approach, they are not fully compact since other grid points, since other grid points (besides the center and adjacent points) are involved.. Chu and Fan [13] proposed a three-point combined compact difference scheme, which can achieve sixth order accuracy for the inner grid points and fifth order accuracy for the boundary grid points.However, their scheme asks to compute the dependent variables and their derivatives together which results in a triple-tridiagonal system with high computational cost for solution.In addition, it has a stability problem that, for certain problems with a large meshsize, the computed solution may be oscillatory [14].Although a finer meshsize may avoid numerical oscillations, the use of fine mesh is contradictory to the motivation of using higher order compact schemes.There are other sixth order finite difference schemes generated similarly [15] [16], but all of them share common weak points such as: 1) derivatives of source term appeared in the right-hand side which require analytical forms or approximations for the derivatives with certain order accuracy; 2) not complete compact schemes which may cause problems near to the boundaries; 3) complicated resulting linear systems which increase the difficulty of choosing effective iterative solvers.
Another category of sixth order compact approximations is based on Richardson extrapolation which is a technique introduced by Lewis Fry Richardson in the early of 20th century [17].In numerical analysis, Richardson extrapolation is a sequence acceleration method used to improve the rate of convergence of a sequence and is also the basis of Romberg integration [18].Although assumptions of smoothness and monotone truncation error convergence in the mesh spacing are involved, using Richardson extrapolation to get higher order accurate solutions is more convenient than developing direct discretizations [19].We can avoid complicated stencils, wider bandwidth matrices, special considerations for near-boundary points, additional stability analysis, etc.In addition, highly efficient solvers for the resulting large sparse linear system can be easily applied in such sixth order methods.Further more, Richardson extrapolation does not require any knowledge of the underlying methodology except the order of accuracy, which guarantees the minimal intervention to the existing computational tools [20].Sun et al. [21] first proposed to use Richardson extrapolation on two fourth order solutions of different grained grids to compute sixth order solutions on the fine grid.Then, Wang [22] developed a multiscale multigrid (MSMG) method as a very efficient solver to compute sixth order solutions for the 2D Poisson equation with Richardson extrapolation.Since the extrapolated sixth order solution is only generated on the coarse grid, the key issue of using Richardson extrapolation in the sixth order solution computation is how to obtain sixth order solutions on the fine grid.The most direct way is to inject sixth order coarse grid solution into the fine grid, which allows a subset of fine grid points to get sixth order solutions.Then, other strategies are required to compute sixth order solutions for the rest of the fine grid points and finally make the entire fine grid solution reach sixth order accuracy.At present, there are three such kind of strategies worthy of analysis and discussion.They are: 1) operator based interpolation; 2) multiple coarse grid computation; 3) complete Richardson extrapolation.
In this paper, our goal is to analyze and compare these Richardson extrapolation-based sixth order methods in computational accuracy.Therefore, we first describe these three Richardson extrapolation-based sixth order methods in Section 2.Then, we give a detailed analysis on their truncation errors in Section 3.After that, we use numerical experiments to verify our theoretical analysis in Section 4. At last, the conclusion and comments are given in Section 5.

Sixth Order Compact Approximations with Richardson Extrapolation
Consider a 2D Poisson equation in the form of where Ω is a rectangular domain with suitable boundary conditions defined on ∂Ω .The solution u and the forcing term ( ) , f x y are assumed to be sufficiently smooth and have necessary continuous partial derivatives.Ω can be discretized with uniform meshsize h in the x and y coordinate directions.The mesh points are ( ) x y with i x ih = and j y jh = .
Assume we have pth order accurate approximate solutions is used to compute the sixth order accurate solution 2 , h i j u  on the coarse grid 2h Ω [22].

Richardson Extrapolation with Operator Based Interpolation
One method using Richardson extrapolation to compute a sixth order accurate solution on the fine grid is to use an operator based interpolation scheme, proposed by Wang and Zhang [22], which iteratively updates the fine grid points in a specific sequence until some convergence condition is satisfied.This process is an iterative refinement procedure and similar to basic iterative methods [5].The operator based interpolation for the 2D Poisson Equation (1) has the formula as [22] ( ) ( ) .

Richardson Extrapolation with Multiple Coarse Grid Computation
The second Richardson extrapolation sixth order computation is to use multiple coarse grids.For a 2D problem, the fine grid can be coarsened into four coarse grids, each of which is composed of one subset of fine grid points from: (even, even) fine grid points, (even, odd) fine grid points, (odd, even) fine grid points, and (odd, odd) fine grid points.The sixth order solution for (even, even) fine grid points comes from Richardson extrapolation.Dai et al. [23] proposed a direct sixth order computation for other three groups of fine grid points based on multiple coarse grids.First, tridiagonal systems are built on the X-odd grid view, which is constructed by (even, even) fine grid points (marked in red) and (odd, even) fine grid points (marked in black) as Figure 1, for computing the sixth order solution of (odd, even) fine grid points.Then, tridiagonal systems are built on the Y-odd grid view, which is constructed by (even, even) fine grid points (marked in red) and (even, odd) fine grid points (marked in blue) as Figure 2, for computing the sixth order solution of (even, odd) fine grid points.At last, for the (odd, odd) fine grid points (marked in green) surrounded by other fine grid points with sixth order solutions (marked in red) as Figure 3,   only one step of operator interpolation is needed to reach the sixth order solution.Computation details can be found in [23].

Completed Richardson Extrapolation
The third approach of using Richardson extrapolation for computing sixth order solutions is Completed Richardson extrapolation.Completed Richardson extrapolation was developed by Roache and Knupp [19] and once used to produce a fourth order accurate solution on the fine grid.It did not use the extrapolated fourth order solution but rather the correction between the second order solution and the fourth order solution in the interpolation process.
Similarly, by using the correction between the fourth order solution and the extrapolated sixth order solution rather than the improved solution itself, we could compute a sixth order solution on the entire fine grid.
Assume , i j u * be the exact solution at node ( ) , i j on the fine grid.With the definition of fourth order accurate solution, we have ( ) , where As are the coefficients of the leading error term.
For the (even, even) fine grid points, we could directly inject extrapolated coarse grid solution to get sixth order solution.The correction between the fourth order solution and the extrapolated sixth order solution, say , h i j c , can be used to approximate the leading error term Then, we could use the correction ( 6) to approximate corrections for other three subsets of fine grid points, and thus compute sixth order solutions for them.
For the (odd, odd) fine grid points, the rotated grid interpolation, as Figure 2(a), is used to approximate the coefficient A of the leading error term is used to obtain the formula for sixth order solution computation as , , , , , ; where ( ) For the (odd, even) and (even, odd) fine grid points, the standard grid interpolation, as Figure 2(b), is used to approximate the coefficient A of the leading error term , even; even, odd.
is used to obtain the formula for sixth order solution computation as , , , , odd, even; even, odd;

Truncation Error Analysis
In this section, we will give an analysis of truncation errors to compare the accuracy of three Richardson extrapolation-based sixth order methods described in Section 2. All of these methods need to use FOC schemes to get the fourth order solutions on fine and coarse grids.We first analyze the truncation error of the FOC schemes.For more general applications, we will derive the truncation error for the FOC scheme with unequal meshsizes [8].

Denote x ∆ and y
∆ to be the meshsizes in the x and y coordinate directions, respectively.The standard second order central difference operators are , By using Taylor series, we have ( ) From Equations ( 11) and ( 12) we can discretize Equation (1) at the grid point By taking two times partial derivatives of x and y on both sides of Equation ( 1), respectively, we have .
Using central difference operators and Taylor series in Equations ( 14) and ( 15) gives By continuously taking partial derivatives of x on both sides of Equation ( 14), we get Similarly, by continuously taking partial derivatives of y on both sides of Equation ( 15), we get .
Substituting Equations ( 18) and (19) in Equation ( 16) gives . 144 360 360 And, substituting Equations ( 20) and (21) in Equation ( 17 Then, we use Equations ( 22) and ( 23) to replace the  x i j y i j i j x i j y i j x i j x i j x i j y i j y i j y i j where , 2 where the coefficients are ( ) ( ) The fourth order truncation error of the FOC scheme ( 25) is ( ) And, the sixth order truncation error of the FOC scheme ( 25) is In a special case with x y h ∆ =∆ = , the FOC scheme has the form as ( ) . 2 The fourth order and sixth order truncation errors of the FOC scheme (28) are ( ) ( ) Now we can take a look at the truncation error after applying Richardson extrapolation.From the definition of the fourth order accurate solutions on the fine and coarse grids, we have For all Richardson extrapolation-based sixth order compact approximations, the truncation error of (even, even) fine grid points is Extrapo τ because the corresponding sixth order solution is injected from the extrapolated solution of the standard coarse grid.For other three subsets of fine grid points, three computational strategies (operator based interpolation, multiple coarse grid computation, and completed Richardson extrapolation) could be used to obtain sixth order solutions.In the following part, the truncation error for these methods are given respectively.

Truncation Error of Operator Based Interpolation
The operator based interpolation scheme is from the 9-point FOC scheme (28), which can be iteratively used to approach sixth order solutions for (odd, odd), (odd, even) and (even, odd) fine grid points.It generates a sixth order error in the form as ( ) ( ) The operator based interpolation Equation ( 4) can be written as .

i j i j op i j i j i j i j i j i j i j i j i j op
.

i j op i j i j i j i j i j i j i j i j i j op
When the iterative process of using operator based interpolation converges, the residual r tends to 0. Then Equation (37) becomes to From Equation (39), we get

Truncation Error of Multiple Coarse Grid Computation
After applying Richardson extrapolation to get sixth order solutions for (even, even) fine grid points, in the multiple coarse grid computation, X-odd grid view and Y-odd grid view are constructed to compute sixth order solutions for (odd, even) and (even, odd) fine grid points, respectively [23].The X-odd grid view composed by (even, even) and (odd, even) indexed fine grid points is a view of unequal meshsize grid with meshsizes h and 2h in the x and y coordinate directions, respectively.The Y-odd grid view composed by (even, even) and (even, odd) indexed fine grid points is a view of unequal meshsize grid with meshsizes 2h and h in the x and y coordinate directions, respectively.The sixth order computations on the X-odd grid view and the Y-odd grid view by using tridiagonal solvers lead to sixth order truncation errors -odd x τ and -odd y τ , respectively.By using the general fourth order truncation error Equation (26) and setting corresponding mesh aspect ratio λ , we have explicit forms of As for the computation of (odd, even) fine grid points on the X-odd grid view with the mesh aspect ratio -odd 1 2 x λ = , the coefficients in Equation ( 25) are set as Assume the truncation error of (odd, even) fine grid points to be mcg α , an equation upon the error of X-odd grid view is generated from Equation (25) with the above coefficients as As for the computation of (even, odd) fine grid points on the Y-odd grid view with the mesh aspect ratio -odd 2 y λ = , the coefficients in Equation ( 25) are set as Assume the truncation error of (even, odd) fine grid points to be mcg β , an equation upon the error of Y-odd grid view is generated from Equation (25) with the above coefficients as which gives The update for (odd, odd) fine grid points uses the operator based interpolation Equation ( 4) and updated sixth order solutions of (even, even), (odd, even) and (even, odd) fine grid points.Assume the truncation error of (odd, odd) fine grid points to be mcg γ , an equation upon the error of the fine grid points is generated by Equation ( ,

Truncation Error of Completed Richardson Extrapolation
The completed Richardson extrapolation uses two kinds of interpolation on the correction Equation ( 6) of (even, even) fine grid points to approximate corrections for other fine grid points.As we all know, interpolation brings error.
Next, we will analyze how much of the error.
The A coefficients in Equations ( 7) and ( 9) can be viewed as a function of u which has the form of ( )

FOC4
A u h τ = .Based on the Taylor series expansion, the ( ) O h term in Equation ( 7) has an explicit form as the ( ) ( ) ; while the ( ) O h term in Equation ( 9) has an explicit form as ( ) ( ) The second order truncation error for Equation (7), which uses rotated grid interpolation to approximate (odd, odd) fine grid points, has the form as ( ) ( ) while the second order truncation error for Equation (9), which uses standard grid interpolation to approximate (odd, even) and (even, odd) fine grid points, has the form as ( ) ( ) We find that RotateInter StandInter 2 τ τ = .
Consider (odd, odd) fine grid points at first.Equation ( 7) can be re-written as The sixth order computation for the (odd, odd) fine grid points is only related to (even, even) fine grid points.For the (even, even) fine grid points, by using definition of fourth order solutions obtained from the FOC scheme, we have 4 even,even even,even even,even FOC6 4
After injecting the extrapolated coarse grid solution into the fine grid, we have 6 even,even even,even Extrapo .
The sixth order computation for the (odd, even) and (even, odd) fine grid points are related to both (even, even) and (odd, odd) fine grid points.For the updated (odd, odd) fine grid points, we have .
We find that the truncation errors of (odd, odd), (odd, even) and (even, odd) fine grid points have the same form as ( ), which is larger than the truncation error of (even, even) fine grid points ( Extrapo τ ) generated from Richardson extrapolation as we expect.It is because another interpolation is involved, i.e., Equation (55) or Equation (60).In summary, these three Richardson extrapolation-based methods are able to compute the sixth order accurate solution on the entire fine grid.For (even, even) fine grid points, all methods use Richardson extrapolation to get the sixth order solution with truncation error Extrapo τ .For other three subsets of fine grid points, different computational strategies are used to obtain sixth order solutions, which add errors of different magnitude on the truncation error Extrapo τ . Table 1 lists the truncation errors of the three different Richardson extrapolation-based sixth order compact computations, respectively.Since the error expressions involve various high-order partial derivatives of u, it is hard to conclude a quantitative relationship.By comparing the coefficients of common items, we predict that the completed Richardson extrapolation method should be more accurate than the other two methods.We think the operator based interpolation method has comparable accuracy to the multiple coarse grid computation method.We cannot estimate which one is more accurate than the other because there exists uncertainty of high order partial differential terms in the explicit truncation errors and it is hard to determine the magnitude and sign of them.on 4h Ω was the zero vector.The V-cycle on 2h Ω and h Ω stopped when the L 2 -norm of the difference of the successive solutions was reduced by a factor of 10 10 .The stopping criteria for the iterative operator based interpolation and Gauss-Seidel procedure was 10 −10 .All reported errors were the maximum absolute errors over the discrete grid of the finest level.
The codes were written in Fortran 77 programming language and run on a PC, which has Intel Core i7-4770 with 3.40 GHz and 16 GB RAM.

Test Problems
Problem 1. ( where the boundary conditions are The parameters are chosen as The analytical solution is

Accuracy and Efficiency Comparison
We refined the grid from 32 N = to 256 N = for both test problems.For convenience, we used three abbreviations to represent three Richardson extrapolation-based sixth order methods to be compared."Op-Six" is short for the sixth order method with Richardson extrapolation and operator based interpolation; "MCG-Six" means the sixth order method with Richardson extrapolation and multiple coarse grid computation; "CR-Six" denotes the sixth order method with completed Richardson extrapolation.
For both test problems, we compare the accuracy and efficiency of three methods by computing maximum errors and accuracy order, and recording the CPU time in seconds.The results of Problem 1 are shown in Table 2 and Table 3 and Figure 3.The results of Problem 2 are shown in Table 4 and Table 5 and     with mixed boundary conditions.The finite difference ghost-cell technique [25] keeps unchanged the symmetry of the stencil and uniform meshsize through adding extra points (ghost points) outside the domain.Another thing we want to point out is Richardson extrapolation can also be used to improve the accuracy of solutions obtained from finite element methods [26], although we used finite difference methods to illustrate the idea in this paper.And, compared to high order spectral methods which ask for simple/regular geometry of the problem domain or special strategies (i.e., spectral elements, domain decomposition, and etc.) for complex/irregular domains, using finite element methods to provide basic solutions and Richardson extrapolation to improve them would be a more flexible way to handle complex geometries appearing in many engineering problems.

Figure 2 .
Figure 2. Illustration of the interpolation strategy in 2D.(a) Rotated grid interpolation scheme; (b) Standard grid interpolation scheme.

Figure 3 .
Figure 3.Comparison of the error and CPU for the Problem 1.
equation Au F = , there is a corresponding residual equationAe r= , where e is the error and r is the residual.According to Equation (36), we could get the operator based interpolation on error as

Figure 4 .
Figure 4. Comparison of the error and CPU for the Problem 2.

Table 2 .
Accuracy comparison results for Problem 1.

Table 3 .
Efficiency comparison results for Problem 1.

Table 4 .
Accuracy comparison results for Problem 2.

Table 5 .
[24]ciency comparison results for Problem 2. are generally consistent with the observation from the theoretical analysis.The completed Richardson extrapolation method usually performs better in accuracy than the operator based interpolation with Richardson extrapolation method and the multiple coarse grid computation with Richardson extrapolation method.As expected, the efficiency of the operator based interpolation with Richardson extrapolation is lower than the other two sixth order methods.The exploration of using Richardson extrapolation on sixth order compact computation for high dimensional problems has already begun.In fact, the operator based interpolation with Richardson extrapolation method and the multiple coarse grid computation with Richardson extrapolation method have been used to solve 3D convection-diffusion equations and show exciting performance[24].The error analysis to 3D cases could be conducted and reported in the future.Since Richardson extrapolation-based sixth order methods require two comparable uniform grids (i.e., fine grid with meshsize h and coarse grid with meshsize 2h), the application of this technique to more practical problems in irregular domain with various boundary conditions (i.e., Dirichlet, Neumann, or mixed conditions) has limitations and is not straight forward.One undergoing study is to use the proposed methodology with finite difference ghost-cell technique to obtain high accuracy high efficiency solutions for Poisson equation accuracy