Forcing cutoffs in spectral density

[gefux]

This discussion follows from Issue #46

Currently the code forces the user to specify a cutoff frequency in the spectrad density of a bosonic bath. This is to make sure that the integration of the function works correctly (see explaination in Issue #46).
Is this a reasonable design or is there a more elegant way?

[piperfw]

This is quite an interesting problem which has definitely made me aware of the perils of (blind) numerical quadrature.

For a simple illustration of the kind of issue that can occur, consider the integration by adaptive quadrature of the Gaussian
$f(x) = \exp(-(x-20)^2)/\sqrt{\pi}$ which is a strongly peaked around $x=20$ and normalised such that $\int_{-\infty}^\infty f(x)\ dx =1$.

from scipy.integrate import quad
f = lambda x: np.exp(-(x-20)**2)/np.sqrt(np.pi)
quad(f, 0, 40)     # 1.0 
quad(f, 0, 400)    # 1.0
quad(f, 0, 4000)   # 0.0
quad(f, 0, np.inf) # 0.0

i.e. already by an upper limit $\sim10^3$ the integrator misses the non-zero feature entirely. For this reason I think we should be careful in advising to simply setting a large (hard) cut-off frequency as that can still lead to issues.

I am yet to find a good discussion of this issue online or in a textbook. The quad function does have a points option (note this can only be used with a finite upper limit) to define 'A sequence of break points in the bounded integration interval where local difficulties of the integrand may occur'. This can be used to correct the above integral, e.g.,

quad(f, 0, 4000, points=(10,30)) # 1.0

which is really just a way of telling scipy an integral with a sensible domain for this function is quad(f, 10, 30). Unfortunately just supplying a single point (e.g. maxmium) of interest does not appear to be reliable!

quad(f, 0, 4000, points=(20))  # 1.0
quad(f, 0, 40000, points=(20)) # 0.5 !