For a given positive measure on a fixed domain, Gaussian
quadrature routines can be defined via their property of integrating an
arbitrary polynomial of degree $2n+1$ exactly using only $n+1$ quadrature
nodes. In the special case of Gauss--Jacobi quadrature, this means that $$\int_{-1}^1 (1+x)^\alpha(1-x)^\beta f(x) dx =
\sum_{j=0}^{n} w_j f(x_j), \quad \alpha, \beta > -1, $$ whenever $f(x)$ is a
polynomial of degree at most $2n+1$. When $f$ is not a polynomial, but a
function analytic in a neighbourhood of $[-1,1]$, the above is not an equality
but an approximation that converges exponentially fast as $n$ is increased.
An undergraduate mathematician taking a numerical
analysis course could tell you that the nodes $x_j$ are roots of the Jacobi
polynomial $P^{\alpha,\beta}_{n+1}(x)$, the degree $n+1$ polynomial orthogonal
with respect to the given weight, and that the quadrature weight at each node
is related to the derivative $P'^{\alpha,\beta}_{n+1}(x_j)$. However, these
values are not generally known in closed form, and we must compute them
numerically... but how?
Traditional approaches involve applying the recurrence
relation satisfied by the orthogonal polynomials, or solving the Jacobi matrix
eigenvalue problem in the algorithm of Golub and Welsch, but these methods are
inherently limited by a minimal complexity of $O(n^2)$. The current
state-of-the-art is the $O(n)$ algorithm developed by Glasier, Liu, and
Rokhlin, which hops from root to root using a Newton iteration evaluated with a
Taylor series defined by the ODE satisfied by $P^{\alpha,\beta}_{n+1}$.
We propose an alternative approach, whereby the function
and derivative evaluations required in the Newton iteration are computed
independently in $O(1)$ operations per point using certain well-known asymptotic
expansions. We shall review these expansions, outline the new algorithm, and
demonstrate improvements in both accuracy and efficiency.
