search

JuliaMath / Cubature.jl

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

h-adaptive algorithm only evaluates the integrand at the interior of the domain? #33

Closed cbruder closed 6 years ago

cbruder commented 6 years ago 
               _
   _       _ _(_)_     |  A fresh approach to technical computing
  (_)     | (_) (_)    |  Documentation: https://docs.julialang.org
   _ _   _| |_  __ _   |  Type "?help" for help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 0.6.0 (2017-06-19 13:05 UTC)
 _/ |\__'_|_|_|\__'_|  |  Official http://julialang.org/ release
|__/                   |  x86_64-pc-linux-gnu

julia> using Cubature
julia> hquadrature(x-> 1/sqrt(1-x^2),0.0,1.0)
(Inf, NaN)

If hquadrature only evaluates the integrand at the interior of the domain I would expect pi/2 in this example.

stevengj commented 6 years ago

The trouble is that it's having difficulty converging to the requested tolerance. It subdivides the integration more and more at the edges until it evaluates the integrand at the endpoint due to roundoff error. If you request a lower tolerance, then it gives the expected result:

julia> Cubature.hquadrature(x-> 1/sqrt(1-x^2),0.0,1.0,reltol=1e-6)
(1.570796265139041, 1.2599034852724116e-6)

The HCubature package uses a similar algorithm with a slightly different error estimation scheme, and can converge to slightly better tolerance on this integrand, but it still runs into trouble if you request too small of an error:

julia> HCubature.hquadrature(x-> 1/sqrt(1-x^2),0.0,1.0,rtol=1e-9)
(1.5707963193638252, 1.360839997999477e-9)

julia> HCubature.hquadrature(x-> 1/sqrt(1-x^2),0.0,1.0,rtol=1e-10)
(Inf, Inf)

If you want to do this kind of singular integral to high accuracy, ultimately you need a specialized quadrature scheme that removes the singularity.

For this particular kind of singularity, a standard technique would be Gauss–Chebyshev quadrature ... those weights can be computed by e.g. https://github.com/ajt60gaibb/FastGaussQuadrature.jl (or you can just use the formula from Wikipedia).

cbruder commented 6 years ago

many thanks!

stevengj commented 6 years ago

A good way to diagnose this is simply to print the integrand evaluations. Define p(f) = x -> let y = f(x); println("f($x) = $y"); y; end and then you can do hquadrature(p(x-> 1/sqrt(1-x^2)),0.0,1.0)