Vladimir Rokhlin

Andrei Osipov, Vladimir Rokhlin

As demonstrated by Slepian et. al. in a sequence of classical papers, 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. However, many of the existing algorithms tend to be relatively slow when $c$ is large (e.g. c>10^4). 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, 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 \cdot \log c)$ operations; the construction of accurate quadrature rules for the integration (and associated interpolation) of bandlimited functions with band limit $c$ requires $O(c)$ operations. All algorithms described in this paper produce results essentially to machine precision. Our results are illustrated via several numerical experiments.

Jhu Heitman, James Bremer, Vladimir Rokhlin

We observe that solutions of a large class of highly oscillatory second order linear ordinary differential equations can be approximated using nonoscillatory phase functions. In addition, we describe numerical experiments which illustrate important implications of this fact. For example, that many special functions of great interest --- such as the Bessel functions $J_ν$ and $Y_ν$ --- can be evaluated accurately using a number of operations which is $O(1)$ in the order $ν$. The present paper is devoted to the development of an analytical apparatus. Numerical aspects of this work will be reported at a later date.

Andrei Osipov, Vladimir Rokhlin

As demonstrated by Slepian et. al. in a sequence of classical papers, prolate spheroidal wave functions (PSWFs) provide a natural and efficient tool for computing with bandlimited functions defined on an interval. As a result, PSWFs are becoming increasing popular in various areas in which such function occur - this includes physics (e.g. wave phenomena, fluid dynamics), engineering (e.g. 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 quadratures, interpolation formulas, etc. Even though, for the last half a century, substantial progress has been made in design of such algorithms, 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 at least $O(c^2)$ operations. Therefore, while the existing algorithms are quite satisfactory for moderate values of $c$ (e.g. $c \leq 10^3$), they tend to be relatively slow when $c$ is large (e.g. $c \geq 10^4$). 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. Also, we perform detailed analysis of the related properties of PSWFs. While the analysis is somewhat involved, 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)$ operations; also, the construction of related accurate quadrature rules requires $O(c)$ operations. Our results are illustrated via several numerical experiments.

Ryan Denlinger, Zydrunas Gimbutas, Leslie Greengard et al.

We present a fast summation method for lattice sums of the type which arise when solving wave scattering problems with periodic boundary conditions. While there are a variety of effective algorithms in the literature for such calculations, the approach presented here is new and leads to a rigorous analysis of Wood's anomalies. These arise when illuminating a grating at specific combinations of the angle of incidence and the frequency of the wave, for which the lattice sums diverge. They were discovered by Wood in 1902 as singularities in the spectral response. The primary tools in our approach are the Euler-Maclaurin formula and a steepest descent argument. The resulting algorithm has super-algebraic convergence and requires only milliseconds of CPU time.

William Leeb, Vladimir Rokhlin

This paper introduces a fast and numerically stable algorithm for the solution of fourth-order linear boundary value problems on an interval. This type of equation arises in a variety of settings in physics and signal processing. Our method reformulates the equation as a collection of second-kind integral equations defined on local subdomains. Each such equation can be stably discretized and solved. The boundary values of these local solutions are matched by solving a banded linear system. The method of deferred corrections is then used to increase the accuracy of the scheme. Deferred corrections requires applying the integral operator to a function on the entire domain, for which we provide an algorithm with linear cost. We illustrate the performance of our method on several numerical examples.

James Bremer, Vladimir Rokhlin

We express a certain complex-valued solution of Legendre's differential equation as the product of an oscillatory exponential function and an integral involving only nonoscillatory elementary functions. By calculating the logarithmic derivative of this solution, we show that Legendre's differential equation admits a nonoscillatory phase function. Moreover, we derive from our expression an asymptotic expansion useful for evaluating Legendre functions of the first and second kinds of large orders, as well as the derivative of the nonoscillatory phase function. Numerical experiments demonstrating the properties of our asymptotic expansion are presented.

Hanwen Zhang, Vladimir Rokhlin

We present a scheme for finding all roots of an analytic function in a square domain in the complex plane. The scheme can be viewed as a generalization of the classical approach to finding roots of a function on the real line, by first approximating it by a polynomial in the Chebyshev basis, followed by diagonalizing the so-called ''colleague matrices''. Our extension of the classical approach is based on several observations that enable the construction of polynomial bases in compact domains that satisfy three-term recurrences and are reasonably well-conditioned. This class of polynomial bases gives rise to ''generalized colleague matrices'', whose eigenvalues are roots of functions expressed in these bases. In this paper, we also introduce a special-purpose QR algorithm for finding the eigenvalues of generalized colleague matrices, which is a straightforward extension of the recently introduced componentwise stable QR algorithm for the classical cases (See [Serkh]). The performance of the schemes is illustrated with several numerical examples.

Jhu Heitman, James Bremer, Vladimir Rokhlin et al.

We introduce a version of the asymptotic expansions for Bessel functions $J_ν(z)$, $Y_ν(z)$ that is valid whenever $|z| > ν$ (which is deep in the Fresnel regime), as opposed to the standard expansions that are applicable only in the Fraunhofer regime (i.e. when $|z| > ν^2$). As expected, in the Fraunhofer regime our asymptotics reduce to the classical ones. The approach is based on the observation that Bessel's equation admits a non-oscillatory phase function, and uses classical formulas to obtain an asymptotic expansion for this function; this in turn leads to both an analytical tool and a numerical scheme for the efficient evaluation of $J_ν(z)$, $Y_ν(z)$, as well as various related quantities. The effectiveness of the technique is demonstrated via several numerical examples. We also observe that the procedure admits far-reaching generalizations to wide classes of second order differential equations, to be reported at a later date.

Vladimir Rokhlin, Mark Tygert

We describe an algorithm that, given any full-rank matrix A having fewer rows than columns, can rapidly compute the orthogonal projection of any vector onto the null space of A, as well as the orthogonal projection onto the row space of A, provided that both A and its adjoint can be applied rapidly to arbitrary vectors. As an intermediate step, the algorithm solves the overdetermined linear least-squares regression involving the adjoint of A (and so can be used for this, too). The basis of the algorithm is an obvious but numerically unstable scheme; suitable use of a preconditioner yields numerical stability. We generate the preconditioner rapidly via a randomized procedure that succeeds with extremely high probability. In many circumstances, the method can accelerate interior-point methods for convex optimization, such as linear programming (Ming Gu, personal communication).