Richard Mikael Slevinsky

Richard Mikael Slevinsky, Hadrien Montanelli, Qiang Du

We present algorithms for solving spatially nonlocal diffusion models on the unit sphere with spectral accuracy in space. Our algorithms are based on the diagonalizability of nonlocal diffusion operators in the basis of spherical harmonics, the computation of their eigenvalues to high relative accuracy using quadrature and asymptotic formulas, and a fast spherical harmonic transform. These techniques also lead to an efficient implementation of high-order exponential integrators for time-dependent models. We apply our method to the nonlocal Poisson, Allen--Cahn and Brusselator equations.

Richard Mikael Slevinsky, Sheehan Olver

We develop a spectral method for solving univariate singular integral equations over unions of intervals by utilizing Chebyshev and ultraspherical polynomials to reformulate the equations as almost-banded infinite-dimensional systems. This is accomplished by utilizing low rank approximations for sparse representations of the bivariate kernels. The resulting system can be solved in ${\cal O}(m^2n)$ operations using an adaptive QR factorization, where $m$ is the bandwidth and $n$ is the optimal number of unknowns needed to resolve the true solution. The complexity is reduced to ${\cal O}(m n)$ operations by pre-caching the QR factorization when the same operator is used for multiple right-hand sides. Stability is proved by showing that the resulting linear operator can be diagonally preconditioned to be a compact perturbation of the identity. Applications considered include the Faraday cage, and acoustic scattering for the Helmholtz and gravity Helmholtz equations, including spectrally accurate numerical evaluation of the far- and near-field solution. The Julia software package SingularIntegralEquations.jl implements our method with a convenient, user-friendly interface.

Richard Mikael Slevinsky

A rapid transformation is derived between spherical harmonic expansions and their analogues in a bivariate Fourier series. The change of basis is described in two steps: firstly, expansions in normalized associated Legendre functions of all orders are converted to those of order zero and one; then, these intermediate expressions are re-expanded in trigonometric form. The first step proceeds with a butterfly factorization of the well-conditioned matrices of connection coefficients. The second step proceeds with fast orthogonal polynomial transforms via hierarchically off-diagonal low-rank matrix decompositions. Total pre-computation requires at best $\mathcal{O}(n^3\log n)$ flops; and, asymptotically optimal execution time of $\mathcal{O}(n^2\log^2 n)$ is rigorously proved via connection to Fourier integral operators.

Richard Mikael Slevinsky

We describe a fast, simple, and stable transform of Chebyshev expansion coefficients to Jacobi expansion coefficients and its inverse based on the numerical evaluation of Jacobi expansions at the Chebyshev--Lobatto points. This is achieved via a decomposition of Hahn's interior asymptotic formula into a small sum of diagonally scaled discrete sine and cosine transforms and the use of stable recurrence relations. It is known that the Clenshaw--Smith algorithm is not uniformly stable on the entire interval of orthogonality. Therefore, Reinsch's modification is extended for Jacobi polynomials and employed near the endpoints to improve numerical stability.

Richard Mikael Slevinsky

We describe a skeletonization of the spherical harmonic connection problem that reduces the storage and pre-computation to superoptimal complexities at the cost of increasing the execution time by the modest multiplicative factor of $\mathcal{O}(\log n)$. One advantage of accelerating the spherical harmonic connection problem over accelerating synthesis and analysis is that neighbouring layers (in steps of two) may be expanded in eachother's bases. The proposed skeletonization maximizes this interconnectivity by overlaying a dyadic partitioning on the connection problem. We derive the symmetric-definite banded generalized eigenvalue problem required to accelerate spherical harmonic transforms. We also include a full analysis of the weighted normalized Jacobi connection problem with applications to fast harmonic polynomial transforms on the disk, triangle, rectangle, deltoid, wedge, and any other geometry with a bivariate analogue of Jacobi polynomials.

Jared Lee Aurentz, Richard Mikael Slevinsky

A mechanism is described to symmetrize the ultraspherical spectral method for self-adjoint problems. The resulting discretizations are symmetric and banded. An algorithm is presented for an adaptive spectral decomposition of self-adjoint operators. Several applications are explored to demonstrate the properties of the symmetrizer and the adaptive spectral decomposition.

Yu Li, Richard Mikael Slevinsky

We present a unified treatment of the Fourier spectra of spherically symmetric nonlocal diffusion operators. We develop numerical and analytical results for the class of kernels with weak algebraic singularity as the distance between source and target tends to $0$. Rapid algorithms are derived for their Fourier spectra with the computation of each eigenvalue independent of all others. The algorithms are trivially parallelizable, capable of leveraging more powerful compute environments, and the accuracy of the eigenvalues is individually controllable. The algorithms include a Maclaurin series and a full divergent asymptotic series valid for any $d$ spatial dimensions. Using Drummond's sequence transformation, we prove linear complexity recurrence relations for degree-graded sequences of numerators and denominators in the rational approximations to the divergent asymptotic series. These relations are important to ensure that the algorithms are efficient, and also increase the numerical stability compared with the conventional algorithm with quadratic complexity.

Julien Molina, Richard Mikael Slevinsky

A rapid algorithm is derived for the Helmholtz--Hodge decomposition on the surface of the sphere in spherical coordinates. The algorithm uncouples modes of spherical harmonics with different absolute order, writes the conversion as barely-overdetermined banded linear systems, and solves them with banded $QR$ decompositions that factor and execute in optimal complexity. Rigorous upper bounds on the $2$-norm relative condition number of the banded linear systems support the observable low error growth with respect to truncation degree.