What would be the computationally faster way to implement this 2D numerical integration?

Question

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

212.96751403808594

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.

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

  294.532276869
  

  while the fixed one took

  224.198880911
  

• 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.