# 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:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![khaledharizb](https://avatars.discourse-cdn.com/v4/letter/k/57b2e6/32.png) [@khaledharizb](https://discourse.julialang.org/u/khaledharizb)\
**Post date:** [March 24, 2023, 9:55am UTC](https://discourse.julialang.org/t/domain-error-with-quadgk-to-derive-caputo-derivative/96550/1 "2023-03-24T09:55:37Z")

</div>

Hi everyone! I’m trying to evaluate caputo derivative using Quadgk. It works unless when `t=3` it gives me Domain error. With Mathematica I got the value `-3.11522`

```julia
using QuadGK, ForwardDiff, SpecialFunctions
f(t)=t^2*sin(t);
g= t-> ForwardDiff.derivative(t->f(t) , t)
F(t) = x-> g(x)/(gamma(1/2)* √(t-x))
caputo(t) = quadgk(F(t), 0, t)[1]
caputo(3)

```

This is the error  
`DomainError with 2.9999999999999787: integrand produced Inf in the interval (2.9999999999999574, 3.0)`

---

<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!

---

<div class="post-metadata">

**Author:** ![khaledharizb](https://avatars.discourse-cdn.com/v4/letter/k/57b2e6/32.png) [@khaledharizb](https://discourse.julialang.org/u/khaledharizb)\
**Post date:** [March 24, 2023, 3:57pm UTC](https://discourse.julialang.org/t/domain-error-with-quadgk-to-derive-caputo-derivative/96550/3 "2023-03-24T15:57:06Z")

</div>

Geat! Thank you @stevengj , I used Gauss quadrature which was very efficient.