How to numerically evaluate this double Integral?

I want to evaluate the following integral:

$$\int_{0}^{60} \ \left(\int_{0}^{2z} 0.5\cdot t \left(\mathrm{erf}(t-a) -1 \right)J_{0}(qt)\mathrm{d}t \right)^2 \mathrm{exp}\left(-\frac{(z-a)^2}{2s^2}\right)\mathrm{d}z$$

Where a, q and s are constants. Without the square outside the first integral, this can be easily evaluated using dblquad from scipy as:

import math
from scipy.special import erf, jv

def h(t, z):
return f(t) * g(z)

def f(t):
return 0.5 * t * (erf(t - a) - 1) * jv(0, q * t)

def g(z):
return math.exp(-((z - a) ** 2)/(2 * (s ** 2)))

if __name__ == '__main__':
a, q, s = 0, 1, 2  # set the constants
result, abserr = dblquad(h, 0, 60, lambda z: 0, lambda z: 2 * z)
print(f'result: {result}, abserr: {abserr}')

I looked about the square of integrals, and as shown here square of integral is equivalent to double integral:

$$\left(\int_a^bf(x)\text{d}x\right)^2 = \int_a^b \int_a^b f(x) f(y) \text{d}x\text{d}y$$

This makes the problem more complicated as now I have to evaluate triple integral which makes evaluation slower. Is there any better way to tackle this problem? Any help would be highly appreciated.

• You could use just quad, defining the function to be integrated with another quad, its upper limit being 2z. May 16 '20 at 7:05
• @AmitHochman Yes, that's a good idea but I am not sure if I would be able to used numba for such a function. I am thinking to increase the speed using jit compilation. May 16 '20 at 9:48
• How is speed an issue? Do you need to evaluate this for various values of a, q, or s? May 16 '20 at 12:56
• @AmitHochman Yes, this function is a part of model that will be used for curve fitting. So, it will be evaluated for an array of q values for many times. May 16 '20 at 18:13
• Then it depends on what is the range of q values. If q is large, q*2*60 >> 1 then the inner integral will be difficult to evaluate numerically because the integrand is highly oscillatory. This case would need careful consideration even for a single evaluation. If q is small, then you could replace the Bessel function with its Taylor series. Then, each term could be integrated analytically and you would have an analytical function of q for the inner integral. That would speed things up if you need many q values. May 16 '20 at 20:18