lebedev-quadrature

Basic info

Lebedev quadrature is a scheme for numerically computing integrals around the unit sphere:

where is the regular sphere, are the Lebedev quadrature points, and are the quadrature weights.

Library features

• Pure C++ library for computing spherical integrals with Lebedev quadrature
• No dependencies – only uses standard library
• Header only, with the option to compile a shared or static library to speed up compilation
• Can be included as a git submodule for easy inclusion in project

Basic usage

Create a lebedev::QuadraturePoints by passing in a lebedev::QuadratureOrder enum (to see which orders are available, see below):

#include <lebedev_quadrature.hpp>

auto quad_order = lebedev::QuadratureOrder::order_590;
auto quad_points = lebedev::QuadraturePoints(quad_order);

From this, one can pass in a function which evaluates the integrand at a point and returns a value:

// define a regular function
double xyz_squared(double x, double y, double z)
{
    return x*x * y*y * z*z;
}

// or a lambda function
auto lambda_func = [](double x, double y, double z) { return x*x * y*y * z*z; };

double quadrature_value = quad_points.evaluate_spherical_integral(xyz_squared);
double lambda_quadrature_value = quad_points.evaluate_spherical_integral(lambda_func);

or a function which takes all of the quadrature points and returns a vector that is the integrand evaluated at all the quadrature points:

using vec = std::vector<double>;

// define a regular function
vec xyz_squared(const vec& x, const vec& y, const vec& z)
{
    vec return_vec(x.size());
    for (std::size_t i = 0; i < x.size(); ++i)
        return_vec[i] = x[i]*x[i] * y[i]*y[i] * z[i]*z[i];

    return return_vec;
}

// or a lambda function
auto lambda_func = [](const vec& x, const vec& y, const vec& z) { 
                        vec return_vec(x.size());
                        for (std::size_t i = 0; i < x.size(); ++i)
                            return_vec[i] = x[i]*x[i] * y[i]*y[i] * z[i]*z[i];

                        return return_vec;
                   };

auto quadrature_value = quad_points.evaluate_spherical_integral(xyz_squared);
auto lambda_quadrature_value = quad_points.evaluate_spherical_integral(lambda_func);

Alternatively, you can get directly get const references the quadrature points and weights (in case you need it for optimization purposes).

auto qx = quad_points.get_x();
auto qy = quad_points.get_y();
auto qz = quad_points.get_z();
auto qw = quad_points.get_weights();

double quadrature_value = 0.0;
for (std::size_t i = 0; i < w.size(); ++i)
    quadrature_value += (x[i]*x[i] * y[i]*y[i] * z[i]*z[i]) * w[i];

quadrature_value *= 4 * M_PI;

Library installation

add_library(lebedev_implementation
    lebedev_implementation.cpp)
target_include_directories(lebedev_implementation
    PRIVATE
    path/to/lebedev_quadrature/src)
target_include_directories(lebedev_implementation
    PUBLIC
    ${CMAKE_CURRENT_SOURCE_DIR})

git submodule add --depth 1 https://github.com/lucasmyers97/lebedev-quadrature.git lebedev-quadrature

Implementation details

For implementation details, see here.

Tests

Lebedev quadrature is designed to exactly integrate polynomials up to some degree which is dependent on the order of quadrature (see below for details). Hence, to test for bugs in this implementation we numerically integrate all monomials of the proper degree and compare that with analytic answers (indeed, there is a way to exactly integrate polynomials on the sphere). See the here.

Sources

This is essentially a rewrite of code by John Burkardt and Dmitri Laikov. It's lovely code and works perfectly, but it's basically just C. I thought it might be nice to rewrite in (somewhat modern) C++ to use things like RAII and objects. I also took the opportunity to try to make the intent of the code a little more explicit in the README files. Finally, I included a function which will calculate the integral, given a std::function object and a lebedev::QuadraturePoints (for ease of operation).

Available quadrature orders

Below is a table which gives the rule number (this is just a way to enumerate the rules), whether it is available in this library, the precision (i.e. the degree of polynomial it can exactly integrate), and the order (i.e. the number of points in the quadrature scheme). If you're interested in doing this programmatically, the get_rule_order(unsigned int rule_number), get_rule_availability(unsigned int rule_number), and get_rule_precision(unsigned int rule_number) functions all do what their names suggest.