# Domain error with QuadGK to derive Caputo derivative

**URL:** <https://discourse.julialang.org/t/domain-error-with-quadgk-to-derive-caputo-derivative/96550>\
**Category:** General Usage\
**Tags:** quadgk, integral\
**Created:** [March 24, 2023, 9:55am UTC](https://discourse.julialang.org/t/domain-error-with-quadgk-to-derive-caputo-derivative/96550 "2023-03-24T09:55:37Z")\
**Posts on this page:** 1\
**Showing post:** 2

<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:** [March 24, 2023, 12:38pm UTC](https://discourse.julialang.org/t/domain-error-with-quadgk-to-derive-caputo-derivative/96550/2 "2023-03-24T12:38:09Z")

</div>

The basic problem is that, because of the (integrable) singularity in the integrand at x=t, it’s not able to converge the integral to the requested tolerance — it defaults to `rtol=sqrt(eps()) ≈ 1.5e-8`, which is too ambitious for singular integrands in `Float64` precision. It keeps refining/subdividing the integration integral near `x=3`, and eventually gets so close that a [cancellation error](https://en.wikipedia.org/wiki/Catastrophic_cancellation) causes it to divide by zero.

If you request a lower tolerance, it works:

```julia
julia> t = 3; quadgk_count(F(t), 0, t; rtol=1e-7)
(-3.115217543805007, 2.4294414106433614e-7, 1275)

```

or if you use greater precision (which Mathematica probably does automatically):

```julia
julia> F(t) = x-> g(x) / (√oftype(x, π) * √(t-x)); # compute gamma(1/2)=√π in precision of x

julia> t = big"3.0"; quadgk_count(F(t), 0, t; rtol=1e-16)
(-3.115217737066816205260802666768569197165490937152478941478554510552125233724725, 2.921873150345359277761508183141050033640978755863127757988065979623519096541608e-16, 4965)

julia> Float64(quadgk_count(F(t), 0, t; rtol=1e-18)[1]) # rounded to ≈ machine precision
-3.1152177370668164

```

Ideally, to compute such integrals accurately and efficiently, you should use a quadrature technique that “knows” the singularity. For example, you could use [Gauss–Jacobi quadrature](https://juliaapproximation.github.io/FastGaussQuadrature.jl/stable/gaussquadrature/#Gauss-Jacobi-quadrature):

```julia
julia> using FastGaussQuadrature

julia> x, w = gaussjacobi(10, -1/2, 0.0); # 10 points gives answer to about machine precision!

julia> t = 3; sum(w .* (x -> g(x)/√π).((x .+ 1) .* (t/2))) * √(t/2)
-3.115217737066818

```

It’s a little more awkward to use because you have to manually transform the interval to (-1,1) by a change of variables (taking out the 1/\sqrt{t-x} term since that is built into the quadrature rule) with an associated Jacobian factor, and it’s not adaptive so you need to decide how many points to use yourself by e.g. repeatedly doubling the order. But it’s fantastically accurate with just a few points, without going to extended precision!

---

_[View the full topic](https://discourse.julialang.org/t/domain-error-with-quadgk-to-derive-caputo-derivative/96550)._