# Numerical integration from 0.0 to Inf

**URL:** <https://discourse.julialang.org/t/numerical-integration-from-0-0-to-inf/40379>\
**Category:** General Usage\
**Created:** [May 28, 2020, 9:37pm UTC](https://discourse.julialang.org/t/numerical-integration-from-0-0-to-inf/40379 "2020-05-28T21:37:41Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Olivier\_Merchiers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivier_merchiers/32/4073_2.png) [@Olivier\_Merchiers](https://discourse.julialang.org/u/Olivier_Merchiers)\
**Post date:** [May 28, 2020, 9:37pm UTC](https://discourse.julialang.org/t/numerical-integration-from-0-0-to-inf/40379/1 "2020-05-28T21:37:41Z")

</div>

Hello everyone,

I was playing around with QuadGK.jl on my function f(x) which has some strong oscillations but a finite value when integrating from 0 to \infty.

I used the transformation described in [this link](https://github.com/stevengj/cubature/blob/master/README.md#infinite-intervals) to solve these integrals.

```julia
using QuadGK
quadgk(g, 0, 1)

```

where g(u) is the transformed function.  
Now I just realized that I could just as well do

```julia
quadgk(f, 0, Inf)

```

where f is my original function f(x).  
And it would give exactly the same result.  
How does this work? Is QuadGK.jl applying the transformation automatically when given Inf to the upper bound?

If indeed some automatic transformation is being performed, does HCubature do a similar transformation on multdimensional integrals?

On a slightly different note, I saw that

```julia
julia> quadgk(x->exp(-x),0.0,Inf)
       
(1.0, 4.5074000326453615e-11)

```

Which indeed corresponds with the analytical solution but then on the other hand:

```julia
julia> quadgk(x->exp(-x),0.0,1e8)
(0.0, 0.0)

```

Some kind of overflow?  
Could someone give me some hints on what is happening?

Many thanks in advance,  
Olivier

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 28, 2020, 9:48pm UTC](https://discourse.julialang.org/t/numerical-integration-from-0-0-to-inf/40379/2 "2020-05-28T21:48:37Z")

</div>

> [@Olivier\_Merchiers](#):
>
> Is QuadGK.jl applying the transformation automatically when given Inf to the upper bound?

Yes. This is mentioned [in the documentation](https://juliamath.github.io/QuadGK.jl/stable/).

> [@Olivier\_Merchiers](#):
>
> If indeed some automatic transformation is being performed, does HCubature do a similar transformation on multdimensional integrals?

Not currently. You have to do any transformation yourself with HCubature.

> [@Olivier\_Merchiers](#):
>
> ```julia
> julia> quadgk(x->exp(-x),0.0,1e8)
> (0.0, 0.0)
> 
> ```
> 
> Some kind of overflow?

In fact, it’s underflow.

The way a quadrature rule like QuadGK works is that it first evaluates your function at some set of points to get an estimated integral and error estimate, and then refines to a denser set of points if needed. This means that you can always “fool” _any_ quadrature routine if your function is exactly zero (or, in fact, matches any low-degree polynomial) on the initial set of points, because then the quadrature rule will think it is done.

That’s what’s happening here, as you can see if you print out the integrand as it is being evaluated:

```julia
julia> quadgk(0, 1e8) do x
           @show x, exp(-x)
           exp(-x)
       end
(x, exp(-x)) = (2.5446043828620757e6, 0.0)
(x, exp(-x)) = (9.745539561713792e7, 0.0)
(x, exp(-x)) = (427231.4439593694, 0.0)
(x, exp(-x)) = (9.957276855604063e7, 0.0)
(x, exp(-x)) = (1.2923440720030276e7, 0.0)
(x, exp(-x)) = (8.707655927996972e7, 0.0)
(x, exp(-x)) = (6.756778832011545e6, 0.0)
(x, exp(-x)) = (9.324322116798845e7, 0.0)
(x, exp(-x)) = (2.970774243113014e7, 0.0)
(x, exp(-x)) = (7.029225756886986e7, 0.0)
(x, exp(-x)) = (2.0695638226615444e7, 0.0)
(x, exp(-x)) = (7.930436177338456e7, 0.0)
(x, exp(-x)) = (5.0e7, 0.0)
(x, exp(-x)) = (3.9610752249605075e7, 0.0)
(x, exp(-x)) = (6.0389247750394925e7, 0.0)
(0.0, 0.0)

```

If your initial interval is [0,10^8], then all of the initial quadrature points are \> 10^6, in which case `exp(-x)` [underflows](https://en.wikipedia.org/wiki/Arithmetic_underflow) to `0.0` in double precision.

---

<div class="post-metadata">

**Author:** ![tbeason](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tbeason/32/15898_2.png) [@tbeason](https://discourse.julialang.org/u/tbeason)\
**Post date:** [May 28, 2020, 9:56pm UTC](https://discourse.julialang.org/t/numerical-integration-from-0-0-to-inf/40379/3 "2020-05-28T21:56:05Z")

</div>

That example was extremely insightful! I think I’ve encountered this before.

---

<div class="post-metadata">

**Author:** ![ericphanson](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ericphanson/32/215186_2.png) [@ericphanson](https://discourse.julialang.org/u/ericphanson)\
**Post date:** [May 28, 2020, 10:09pm UTC](https://discourse.julialang.org/t/numerical-integration-from-0-0-to-inf/40379/4 "2020-05-28T22:09:39Z")

</div>

Interesting, when I read the [`quadgk` docs](https://juliamath.github.io/QuadGK.jl/latest/#QuadGK.quadgk),

> Returns a pair `(I,E)` of the estimated integral `I` and an estimated upper bound on the absolute error `E` . If `maxevals` is not exceeded then `E <= max(atol, rtol*norm(I))` will hold.

I thought `E` was an estimate only in that if `maxevals` is exceeded then all bets are off, but otherwise it should be an upper bound. I see instead that it’s not necessarily an upper bound on the error even if `maxevals` is not exceeded. With a careful second reading though, I see actually nothing is promised about `E` truly being an upper bound in any case. I wonder if the docs could be clarified somehow?

---

<div class="post-metadata">

**Author:** ![Olivier\_Merchiers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/olivier_merchiers/32/4073_2.png) [@Olivier\_Merchiers](https://discourse.julialang.org/u/Olivier_Merchiers)\
**Post date:** [May 29, 2020, 9:43am UTC](https://discourse.julialang.org/t/numerical-integration-from-0-0-to-inf/40379/5 "2020-05-29T09:43:30Z")

</div>

Thank you very much Steven Johnson for this very nice explanation!

> Yes. This is mentioned [in the documentation](https://juliamath.github.io/QuadGK.jl/stable/).

Yes indeed! I missed it. Sorry for that.

> If your initial interval is [0,10^8] , then all of the initial quadrature points are \>10^6

Indeed. When I try a smaller upperbound such as 10^5, I see that quadgk computes more points over the interval than using 10^8 which I find surprising since the interval is smaller.  
Is it because he detects many values are strictly zero and then stops earlier than he should?

Thanks again for the great explanation.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 29, 2020, 12:00pm UTC](https://discourse.julialang.org/t/numerical-integration-from-0-0-to-inf/40379/6 "2020-05-29T12:00:03Z")

</div>

> [@Olivier\_Merchiers](#):
>
> Is it because he detects many values are strictly zero and then stops earlier than he should?

Yes. As far as it can tell, for the interval from [0,10^8] case your function is identically zero and so `quadgk` can stop immediately.