fast double integration over a sphere

Question

I need to perform 5 numerical integrations over the surface of a unit sphere.

I have tried with dblquad but this is taking quite some time.

Here is a MRE:

import numpy as np
from scipy import integrate

def foo():

    f1 = lambda theta, phi: np.cos(theta)**2
    f2 = lambda theta, phi: f1(theta, phi)*np.cos(theta)**2*np.sin(theta)  
    
    ret = np.zeros(5)
    
    ret[0] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]
    ret[1] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]
    ret[2] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]
    ret[3] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]
    ret[4] = integrate.dblquad(f2, 0, 2*np.pi, 0, np.pi)[0]

    return ret
    
def main():
    tst = np.zeros((1000, 5))
    for i in range(1000):
        tst[i] = foo()
        
    
main()

In my machine when I do %timeit main() I get 21.4 s ± 126 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

This code structure is taking a long like to perform the computations, any way to improve this?

I'm assuming these are just example functions to show that it takes a while. Correct? — jared, Aug 23 '23 at 23:25
See https://stackoverflow.com/questions/71244504/reducing-redundancy-for-calculating-large-number-of-integrals-numerically/71245570#71245570. It provides a fast Numba implementation (thanks to `LowLevelCallable`) for `quad` which is similar to `dblquad`. You can also use Cython. Simple precision might also be faster if you do not care about having precise results. — Jérôme Richard, Aug 23 '23 at 23:26
Besides, you can compute `np.cos(theta) ** 2` once (otherwise it will be computed twice, certainly even in Numba). In fact, you can just compute directly `np.cos(theta) ** 4`. Actually you can even precompute `y = sin(x)` and then compute `y-2*y**3+y**5`. This means only 1 trigonometric function is computed per call rather than 3 initially. Trigonometric functions are pretty expensive so avoiding them is better. — Jérôme Richard, Aug 23 '23 at 23:38

Andrej Kesely · Answer 1 · 2023-08-23T23:13:10.310

You can make the computation faster by using numba and jitting the f1/f2 functions:

import numpy as np
from numba import njit
from scipy import integrate


@njit
def f1(theta, phi):
    return np.cos(theta) ** 2


@njit
def f2(theta, phi):
    return f1(theta, phi) * np.cos(theta) ** 2 * np.sin(theta)


def foo():
    ret = np.zeros(5)

    ret[0] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]
    ret[1] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]
    ret[2] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]
    ret[3] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]
    ret[4] = integrate.dblquad(f2, 0, 2 * np.pi, 0, np.pi)[0]

    return ret


def main():
    tst = np.zeros((1000, 5))
    for i in range(1000):
        tst[i] = foo()


main()

Running it with time on my machine (AMD 5700x/Python 3.11):

andrej@MyPC:/app$ time python3 script.py

real    0m2,204s
user    0m0,013s
sys     0m0,004s

For comparison, without the @njit decorator:

andrej@MyPC:/app$ time python3 script.py

real    0m13,505s
user    0m0,014s
sys     0m0,004s

This is perfect. I should look more into these libraries! With this version of the code I get: `1.6 s ± 24.3 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)` which is ~13.5x faster — tommy J, Aug 24 '23 at 08:09

fast double integration over a sphere

1 Answers1