2024 CMS Summer Meeting

Saskatchewan, May 30 - June 3, 2024

Numerical Methods for and with Special Functions
Org: James Bremer (University of Toronto), Timon Gutleb (University of British Columbia) and Richard Slevinsky (University of Manitoba)
[PDF]

CADE BALLEW, University of Washington Numerical solutions of Riemann--Hilbert problems on disjoint intervals [PDF]

We present a general approach to numerically compute the solutions of Riemann--Hilbert problems with jump conditions supported on disjoint intervals. Applied to the Fokas--Its--Kitaev Riemann--Hilbert problem, this enables the computation of Chebyshev-like polynomials on multiple intervals, requiring only $\mathrm{O}(N)$ arithmetic operations to compute the first $N$ recurrence coefficients. Moreover, expansions in these orthogonal polynomials yield a novel iterative method for solving indefinite linear systems and computing matrix functions. This method applies in settings where classical polynomial approximations behave poorly and are therefore not applicable. We also discuss an application to the computation of finite-genus and soliton gas solutions of the Korteweg–de Vries equation.

THOMAS BOTHNER, University of Bristol Universality for random matrices with an edge spectrum singularity [PDF]

We study invariant random matrix ensembles defined on complex Hermitian matrices with a single root type singularity and one-cut regular density of states. Assuming that the singularity lies within the soft edge boundary layer we compute asymptotics of the model's generating functional by using Riemann-Hilbert problems for orthogonal polynomials and integrable operators. This extends an old result by Forrester and Witte and is based on ongoing joint work with Toby Shepherd (Bristol).

JAMES BREMER, University of Toronto Frequency-indepedent solvers for linear ODEs [PDF]

I will discuss a class of solvers for linear scalar ordinary differential equations which run in time bounded independent of frequency. They operate by producing exponential representations of a basis in the space of solutions of the equation. These exponential representations can be used to rapidly evaluate any desired solution of the differential equation at any point in the solution domain with accuracy on the order of the condition number of the problem. I will also present a theorem which bounds the complexity of the exponential representations as a function of a measure of the complexity of the equation's coefficients.

I will discuss applications of this work to the numerical evaluation of special functions and the numerical solution of Sturm-Liouville problems.

AMPARO GIL, Universidad de Cantabria Computation and inversion of some cumulative distribution functions [PDF]

Some special functions hold particular importance in Applied Probability and Statistics. Notably, the incomplete gamma and beta functions serve as (with normalization factors) the cumulative central gamma and beta distribution functions, respectively. Additionally, the corresponding noncentral distributions—like the Marcum-Q function and the cumulative noncentral beta distribution function—play significant roles across various applications. These functions' inversion proves valuable in hypothesis testing and random sample generation following the respective probability density functions.

In this talk we describe developments in the asymptotic and numerical computation and inversion of the beta cumulative distribution function (both central and non-central). The effectiveness of these methods will be demonstrated through numerical examples.

Joint work with Javier Segura and Nico M. Temme.

TIMON S. GUTLEB, University of British Columbia A frame approach for equations involving the fractional Laplacian [PDF]

Exceptionally elegant hypergeometric formulae exist for the fractional Laplacian operator applied to weighted classical orthogonal polynomials. We utilize these results to construct a spectral method based on frame properties for solving equations involving the fractional Laplacian on whole space.

MOHAMMAD HAMDAN, University of New Brunswick Polynomials of the Higher Derivatives of the Nield-Kuznetsov Integral Function [PDF]

M.H. Hamdan and T.L. Alderson Abstract Airy’s functions of the first and second kinds represent two linearly independent solutions to the homogeneous Airy’s equation. Airy’s polynomials arise when one considers higher derivatives of Airy’s functions or derivatives of their products. Airy’s functions and associated polynomials are of importance in the study of circuit theory, systems theory and signal processing and arise in solutions to Stark, Schrodinger, and Tricomi’s equations. Many differential equations in quantum theory can be reduced to Airy’s equation by an appropriate change of variables, thus adding to the importance of studies of Airy’s and related functions.

The inhomogeneous Airy’s equation has been shown to have a particular solution expressible in terms of the Nield-Kuznetsov integral function, defined in terms of Airy’s functions and their integrals. Associated with higher derivatives of the Nield-Kuznetsov function are three sets of polynomials of non-equal degrees that are functions of the order of the derivative involved. Two of the arising sets of polynomials are Airy’s polynomials, while the third set arises from coefficients of the Wronskian of Airy’s functions. This Wronskian appears in the Nield-Kuznetsov function and its derivatives. Our objective are: 1) Analyze the Nield-Kuznetsov polynomials and express them in terms of the Nield-Kuznetsov first derivative. 2) Implement the resulting polynomials in expressing the Nield-Kuznetsov function in terms of Bessel functions. 3) Derive ascending series and asymptotic series expressions of the Nield-Kuznetsov function in terms of arising polynomials and use them in computations of these functions.

CECILE PIRET, Michigan Technological University Computing generalized hypergeometric functions in the complex plane using an end-corrected trapezoidal rule [PDF]

Generalized hypergeometric functions $_pF_q$ are ubiquitous in the scientific and engineering fields, for which their accurate evaluation is essential. Although a myriad of algorithms exists for their evaluation in the complex plane, most commonly with low parameters $p$ and $q$, no general framework exists that is adaptable to low and high $p$ and $q$ values, wide ranges of coefficients, and small to large evaluation domains. This results in an intricate patchwork of algorithms with sometimes radically different orders of accuracy and computational cost. We introduce a high order method ($>20^{th}$ order) which addresses this issue by its wide ranging validity. It is based on the Euler's integral transformation of the hypergeometric functions formula, and therefore shares its parameters ($p\le q+1$) and coefficient constraints. The method is based on an end-corrected trapezoidal rule applied to singular integrals, first introduced in [1] in the context of fractional derivatives.

[1] B. Fornberg and C. Piret, Computation of Fractional Derivatives of Analytic Functions. J Sci Comput 96, 79 (2023).

DIEGO RUIZ-ANTOLÍN, Universidad de Cantabria Asymptotic and numerical approximations to the zeros of parabolic cylinder functions [PDF]

The zeros of parabolic cylinder functions find applications in different areas of science and engineering. For example, they are needed in the design and optimization of waveguides, including the determination of cutoff frequencies and propagation characteristics of different modes. In this talk, uniform asymptotic approximations to the zeros of the parabolic cylinder function $U(a,z)$ involving certain combinations of the zeros of Airy functions are discussed. The accuracy of the expansions is tested using a numerical implementation of a method for finding the complex zeros of solutions of second order ODEs described in [2]. For the numerical algorithm, we use the recent results obtained in [1].\newline

[1] T.M. Dunster, A. Gil, J. Segura. Computation of parabolic cylinder functions having complex argument. Appl. Numer. Math. 197 (2024), 230-242. \newline [2] J. Segura. Computing the complex zeros of special functions. Numer. Math. 124 (4) (2013) 723-752.

JAVIER SEGURA, Universidad de Cantabria Computation of classical Gaussian quadratures and associated barycentric interpolation [PDF]

Algorithms for computing the three classical Gaussian rules based on asymptotic methods and globally convergent iterative methods are presented. The Gauss-Radau and Gauss-Lobatto variants are also considered, alongside the computation of barycentric weights for Lagrange interpolation. The asymptotic and iterative algorithms offer distinct advantages: asymptotic methods are highly accurate for large degrees, while iterative methods are generally faster and valid for a broader range of parameters. The combination of both methods provides the fastest and most accurate double precision methods to date, with an extended range of validity compared to previous methods. We also discuss how the iterative methods can be used for arbitrary precision computations, and this is illustrated with some Maple implementations of the algorithms.

RICHARD M. SLEVINSKY, University of Manitoba Fast and stable rational approximation of generalized hypergeometric functions [PDF]

Rational approximations of generalized hypergeometric functions ${}_pF_q$ of type $(n+k,k)$ are constructed by the Drummond and factorial Levin-type sequence transformations. We derive recurrence relations for these rational approximations that require $\mathcal{O}[\max\{p,q\}(n+k)]$ flops. These recurrence relations come in two forms: for the successive numerators and denominators; and, for an auxiliary rational sequence and the rational approximations themselves. Numerical evidence suggests that these recurrence relations are much more stable than the original formulae for the Drummond and factorial Levin-type sequence transformations. Theoretical results on the placement of the poles of both transformations confirm the superiority of factorial Levin-type transformation over the Drummond transformation.

TOM TROGDON, University of Washington Some old and new perspectives on the convergence of spectral methods [PDF]

This talk will concern the convergence theory for Fourier and Chebyshev (ultraspherical) spectral methods for operator equations. The classical convergence theory typically succeeds by showing that the operator under consideration is relatively compact with respect to an operator that is sufficiently simple. In this vein, we discuss the results of G. M. Vainikko [Krasnosel'skii et al., 1972] and a modern reimplementation of the ideas for the convergence of the Fourier-Floquet-Hill method. The ideas are also applicable to the Riemann--Hilbert (Wiener--Hopf) problem on the circle. We then consider the convergence of Chebyshev collocation methods for boundary-value problems and use yet another result of G. M. Vainikko [Krasnosel'skii et al., 1972] to establish convergence of the rectangular collocation method [Driscoll and Hale, 2016], for a special class of boundary conditions. Lastly, building on these ideas, and the work of Olver and Townsend, we develop an ultraspherical collocation method for boundary-value problems that is provably convergent for all (reasonable) boundary conditions.

MOHAN ZHAO, University of Toronto The Approximation of Singular Functions by Series of Non-integer Powers [PDF]

In this talk, we describe an algorithm for approximating functions of the form $f(x) = \langle \sigma(\mu),x^\mu \rangle$ over the interval $[0,1]$, where $\sigma(\mu)$ is some distribution supported on $[a,b]$, with $0<a<b<\infty$. Given a desired accuracy and the values of $a$ and $b$, our method determines a priori a collection of non-integer powers, so that functions of this form are approximated by expansions in these powers, and a set of collocation points, such that the expansion coefficients can be found by collocating a given function at these points. Our method has a small uniform approximation error which is proportional to the desired accuracy multiplied by some small constants, and the number of singular powers and collocation points grows logarithmically with the desired accuracy. This method has applications to the solution of partial differential equations on domains with corners.

top of page