Handling an overflow example in QuadGK

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)

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

julia> exp(-2.5 * 1000) # underflows

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 ythis 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.


Might be useful also in QuadGK.jl docs?

1 Like

I’m not sure; it’s not really about integration at all, since it’s just a matter of calculating the integrand accurately.

Is there a package that provides macros to automate such transformations? I guess it’ll have to be a symbolic algebra package

I’ve never heard of such a thing (in any language); numerical analysis has thus far been difficult to automate.

Herbie does a good job some of the time, but it’s very far from being reliable