Using scipy.integrate.quad to perform 3D integral

Question

Motivation for the question

I'm trying to integrate a function f(x,y,z) over all space.

I have tried using scipy.integrate.tplquad & scipy.integrate.nquad for the integration, but both methods return the integral as 0 (when the integral should be finite). This is because, as the volume of integration increases, the region where the integrand is non-zero gets sampled less and less. The integral 'misses' this region of space. However, scipy.integrate.quad does seem to be able to cope with integrals from [-infinity, infinity] by performing a change of variables...

Question

Is it possible to use scipy.integrate.quad 3 times to perform a triple integral. The code I have in mind would look something like the following:

x_integral = quad(f, -np.inf, np.inf)
y_integral = quad(x_integral, -np.inf, np.inf)
z_integral = quad(y_integral, -np.inf, np.inf)

where f is the function f(x, y, z), x_integral should integrate from x = [- infinity, infinity], y_integral should integrate from y = [- infinity, infinity], and z_integral should integrate from z = [- infinity, infinity]. I am aware that quad wants to return a float, and so does not like integrating a function f(x, y, z) over x to return a function of y and z (as the x_integral = ... line from the code above is attempting to do). Is there a way of implementing the code above?

Thanks

Answer 1

Here is an example with nested call to quad performing the integration giving 1/8th of the sphere volume:

import numpy as np
from scipy.integrate import quad

def fz(x, y):
    return quad( lambda z:1, 0, np.sqrt(x**2+y**2) )[0]

def fy(x):
    return quad( fz, 0, np.sqrt(1-x**2), args=(x, ) )[0]

def fx():
    return quad( fy, 0, 1 )[0]

fx()
>>> 0.5235987755981053

4/3*np.pi/8
>>> 0.5235987755982988

Answer 2

I'm trying to integrate a function f(x,y,z) over all space.

First of all you'll have to ask yourself why the integral should converge at all. Does it have a factor exp(-r) or exp(-r^2)? In both of these cases, quadpy (a project of mine has something for you), e.g.,

import quadpy

scheme = quadpy.e3r2.stroud_secrest_10a()
val = scheme.integrate(lambda x: x[0]**2)
print(val)

2.784163998415853