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.

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

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

12 Likes

This would be very helpful. Thanks!