Fully adaptive algorithms for multivariate integral equations using the non-standard form and multiwavelets with applications to the Poisson and bound-state Helmholtz kernels in three dimensions

Abstract

We have developed and implemented a general formalism for fast numerical solution of time-independent linear partial differential equations as well as integral equations through the application of numerically separable integral operators in d ≥ 1 dimensions using the non-standard (NS) form. The proposed formalism is universal, compact and oriented towards the practical implementation into a working code using multiwavelets. The formalism is applied to the case of Poisson and bound-state Helmholtz operators in d = 3. Our algorithms are fully adaptive in the sense that the grid supporting each function is obtained on the fly while the function is being computed. In particular, when the function g = O  f is obtained by applying an integral operator O, the corresponding grid is not obtained by transferring the grid from the input function f. This aspect has significant implications that will be discussed in the numerical section. The operator kernels are represented in a separated form with finite but arbitrary precision using Gaussian functions. Such a representation combined with the NS form allows us to build a sparse, banded representation of Green’s operator kernel. We have implemented a code for the application of such operators in a separated NS form to a multivariate function in a finite but, in principle, arbitrary number of dimensions. The error of the method is controlled, while the low complexity of the numerical algorithm is kept. The implemented code explicitly computes all the 2^2d components of the d-dimensional operator. Our algorithms are described in detail in the paper through pseudo-code examples. The final goal of our work is to be able to apply this method to build a fast and accurate Kohn–Sham solver for density functional theory.

Keywords:

• non-standard form
• multiwavelets
• Green’s operator
• separation of kernels
• sparse grids
• sparse operators
• truncation of functions and operators
• operator bandwidths
• guaranteed finite arbitrary precision
• scaling of computational cost with system size

Acknowledgements

We acknowledge the support received from the Norwegian Research Council through a Center of Excellence Grant (Grant No. 179568/V30), a YFF Grant to KR (Grant No 162746/V00), as well as through a grant of computer time from the Norwegian Supercomputing Program. It is a privilege for us to dedicate this paper to our long-term mentor, collaborator and friend Prof. Trygve Helgaker on the occasion of his 60th birthday.

Notes

1. A couple of examples of useful alternative tests: (1) for a very narrow Gaussian, the initial refinement can be ‘forced’ if the node contains the Gaussian centre even if the quadrature is still not able to reveal it; and (2) for a homogeneous kernel, refinement is only necessary in one direction, the other directions being obtained by symmetry.

2. In practice, this has been computed by projecting g on an adaptive grid with a 100 times tighter accuracy than the one for the operator application and taking the difference between the two representations. In the Supporting Information (available with the online version of this article), an analysis of the pointwise error computed the same way can also be found.