How SciPy’s quad Uses QUADPACK qagie to Integrate Over Infinite Intervals via [0,1] Reparameterization

Learn how SciPy quad and QUADPACK qagie map infinite integrals to [0,1] to improve stability, use points, and understand convergence with clear Python examples.

Integrating functions on infinite intervals looks deceptively simple when you call quad in SciPy, but there is precise machinery underneath. Understanding the coordinate mapping used by QUADPACK’s qagie routine helps to reason about numerical stability, convergence, and how to work with features like points on a finite domain.

Problem setup

Here’s a minimal example that integrates the standard normal PDF from 0 to infinity using SciPy. The goal is to see how quad handles infinite bounds.

import scipy
import numpy as np


def pdf_normal(z):
    return np.exp(-(z**2) / 2) / np.sqrt(2 * np.pi)


val, err, dbg = scipy.integrate.quad(pdf_normal, 0.0, np.inf, full_output=True)
print(val)

The SciPy docs state that infinite bounds are handled by a coordinate transformation via QUADPACK’s qagie. In their words:

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.

The interesting part is the exact mapping to the finite interval [0, 1] and what it implies for the integrand shape and numerical behavior.

What actually happens under the hood

The mapping used by qagie is explicit. For a semi-infinite integral with lower limit A, the variable substitution is T in (0, 1] mapped to X in [A, ∞) by X = A + (1 − T)/T. This transformation turns the original integral into one over the unit interval with a modified integrand. For two-sided infinity, the same idea is applied symmetrically by evaluating the underlying function at both X and −X for each T.

The details are documented in: Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. Pages 64–65.

Single infinite limit: [A, ∞)

Use the substitution X = A + (1 − T)/T with T in (0, 1]. When T → 1, X → A. When T → 0+, X → ∞. Differentiating gives dX/dT = −T^(−2). Reversing limits removes the minus sign, yielding an integral over T from 0 to 1 with integrand F(A + (1 − T)/T) · T^(−2). This is the core of how qagie maps a half-infinite integral onto [0, 1].

Two infinite limits: (−∞, ∞)

For integrals over the full real line, the routine applies the same mapping idea but evaluates both tails symmetrically. Each T in (0, 1] corresponds to X = (1 − T)/T and the integrand is sampled as F(X) and F(−X), combined as F(X) + F(−X), again with the T^(−2) weight.

Solution in practice: manual reparametrization

You can apply the same substitution yourself, integrate over [0, 1], and observe how the transformed integrand behaves. This can be useful for diagnostics and for using options like points that are otherwise incompatible with infinite bounds.

import scipy
import numpy as np
import matplotlib.pyplot as plt


def pdf_normal(z):
    return np.exp(-(z**2) / 2) / np.sqrt(2 * np.pi)


def reparam_semi_infinite(fn, a):
    """Map T in (0, 1] to X in [a, +inf) via X = a + (1 - T)/T."""
    def wrapped(t):
        return fn(a + (1 - t) / t) * (t ** -2)
    return wrapped


def reparam_full_line(fn):
    """Map T in (0, 1] to X in (-inf, +inf) using symmetry."""
    def wrapped(t):
        x = (1 - t) / t
        return (fn(x) + fn(-x)) * (t ** -2)
    return wrapped


# Semi-infinite case: integrate normal PDF from 0 to +inf via T in [eps, 1]
pdf_mapped = reparam_semi_infinite(pdf_normal, a=0.0)

# The transformed integrand is singular at T=0 (corresponds to X = +inf),
# so start a bit away from 0.
t0 = 1e-8
res, err, meta = scipy.integrate.quad(pdf_mapped, t0, 1.0, full_output=True)
print(res)

# Visualize the original and transformed integrands
xs1 = np.linspace(0.0, 5.0, 100)
plt.title("Normal PDF in X domain")
plt.plot(xs1, pdf_normal(xs1))
plt.show()

xs2 = np.linspace(t0, 1.0, 100)
plt.title("Transformed integrand in T domain")
plt.plot(xs2, pdf_mapped(xs2))
plt.show()

The integral in T is not defined at T = 0, which corresponds to X = ∞ under the substitution; that’s why the lower limit for T needs to start at a small positive value. With this setup you integrate on [t0, 1] and recover exactly the same value as the original half-infinite integral.

If you need to specify breakpoints for the integrand, the points argument can now be used because the domain is [0, 1]. Breakpoints defined in the original X-domain can be transported to T via t = 1 / (1 − a + x), where a is the lower limit of the original integral.

Why this matters

Knowing the mapping clarifies three practical aspects. First, any singularity at T = 0 is an artifact of the transformation and corresponds to behavior at infinity in the original variable; that helps explain integrator difficulty or the need to avoid T = 0 exactly. Second, inspecting the transformed integrand can reveal whether the integrator has to fight a rapidly varying shape near the endpoints. Third, once the integral lives on [0, 1], features like points become available, which can improve reliability around sharp changes after you relocate breakpoints using the given formula.

Takeaway

quad handles infinite ranges by reparameterizing onto [0, 1]. For [A, ∞) it uses X = A + (1 − T)/T with the T^(−2) weight; for (−∞, ∞) it evaluates both tails symmetrically as F(X) + F(−X) with the same weight. You can reproduce the transformation manually to visualize the integrand, to place breakpoints using t = 1/(1 − a + x), and to better understand where numerical challenges may arise. When you see trouble near the T endpoints, remember it reflects behavior at infinity in the original variable and adjust your integration limits or diagnostics accordingly.