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.integrate import dblquad
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.
