Larger error on specific points

[putianyi889]

I was using this package to benchmark another algorithm I'm developing and I found something weird.

I evaluated mittleff(0.125, -x^0.125) where x in 0:0.01:2. The function gives 1e-15 accuracy uniformly except for some points where the errors suddenly jump by magnitudes.

I wrote up a naive implementation using the infinite series and the results show that the errors come from this package rather than my algorithm.

function mittleff_series(α,β,z,n)
    ret=BigFloat(0);
    for k=0:n
        ret += z^k/gamma(β+α*k);
    end
    return ret
end

first make sure this method converges:

julia> mittleff_series(BigFloat(0.125),1,-1,500)
0.4819520815350483462338668282319733966720057670467247536146370806425991075905526

julia> mittleff_series(BigFloat(0.125),1,-1,1000)
0.4819520815350483462338668282319733966720057670467247536146370806425991075905526

then compare with MittagLeffler.jl:

julia> mittleff(0.125,1,-1)
0.4819520815350949

you can see the last 3 digits are wrong. Surely this is still pretty accurate, but it's unexpectedly worse than what we would expect from other points:

julia> xgrid=-(0:0.01:2).^0.125;

julia> y1=mittleff.(0.125,1,xgrid);

julia> y2=mittleff_series.(BigFloat(0.125),1,xgrid,1000);

julia> plot(0:0.01:2,log10.(abs.(y1-y2)))

You can see the two spikes from the plot.

[jlapeyre]

Thanks. This is bad behavior. mittleff is rather complicated. In this case it dispatches to quadgk eventually. It will take a bit of work to figure this out.

[jlapeyre]

By changing the order of the integration scheme here
https://github.com/jlapeyre/MittagLeffler.jl/blob/750bcc2c1e9376fb78ead2becc31ca40122190ba/src/mittag_leffler.jl#L24
I can decrease the error in you example. But, it increases elsewhere. Choosing order=5 minimizes the maximum over the errors in your example range above for order between 1 and 9.

EDIT: with the following

function test_quadgk()
    xgrid=-(0:0.01:2).^0.125;
    y1=mittleff.(0.125,1,xgrid);
    y2=mittleff_series.(BigFloat(0.125),1,xgrid,1000);
    diffs = abs.(y1 .- y2)
    return extrema(diffs)
end

and order=17, I find test_quadgk() gives about (0.0, 1.10e-15)

[jlapeyre]

This commit jlapeyre@352c2a0 decreases the error in your examples.