High order solutions to self adjoint elliptic equations by reduction to the constant coefficient case

We develop a solver for nonseparable, self adjoint elliptic equations with a variable coefficient. If the coefficient is the square of a harmonic function, a transformation of the dependent variable, results in a constant coefficient Poisson equation. A highly accurate, fast, Fourier-spectral algorithm can solve this equation. When the square root of the coefficient is not harmonic, we approximate it by a harmonic function. A small number of correction steps are then required to achieve high accuracy. The procedure is particularly efficient when the approximation error is small. For functions with large Gaussian curvature, harmonic functions do not give optimal approximations. We investigate sub (and super) harmonic functions which are themselves solutions to constant coefficient elliptic equations. These functions give much smaller approximation errors. For a given function the approximation error becomes smaller as the size of the domain decreases. A highly parallelizable hierarchical procedure allows a decomposition into small sub-domains where the solution is efficiently computed this step is followed by a hierarchical matching procedure to reconstruct the global solution. Numerical experiments illustrate the high accuracy of the approach even at very coarse resolutions.