A least square extrapolation method for improving solution accuracy of PDE computations

Richardson extrapolation (RE) is based on a very simple and elegant mathematical idea that has been successful in several areas of numerical analysis such as quadrature or time integration of ODEs. In theory, RE can be used also on PDE approximations when the convergence order of a discrete solution is clearly known. But in practice, the order of a numerical method often depends on space location and is not accurately satisfied on different levels of grids used in the extrapolation formula. We propose in this paper a more robust and numerically efficient method based on the idea of finding automatically the order of a method as the solution of a least square minimization problem on the residual. We introduce a two-level and three-level least square extrapolation method that works on nonmatching embedded grid solutions via spline interpolation. Our least square extrapolation method is a post-processing of data produced by existing PDE codes, that is easy to implement and can be a better tool than RE for code verification. It can be also used to make a cascade of computation more numerically efficient. We can establish a consistent linear combination of coarser grid solutions to produce a better approximation of the PDE solution at a much lower cost than direct computation on a finer grid. To illustrate the performance of the method, examples including two-dimensional turning point problem with sharp transition layer and the Navier-Stokes flow inside a lid-driven cavity are adopted.