diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp index f07427e..4238886 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev1.hpp @@ -77,8 +77,17 @@ struct quadrature_traits> { auto wi = pi_ov_npts; // However, as we're integrating f(x) not \frac{f(x)}{\sqrt{1 - - // x^2}}, we must factor the \sqrt{1-x^2} into the weight - wi *= std::sqrt(1.0 - xi * xi); + // x^2}}, we must factor the \sqrt{1-x^2} into the weight. + // Since x_i = cos(t_i) with t_i in (0,pi), sqrt(1-x^2) is sin(t_i) + // exactly; forming 1 - x*x instead cancels badly as the nodes approach + // +-1, losing accuracy as O(n^2). + // + // t_i lies in (pi/2,pi) here, where sin(t) ~ pi - t, so the O(ulp(pi)) + // absolute rounding of t_i would show up as an O(n eps) *relative* + // error in sin(t_i), and the reflection below copies it to both ends. + // The complement is exact in the index, t_i = pi - (2 idx + 1) h, so + // evaluate the sine at that small angle instead. + wi *= std::sin((2.0 * idx + 1.0) * pi_ov_2npts); // Store into memory points[idx] = xi; diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp index ef64afe..335e6ac 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev2.hpp @@ -78,8 +78,14 @@ struct quadrature_traits> { const auto xi = std::cos(ti); // The quadrature weight, transformed to unit weight factor in - // [-1,1] is given by (see comments above for explanation) - const auto wi = pi_ov_npts_p_1 * std::sqrt(1.0 - xi * xi); + // [-1,1] is given by (see comments above for explanation). Since + // x_i = cos(t_i) with t_i in (0,pi), sqrt(1-x^2) is sin(t_i) exactly; + // forming 1 - x*x instead cancels badly as the nodes approach +-1. + // + // As in the T1 rule, t_i lies in (pi/2,pi), where the absolute rounding + // of t_i is an O(n eps) relative error in sin(t_i). Use the complement, + // which is exact in the index: t_i = pi - (idx + 1) h. + const auto wi = pi_ov_npts_p_1 * std::sin((idx + 1.0) * pi_ov_npts_p_1); // Store into memory points[idx] = xi; diff --git a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp index bb242d8..3d3b35c 100644 --- a/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp +++ b/include/integratorxx/quadratures/primitive/gausschebyshev3.hpp @@ -72,11 +72,21 @@ struct quadrature_traits> { const auto ti = 0.5 * (2 * i - 1) * pi_ov_2n_p_1; const auto cti = std::cos(ti); const auto xi = cti * cti; // cos^2(t) - auto wi = 2.0 * pi_ov_2n_p_1 * xi; - // However, since we want the rule with a unit weight factor, we - // divide the weights by sqrt(x/(1-x)). - wi *= std::sqrt((1.0 - xi) / xi); + // The standard weight is 2h x_i, and since we want the rule with a unit + // weight factor we divide by sqrt(x/(1-x)). With x_i = cos^2(t_i) and + // t_i in (0,pi/2) the whole product collapses exactly: + // + // 2h cos^2(t) sqrt((1-cos^2 t)/cos^2 t) = 2h cos(t) sin(t) = h sin(2t) + // + // Forming 1 - x instead cancels badly as the nodes approach 1. + // + // 2 t_i is m h with m = 2i - 1 an exact integer, and for m h past pi/2 + // the sine of the supplement is the relatively accurate one, since + // sin(s) ~ pi - s there. Reflect about pi/2 using the index alone. + const size_t m = 2 * i - 1; + const size_t k = (m <= npts) ? m : (2 * npts + 1 - m); + const auto wi = pi_ov_2n_p_1 * std::sin(k * pi_ov_2n_p_1); // Copy to storage points[idx] = xi; diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index a7cc128..fcc5591 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -184,6 +184,60 @@ TEST_CASE( "Gauss-Chebyshev T3 Quadratures", "[1d-quad]" ) { } +TEST_CASE( "Gauss-Chebyshev extremal weights at large N", "[1d-quad]" ) { + + // The extremal weight of each Chebyshev rule is sin(t) at an angle that the + // generators reach as pi - t. Rounding t to the nearest double is an + // O(ulp(pi)) absolute error there, which is an O(N eps) *relative* error in + // the weight: ~1e-13 at N=800, growing with N. Evaluating the sine at the + // complement, which is exact in the loop index, keeps it at O(eps). + // + // The reference is the Taylor series of sin about 0. At these angles + // (h < 2e-3) the first omitted term is O(h^7) ~ 1e-22 relative, so the + // series is good to rounding and is independent of std::sin. + auto sin_series = [](double h) { + const double h2 = h * h; + return h * (1.0 - h2 / 6.0 * (1.0 - h2 / 20.0 * (1.0 - h2 / 42.0))); + }; + + // An absolute 10*eps check cannot see this: the weights are themselves + // O(1/N^2) here, so the error hides below the tolerance. + const double rel_tol = 1e-14; + const size_t npts = 800; + + SECTION("T1") { + // w_0 = (pi/N) sin(h), h = pi/2N, at the node nearest -1 + const double h = M_PI / (2.0 * npts); + const double ref = (M_PI / npts) * sin_series(h); + const GaussChebyshev1 q(npts); + const auto& wgts = q.weights(); + CHECK_THAT( wgts.front(), Catch::Matchers::WithinRel(ref, rel_tol) ); + CHECK_THAT( wgts.back(), Catch::Matchers::WithinRel(ref, rel_tol) ); + } + + SECTION("T2") { + // w_0 = h sin(h), h = pi/(N+1) + const double h = M_PI / (npts + 1); + const double ref = h * sin_series(h); + const GaussChebyshev2 q(npts); + const auto& wgts = q.weights(); + CHECK_THAT( wgts.front(), Catch::Matchers::WithinRel(ref, rel_tol) ); + CHECK_THAT( wgts.back(), Catch::Matchers::WithinRel(ref, rel_tol) ); + } + + SECTION("T3") { + // front() is the node nearest 0: there t = pi/2 - h exactly, so + // x = sin^2(h) and the weight is h sin(2h). back() is the node nearest + // 1, with weight h sin(h). + const double h = M_PI / (2.0 * npts + 1); + const GaussChebyshev3 q(npts); + const auto& wgts = q.weights(); + CHECK_THAT( wgts.front(), Catch::Matchers::WithinRel(h * sin_series(2.0*h), rel_tol) ); + CHECK_THAT( wgts.back(), Catch::Matchers::WithinRel(h * sin_series(h), rel_tol) ); + } + +} + TEST_CASE( "Euler-Maclaurin Quadratures by Murray, Handy, and Laming", "[1d-quad]" ) { IntegratorXX::MurrayHandyLaming quad(150); const auto msg = "Euler-Maclaurin N = " + std::to_string(quad.npts());