Lebedev quadrature

In numerical analysis, Lebedev quadrature, named after Vyacheslav Ivanovich Lebedev, is an approximation to the surface integral of a function over a three-dimensional sphere. The grid is constructed so to have octahedral rotation and inversion symmetry. The number and location of the grid points together with a corresponding set of integration weights are determined by enforcing the exact integration of polynomials (or equivalently, spherical harmonics) up to a given order, leading to a sequence of increasingly dense grids analogous to the one-dimensional Gauss-Legendre scheme.

The Lebedev grid is often employed in the numerical evaluation of volume integrals in the spherical coordinate system, where it is combined with a one-dimensional integration scheme for the radial coordinate. Applications of the grid are found in fields such as computational chemistry and neutron transport.^[1][2]

The surface integral of a function over the unit sphere,

${\displaystyle I[f]=\int \mathrm {d} \Omega \ f(\Omega )=\int _{0}^{\pi }\sin(\theta )\mathrm {d} \theta \int _{0}^{2\pi }\mathrm {d} \varphi \ f(\theta ,\varphi ),}$

is approximated in the Lebedev scheme as

${\displaystyle {\tilde {I}}[f]=4\pi \sum _{i=1}^{N}\ w_{i}f(\theta _{i},\varphi _{i}),}$

where the particular grid points and grid weights are to be determined. The use of a single sum, rather than two one dimensional schemes from discretizing the θ and φ integrals individually, leads to more efficient procedure: fewer total grid points are required to obtain similar accuracy. A competing factor is the computational speedup available when using the direct product of two one-dimensional grids. Despite this, the Lebedev grid still outperforms product grids.^[3] However, the use of two one-dimensional integration better allows for fine tuning of the grids, and simplifies the use of any symmetry of the integrand to remove symmetry equivalent grid points.

The Lebedev grid points are constructed so as to lie on the surface of the three-dimensional unit sphere and to be invariant under the octahedral rotation group with inversion.^[4] For any point on the sphere, there are either five, seven, eleven, twenty-three, or forty-seven equivalent points with respect to the octahedral group, all of which are included in the grid. Further, all points equivalent under the rotational and inversion group share the same weights. The smallest such set of points is constructed from all six permutations of (±1, 0, 0) (collectively denoted as a¹), leading to an integration scheme

${\displaystyle {\tilde {I}}_{6}[f]=A_{1}\sum _{i=1}^{6}f(a_{i}^{1}),}$

Distinct classes of grid points

	Typical element	Constraint	Number of points
${\displaystyle a^{1}\,}$	${\displaystyle (1,0,0)\,}$		6
${\displaystyle a^{2}\,}$	${\displaystyle {\frac {1}{\sqrt {2}}}(1,1,0)}$		12
${\displaystyle a^{3}\,}$	${\displaystyle {\frac {1}{\sqrt {3}}}(1,1,1)}$		8
${\displaystyle b^{k}\,}$	${\displaystyle (l_{k},l_{k},m_{k})\,}$	${\displaystyle 2l_{k}^{2}+m_{k}^{2}=1}$	24
${\displaystyle c^{k}\,}$	${\displaystyle (p_{k},q_{k},0)\,}$	${\displaystyle p_{k}^{2}+q_{k}^{2}=1}$	24
${\displaystyle d^{k}\,}$	${\displaystyle (r_{k},S_{k},W_{k})\,}$	${\displaystyle r_{k}^{2}+S_{k}^{2}+W_{k}^{2}=1}$	48

where the grid weight is A₁. Geometrically these points correspond to the vertices of a regular octahedron when aligned with the Cartesian axes. Two more sets of points, corresponding to the centers and vertices of the octahedron, are all eight uncorrelated permutations of ${\displaystyle 1/{\sqrt {3}}(\pm 1,\pm 1,\pm 1)}$ (denoted as a³), and all twelve permutations of ${\displaystyle 1/{\sqrt {2}}(\pm 1,\pm 1,0)}$ (denoted as a²). This selection of grid points gives rise to the scheme

${\displaystyle {\tilde {I}}_{26}[f]=A_{1}\sum _{i=1}^{6}f(a_{i}^{1})+A_{2}\sum _{i=1}^{12}f(a_{i}^{2})+A_{3}\sum _{i=1}^{8}f(a_{i}^{3}),}$

where A₁, A₂, and A₃ are the weight functions that still need to be determined. Three further types of points can be employed as shown in the table. Each of these types of classes can contribute more than one set of points to the grid. In complete generality, the Lebedev scheme is

${\displaystyle {\begin{aligned}{\tilde {I}}_{N}[f]=A_{1}\sum _{i=1}^{6}f(a_{i}^{1})&+A_{2}\sum _{i=1}^{12}f(a_{i}^{2})+A_{3}\sum _{i=1}^{8}f(a_{i}^{3})\\&+\sum _{k=1}^{N_{1}}B_{k}\sum _{i=1}^{24}f(b_{i}^{k})+\sum _{k=1}^{N_{2}}C_{k}\sum _{i=1}^{24}f(c_{i}^{k})+\sum _{k=1}^{N_{3}}D_{k}\sum _{i=1}^{48}f(d_{i}^{k}),\end{aligned}}}$

where the total number of points, N, is

${\displaystyle N=26+24(N_{1}+N_{2})+48N_{3}.\,}$

but in some cases A₂ is set to zero and the number of points is

${\displaystyle N=14+24(N_{1}+N_{2})+48N_{3}.\,}$

The determination of the grid weights is achieved by enforcing the scheme to integrate exactly all polynomials up to a given order. On the unit sphere, this is equivalent to integrating all spherical harmonics up to the same order. This problem is simplified by a theorem of Sergei Lvovich Sobolev implying that this condition need be imposed only on those polynomials which are invariant under the octahedral rotation group with inversion.^[5] Enforcing these conditions leads to a set of nonlinear equations which have been solved and tabulated up to order 131 in the polynomial.^{[4][6][7][8][9][10]}

• Fortran code for evaluating Lebedev grid points and weights
• Python codes: quadpy and CasperBeentjes
• [1] Downloadable grid points