Quadpy should have dblquad.

[shankardutt]

Hi, I am using scipy.integrate.dblquad to integrate the function as shown below. Obviously, it is quite slow. Is there any way to implement the same in quadpy so as to achieve better performance.

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)

Total time is

212.96751403808594

[nschloe]

As mentioned on stackoverflow, https://stackoverflow.com/questions/61026523/what-would-be-the-computationally-faster-way-to-implement-this-2d-numerical-inte?noredirect=1#comment107965006_61026523, feel free to submit a PR.