... or not. I converted it to C and, well, no luck so far. Here's what it generated for n=15 and machine epsilon = 2e-17, just using IEEE doubles to see if I'm on the right track (abscissas on the left, weights on the right):
Code: Select all
n=15
0.9879925180204854884281 0.0307532419961148821563
0.9372733924007059513883 0.0703660474881080549858
0.8482065834104272061822 0.1071592204671719394948
0.7244177313601700696211 0.1395706779261543240001
0.5709721726085388304739 0.1662692058169939757217
0.3941513470775633853904 0.1861610000155621003071
0.2011940939974345421426 0.1984314853271115786093
-0.0000000000000000000000 0.2025782419255612865072
-0.2011940939974345421426 0.1984314853271115786093
-0.3941513470775633853904 0.1861610000155621003071
-0.5709721726085388304739 0.1662692058169939757217
-0.7244177313601700696211 0.1395706779261543240001
-0.8482065834104272061822 0.1071592204671719394948
-0.9372733924007059513883 0.0703660474881080549858
-0.9879925180204854884281 0.0307532419961148821563
Not even close; the constants from the GSL library, from qk15:
Code: Select all
static const double xgk[8] = /* abscissae of the 15-point kronrod rule */
{
0.991455371120812639206854697526329,
0.949107912342758524526189684047851,
0.864864423359769072789712788640926,
0.741531185599394439863864773280788,
0.586087235467691130294144838258730,
0.405845151377397166906606412076961,
0.207784955007898467600689403773245,
0.000000000000000000000000000000000
};
I have this bad feeling I'm massively misunderstanding what they're trying to do here, or else I've somehow got the computations wrong because the interval is different ( over [-1,1] for Gauss-Legendre, vs. ?? for Gauss-Kronrod).
From what I can tell, the main point of using Gauss-Kronrod (or any of the Gauss-derived quadrature methods) is that it's a
lot faster than the Simpson method, plus it comes with its own error estimate built in. (Kronrod's stroke of genius was designing the sample points such that the abscissas and weights would overlap between the sample and the error estimate.) I may change my mind and try to do this using Simpson's method just because I can at least get the darn thing working!
I may ask the guys who built the qd library if they've ever attempted to do this or if anyone they were aware of had done so.
