How do I choose between QuadGK and Cubature when I do singular integral of a complex-valued function?

Since this is a 1d integral, I would use QuadGK:

julia> quadgk(x -> 1 / (x + im*10^(-6)), -1, 1, rtol=1e-8)
(-3.3306690738754696e-16 - 3.141590653589794im, 1.4895084794082606e-9)

Note that it correctly finds that the integral is close to -i\pi. (In fact, the exact answer is -2im*atan(1e6), which it’s finding to nearly 16 digits. Because QuadGK is h-adaptive, it will place more quadrature points close to the singularity. (But it still needs 1245 evaluation points, as you can see by wrapping the integrand in counter(f) = x -> (global count += 1; f(x)) and initializing a global count = 0.)

Don’t use Cubature.jl for this — since Cubature.jl is based on a C library, it doesn’t directly support any return type other than Float64 or Vector{Float64}, so to integrate complex-valued integrands you would have use a vector of the real and imaginary parts. HCubature.jl is the native-Julia analogue of Cubature.jl, and supports complex-valued integrands directly. In 1d, HCubature’s algorithm is the same as the default one in QuadGK, but QuadGK’s is more optimized for the 1d case (and also provides the option of using higher-order rules).

That being said, if you are doing lots of near-singular integrals like this, especially in higher dimensions, you should consider using a semi-analytical singularity-subtraction procedure to handle the near-singular part analytically and only use quadrature for the rest. This is a common procedure in integral-equation methods, for example, where the integrands often have integrable singularities.

See also this thread if you are interested in Cauchy principal values: Numerical integration of cauchy principal value

10 Likes