# Handling an overflow example in QuadGK
**URL:** https://discourse.julialang.org/t/handling-an-overflow-example-in-quadgk/101657
**Category:** Numerics
**Tags:** integral
**Created:** [July 15, 2023, 1:57pm UTC](https://discourse.julialang.org/t/handling-an-overflow-example-in-quadgk/101657 "2023-07-15T13:57:20Z")
**Posts on this page:** 6
**Page:** 1
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [July 15, 2023, 1:57pm UTC](https://discourse.julialang.org/t/handling-an-overflow-example-in-quadgk/101657/1 "2023-07-15T13:57:20Z")
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](https://github.com/JuliaMath/QuadGK.jl), it gives a `NaN` error:
```julia
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
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](https://en.wikipedia.org/wiki/Double-precision_floating-point_format), due to [underflow and overflow](https://stackoverflow.com/questions/40082459/what-is-overflow-and-underflow-in-floating-point). 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
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](https://math.stackexchange.com/questions/3383316/what-is-the-idea-behind-expm1-to-avoid-cancellation-error) to compute it more accurately. So, I would do something like:
```julia
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.
---
### Author:  [@giordano](https://discourse.julialang.org/u/giordano)
#### Post date: [July 15, 2023, 7:42pm UTC](https://discourse.julialang.org/t/handling-an-overflow-example-in-quadgk/101657/2 "2023-07-15T19:42:08Z")
> [@stevengj](#):
>
> I decided to post my answer here since it might be of more general interest
Might be useful also in QuadGK.jl docs?
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [July 16, 2023, 4:06am UTC](https://discourse.julialang.org/t/handling-an-overflow-example-in-quadgk/101657/3 "2023-07-16T04:06:13Z")
> [@giordano](#):
>
> Might be useful also in QuadGK.jl docs?
I’m not sure; it’s not really about integration at all, since it’s just a matter of calculating the _integrand_ accurately.
---
### Author:  [@jishnub](https://discourse.julialang.org/u/jishnub)
#### Post date: [July 16, 2023, 12:42pm UTC](https://discourse.julialang.org/t/handling-an-overflow-example-in-quadgk/101657/4 "2023-07-16T12:42:58Z")
Is there a package that provides macros to automate such transformations? I guess it’ll have to be a symbolic algebra package
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [July 16, 2023, 12:52pm UTC](https://discourse.julialang.org/t/handling-an-overflow-example-in-quadgk/101657/5 "2023-07-16T12:52:52Z")
> [@jishnub](#):
>
> Is there a package that provides macros to automate such transformations?
I’ve never heard of such a thing (in any language); numerical analysis has thus far been difficult to automate.
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [July 16, 2023, 2:06pm UTC](https://discourse.julialang.org/t/handling-an-overflow-example-in-quadgk/101657/6 "2023-07-16T14:06:11Z")
Herbie does a good job some of the time, but it’s very far from being reliable