diff --git a/include/integratorxx/quadratures/gausscheby1.hpp b/include/integratorxx/quadratures/gausscheby1.hpp index c9b28b7..20963b9 100644 --- a/include/integratorxx/quadratures/gausscheby1.hpp +++ b/include/integratorxx/quadratures/gausscheby1.hpp @@ -1,6 +1,7 @@ #pragma once #include +#include namespace IntegratorXX { @@ -64,16 +65,21 @@ struct quadrature_traits> { weight_container weights(npts); const weight_type pi_ov_npts = M_PI / npts; - const weight_type two_npts_x_pi = 2 * npts * M_PI; + const weight_type pi_ov_2npts = pi_ov_npts / 2; for (size_t i = 0; i < npts; ++i) { + const auto ti = (2.0 * (i + 1) - 1.) * pi_ov_2npts; // The standard nodes and weights are given by - points[i] = std::cos((2.0 * (i + 1) - 1.) / two_npts_x_pi); - weights[i] = pi_ov_npts; - + const auto xi = std::cos(ti); + 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 - weights[i] *= std::sqrt(1. - std::pow(points[i],2)); + wi *= std::sqrt(1. - xi*xi); + + // Store into memory + points[i] = xi; + weights[i] = wi; } return std::make_tuple(points, weights); diff --git a/include/integratorxx/quadratures/gausscheby2.hpp b/include/integratorxx/quadratures/gausscheby2.hpp index aaefbe2..8121112 100644 --- a/include/integratorxx/quadratures/gausscheby2.hpp +++ b/include/integratorxx/quadratures/gausscheby2.hpp @@ -65,14 +65,19 @@ struct quadrature_traits> { weight_container weights(npts); point_container points(npts); for (size_t i = 0; i < npts; ++i) { - // The standard nodes and weights are given by - points[i] = std::cos((i + 1) * pi_ov_npts_p_1); - weights[i] = - pi_ov_npts_p_1 * std::pow(std::sin((i + 1) * pi_ov_npts_p_1), 2); + const auto ti = (i + 1) * pi_ov_npts_p_1; + // The standard nodes are given by + const auto xi = std::cos(ti); + // Avoid float point cancellations - form modified weight directly in following statements // However, since we want the rule with a unit weight factor, we // divide the weights by sqrt(1-x^2). - weights[i] /= std::sqrt(1.0 - std::pow(points[i], 2)); + // PI / (n+1) * sin^2(x) / sqrt(1 - cos^2(x)) = PI / (n+1) * sqrt(1 - cos^2(x)) + const auto wi = pi_ov_npts_p_1 * std::sqrt(1.0 - xi*xi); + + // Store into memory + points[i] = xi; + weights[i] = wi; } return std::make_tuple(points, weights); diff --git a/include/integratorxx/quadratures/gausscheby2_mod.hpp b/include/integratorxx/quadratures/gausscheby2_mod.hpp index d793cfe..c66a895 100644 --- a/include/integratorxx/quadratures/gausscheby2_mod.hpp +++ b/include/integratorxx/quadratures/gausscheby2_mod.hpp @@ -68,9 +68,10 @@ struct quadrature_traits> { point_container points(npts); weight_container weights(npts); for (size_t i = 1; i <= npts; ++i) { - auto sine = sin(i * M_PI * oonpp); - auto sinesq = sine * sine; - auto cosine = cos(i * M_PI * oonpp); + const auto ti = i * M_PI * oonpp; + const auto sine = std::sin(ti); + const auto sinesq = sine * sine; + const auto cosine = std::cos(ti); points[i - 1] = 1.0 - 2.0 * i * oonpp + M_2_PI * (1.0 + 2.0 / 3.0 * sinesq) * cosine * sine; diff --git a/include/integratorxx/quadratures/gausscheby3.hpp b/include/integratorxx/quadratures/gausscheby3.hpp index 49137c7..18f053a 100644 --- a/include/integratorxx/quadratures/gausscheby3.hpp +++ b/include/integratorxx/quadratures/gausscheby3.hpp @@ -62,18 +62,20 @@ struct quadrature_traits> { weight_container weights(npts); point_container points(npts); - for (size_t i = 1; i <= npts; ++i) { + for (size_t i = 0; i < npts; ++i) { + const auto ti = 0.5 * (2*(i+1) - 1) * pi_ov_2n_p_1; + const auto cti = std::cos(ti); // The standard nodes and weights are given by - points[i] = std::cos(0.5 * (2*i - 1) * pi_ov_2n_p_1); - weights[i] = 2.0*pi_ov_2n_p_1 * points[i]; + 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((1-x)/x). - weights[i] /= std::sqrt((1.0 - points[i])/points[i]); + // divide the weights by sqrt(x/(1-x)). + wi *= std::sqrt((1.0 - xi)/xi); // Finally, convert from [0,1] to [-1,1] - points[i] = 2*points[i] - 1.0; - weights[i] *= 2.0; + points[i] = 2*xi - 1.0; + weights[i] = 2.0 * wi; } return std::make_tuple(points, weights); diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 310ffbe..159d995 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -141,8 +141,9 @@ TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) { auto f = [=]( double x ){ return gaussian(x); }; double res = 0.; - for( auto i = 0; i < quad.npts(); ++i ) + for( auto i = 0; i < quad.npts(); ++i ) { res += wgt[i] * f(pts[i]); + } CHECK( res == Catch::Approx(ref_gaussian_int(-1.,1.)) ); } @@ -167,6 +168,41 @@ TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) { } +TEST_CASE( "Gauss-Chebyshev Quadratures", "[1d-quad]") { + + constexpr unsigned order = 200; + auto integrate = [&](auto& quad) { + const auto& pts = quad.points(); + const auto& wgt = quad.weights(); + + auto f = [=]( double x ){ return gaussian(x); }; + + double res = 0.; + for( auto i = 0; i < quad.npts(); ++i ) { + res += wgt[i] * f(pts[i]); + } + + CHECK( res == Catch::Approx(ref_gaussian_int(-1.,1.)) ); + }; + + SECTION("First Kind") { + IntegratorXX::GaussChebyshev1 quad(order, -1., 1.); + integrate(quad); + } + SECTION("Second Kind") { + IntegratorXX::GaussChebyshev2 quad(order, -1., 1.); + integrate(quad); + } + SECTION("Second Kind (Modified)") { + IntegratorXX::GaussChebyshev2Modified quad(order, -1., 1.); + integrate(quad); + } + SECTION("Third Kind") { + IntegratorXX::GaussChebyshev3 quad(order, -1., 1.); + integrate(quad); + } +} + TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) { IntegratorXX::MurrayHandyLaming quad( 150 ); @@ -195,7 +231,6 @@ TEST_CASE( "Ahlrichs Quadratures", "[1d-quad]" ) { double res = 0.; for( auto i = 0; i < quad.npts(); ++i ) { - //std::cout << wgt[i] << ", " << pts[i] << std::endl; res += wgt[i] * f(pts[i]); }