Domain error with QuadGK to derive Caputo derivative

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

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)

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 causes it to divide by zero.

If you request a lower tolerance, it works:

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

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!

2 Likes

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