QuadGK returns NaNs

Hello everyone :slight_smile:

Would someone know why

using QuadGK
n = 538
_, w, _ = kronrod(n,1e-11,40); @show w

returns a vector of NaN weights at n=538 and not before?

This is a bug: spurious underflow in kronrod for large n · Issue #93 · JuliaMath/QuadGK.jl · GitHub

I should have a fix pushed shortly.

That being said, you will run into other floating-point problems for large n, unless you go to higher precision; Laurie’s algorithm for the Gauss–Kronrod weights is not really designed to scale to n more than a few hundred, I think.

Adding link to associated discussion on StackOverflow: integration - Optimal quadrature rule for heavy tail measure - Computational Science Stack Exchange

Summary: splitting the ray with an endpoint singularity into an interval with endpoint singularity and a ray may lead to more efficient quadratures

Since the measure of interest is a power law, using Gauss quadrature with a Jacobi weight function should perform well for an interval with the endpoint singularity. Computing those quadrature rules is a feature built into QuadGK, and I wrote a wrapper specifically for Jacobi weights called QuadJGK. It can be used to compute a 10-point Gauss rule for polynomials orthogonal w.r.t the weight (1-x)^0.0 * (1+x)^-0.8 for the interval [-1, 1] like this

using QuadJGK
x, w = QuadJGK.jacobigauss(JacobiSpace(0.0, -0.8), 10)

That package is still a work in progress, and it does allow h-adaptive quadrature when Kronrod rules exist for the Jacobi weight, but they do not always exist, according to this paper.

1 Like