Integrate an indicator function in python with quad

Question

I am integrating an indicator function in python with scipy.integrate.quad. However, the value I derived seems to be incorrect.

import numpy as np
from scipy.integrate import quad

def indac(x, xc, rad):
    if xc - rad <= x <= xc + rad:
        return 1
    else:
        return 0

phi = lambda ii, x: np.sin(ii * x)


xc = 0.1586663
rad = 0.01 * np.pi

result, _ = quad(lambda x: phi(1, x) * indac(x, xc, rad), 0., np.pi)
print(result)  # 0.0

a, b = xc - rad, xc + rad
result, _ = quad(lambda x: phi(1, x) * indac(x, xc, rad), a, b)
print(result)  # 0.009925887836572549

This gives me zero. However, the correct value can be acquired if I integrate from a to b as commented in the code.

It seems it is related to the quadrature method: similarly found in the link.

Any hints will be really appreciated!

Answer 1

Note that quad does not evaluate the integrand at all possible points or even many points. As an adaptive quadrature method, it approximates the integral and estimates its own error, and (for efficiency) it terminates when its error estimate falls below the requested tolerances. Its error estimation strategy is excellent for many well-behaved functions, but yours does not satisfy its assumptions. After only 21 function evaluations that return exactly 0.0, its error estimate is zero, so it terminates with an integral of 0.0.

With the tighter interval, it sees the points at which the integrand is 1, so it returns the correct result.

Consider qmc_quad which samples the space at as many quasi-random points as you please:

import numpy as np
from scipy import integrate

def indac(x, xc, rad):
    return (xc - rad <= x) & (x <= xc + rad)

phi = lambda ii, x: np.sin(ii * x)

xc = 0.1586663
rad = 0.01 * np.pi

# The integrand callable needs to be vectorized to evaluate
# the integrand at `n_points` points in a single call. 
# Increase `n_points` for more accurate results.
res = integrate.qmc_quad(lambda x: phi(1, x) * indac(x, xc, rad), 
                         0., np.pi, n_points=10000)

print(res)
# QMCQuadResult(integral=0.009904273812591187, standard_error=1.5619537172522532e-05