Error due to Floating Point inaccuracies (Documentation)

[shivin9]

The sum of calculated stencil doesn't come out to be exactly 0 with the error increasing as we increase the order of approximation.
For a stencil for 4th order derivative and 10th order approximation, we have

julia> sum(A.stencil_coefs)
-1.4052778038453617e-14

This is not a problem with Float64 as tested with BigFloat also.
So the solution is to adjust the middle coefficient of the stencil to set the overall sum = 0.
Reference:
SO Post: https://scicomp.stackexchange.com/questions/11249/numerical-derivative-and-finite-difference-coefficients-any-update-of-the-fornb
Improved Foenberg paper: http://epubs.siam.org/doi/pdf/10.1137/S0036144596322507

julia> A.stencil_coefs[7] -= sum(A.stencil_coefs)
1.371740740740740209188154585806529439651058055460453033447265625000000000000000e+01

julia> sum(A.stencil_coefs)
0.000000000000000000000000000000000000000000000000000000000000000000000000000000

julia> res[500]
2.399999481622449663473228831378492031944915652275085449218750000000000000000000e+01

which is what we want.

[dextorious]

Good catch. In general, I've never seen a practical case (like a PDE discretization) where this numerical issue would matter (say, by significantly impacting the numerical stability), but the fix is cheap and there's no reason not to do it.

[shivin9]

[Reminder to self]: Use @code_warntype