High Order of Accuracy for Poisson Equation Obtained by Grouping of Repeated Richardson Extrapolation with Fourth Order Schemes

Luciano Pereira da Silva,Bruno Benato Rutyna,Aline Roberta Santos Righi and Marcio Augusto Villela Pinto

1Graduate Program in Numerical Methods in Engineering,Federal University of Paran´a,Curitiba,81531-980,Brazil

2Graduate Program in Space Sciences and Technologies,Aeronautics Institute of Technology,S˜ao Jos´e dos Campos,12228-900,Brazil

3Department of Mechanical Engineering,Federal University of Paran´a,Curitiba,81531-980,Brazil

ABSTRACT In this article,we improve the order of precision of the two-dimensional Poisson equation by combining extrapolation techniques with high order schemes.The high order solutions obtained traditionally generate non-sparse matrices and the calculation time is very high.We can obtain sparse matrices by applying compact schemes.In this article,we compare compact and exponential finite difference schemes of fourth order.The numerical solutions are calculated in quadruple precision (Real ⋆16 or extended precision)in FORTRAN language,and iteratively obtained until reaching the round-offerror magnitude around 1.0E−32.This procedure is performed to ensure that there is no iteration error.The Repeated Richardson Extrapolation (RRE)method combines numerical solutions in different grids,determining higher orders of accuracy.The main contribution of this work is based on a process that initializes with fourth order solutions combining with RRE in order to find solutions of sixth,eighth,and tenth order of precision.The multigrid Full Approximation Scheme(FAS)is also applied to accelerate the convergence and obtain the numerical solutions on the fine grids.

KEYWORDS Tenth order accuracy; RRE; compact scheme; exponential scheme; multigrid; finite difference

1 Introduction

One of the main goals of researchers is to obtain improved numerical methods capable of solving problems that hinder technological progress,as well as those related to Engineering.In particular,one of these problems is characterized by the determination of the numerical solution of the Navier-Stokes (NS)equations,which require a large computational effort when determining the numerical solution of the continuity equation,represented by the Poisson equation.The simplified two-dimensional Poisson equation is given by

along with boundary conditions relating to theΩdomain of the problem.Some studies,including Poisson equation,present solutions using compact fourth and sixth order methods with direct methods to solve the system of equations [1].In addition,they have found results including higher order of accuracy on non-uniform grids to determine critical grids that violate [2] stability conditions.In [3] the Eq.(1)was solved numerically by combining the standard multigrid method and Repeated Richardson Extrapolations (RRE),where solutions of sixth order of accuracy were obtained.In the work of [4] an unconventional compact nine points form was introduced,named as Exponential Difference Scheme based on the Finite Difference Method (FDM).

This work includes higher order methods for discretizing Eq.(1)based on FDM [5].These methods belong to a group of higher order schemes,where the combination of FDM approximations results in new methods that consider a larger number of neighboring points.They are:the Compact Finite Difference Scheme of fourth order (Compact−4)[3]; and the Exponential Finite Difference Scheme of fourth order (Exponential−4)[4].In addition,the higher order methods in conjunction with the standard multigrid Full Approximation Scheme (FAS)[6] provide faster convergence to the numerical solution.One can also determine significantly better results when applying a fourth order of accuracy scheme with nine points,due to the decreased discretization error.To determine numerical solutions with higher order of accuracy,the post-processing method known as RRE can be used [7–9].

The main goal of this work is to obtain successive accurate solutions for Eq.(1)in an efficient way without increasing the computational time.For this,we use RRE that combines nine numerical solutions using quadruple precision (Real ∗16 or extended precision)in order to obtain tenth order precision solutions.The grids have a refinement ratioq=2 from 9×9 to 2049×2049.The numerical solutions included in the present work were obtained using FORTRAN language.

This work is subdivided as follows:In Section 2 we present the proposed mathematical and numerical model; In Section 3 the applied fourth order methods; In Section 4 we describe the RRE method used to obtain the tenth order of accuracy of the numerical solution; In Section 5 we show the standard multigrid method; In Section 6 the numerical tests and results,and,in Section 7 the conclusions of the work.

2 Mathematical and Numerical Models

The main goal is to solve the two-dimensional Poisson equation,represented by (1),and to achieve tenth order of accuracy numerical solutions.For this,the first step is to find the best fourth order scheme applying the standad FAS,and from that apply the RRE.Therefore,there is no problem in taking a simple domain with a discretization in simple grids,because the standard multigrid could be applied in non-orthogonal grids or unstructured grids [7].As the RRE is a post-processing [10],its application is independent of how data were collected.

This section describes the model of the equation,its domain,boundary conditions and analytical solution that will be used by means of comparison and verification.

Given the equation

with boundary conditions obtained by the analytical solution,the unit square domain was discretized applying uniform grids withh:=1/(Nx−1)andk:=1/(Ny−1),forxandydirections,respectively.Here,NxandNyrepresent the number of grid points in coordinate directionsxandy,respectively.Therefore,these grid points are defined by(xi,yj),wherexi:=ihandyj:=jk,for 0 ≤i≤Nx−1 and 0 ≤j≤Ny−1.

Considering an uniform grid by direction with elements sizehandkfor each of the coordinate directionsxandy,respectively.The approximations for the partial derivarives are given by Central Difference Scheme of second order (CDS−2)[5].

and

Replacing Eqs.(3)and (4)in Eq.(1)the equation is obtained from which the estimate Eq.(5)is derived

The approximations given by Eqs.(3)and (4)result in the asymptotic orderpL=2 and true orderspm=2,4,6,...in Eq.(5)[11].We describe how to increase the order of this discretization method for two known high order methods in the next section.

3 High Order Methods for Poisson Equation

3.1 Compact Finite Difference Scheme

Given Eq.(1),whereUis the solution of Eq.(5)through a second order approximation.Therefore,in order to determine approximations for fourth order derivatives and obtain the known compact fourth order method,it is necessary to expand the following terms:

and

Then,replacing Eqs.(6)and (7)in Eqs.(3)and (4),respectively.Furthermore,considering aCDS−2 approximation for the second order partial derivative off(x,y),we have

and

For∂4u/(∂x2∂y2)|i,j,the approximation is

Next,Eq.(10)is replaced in Eqs.(8)and (9).After that,summing up both resulting equations leads to

Later,consideringk=hand rearranging Eq.(11),we have

where

According to [12],the local truncation error is given by

Finally,by Eq.(14)it is possible to stablish the asymptotic,pL=4,and true orders,pm=4,6,8,...,respectively.

After some algebraic manipulations,the stencil model is

Therefore,Eq.(15)is the stencil of the approximation byCompact−4 for the two-dimensional Poisson equation in an unit square domain applying uniform grid elements based on approximations of [3,13,14].

3.2 Exponential Finite Difference Scheme

According to [4,15,16],the approximationExponential−4 is accomplished from the following formulation:

Then,assuming thatk=h,by substituting the approximations from Eq.(11)on the left of the equality in Eq.(16),we get

where

Furthermore,the laplacian of Eq.(18)can be approximated,according to [5],byCDS−2.Thus,it results in

then,

According to [4],the constants area0=4,a1=1,a2=−20 andb0=6.Therefore,Eq.(17)can be rewrite as

The local truncation error is given by [4] as

Consequently,also according to [4],combining Eqs.(21)and (22)specifies the asymptotic,pL=4,and true orders,pm=4,6,8,...,respectively.

The stencil of the approximation for the Poisson equation is given by

Therefore,Eqs.(23)and (24)stablish the approximation forExponential−4 for the twodimensional Poisson equation with uniform grid elements based on the approximations of [4].

Both for Eqs.(12)and (13)of Compact Scheme (with stencil given by Eqs.(15)and (21)of Exponential Scheme (with stencils given by Eqs.(23)and (24)),a system of the typeAU=f is generated.This is a linear system with a sparse and large matrix,where iterative methods are the most suitable [5].

In the present work,we chose the nine points Gauss-Seidel solver (method for solving the linear system)in a lexicographic order for both discretizations in order to generate numerical solutions of such linear equations systems,whose matrices of coefficients are 9-diagonal diagonal given by algebraic manipulation of Eqs.(12)and (13)for Compact Scheme and algebraic manipulation of Eq.(21)for Exponential Scheme.

This Gauss-Seidel will be used as a smoother in the geometric multigrid method to determine the numerical solutions of the linear systems (see Section 4).

4 Repeated Richardson Extrapolation

For numerical analysis,Richardson Extrapolation (RE)[5] is an efficient technique of accelerating terms of a sequence that can be applied to reduce the discretization error of nodal numerical solutions.According to the Eq.(25),when combining two nodal solutions whith a refinement ratioqand order of accuracypof the numerical method it is possible to determineφ∞,which is an estimate for the analytical solution.

Then,we have the equation

whereφ2andφ1are solutions at the fine and coarse grid,respectively.

By applying RE repeatedly,we obtain the technique known as Repeated Richardson Extrapolation (RRE)[7,8,10,17].The RRE is a post-processing technique used to decrease the discretization error of numerical solutions without increasing the number of terms in the finite difference equation obtained by the Taylor series,it depends only on the collected data.The RRE combines different grid solutions and for each extrapolation,the obtained solution presents a much smaller discretization error.This magnitude obeys the true orders in the finite difference equation used to approximate the derivatives of the differential equation.For this case,models of centered fourth orders were used based onCDS−2 approximations,thus,the true orders arepm={4,6,8,...}[17,18].The expression for RRE is given by

whereφis the numerical solution for the variable of interest (see Section 6),gis the grid where the numerical solution was obtained,mis the number of Repeated Richardson Extrapolations,pmare the true orders for the equation of the discretization error andq=hg−1/hgis the refinement ratio.In this work we useq=2.

It is possible to approximate the asymptotic order of accuracy of each numerical method by the effective and apparent orders [11],respectively.The basis for the approximation is the numerical error (difference between analytical and numerical solutions)or by estimating the error of the numerical solution.The effective order (based on the numerical error)is given by [18].

whereΦis the analytical solution andφis the numerical solution of the variable of interest.In addition,hg−1andhgrepresent the course and fine grids,respectively.

The apparent order (based on its estimation of the numerical solution error)is represented by [18]whereqis the refinement ratio,φis the numerical solution of the variable of interest andhg−2,hg−1andhgrepresents the coarse grid,intermediate and fine,respectively.

To execute the RRE,pm=pEorpm=pUcan be used in Eq.(26),which are the orders of accuracy of the numerical scheme,calculated according to Eqs.(27)and (28),respectively.The accuracy orderspEandpUare equivalent for the RRE.However,in order to usepE,the analytical solution must be known.In addition,to calculatepU,three nodal solutions are needed,while to calculatepEonly two nodal solutions are sufficient.To show how the RRE is calculated,pm=pUis chosen,because in this case it is not necessary to know the analytical solution of the differential equation,and this makes the technique more robust.

The Algorithm 1 describes the steps to implement the RRE usingpU,but can be easily adapted to usepE.In Fig.1,the representation of one nodal numerical solution for the same point is shown in different grids and levels of extrapolation,whereφg=G=9,m=4is the solution of high order of accuracy associated with theg=9 grid (note that there are only nine grids in the representation,sog=G=9)after four Richardson extrapolations (m=4)that are featured in the Fig.1.It is important to note that at levelm=0 there is no extrapolation in the numerical solutions.The solutionφ9,4represents the output of Algorithm 1 forφ1...9,0,G=9 andq=2 as input data.

Algorithm 1:Repeated Richardson Extrapolation (RRE)Input:solutions φ1...G,0,number of grids (G)and refinement ratio (q)Output φ with high order 1 Inicialization: φ1,0=φ1,φ2,0=φ2,φ3,0=φ3,...,φG,0=φG 2 for g=3,...,G do 3for m=0,...,[integer part(g −1)/2] do 4if m<=[integer part(g −3)/2] then 5(pU)g,m=log(|φg−1,m−φg−2,m||φg,m−φg−1,m|)log(q)6end 7if m>=1 then 8φg,m=φg,m−1+φg,m−1 −φg−1,m−1 q(pU)g,m−1 −1 9end 10end 11 end 12 end

Figure 1:Schematic representation of RRE for g=9 and m=4

5 Multigrid Method(MG)

In general,the decrease in the discretization error is associated with the grid size,in other words,as much as the grid is refined,lower is the discretization error.However,highly refined grids result in large scale linear equation systems,which requires very high computational effort to achieve a numerical solution.Given that,in order to minimize the computational time used with highly refined grids,we adopted the multigrid method,which accelerates the convergence of standard methods using coarser grids as a support base.

Some classical iterative methods can reduce very quickly oscillatory modes errors of the finest grid leaving only smooth modes,where the method loses the effectiveness.Later,one should transfer information to the coarser grids (using restriction operators),where the smooth modes become more oscillatory,and then,the method becomes more effective.When reaching the possible or desired coarser grid,one should return with the information to the original finest grid of the problem (using prolongation operators).It speeds up the convergence of the method [6,19,20].

Depending on the type of information that is transferred between grids,one can use the Correction Scheme (CS)where only the residual is transferred,or the Full Approximation Scheme(FAS),where both residual and solution are transferred.In the present work,we chose the FAS method [6,19].A cycle is the way that different grids are traversed.There are different types of cycles,for example,V-,W-,F-cycle [21,22] and we choose to appply the V-cycle because its lower computational coast [6,19].The V-cycle is defined asV(ν1,ν2),whereν1andν2are the numbers of pre and post-smoothing the multigrid cycle,in other words,the number of smoothing in restriction and prolongation process,respectively.Coarsening ratio (re)is the ratio of the distance between the consecutive points of the grid from de fine grid (h)and the coarse grid (H).In this work,we usere=2,i.e.,H=2h[6,19].

In the present work,we adopt the standard geometric multigrid [23,24],that is usually recommended to orthogonal structured and non-orthogonal structured grids.However,for unstructured grids,the algebric multigrid can be applied [7,25,26].

The process of restriction (transfer information from the fine grid to the coarser grid)is carried out by the full weighting operator,defined by [6,19]

The process of prolongation (transfer information from the coarse grid to the finer grid)is carried out by the bilinear interpolation operator,defined by [6,19]

6 Numerical Results

Approximations ofAU=f were calculated in FORTRAN using quadruple precision,whereAN×Nis the matrix of coefficients,UN×1is the vector of variables whose solution represents the temperature and fN×1=f(xi,yj)=fi,jis the source term in lexicographic order,given byf(x,y)=−cos(x+y)+cos(x−y).The boundary conditions obtained by the analytical solutionu(x,y)=cos(x)cos(y).In addition,Rit=f −AUitis the defect of every each iteration (it).

The variables of interests are:the infinity norm (l∞)of the error based on the temperatureUi,j,the apparent and effective orders of the error based on the temperature at the midpoint,U(1/2,1/2),also,its defect |R(1/2)|.By that,one can analyze the decreasing in the module of the defect in the central point (|R(1/2)|)vs.the number of V-cycles (V(2,1))by its dimensionless form represented by |Rit|/|R0|.To assure that the numerical solution admit only discretization errors,the absolute value ofrit≈1×10−32was calculated until it achieves its machine round-off error for the gridsN=9×9,17×17,33×33,...,4097×4097.The numerical solutions for problems are obtained using computational programs written by the authors in Fortran 95 and using the Microsoft®Visual Studio®2008 compiler v.9.0.21022.8 RTM with quadruple precision and are run on a computer with a 3.40 GHz Intel®Core™i7-6700 processor with 16 GB of RAM hosting 64-bit Windows®10.

The multigrid method FAS uses V-cycles,V(ν1,ν2)=V(2,1)with the maximum possible value of levels of grids.Operator of restrictions and prolongations are defined by Eqs.(29)and (30).

Our main contribution is,first,to find the best fourth order method,and then,combine it with the RRE.In this case,using the multigrid method in order to solve the problem in refined grids.Therefore,theCDS−2 is used only for purpose of comparison and verification of the computational code.

6.1 Verification of Truncation Error

Considering all the grids studied in this work,for keeping the same standard,we chose grids withN= 2049×2049 eN= 4097×4097 to represent all the ones adopted.It is possible to observe in Fig.2 (forN=2049×2049 eN=4097×4097)that iteration errors for theCDS−2,Compact−4 eExponential−4 methods decrease approximately one order of magnitude for everyV(2,1)-cycle.Moreover,for both fourth order methods,the numbers of cycles were lower than 30 to reach the round-off error,while the second order method needed a slightly more than 35V(2,1)-cycles.

Figure 2:Iteration error vs. number of V-cycle for (a)2049×2049 and (b)4097×4097

Figure 3:Effective order pE of accuracy based on l∞of error

Fig.3 represents the effective order of the infinity norm (l∞)based on the error of the variableU(xi,yj).The Fig.4 represents the effective order of the numerical solution error at the midpoint and Fig.5 represents the apparent order of the numerical solution error at the midpoint.

Figure 4:Effective order pE of accuracy in approximating U(1/2,1/2)

Figure 5:Apparent order pU of accuracy in approximating U(1/2,1/2)

One can observe in these figures that even forpEas forpU,also for any other variable of interest that was studied in this paper,the second order method and the fourth order methods converge to their asymptotic orders predicted in Sections 2 and 3.

6.2 Analysis of Convergence

In order to analyze the results when approximatingU(1/2,1/2),one can verify the amount of V-cycles needed to obtain solutions that respect the stop criteria of 1E−10.Tab.1 shows the comparison between the number of cyclesV(2,1)that each method requires to obtain solutions.

Table 1:Number of cycles V(2,1) to reach the tolerance 1E-10

As one can see in Tab.1,both theCompact−4 and theExponential−4 methods obtained exactly the same amount of V-cycles.However,a slightly less number of cycles than for theCDS−2,which indicates that fourth orders methods have a slightly better convergence factor.

Tab.2 shows another analysis only for fourth order methods,considering time computational needed to obtain the same numerical solutions.

Table 2:CPU time(s)for fourth order methods

As one can see by Tab.2,for both methods of fourth order hold very similar computational time.Therefore,it is not possible yet to determine which one is the more efficient.

Given that,Tab.3 shows the complexity order of theCompact−4 andExponential−4 methods using multigrid FAS through a geometric curve fitting [5],represented by equationtcpu(N)=cNp.It is known that,according to [6,19],theoretically thepfactor must be close to the unit.

Table 3:Coefficient for the geometric curve fitting

Finally,after applying the curve fitting it is possible to verify that theExponential−4 method obtained a more adequate complexity factor when obtaining the numerical solutions for the variable of interestU(1/2,1/2),which we defined before in this paper.Note that theExponential−4 method is not part of the classical compact schemes,so we believe that it is possible to extend its application to obtain numerical solutions for more complex equations,such as the Poisson equations with non-constant or anisotropic coefficients.

6.3 Effect of the Number of Extrapolations

The results present in Tab.3 show that the complexity factorpof theExponential−4 method is closer to unity,agreeing with the literature [6].Therefore,from now on we choose theExponential−4 method to determine the numerical solutions for the variableU(1/2,1/2)and generate new extrapolated solutions with sixth,eighth and tenth orders of accuracy.In order to do repeated extrapolations using the RRE,the analytical solution was not used,for instance,all the approximations were based on the apparent orderpU.For that,one can remember that the true orders of theExponential−4 method arepm=4,6,8,...,and asymptotic order [17,18] given bypL=4.

It is possible to observe through Fig.6 the behavior of the true and apparent orders in the numerical solutions:without performing any extrapolations (m=0,p0andUg,0).When generating the first (m=1,p1andUg,1),second (m=2,p2andUg,2),third (m=3,p3andUg,2)and fourth(m=4,p4andUg,4)Richardson extrapolations.

Figure 6:Repeated Richardson extrapolation for U(1/2,1/2)

Observing Fig.6,we concluded that numerical solutions are found with several precision orders.

The main contribution of this article is to show the possiblity to obtain a new higher accuracy solution without the require to generate a new approximation involving many neighbor points.For this,we start with high-order solution and then apply the post-processing RRE method.This would imply in a large-scale linear systems,so that a computational limitation would happen.

In the Fig.6 we can observe the influence of the round-off error on the discretization error,which degenerated the orders.Note that theUg,4solution is affected by this type of error.With this,we believe that only solutions of tenth orderUg,3can be considered,excluding the degenerated part.

7 Conclusions

In the present paper,our results were based on a comparison betweenCompact−4 andExponential−4 methods with the goal of finding the best efficiency of the fourth order methods when decreasing the discretization error.In addition,without increasing the computational time when solving linear systems of large scale resulting of the discretization of the two-dimensional Poisson equation.

The fourth order methods achieved similar results when considering the number ofV(2,1)cycles when reaching the convergence,which indicates that they have a very close convergence factor.

From the curve fitting complexity analysis,theExponential−4 method obtained a better factor (smaller and closer to the unit)when compared to theCompact−4 method.However,both factors are close to the theoretical expected in the FAS multigrid case.

Given that,using theExponential−4 method and combining numerical solutions for theU(1/2,1/2)variable in the different grids using the RRE,we obtained new solutions with sixth,eighth and tenth orders of accuracy,as expected by the true orders.

However,it is possible to obtain solutions with superior orders by successive applications of RRE.For that,it is only needed enough different solutions in gridsg−1 to generatemrepeated extrapolations,i.e.,m=g−1.

In practice,our contribution is related to the significant reduction of computational costs during the execution of Engineering projects.It is common that projects are based on numerical solutions of problems that model,for example,heat diffusion.Thus,the strategy presented in this work allows the study of more physical phenomena,ensuring project cost reduction and high reliability.

This work is related to the exclusive approach of the Laplacian operator.However,we believe that extending this technique to the advection-diffusion operator,for example,is not difficult.

Acknowledgement:We would like to acknowledge the graduate program of Numerical Methods in Engineering—PPGMNE from the Federal University of Paraná,Curitiba—PR,Brazil.This study was financed in party by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brasil (CAPES)—Finance Code 001.

Funding Statement:The authors received no specific funding for this study.

Conflicts of Interest:The authors declare that they have no conflicts of interest to report regarding the present study.