search

JuliaMath / Cubature.jl

One- and multi-dimensional adaptive integration routines for the Julia language
Other
126 stars 21 forks source link

Improper integration problem. #37

Closed jmarcellopereira closed 6 years ago

jmarcellopereira commented 6 years ago 
using QuadGk
value , error = quadgk(f9, 1, 1e20)
(0.785398163397316, 2.03069694813391e-9)

using Cubature 
value , error = hquadrature(f9, 1, 1e20)
(1.5707963267948977, 3.9302603196994814e-10)
stevengj commented 6 years ago

Putting in a huge number for the integration limit is a terrible way to do improper integrals numerically, because any quadrature scheme will initially sample the integrand extremely coarsely over the interval and will probably get fooled, e.g. if the integrand decays rapidly. (In this case, both quadgk and hquadrature are using the same Gauss–Kronrod quadrature schemes IIRC, but I think they default to different orders.)

A better way to do improper integrals is to do a coordinate transformation to a finite interval, as explained here: https://github.com/stevengj/cubature#infinite-intervals … if you do quadgk(f9, 1, Inf), the quadgk routine will do this coordinate transformation for you automatically.

(There are other fancier methods, e.g. Gauss–Laguerre quadrature, but they depend upon you knowing the asymptotic decay rate of your integrand.)