I am interested in doing a 2D numerical integration. Right now I am using the scipy.integrate.dblquad but it is very slow. Please see the code below. My need is to evaluate this integral 100s of times with completely different parameters. Hence I want to make the processing as fast and efficient as possible. The code is:

import numpy as np
from scipy import integrate
from scipy.special import erf
from scipy.special import j0
import time

q = np.linspace(0.03, 1.0, 1000)

start = time.time()

def f(q, z, t):
    return t * 0.5 * (erf((t - z) / 3) - 1) * j0(q * t) * (1 / (np.sqrt(2 * np.pi) * 2)) * np.exp(
        -0.5 * ((z - 40) / 2) ** 2)

y = np.empty([len(q)])
for n in range(len(q)):
    y[n] = integrate.dblquad(lambda t, z: f(q[n], z, t), 0, 50, lambda z: 10, lambda z: 60)[0]

end = time.time()
print(end - start)

Time taken is


This is too much. Please suggest a better way to achieve what I want to do. I have read quadpy can do this job better and very faster but I have no idea how to implement the same. Also, I tried to use cython-prange but scipy doesn't work without gil. I tried numba but again it didn't work for scipy. Please help.


1 Answer 1

  • A very low hanging fruit, namely lifting the constant factor

    0.5 * (1 / (np.sqrt(2 * np.pi) * 2))

    out from the integral, reduces time by a quarter. On my system the original code took


    while the fixed one took

  • Another fruit, a bit above, is to abandon dblquad whatsoever. It is just a wrong tool here.

    Notice that the only dependency on q is in j0(q * t), and that it does not depend on z either. From the mathematical viewpoint it means that an integral of erf(...) * exp(...) over dz can be tabulated once as a function of t, say F(t) for the lack of better name, which then can be fed into the final integration as t * j0(q*t) * F(t).

    Of course you'd need a lambda to interpolate tabulated values, and to manually take care of precision, and maybe do something else that dblquad does under the hood. Nevertheless, expect a thousand-fold speedup.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.