A simple example is
f(x) = 1 / (x + im*10^(-6))
and calculate the integral of f(x) from -1 to 1. How do I write the Julia code?
Note that I don’t want to substitute f(x) with a delta function here.
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
For example, suppose that you integrating:
I = \int_a^b \frac{g(x)}{x - i\alpha} dx
for a small 0 < \alpha \ll 1. For \alpha \to 0^+, it approaches i\pi g(0) if a = -b, but for small \alpha > 0 you have to numerically integrate (for a general function g(x)) a function with a sharp spike at x=0, which will require a large number of quadrature points. But you can subtract out the singularity analytically:
I = \int_a^b \left[ \frac{g(x)-g(0)}{x - i\alpha} + \frac{g(0)}{x - i\alpha} \right] dx \\
= \int_a^b \frac{g(x)-g(0)}{x - i\alpha}dx + \underbrace{g(0) \left[\frac{1}{2}\log(x^2 + \alpha^2) + i\tan^{-1}(x/\alpha) \right]_a^b}_{I_0}
and then you only need to numerically integrate I - I_0, which has the spike subtracted. A little caution is required in the tolerance of this numerical integral because you only need a small relative error compared to I_0, not compared to the remainder integral, so you should pass an absolute tolerance like atol = rtol*I₀ if you want a certain relative tolerance rtol in the overall integral I. Or an even simpler trick is just to add I_0 / (b-a) to the integrand of your numerical integral.
Another useful trick is to tell quadgk to compute \int_a^0 + \int_0^b, i.e. give it x=0 as an explicit endpoint, since we know the integrand is badly behaved there (and may even be Inf or NaN if \alpha = 0). This ensures that quadgk never evaluates the integrand exactly at x=0, and also can help it adaptively subdivide the domain more efficiently.
In code:
using QuadGK
function int_slow(g, α, a, b; kws...)
if a < 0 < b
# put an explicit endpoint at x=0 since we know it is badly behaved there
return quadgk(x -> g(x) / (x - im*α), a, 0, b; kws...)
else
return quadgk(x -> g(x) / (x - im*α), a, b; kws...)
end
end
function int_fast(g, α, a, b; kws...)
g₀ = g(0)
denom_int(x) = log(x^2 + α^2)/2 + im * atan(x/α)
I₀ = g₀ * (denom_int(b) - denom_int(a))
if a < 0 < b
# put an explicit endpoint at x=0 since we know it is badly behaved there
return quadgk(x -> I₀/(b-a) + (g(x) - g₀) / (x - im*α), a, 0, b; kws...)
else
return quadgk(x -> I₀/(b-a) + (g(x) - g₀) / (x - im*α), a, b; kws...)
end
end
If you do int_slow(cos, 1e-6, -1, 1) and int_fast(cos, 1e-6, -1, 1) (with the default rtol ≈ 1e-8), they agree to about 13 digits, but the slow brute-force method requires 1230 function evaluations while the fast singularity-subtracted method requires only 31 function evaluations.
As an added bonus, int_fast works even for α = 0, where it gives you i\pi g(0) (for 0 \in (a,b)) plus the Cauchy principal part.
13 Likes
I just added this example (along with several others) to the QuadGK manual.
5 Likes
This would be very helpful. Thanks!