[Python] What coordinate transform does scipy.integrate.quad use for infinite bounds?

In SciPy's quad function, it is possible to integrate over an integral with infinite bounds.

import scipy
import numpy as np

def normal_distribution_pdf(x):
return np.exp(-x**2 / 2) / np.sqrt(2 * np.pi)

result, eps, info = scipy.integrate.quad(normal_distribution_pdf, 0, np.inf, full_output=True)
print(result)


According to the SciPy docs, infinite bounds are handled using a coordinate transformation, using the QUADPACK qagie subroutine.

qagie

handles integration over infinite intervals. The infinite range is mapped onto a finite interval and subsequently the same strategy as in QAGS is applied.
​

But it is vague about how the infinite range is mapped onto the finite range. Is it something like log(x)? Or another function?

I also examined the QUADPACK documentation. It provides the destination range, [0, 1], but is also vague about how it happens.

Another question discusses what happens when a sin or cos weighted infinite integral is computed, in which case qawfe and Fourier integrals are used. However, I am specifically asking about unweighted integrals where weight=None.

I am interested in knowing the specifics of how an infinite range is mapped to the range [0, 1], because I am interested in understanding the limitations of quad().

Continue reading...