diff --git a/include/integratorxx/quadratures/gausscheby1.hpp b/include/integratorxx/quadratures/gausscheby1.hpp index 7647010..c9b28b7 100644 --- a/include/integratorxx/quadratures/gausscheby1.hpp +++ b/include/integratorxx/quadratures/gausscheby1.hpp @@ -16,6 +16,11 @@ namespace IntegratorXX { * \f{eqnarray*}{ x_{i} = & \cos \left( \displaystyle \frac {2i-1} {2n} \pi * \right) \\ w_{i} = & \displaystyle \frac \pi n \f} * + * Reference: + * Abramowitz and Stegun, Handbook of Mathematical Functions with + * Formulas, Graphs, and Mathematical Tables, Tenth Printing, + * December 1972, p. 889 + * * To transform the rule to the form \f$ \int_{-1}^{1} f(x) {\rm d}x \f$, the * weights have been scaled by \f$ w_i \to w_i \sqrt{1-x_{i}^2} \f$. * diff --git a/include/integratorxx/quadratures/gausscheby2.hpp b/include/integratorxx/quadratures/gausscheby2.hpp index 6b47c20..aaefbe2 100644 --- a/include/integratorxx/quadratures/gausscheby2.hpp +++ b/include/integratorxx/quadratures/gausscheby2.hpp @@ -17,6 +17,11 @@ namespace IntegratorXX { * \right) \\ w_{i} = & \displaystyle \frac \pi {n+1} \sin^2 \left( \frac {i} * {n+1} \pi \right) \f} * + * Reference: + * Abramowitz and Stegun, Handbook of Mathematical Functions with + * Formulas, Graphs, and Mathematical Tables, Tenth Printing, + * December 1972, p. 889 + * * To transform the rule to the form \f$ \int_{-1}^{1} f(x) {\rm d}x \f$, the * weights have been scaled by \f$ w_i \to w_i \sqrt{1-x_{i}^2} \f$. * diff --git a/include/integratorxx/quadratures/gausscheby3.hpp b/include/integratorxx/quadratures/gausscheby3.hpp new file mode 100644 index 0000000..49137c7 --- /dev/null +++ b/include/integratorxx/quadratures/gausscheby3.hpp @@ -0,0 +1,84 @@ +#pragma once + +#include + +namespace IntegratorXX { + +/** + * @brief Implementation of Gauss-Chebyshev quadrature of the third kind. + * + * This quadrature is originally derived for integrals of the type + * + * \f$ \displaystyle \int_{0}^{1} f(x) \sqrt{\frac {x} {1-x}} {\rm d}x \f$. + * + * The original nodes and weights from the rule are given by + * + * \f{eqnarray*}{ x_{i} = & \cos^2 \displaystyle \frac {(2i-1) \pi} {2(2n+1)} \\ w_{i} = & \displaystyle \frac {2 \pi} {2n+1} x_i \f} + * + * Reference: + * Abramowitz and Stegun, Handbook of Mathematical Functions with + * Formulas, Graphs, and Mathematical Tables, Tenth Printing, + * December 1972, p. 889 + * + * To transform the rule to the form \f$ \int_{0}^{1} f(x) {\rm d}x \f$, the + * weights have been scaled by \f$ w_i \to w_i \sqrt{\frac {1-x_{i}} {x_i}} \f$. + * + * @tparam PointType Type describing the quadrature points + * @tparam WeightType Type describing the quadrature weights + */ + +template +class GaussChebyshev3 + : public Quadrature> { + + using base_type = Quadrature>; + +public: + using point_type = typename base_type::point_type; + using weight_type = typename base_type::weight_type; + using point_container = typename base_type::point_container; + using weight_container = typename base_type::weight_container; + + GaussChebyshev3(size_t npts, point_type lo, point_type up) + : base_type(npts, lo, up) {} + + GaussChebyshev3(const GaussChebyshev3 &) = default; + GaussChebyshev3(GaussChebyshev3 &&) noexcept = default; +}; + +template +struct quadrature_traits> { + + using point_type = PointType; + using weight_type = WeightType; + + using point_container = std::vector; + using weight_container = std::vector; + + inline static std::tuple + generate(size_t npts, point_type lo, point_type up) { + + const weight_type pi_ov_2n_p_1 = M_PI / (2*npts + 1); + + weight_container weights(npts); + point_container points(npts); + for (size_t i = 1; i <= npts; ++i) { + // 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]; + + // 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]); + + // Finally, convert from [0,1] to [-1,1] + points[i] = 2*points[i] - 1.0; + weights[i] *= 2.0; + } + + return std::make_tuple(points, weights); + } +}; + + +} // namespace IntegratorXX diff --git a/test/1d_quadratures.cxx b/test/1d_quadratures.cxx index 9c7196f..310ffbe 100644 --- a/test/1d_quadratures.cxx +++ b/test/1d_quadratures.cxx @@ -7,6 +7,7 @@ #include #include #include +#include #include #include