FractionalCalculus.jl icon indicating copy to clipboard operation
FractionalCalculus.jl copied to clipboard

Numerical error too high

Open AndresCenteno opened this issue 11 months ago • 4 comments

I'm doing the most simple Caputo derivative $$^C_0 D_t^\alpha t = \frac{t^{1-\alpha}}{\Gamma(2-\alpha)}$$ and I get a relative error of $0.024$ for a time step of $10^{-9}$, is this supposed to be working properly?

using FractionalCalculus
using SpecialFunctions: gamma

T = 0.6; α = 0.9; Δt = 1e-9
numerical_caputo_der = fracdiff(t->t,α,0,T,Δt,CaputoDirect())[1]
theoretical_caputo_der = T^(1-α)/gamma(2-α)
abs(numerical_caputo_der-theoretical_caputo_der)/abs(theoretical_caputo_der) # 0.0239

AndresCenteno avatar Mar 22 '24 11:03 AndresCenteno

What is the version of package you are using for the computing? I just tested on master branch: Numerical:

julia> numerical_caputo_der = fracdiff(t->t,α,T,0.001,CaputoDiethelm())
0.9987906107845391

Analytical:

theoretical_caputo_der = T^(1-α)/gamma(2-α)
0.9987906107845398

Error:

julia> abs(numerical_caputo_der-theoretical_caputo_der)/abs(theoretical_caputo_der)
6.669404053086289e-16

ErikQQY avatar Mar 22 '24 12:03 ErikQQY

Also, using such a small step size is bad

ErikQQY avatar Mar 22 '24 12:03 ErikQQY

FractionalCalculus v0.2.10 I tried with a bunch of stepsizes! But the CaputoDiethelm works for me, so... it's ok now I saw in the docs that you use CaputoTrap() for example, but it throws me an error saying it is not defined

AndresCenteno avatar Mar 22 '24 14:03 AndresCenteno

Just tagged a new major release, which should fix the "not defined" errors

ErikQQY avatar Mar 23 '24 08:03 ErikQQY