SciFracX / FractionalCalculus.jl

FractionalCalculus.jl: A Julia package for high performance, comprehensive and high precision numerical fractional calculus computing.
http://scifracx.org/FractionalCalculus.jl/dev/
MIT License
35 stars 6 forks source link

Numerical error too high #21

Open AndresCenteno opened 6 months ago

AndresCenteno commented 6 months ago

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
ErikQQY commented 6 months ago

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 commented 6 months ago

Also, using such a small step size is bad

AndresCenteno commented 6 months ago

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

ErikQQY commented 6 months ago

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