Someone asked me this question over email, and I decided to post my answer here since it might be of more general interest.
Suppose that we are trying to compute \int_0^\infty (e^y - 1) \frac{e^{-2.5y}}{y^{1.2}} dy. This integral should be well defined since it decays exponentially as y \to \infty and blows up as 1/y^{0.2} (an integrable singularity) as y \to 0. However, if we naively plug this into QuadGK, it gives a NaN error:
julia> quadgk(y -> (exp(y) - 1) * exp(-2.5y) / y^1.2, 0, Inf)
ERROR: DomainError with 0.875:
integrand produced NaN in the interval (0.75, 1.0)
Is this a bug in QuadGK? No. In fact, if we evaluate this integrand for large y, it really does give NaN:
julia> f(y) = (exp(y)-1) * exp(-2.5y) / y^1.2
f (generic function with 1 method)
julia> f(1000)
NaN
Is this a bug in Julia? No. In fact, exactly the same thing will happen in any language (Python, Matlab, R, etc.) that uses the usual double-precision floating-point (Float64) arithmetic, due to underflow and overflow. Let’s look at this example more closely. The problem is that exp(1000) ≈ 1.97e434 is larger than the largest representable Float64 value, so it overflows to Inf, while exp(-2.5 * 1000) ≈ 1.84e-1086 is smaller than the smallest representable magnitude so it underflows to 0.0:
julia> exp(1000) # overflows
Inf
julia> exp(-2.5 * 1000) # underflows
0.0
and Inf * 0.0 gives NaN — there is no way to evaluate it to a number because the computer has lost track of how large and small the magnitudes are.
How can we avoid this? The usual solution is to refactor your expression to combine the values analytically before they over/underflow. In this case, we want to combine e^y e^{-2.5y} = e^{-1.5y}.
However, we should also be careful at y \to 0, since if you compute e^y - 1 naively by exp(y) - 1, it is subject to severe cancellation error for small y — this is why most numerics libraries provide an expm1 function to compute it more accurately. So, I would do something like:
julia> f(y) = (y < 10 ? expm1(y) * exp(-2.5y) : (exp(-1.5y) - exp(-2.5y))) / y^1.2
f (generic function with 1 method)
julia> quadgk(f, 0, Inf) # default tolerance rtol = √ε ≈ 1e-8
(0.6790524782732376, 9.565804913317641e-9)
julia> quadgk(f, 0, Inf, rtol=1e-14)
(0.6790524809899963, 5.8288819115959285e-15)
which uses expm1 for sufficiently small arguments and explicitly combines the exp factors for large arguments. The quadgk function now works well even if you require a very small error tolerance close to machine precision.