Numerical integration from 0.0 to Inf

General Usage
1

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 to solve these integrals.

using QuadGK
quadgk(g, 0, 1)

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

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

2

Yes. This is mentioned in the documentation.

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

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> 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 to 0.0 in double precision.

6 Likes
3

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

4

Interesting, when I read the quadgk docs,

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?

5

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

Yes. This is mentioned in the documentation .

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.

6

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.