Higher Order Differences

Most textbook discussions of relaxation methods for solving partial differential equations use the familiar second order form for the Laplacian. Finite difference Laplacian representations result from Taylor series (i.e. polynomial) expansions of a function centered on the grid point x. We give the prescription for the second order formula and then it is apparent how to proceed to higher orders. Expand the function in the positive and negative directions to fourth order:

$$f(x + h) = f(x) + f'(x)h + \frac{1}{2} f''(x)h^2 + \frac{1}{6} f'''(x)h^3 + O(h^4)$$ (24)

and

$$f(x - h) = f(x) - f'(x)h + \frac{1}{2} f''(x)h^2 - \frac{1}{6} f'''(x)h^3 + O(h^4)$$ (25)

where f(x-h),f(x), and f(x+h) are the known function values on the grid points. The variable that we want to estimate is f''(x). Add the two equations, then solve for f''(x):

$$f''(x) = \frac{1}{h^2} [f(x - h) -2f(x) +f(x + h)] + O(h^2)$$ (26)

The same technique can be used to generate the 4th order formula for the Laplacian by including extra terms on either side of the 3pt formula giving a 5pt formula:

$$f''(x) = \frac{1}{12h^2} [-f(x-2h)+16f(x-h)-30f(x)+16f(x+h)-f(x+2h)] +O(h^4)$$ (27)

In two dimensions, the 4th order formula just contains the two `strips' along the x and y axes ($\nabla^2$ is a direct product), i.e. a total of 9 terms. In three dimensions, the total number of terms in the Laplacian at a given point are 7, 13, 19, etc. for 2nd, 4th, 6th ... order formulas. So for a minimal increase in computational cost, significantly more accurate formulas are obtained. However, there is a tradeoff in that the higher order forms can `span' a large portion of the grid for the coarser scales, and problems may be encountered near boundaries since so many terms must be fixed there (or else the function must be interpolated significantly out beyond the fixed boundary). This last factor must be considered in developing viable parallel algorithms. The best strategy thus entails a trade-off between higher order accuracy and computational feasibility. High order multigrid algorithms are a relatively new field, and much remains to be worked out in designing the optimal solver for a given problem.

An elegant technique which is a general method for generating high order formulas of many sorts has been developed by Hamming [26]. There the coefficients are obtained by computing the `moments' for the operation of interest (differentiation, integration, etc.). That is, let the operator of interest act on x<sup>n</sup> ; n=0,1,2,... with n going up to the number of grid points to be used in the formula. This will generate a vector to be multiplied by a `universal' matrix which is independent of the particular operation. These matrices are tabulated up to high orders in Hamming's book. The matrix times vector multiplication yields directly the coefficients in the expansion.