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

**URL:** <https://discourse.julialang.org/t/how-do-i-choose-between-quadgk-and-cubature-when-i-do-singular-integral-of-a-complex-valued-function/94198>\
**Category:** General Usage\
**Tags:** question\
**Created:** [February 7, 2023, 11:56am UTC](https://discourse.julialang.org/t/how-do-i-choose-between-quadgk-and-cubature-when-i-do-singular-integral-of-a-complex-valued-function/94198 "2023-02-07T11:56:08Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![f1974939505](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/f1974939505/32/44543_2.png) [@f1974939505](https://discourse.julialang.org/u/f1974939505)\
**Post date:** [February 7, 2023, 11:56am UTC](https://discourse.julialang.org/t/how-do-i-choose-between-quadgk-and-cubature-when-i-do-singular-integral-of-a-complex-valued-function/94198/1 "2023-02-07T11:56:08Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 7, 2023, 1:27pm UTC](https://discourse.julialang.org/t/how-do-i-choose-between-quadgk-and-cubature-when-i-do-singular-integral-of-a-complex-valued-function/94198/2 "2023-02-07T13:27:59Z")

</div>

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

```julia
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](https://github.com/JuliaMath/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](https://math.mit.edu/~stevenj/papers/ReidWhJo14.pdf), 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](https://discourse.julialang.org/t/numerical-integration-of-cauchy-principal-value/36059)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 7, 2023, 4:29pm UTC](https://discourse.julialang.org/t/how-do-i-choose-between-quadgk-and-cubature-when-i-do-singular-integral-of-a-complex-valued-function/94198/3 "2023-02-07T16:29:31Z")

</div>

> [@stevengj](#):
>
> 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](https://math.mit.edu/~stevenj/papers/ReidWhJo14.pdf), for example, where the integrands often have integrable singularities.

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:

```julia
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.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [February 8, 2023, 1:24am UTC](https://discourse.julialang.org/t/how-do-i-choose-between-quadgk-and-cubature-when-i-do-singular-integral-of-a-complex-valued-function/94198/4 "2023-02-08T01:24:16Z")

</div>

> [@stevengj](#):
>
> For example, …

I just added this example (along with several others) to the [QuadGK manual](https://juliamath.github.io/QuadGK.jl/stable/quadgk-examples/#Nearly-singular-integrands).

---

<div class="post-metadata">

**Author:** ![f1974939505](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/f1974939505/32/44543_2.png) [@f1974939505](https://discourse.julialang.org/u/f1974939505)\
**Post date:** [February 11, 2023, 1:10am UTC](https://discourse.julialang.org/t/how-do-i-choose-between-quadgk-and-cubature-when-i-do-singular-integral-of-a-complex-valued-function/94198/5 "2023-02-11T01:10:23Z")

</div>

This would be very helpful. Thanks!