Lebedev.jl

A julia package to compute the Lebedev quadrature rule over the surface of a three-dimensional unit sphere.

This formula allows to approximate a surface integral as follows

where the number of points in the sum depends on the quadrature order chosen. The coefficients generated by the Lebedev rule allow to calculate integrals of polynomials

with a relative accuracy in the order of 5⋅10⁻¹⁴. The points created by the rule have octahedral rotation (Oₕ point group) and inversion symmetry, while the coefficients are normalized such that they sum up to one.

The Lebedev.jl package is tested with Julia 1.5, 1.6 and 1.7 on Linux, MacOS and Windows. Previous versions of Julia starting at 1.0 might be supported as well, but are not tested; use it at your own risk.

To install the latest version of the package, simply enter into the Julia package manager by typing ] in the REPL and issue the command

pkg> add Lebedev

That's it! You're good to go.

After installing the package, type

julia> using Lebedev

in the REPL or add it to your script.

The Lebedev.jl package provides essentially two functions to generate quadrature points and weights:

lebedev_by_points(n::Integer) -> x,y,z,w
lebedev_by_order(n::Integer) -> x,y,z,w

Both of them return four one-dimensional arrays of the same length, the first three containing the x, y and z Cartesian coordinates of the quadrature points (which lie on the unit sphere) and the last one containing the associated weights w. The lebedev_by_points function returns the points and weights corresponding to the n-point Lebedev rule, while the lebedev_by_order function returns the points and weights ensuring "exact" integration for a polynomial of order up to n.

Lebedev determined 65 quadrature rules ranging from order 3 up to order 131, increasing two by two (hence only odd-numbered orders are known), however, this package only provides a subset of them. Use the function isavailable(n::Integer) to know if the quadrature rule for a given order n is available and the function availablerules() to print out all available orders and corresponding number of points. If you need to access the list of available orders or points, use getavailablerules() (returns a order⇒points dictionary), getavailableorders(), and getavailablepoints(); the latters return lists of integers sorted in ascending order.

julia> @time lebedev_by_points(2702);
  0.000053 seconds (8 allocations: 84.812 KiB)

julia> @time lebedev_by_order(89);
  0.000035 seconds (8 allocations: 84.812 KiB)

julia> f(x,y,z) = x^2 * y^4 * z^6
f (generic function with 1 method)

# we need at least a rule of order n = 2+4+6 = 12
julia> isavailable(12)
false

# only odd-numbered rules are implemented, check if n = 13 is available
julia> isavailable(13)
true

julia> @time x,y,z,w = lebedev_by_order(13);
  0.000007 seconds (4 allocations: 2.625 KiB)

# integrates f(x,y,z) = x²y⁴z⁶ on the unit sphere
julia> @time 4 * pi * dot(w,f.(x,y,z))
  0.000063 seconds (5 allocations: 768 bytes)
0.0041846055991871965

Vyacheslav Lebedev, Dmitri Laikov, “A quadrature formula for the sphere of the 131st algebraic order of accuracy”, Doklady Mathematics, 59 (3), 477-481 (1999).

The Lebedev.jl package is released under the GNU General Public License, version 3.0. This implementation is based on a C source code developed by Dmitri Laikov and John Burkardt, which can be found at https://people.sc.fsu.edu/~jburkardt/c_src/sphere_lebedev_rule/sphere_lebedev_rule.html.