Hello everyone 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?

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