How do I integrate a nearly Dirac-delta-like function with scipy.integrate.quad?

Question

I want to integrate a function from 0 to +inf, figured as below.

It goes to 0 rapidly outside a small interval centered at x=1. The positions of the two maxima are unknown. The positions of the minima is known.

Answer 1

If you know the width of the pulse, but don't know its location you will have to scan your domain using steps 50% of the width.

You can get quad to this this for you using the points argument. However, quad has also a limit of how many times to sample a function, so increase that accordingly.

Here an example using a standard quad call and a call with the recommended parameters

The r1[0] vanishes because quad don't realize that there those pulses exist, r1[0] decreases linearly with the width of the pulse, however the it becomes slower as the pulses becomes narrower.

from scipy.integrate import quad
import numpy as np
w = 1
fn = lambda x: sum(np.exp(-(x - cx)**2/(2*w**2)) for cx in [1, 3**0.5, 2*np.pi])
for w in [1, 0.1, 0.01, 0.001]:
    r1 = quad(fn, -10, 10)
    r2 = quad(fn, 0, 10, 
              points=np.arange(-10, 10, w/2), 
              limit=max(100, int(4*20/w)))
    print(w, r1, r2)