Skip to content

Clean up the Gauss-Legendre quadrature code #35

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
wavefunction91 merged 1 commit into from
Jul 9, 2023
Merged

Clean up the Gauss-Legendre quadrature code #35

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
138 changes: 79 additions & 59 deletions include/integratorxx/quadratures/gausslegendre.hpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@ namespace IntegratorXX {


template <typename PointType, typename WeightType>
class GaussLegendre :
class GaussLegendre :
public Quadrature<GaussLegendre<PointType,WeightType>> {

using base_type = Quadrature<GaussLegendre<PointType,WeightType>>;
Expand All @@ -18,18 +18,48 @@ class GaussLegendre :
using weight_type = typename base_type::weight_type;
using point_container = typename base_type::point_container;
using weight_container = typename base_type::weight_container;
GaussLegendre(size_t npts, point_type lo, point_type up):
base_type( npts, lo, up ) { }

GaussLegendre(size_t npts):
base_type( npts ) { }

GaussLegendre( const GaussLegendre& ) = default;
GaussLegendre( GaussLegendre&& ) noexcept = default;
};





/**
* Evaluates the Legendre polynomial at z with the recurrence
* relation \f$(n+1) P_{n+1}(x) = (2n+1) x P_{n}(x) - n P_{n-1}(x)
* \f$, and its derivative with the relation \f$ \frac {x^2-1} {n}
* \frac {\rm d} {{\rm d}x} P_{n}(x) = x P_{n}(x) - P_{n-1}(x) \f$.
*
* Inputs:
* x - x coordinate
* n - order of polynomial
*
* Outputs: pair (p_n, dp_n)
* p_n - value of the Legendre polynomial, P_{n}(x)
* dp_n - value of the derivative of the Legendre polynomial, dP_{n}/dx
*/
template <typename PointType> inline std::pair<PointType,PointType> eval_Pn(PointType x, size_t n) {
// Evaluate P_{n}(x) by iteration starting from P_0(x) = 1
PointType p_n = 1.0;
// Values of P_{n-1}(x) and P_{n-2}(x) used in the recursions
PointType p_n_m1(0.0), p_n_m2(0.0);

// Recurrence can be rewritten as
// P_{n}(x) = [(2n-1) x P_{n-1}(x) - (n-1) P_{n-2}(x)] / n
for( size_t m = 1; m <= n; ++m ){
p_n_m2 = p_n_m1;
p_n_m1 = p_n;
p_n = ((2.0 * m - 1.0) * x * p_n_m1 - (m-1.0)*p_n_m2)/m;
}

// Compute the derivative
// dP_{n}(x) = (x P_{n}(x) - P_{n-1}(x)) * n/(x^2-1)
PointType dp_n = (x * p_n - p_n_m1) * (n / (x * x - 1.0));

return std::make_pair( p_n, dp_n );
}


template <typename PointType, typename WeightType>
Expand All @@ -44,63 +74,53 @@ struct quadrature_traits<
using weight_container = std::vector< weight_type >;

inline static std::tuple<point_container,weight_container>
generate( size_t npts, point_type lo, point_type up ) {

assert( npts % 2 == 0 );
assert( lo != -std::numeric_limits<double>::infinity() );
assert( up != std::numeric_limits<double>::infinity() );

generate( size_t npts ) {
point_container points( npts );
weight_container weights( npts );

// Absolute precision for the nodes
const double eps=3.0e-11;

// Since the rules are symmetric around the origin, we only need
// to compute one half of the points
const size_t mid = (npts + 1) / 2;
const double eps = 3.e-11; // Convergence tolerance

for( size_t i = 1; i <= mid; ++i ) {

const point_type ip = i;

point_type z = std::cos( M_PI * (ip - 0.25) / (npts + 0.5));
point_type pp(0), z1(0);

// Iteratively determine the i-th root
while(std::abs(z-z1) > eps) {
point_type p1(1.), p2(0.);

// Loop over the recurrence relation to evaluate the
// Legendre polynomial at position z
for( size_t j = 1; j <= npts; ++j ){
const point_type jp = j;

point_type p3 = p2;
p2 = p1;
p1 = ( (2. * jp - 1.) * z * p2 - (jp - 1.) * p3) / jp;
} // end j for

// p1 is now the desired Legrendre polynomial. We next compute
// pp, its derivative, by a standard relation involving also p2,
// the polynomial of one lower order
pp = npts * (z * p1 - p2) / (z * z - 1.);
z1 = z;
z = z1 - p1 / pp;
} // end while

point_type pt = z;
weight_type wgt = 2. / (1. - z*z) / pp / pp;

// Transform points and populate arrays
std::tie(pt,wgt) = transform_minus_one_to_one(lo,up,pt,wgt);
points[i-1] = pt;
weights[i-1] = wgt;

// Reflect the points
points[npts-i] = (lo + up) - pt;
weights[npts-i] = wgt;

}; // Loop over points
for( size_t idx = 0; idx < mid; ++idx ) {
// Index
const point_type i = idx+1;

// Standard initial guess for location of i:th root in [-1, 1]
point_type z = std::cos( M_PI * (i - 0.25) / (npts + 0.5));
// Old value of root
point_type z_old = -2.0;
// Values of P_n(x) and dP_n/dx
point_type p_n, dp_n;

// Solve the root to eps absolute precision
while(std::abs(z-z_old) > eps) {
// Evaluate the Legendre polynomial at z and its derivative
auto pn_eval = eval_Pn(z,npts);
p_n = pn_eval.first;
dp_n = pn_eval.second;

// Newton update for root
z_old = z;
z = z_old - p_n / dp_n;
} // end while

// Quadrature node is z
point_type pt = z;
// Quadrature weight is 2 / [ (1-x^2) dPn(x)^2]
weight_type wgt = 2. / ((1. - z*z) * dp_n*dp_n);

// Store the symmetric points
points[idx] = -pt;
weights[idx] = wgt;

points[npts-1-idx] = pt;
weights[npts-1-idx] = wgt;
} // Loop over points

return std::make_tuple( points, weights );

}

};
Expand Down
44 changes: 14 additions & 30 deletions test/1d_quadratures.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -130,42 +130,26 @@ constexpr T real_spherical_harmonics( int64_t l, int64_t m, Args&&... args ) {

TEST_CASE( "Gauss-Legendre Quadratures", "[1d-quad]" ) {

constexpr unsigned order = 10;
for(unsigned order=10;order<14;order++) {
std::ostringstream oss;
oss << "order " << order;
SECTION( oss.str()) {

SECTION( "untransformed bounds" ) {
IntegratorXX::GaussLegendre<double,double> quad( order, -1, 1 );
IntegratorXX::GaussLegendre<double,double> quad( order );

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.)) );
}
const auto& pts = quad.points();
const auto& wgt = quad.weights();

auto f = [=]( double x ){ return gaussian(x); };

SECTION( "transformed bounds" ) {
const double lo = 0.;
const double up = 4.;
IntegratorXX::GaussLegendre<double,double> quad( 100, lo, up );

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]);
double res = 0.;
for( auto i = 0; i < quad.npts(); ++i ) {
res += wgt[i] * f(pts[i]);
}

CHECK( res == Catch::Approx(ref_gaussian_int(lo,up)) );
CHECK( res == Catch::Approx(ref_gaussian_int(-1.,1.)) );
}
}

}

TEST_CASE( "Gauss-Chebyshev Quadratures", "[1d-quad]") {
Expand Down
6 changes: 3 additions & 3 deletions test/quadrature_manipulation.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -7,19 +7,19 @@ TEST_CASE( "Constructor", "[quad]" ) {
using base_type = IntegratorXX::Quadrature<quad_type>;

SECTION("Basic") {
quad_type q( 10, -1, 1 );
quad_type q( 10 );
}

SECTION("Copy To Base") {
quad_type q( 10, -1, 1 );
quad_type q( 10 );
base_type b(q);

CHECK( q.points() == b.points() );
CHECK( q.weights() == b.weights() );
}

SECTION("Move to Base") {
quad_type q( 10, -1, 1 );
quad_type q( 10 );
std::vector<double> pts = q.points();
std::vector<double> wgt = q.weights();

Expand Down