python で Lebedev quadrature を使う。

注意：情報を更新しました。
koideforest.hatenadiary.com

注意：以下の内容は古いので、テストコードはそのままでは通りません。

以前にLebedev求積法について紹介した。
koideforest.hatenadiary.com

最近ではpythonでもLebedev求積法が使えるmoduleがあるらしい。
quadpy · PyPI

とりあえずテストコードを書いた。

import numpy as np
import quadpy

def f_cart( xyz ):
    # np.shape( xyz ) -> ( 3, number of points )
    # n_points = np.shape( xyz )[1]
    return np.ones( n_points )
    # DO NOT WRITE "return 1"

def f_sphere( phi, theta ):
    # phi: azimuth
    # theta: polar
    # np.shape( phi ) = np.shape( theta ) = ( number of points )
    n_points = np.size( theta )
    return np.ones( n_points )
    # DO NOT WRITE "return 1"

scheme = quadpy.sphere.lebedev_019()

# \int^\pi_0 d\theta \sin\theta \int^{2\pi}_0 d\phi = 4 \pi

xyz_center = ( 0., 0., 0. )
radius = 1.
val_cart = scheme.integrate( f_cart , xyz_center, radius )

val_sphere  = scheme.integrate_spherical( f_sphere )

print( val_cart, val_sphere, 4*np.pi )

以前に紹介したFortranのコードとは異なり、こちらは内部で を掛けてくれるため、被積分関数にだけ気を付ければ良い。

注意するべきは、被積分関数の仕様である。
quadpyの中身の一部を抜粋すると以下のようになる。

class SphereScheme:
    def __init__(self, name, weights, points, azimuthal_polar, degree, citation):
        self.name = name
        self.degree = degree
        self.citation = citation

        if weights.dtype == numpy.float64:
            self.weights = weights
        else:
            assert weights.dtype in [numpy.dtype("O"), numpy.int64]
            self.weights = weights.astype(numpy.float64)
            self.weights_symbolic = weights

        if points.dtype == numpy.float64:
            self.points = points
        else:
            assert points.dtype in [numpy.dtype("O"), numpy.int64]
            self.points = points.astype(numpy.float64)
            self.points_symbolic = points

        if azimuthal_polar.dtype == numpy.float64:
            self.azimuthal_polar = azimuthal_polar
        else:
            assert azimuthal_polar.dtype in [numpy.dtype("O"), numpy.int64]
            self.azimuthal_polar = azimuthal_polar.astype(numpy.float64)
            self.azimuthal_polar_symbolic = azimuthal_polar
        return

    def integrate(self, f, center, radius, dot=numpy.dot):
        """Quadrature where `f` is defined in Cartesian coordinates.
        """
        center = numpy.array(center)
        rr = numpy.multiply.outer(radius, self.points)
        rr = numpy.swapaxes(rr, 0, -2)
        ff = numpy.array(f((rr + center).T))
        return area(radius) * dot(ff, self.weights)

    def integrate_spherical(self, f, dot=numpy.dot):
        """Quadrature where `f` is a function of the spherical coordinates
        `azimuthal` and `polar` (in this order).
        """
        flt = numpy.vectorize(float)
        rr = numpy.swapaxes(flt(self.azimuthal_polar), 0, -2)
        ff = numpy.array(f(rr.T[0], rr.T[1]))
        return area(1.0) * dot(ff, flt(self.weights))

def area(radius):
    return 4 * numpy.pi * numpy.array(radius) ** 2

重要なのは "dot(ff, self.weights)" であり、和を内積で処理しているため、"ff" は "numpy.ndarray" の型でないといけない。
つまり、定義した被積分関数はスカラーでなく配列を返すように工夫する必要がある。
上で挙げたテストコードで、被積分関数が "1" だからといって単純に "1" のみを返すようにすると、"dot(1, self.weights) -> self.weights" となるため、積分値（スカラー）ではなく（各点の重みを表す）配列が返って来てしまう。
また、直交座標を用いる "integrate" では、"ff" に渡す座標の配列が "self.points" に対して転置を取っているため、配列の順番に気を付ける必要がある。

余談：
”numpy.swapaxes(rr, 0, -2)” で何をやっているのか不明だが、とりあえず気にしなくても動くっぽい。
swapaxes
numpy.swapaxes — NumPy v1.15 Manual