How to stabilize scipy.integrate.quad for improper integrals in the Thomas–Hopfield model

Stabilize SciPy quad on the Thomas–Hopfield model: guide adaptive quadrature with points, set epsabs=0 and tight epsrel, and split the 0→∞ integral properly.

When you turn a physics model with two improper integrals into code, small numerical details can make or break your plot. That is exactly what happens with the Thomas–Hopfield model when it is evaluated with scipy.integrate.quad: the output appears to depend too much on the upper integration bound, and vectorizing for plotting adds to the confusion. The good news is that this behavior follows directly from how adaptive quadrature works and can be fixed in a principled way.

Reproducing the issue

The model integrates two kernels over r from 0 to infinity and composes them into a time-dependent signal. Directly integrating to infinity is problematic, so a finite cutoff r_max is used. The sensitivity to r_max is the symptom to focus on.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
# Recombination rate (decays exponentially with distance)
def rate_fn(r, w0, r0):
    return w0 * np.exp((-2 * r) / r0)
# First kernel of the model
def kern_a(r, w0, r0, t):
    return rate_fn(r, w0, r0) * np.exp(-rate_fn(r, w0, r0) * t) * r**2
# Second kernel of the model (used inside an exponent)
def kern_b(r, w0, r0, t):
    return (np.exp(-rate_fn(r, w0, r0) * t) - 1) * r**2
# Thomas-Hopfield model composed from the two integrals
def hopfield_fn(t, w0, r0, dens, amp, bias):
    r_cap = 250e-7  # pragmatic cutoff in place of infinity
    val1 = quad(kern_a, 0, r_cap, args=(w0, r0, t))[0]
    val2 = quad(kern_b, 0, r_cap, args=(w0, r0, t))[0]
    return amp * ((4 * np.pi * dens * val1) * np.exp(4 * np.pi * dens * val2)) + bias
# Vectorized wrapper to accept an array of t
hopfield_vec = np.vectorize(hopfield_fn)
# Example parameters
w0 = 1e6
r0 = 1.2e-9
dens = 3e-19
amp = 1
bias = 0
# Theoretical t=0 value from literature formula
teor_t0 = dens * np.pi * w0 * r0**3
calc_t0 = hopfield_fn(0, w0, r0, dens, amp, bias)
print('Theoretical value for t=0:', teor_t0)
print('Calculated value for t=0:', calc_t0)
# Plot on a log-log grid
T = np.linspace(0, 100000000, 100000)
Y = hopfield_vec(T, w0, r0, dens, amp, bias)
plt.plot(T, Y, color='r', linestyle='-', linewidth=0.7)
plt.yscale('log')
plt.xscale('log')
plt.show()

Changing r_cap stretches the curve, and the match at t=0 can get better with a smaller r_cap than with a larger one. This is not random; it is the expected behavior of an adaptive integrator that has not been told where the difficult region of the integrand is.

What quad actually does and why this matters

As someone using scipy.integrate.quad for the first time, I am asking myself how quad is computing the result (in which order it takes in the array, computes the integrals and returns the values).

quad uses a Gauss–Kronrod adaptive scheme. It samples the integrand at a set of nodes, builds both a low-order and a high-order estimate over the current sub-interval, compares them to estimate local error, and refines only where the discrepancy indicates something interesting is happening. If none of the sampled nodes land near the peak or rapid variation region, the algorithm can mistakenly conclude that the function is flat there and move on.

A minimal example makes this concrete. Integrating a standard normal pdf centered at 0 over a wide symmetric interval returns 1 exactly as you would expect.

import scipy
import numpy as np
pdf0 = scipy.stats.norm(loc=0, scale=1)
val, err = scipy.integrate.quad(pdf0.pdf, -1000, 1000)
print('result', val, 'error', err)

Now shift the same distribution to be centered at 100 and integrate over the same bounds. Theoretically, nothing changes. In practice, if the adaptive sampler never probes near x=100, it concludes the function is negligible everywhere and returns zero.

pdf100 = scipy.stats.norm(loc=100, scale=1)
val, err = scipy.integrate.quad(pdf100.pdf, -1000, 1000)
print('result', val, 'error', err)

To resolve this, you must tell quad where the function has structure. The points argument does exactly that by forcing subdivision around specific coordinates.

val, err = scipy.integrate.quad(pdf100.pdf, -1000, 1000, points=[100])
print('result', val, 'error', err)

Back to the Thomas–Hopfield model, the integrands have scale-sensitive features around small r, roughly near 1e-9 to 1e-7 in the shared example. If the integrator is not nudged to look there, results depend heavily on arbitrary bounds. There is a second ingredient here too. The default absolute tolerance for quad is meant for order-one magnitudes. With extremely small values, the default epsabs is much too large. Setting epsabs to zero and using a relative tolerance gives the integrator the right target accuracy.

A robust fix for the finite bound case

This version explicitly tells quad where the action is and tightens tolerances. The only change in behavior is improved accuracy and stability.

def hopfield_fn(t, w0, r0, dens, amp, bias):
    r_cap = 1.0
    anchors = np.array([1e-9, 1e-8, 1e-7])
    val1 = quad(kern_a, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    val2 = quad(kern_b, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    return amp * ((4 * np.pi * dens * val1) * np.exp(4 * np.pi * dens * val2)) + bias

This immediately removes the r_cap sensitivity within the finite interval and aligns the t=0 value with the theoretical expression when the integral is accurate over the domain of interest.

Extending to integrals over [0, ∞)

The points argument is incompatible with infinite bounds in the underlying library. That is easy to work around by splitting the integral into a finite segment that uses points and a tail segment that runs to infinity without them. In practice, choose a cutoff where the tail is negligible or at least smooth.

def hopfield_fn(t, w0, r0, dens, amp, bias):
    r_cap = 1.0
    anchors = np.array([1e-9, 1e-8, 1e-7])
    a_finite = quad(kern_a, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    a_tail   = quad(kern_a, r_cap, np.inf, args=(w0, r0, t), epsabs=0, epsrel=1.49e-8)[0]
    b_finite = quad(kern_b, 0, r_cap, args=(w0, r0, t), points=anchors, epsabs=0, epsrel=1.49e-8)[0]
    b_tail   = quad(kern_b, r_cap, np.inf, args=(w0, r0, t), epsabs=0, epsrel=1.49e-8)[0]
    a_val = a_finite + a_tail
    b_val = b_finite + b_tail
    return amp * ((4 * np.pi * dens * a_val) * np.exp(4 * np.pi * dens * b_val)) + bias

With this change, the reported theoretical check at t=0 matches over [0, ∞) as follows.

Theoretical value for t=0: 1.628601631620949e-39
Calculated value for t=0: 1.628601631620949e-39

About np.vectorize and alternatives

np.vectorize is not the source of the integration issue. It is a convenience wrapper that effectively loops over t and collects results. A direct Python loop is equivalent and can be easier to debug.

grid_t = np.linspace(0, 1e8, 1000)
vals = np.array([hopfield_fn(tt, w0, r0, dens, amp, bias) for tt in grid_t])

quad itself is a scalar integrator: its integrand is evaluated at scalar r, so passing t one at a time is the right way to drive it.

Optional: automatically finding “interesting” r

In situations where you cannot easily pre-guess good anchor points, one programmatic approach is to find nontrivial stationary points of the integrands and use them to segment the domain. The following program differentiates each kernel symbolically, locates a root of the first derivative away from r=0, and uses that point to place anchors. It then feeds those anchors to quad with tightened tolerances. This is more involved, and is shown here as an automated alternative.

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
import sympy
import mpmath
import scipy
import math
from collections import defaultdict
mpmath.mp.dps = 50
def rate_fn(r, w0, r0):
    return w0 * np.exp((-2 * r) / r0)
def kern_a(r, w0, r0, t):
    return rate_fn(r, w0, r0) * np.exp(-rate_fn(r, w0, r0) * t) * r**2
def kern_b(r, w0, r0, t):
    return (np.exp(-rate_fn(r, w0, r0) * t) - 1) * r**2
w0_symb = sympy.Symbol('w0')
r0_symb = sympy.Symbol('r0')
t_symb = sympy.Symbol('t')
r_symb = sympy.Symbol('r')
w0c = 1e6
r0c = 1.2e-9
def rate_expr():
    return w0_symb * sympy.exp((-2 * r_symb) / r0_symb)
def kern_a_expr():
    return rate_expr() * sympy.exp(-rate_expr() * t_symb) * r_symb**2
def kern_b_expr():
    return (sympy.exp(-rate_expr() * t_symb) - 1) * r_symb**2
def expr_derivative(expr, w0, r0, n=1):
    expr = expr.subs(w0_symb, w0)
    expr = expr.subs(r0_symb, r0)
    deriv = sympy.diff(expr, r_symb, n)
    deriv = sympy.simplify(deriv)
    return sympy.lambdify((r_symb, t_symb), deriv, modules=['mpmath'])
def expand_bracket(func, cond, args, start, mult, maxiter):
    x = start
    if cond(func(x, *args)):
        raise Exception('invalid start point')
    for _ in range(maxiter):
        last_x = x
        x *= mult
        if cond(func(x, *args)):
            return last_x, x
    raise Exception(f'Unable to find starting area for {func} with {args}')
def safe_underflow_wrapper(func):
    def inner(x, *args):
        y_mpf = func(x, *args)
        y_f = float(y_mpf)
        if y_f == 0 and y_mpf != 0:
            y_f = float(mpmath.sign(y_mpf)) * math.ulp(0)
        return y_f
    return inner
def root_after(function, derivative, t, lower_bound):
    t = np.clip(t, 1e-30, np.inf)
    if derivative(lower_bound, t) > 0:
        cond = lambda x: x <= 0
    else:
        cond = lambda x: x > 0
    lo, hi = expand_bracket(derivative, cond, (t,), lower_bound, 5, 100)
    rng = hi - lo
    xtol = rng * 1e-3
    res = scipy.optimize.bisect(
        safe_underflow_wrapper(derivative), lo, hi, args=(t,), xtol=xtol
    )
    return res
ka_d0 = expr_derivative(kern_a_expr(), w0c, r0c, n=0)
ka_d1 = expr_derivative(kern_a_expr(), w0c, r0c, n=1)
kb_d0 = expr_derivative(kern_b_expr(), w0c, r0c, n=0)
kb_d1 = expr_derivative(kern_b_expr(), w0c, r0c, n=1)
last_root_a = None
last_root_b = None
def pick_points_a(t):
    r_star = root_after(ka_d0, ka_d1, t, lower_bound=1e-15)
    return [r_star * 0.01, r_star * 0.9, r_star, r_star * 1.3, r_star * 20]
def pick_points_b(t):
    r_star = root_after(kb_d0, kb_d1, t, lower_bound=1e-15)
    return [r_star * 0.01, r_star * 0.03, r_star, r_star * 1.4, r_star * 20]
def hopfield_fn_auto(t, w0, r0, dens, amp, bias):
    r_cap = 1.0
    pts_a = pick_points_a(t)
    pts_b = pick_points_b(t)
    a_val, _ = quad(kern_a, 0, r_cap, args=(w0, r0, t), points=pts_a, epsabs=0, epsrel=1.49e-8)
    b_val, _ = quad(kern_b, 0, r_cap, args=(w0, r0, t), points=pts_b, epsabs=0, epsrel=1.49e-8)
    return amp * ((4 * np.pi * dens * a_val) * np.exp(4 * np.pi * dens * b_val)) + bias
dens = 3e-19
amp = 1
bias = 0
teor_t0 = dens * np.pi * w0c * r0c**3
calc_t0 = hopfield_fn_auto(0, w0c, r0c, dens, amp, bias)
print('Theoretical value for t=0:', teor_t0)
print('Calculated value for t=0:', calc_t0)
T = np.linspace(0, 100000000, 1000)
Y = np.vectorize(hopfield_fn_auto)(T, w0c, r0c, dens, amp, bias)
plt.plot(T, Y, color='r', linestyle='-', linewidth=0.7)
plt.yscale('log')
plt.xscale('log')
plt.show()

Why this is worth knowing

Adaptive quadrature is powerful but not omniscient. If you do not guide the sampler toward scale-sensitive regions or if tolerances are mismatched to the magnitudes you compute, seemingly arbitrary dependencies on cutoff parameters appear. Understanding points and tolerances lets you keep the physics and remove the artifacts.

Practical wrap-up

When integrating kernels like those in the Thomas–Hopfield model, make the integrator look where the structure is by supplying points. Tighten epsabs for very small magnitudes by setting it to zero and rely on a relative tolerance. If you need the true improper integral from 0 to infinity, split it into a finite part with points and an infinite tail without points. Using np.vectorize is fine; a plain Python loop is equally valid and sometimes more transparent. Finally, when you compute expressions like exp(x) - 1 for small x in similar contexts, consider np.expm1(x) to avoid catastrophic cancellation.

With these adjustments, your plots become stable, your t=0 check lines up with the theoretical value, and you can move on to fitting the model with confidence rather than fighting the integrator.