On the evaluation of prolate spheroidal wave functions and associated quadrature rules

As demonstrated by Slepian et al. in a sequence of classical papers (see Slepian (1983) [33], Slepian and Pollak (1961) [34], Landau and Pollak (1961) [18], Slepian and Pollak (1964) [35], Slepian (1965) [36]), prolate spheroidal wave functions (PSWFs) provide a natural and efficient tool for computing with bandlimited functions defined on an interval. Recently, PSWFs have been becoming increasingly popular in various areas in which such functions occur – this includes physics (e.g. wave phenomena, fluid dynamics), engineering (signal processing, filter design), etc.To use PSWFs as a computational tool, one needs fast and accurate numerical algorithms for the evaluation of PSWFs and related quantities, as well as for the construction of corresponding quadrature rules, interpolation formulas, etc. During the last 15 years, substantial progress has been made in the design of such algorithms – see, for example, Xiao et al. (2001) [40] (see also Bowkamp (1947) [6], Slepian and Pollak (1961) [34], Landau and Pollak (1961) [18], Slepian and Pollak (1964) [35] for some classical results).The complexity of many of the existing algorithms, however, is at least quadratic in the band limit c. For example, the evaluation of the n  th eigenvalue of the prolate integral operator requires O(c2+n2)O(c2+n2) operations (see e.g. Xiao et al. (2001) [40]); the construction of accurate quadrature rules for the integration (and associated interpolation) of bandlimited functions with band limit c   requires O(c3)O(c3) operations (see e.g. Cheng et al. (1999) [8]). Therefore, while the existing algorithms are satisfactory for moderate values of c   (e.g. c⩽103c⩽103), they tend to be relatively slow when c   is large (e.g. c⩾104c⩾104).In this paper, we describe several numerical algorithms for the evaluation of PSWFs and related quantities, and design a class of PSWF-based quadratures for the integration of bandlimited functions. While the analysis is somewhat involved and will be published separately (currently, it can be found in Osipov and Rokhlin (2012) [27]), the resulting numerical algorithms are quite simple and efficient in practice. For example, the evaluation of the n  th eigenvalue of the prolate integral operator requires O(n+c⋅logc) operations; the construction of accurate quadrature rules for the integration (and associated interpolation) of bandlimited functions with band limit c   requires O(c)O(c) operations. All algorithms described in this paper produce results essentially to machine precision. Our results are illustrated via several numerical experiments.