Fast potential theory. II. layer potentials and discrete sums

We present three new families of fast algorithms for classical potential theory, based on Ewald summation and fast transforms of Gaussians and Fourier series. Ewald summation separates the Green function for a cube into a high-frequency localized part and a rapidly-converging Fourier series. Each part can then be evaluated efficiently with appropriate fast transform algorithms. Our algorithms are naturally suited to the use of Green functions with boundary conditions imposed on the boundary of a cube, rather than free-space Green functions. Our first algorithm evaluates classical layer potentials on the boundary of a d-dimensional domain, with d equal to two or three. The quadrature error is O(hm) +ε, where h is the mesh size on the boundary and m is the order of quadrature used. The algorithm evaluates the discretized potential using N elements at O(N) points in O(N log N) arithmetic operations. The constant in O(N log N) depends logarithmically on the desired error tolerance. Our second scheme evaluates a layer potential on the domain itself, with the same accuracy. It produces Md values using N boundary elements in O((N + Md) log M) arithmetic operations. Our third method evaluates a discrete sum of values of the Green function, of the type which occur in particle methods. It attains error ε{lunate} at a cost O(Na log N), where a = 2 (1 + D d) and D is the Hausdorff dimension of the set where the sources concentrate in the limit N √ ∞. Thus it is O(N log N) when the sources do not cluster too much and close to O(N log N) in the important practical case when the points are uniformly distributed over a hypersurface. We also sketch an O(N log N) algorithm based on special functions. Two-dimensional numerical results are presented for all three algorithms. Layer potentials are evaluated to second-order accuracy, in times which exhibit considerable speedups even over a reasonably sophisticated direct calculation. Discrete sum calculations are speeded up astronomically; our algorithm reduces the CPU time required for a calculation with 40,000 points from six months to one hour. © 1992