Gauss–Kronrod points and weights?

[stevengj]

It would be nice to have support for generating Gauss–Kronrod points and weights for the weight functions in FastGaussQuadrature, at least when they exist, in order to have an error estimate from the resulting quadrature rule.

JuliaMath/QuadGK.jl/pull/83 should land shortly, and provides an implementation of an O(n²) algorithm by Laurie (1997) that takes as input a SymTridiagonal Jacobi matrix J for a weight function, and outputs a "Kronrod–Jacobi matrix" QuadGK.kronrodjacobi(J, n) whose eigenvalues/eigenvectors yield an order-n Gauss–Kronrod rule (exactly as an ordinary Jacobi matrix yields the Gaussian quadrature rule). Arbitrary-precision arithmetic is supported.

The two caveats are (a) it's O(n²), and (b) a Gauss–Kronrod rule may not always exist for all weight functions. (If real points and weights don't exist, kronrodjacobi throws an error.) Even if the Kronrod rule exists, it may include quadrature points that lie outside of the original interval.

[dlfivefifty]

I thought Gauss–Kronrod tends to be low order eg n = 15. So does point (a) really matter?

[stevengj]

h-adaptive quadrature schemes based on Gauss–Kronrod typically use relatively low order. That's why I don't worry about the O(n²) in QuadGK. For p-adaptive quadrature, one does often go to high order, but in this case you typically want a whole family of nested rules, not just two — either Gauss–Patterson or Clenshaw–Curtis, typically.

In principle, however, you can use Gauss–Kronrod at any order if you want an error estimate from your quadrature scheme.

Since "FastGaussQuadrature.jl" was originally about O(n) schemes, I mainly wanted to warn you about the O(n²) complexity of Laurie's algorithm in case you felt like it breaks with the spirit of this package.

[dlfivefifty]

Sounds more like a research question (not just a coding one) in the first instance. There would have to be a very strong application for this to be pursued though...

[stevengj]

The application is anyone who wants to use a quadrature rule and have an error estimate - it’s not that mysterious. Personally, I find Gaussian quadrature rules hard to use in practice because they lack error estimates.

But if you only want O(n) algorithms in this package, I agree that’s a research problem for Kronrod.

[dlfivefifty]

I don't think we are not religious about O(n) so it could live here.

I see your point about error estimate.