Add Unit Tests for Gauss-Chebyshev Quadratures

include/integratorxx/quadratures/gausscheby1.hpp

@@ -1,6 +1,7 @@
 #pragma once
 
 #include <integratorxx/quadrature.hpp>
+#include <iostream>
 
 namespace IntegratorXX {
 
@@ -64,16 +65,21 @@ struct quadrature_traits<GaussChebyshev1<PointType, WeightType>> {
     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);

include/integratorxx/quadratures/gausscheby2.hpp

@@ -65,14 +65,19 @@ struct quadrature_traits<GaussChebyshev2<PointType, WeightType>> {
     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);

include/integratorxx/quadratures/gausscheby2_mod.hpp

@@ -68,9 +68,10 @@ struct quadrature_traits<GaussChebyshev2Modified<PointType, WeightType>> {
     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;

include/integratorxx/quadratures/gausscheby3.hpp

@@ -62,18 +62,20 @@ struct quadrature_traits<GaussChebyshev3<PointType, WeightType>> {
 
     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);

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<double, double> quad(order, -1., 1.);
+    integrate(quad);
+  }
+  SECTION("Second Kind") {
+    IntegratorXX::GaussChebyshev2<double, double> quad(order, -1., 1.);
+    integrate(quad);
+  }
+  SECTION("Second Kind (Modified)") {
+    IntegratorXX::GaussChebyshev2Modified<double, double> quad(order, -1., 1.);
+    integrate(quad);
+  }
+  SECTION("Third Kind") {
+    IntegratorXX::GaussChebyshev3<double, double> quad(order, -1., 1.);
+    integrate(quad);
+  }
+}
+
 TEST_CASE( "Euler-Maclaurin Quadratures", "[1d-quad]" ) {
 
   IntegratorXX::MurrayHandyLaming<double,double> 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]);
   }