We review recent advances in algorithms for quadrature, transforms, differential equations and singular integral equations using orthogonal polynomials. Quadrature based on asymptotics has facilitated optimal complexity quadrature rules, allowing for efficient computation of quadrature rules with millions of nodes. Transforms based on rank structures in change-of-basis operators allow for quasi-optimal complexity, including in multivariate settings such as on triangles and for spherical harmonics. Ordinary and partial differential equations can be solved via sparse linear algebra when set up using orthogonal polynomials as a basis, provided that care is taken with the weights of orthogonality. A similar idea, together with low-rank approximation, gives an efficient method for solving singular integral equations. These techniques can be combined to produce high-performance codes for a wide range of problems that appear in applications.

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 O(m 2 n) 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 O(mn) 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.

Received on ; revised on ]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.

A quantum anharmonic oscillator is defined by the Hamiltonian H = − d 2 dx 2 + V (x), where the potential is given byUsing the Sinc collocation method combined with the double exponential transformation, we develop a method to efficiently compute highly accurate approximations of energy eigenvalues for anharmonic oscillators. Convergence properties of the proposed method are presented. Using the principle of minimal sensitivity, we introduce an alternate expression for the mesh size for the Sinc collocation method which improves considerably the accuracy in computing eigenvalues for potentials with multiple wells.We apply our method to a number of potentials including potentials with multiple wells. The numerical results section clearly illustrates the high efficiency and accuracy of the proposed method. All our codes are written using the programming language Julia and are available upon request.

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 O(n 3 log n) flops; and, asymptotically optimal execution time of O(n 2 log 2 n) is rigorously proved via connection to Fourier integral operators.

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.

Abstract. Sturm-Liouville problems are abundant in the numerical treatment of scientific and engineering problems. In the present contribution, we present an efficient and highly accurate method for computing eigenvalues of singular Sturm-Liouville boundary value problems. The proposed method uses the double exponential formula coupled with Sinc collocation method. This method produces a symmetric positive-definite generalized eigenvalue system and has exponential convergence rate. Numerical examples are presented and comparisons with single exponential Sinc collocation method clearly illustrate the advantage of using the double exponential formula.

Abstract. We investigate the use of conformal maps for the acceleration of convergence of the trapezoidal rule and Sinc numerical methods. The conformal map is a polynomial adjustment to the sinh map, and allows the treatment of a finite number of singularities in the complex plane. In the case where locations are unknown, the so-called Sinc-Padé approximants are used to provide approximate results. This adaptive method is shown to have almost the same convergence properties. We use the conformal maps to generate high accuracy solutions to several challenging integrals, nonlinear waves, and multidimensional integrals. 1. Introduction. The trapezoidal rule is one of the most well-known methods in numerical integration. While the composite rule has geometric convergence for periodic functions, in other cases it has been used as the starting point of effective methods, such as Richardson extrapolation [36] and Romberg integration [37]. The geometric convergence breaks down with endpoint singularities, and this issue inspired a different approach to improve on the composite rule. From the Euler-Maclaurin summation formula, it was noted that some form of exponential convergence can be obtained for integrands which vanish at the endpoints, suggesting that undergoing a variable transformation may well induce this convergence [38,40,48]. After this observation, the race was on to determine exactly which variable transformation, and therefore which decay rate, is optimal. Numerical experiments showed the exceptional promise of rules such as the tanh substitution [9], the erf substitution [49], the IMT rule [19], and the tanh-sinh substitution [50], among others [22]. But exactly which one is optimal, and in which setting?Using a functional analysis approach, this question was beautifully answered by establishing the optimality of a double exponential endpoint decay rate for the trapezoidal rule on the real line for approximating analytic integrands [44]. The domain of analyticity is described in terms of a strip of maximal width π centred on the real axis in the complex plane. This optimality also prescribed the optimal step size and a near-linear convergence rate O(e −kN/ log N ), where N is the number of sample points and k is a constant proportional to the strip width.The results allowed for displays of strong performance for integrals with integrable endpoint singularities without changing the rule in any way [23][24][25]47]. The double exponential transformation was also adapted to Fourier and general oscillatory integrals in [34,35]. Recognizing the trapezoidal rule as the integration of a Sinc expansion of the integrand, the double exponential advocates adapted their analysis to Sinc approximations [46], and also to all the numerical methods therewith derived, such as Sinc-Galerkin and Sinc-collocation methods [17,28,45] for initial and boundary