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?
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?
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.