problems with finite difference

[JianghuiDu]

I think there is a problem with the central difference code. Here is a simply example:

julia> Calculus.derivative(x -> x/(x+ 1.4424183196362515e-9),2e-8,:central)
-39.33717713979761

julia> Calculus.derivative(x -> x/(x+ 1.4424183196362515e-9),2e-8,:forward)
1.8509290259770826e6

We know that the analytical derivative is

1.4424183196362515e-9/(x+1.4424183196362515e-9)

which at x=2e-8 is 3.137210795286552e6, similar to the forward difference result (at least on the same order of magnitude) but there is a 5 orders of magnitude difference compared to the central difference result, which is supposed to be 2nd order accurate! I looked through the source code, and I think the problem is with the macro @centralrule

julia> @centralrule 2e-8 epsilon
6.0554544523933395e-6

in comparison

julia> @forwardrule 2e-8 epsilon
1.4901161193847656e-8

The stepping given by @centralrule is clearly unacceptably large (2 orders of magnitude greater than x!. I don't quite understand why the stepping in the @centralrule is taken as the cubic square root of eps() while the @forwardrule uses the square root of eps()?

macro forwardrule(x, e)
    x, e = esc(x), esc(e)
    quote
        $e = sqrt(eps(eltype($x))) * max(one(eltype($x)), abs($x))
    end
end

macro centralrule(x, e)
    x, e = esc(x), esc(e)
    quote
        $e = cbrt(eps(eltype($x))) * max(one(eltype($x)), abs($x))
    end
end

In this case an appropriate stepping show be smaller than the order of 1e-10, which nether of these two methods gives. Why not just use eps() as the step size?

[johnmyleswhite]

Are you sure there's a soluble mistake here? Reasoning based on this specific example doesn't seem like it will provide you with a totally reliable debugging strategy. I'd suggest doing some quick reading of the literature on epsilon-selection rules to understand how people reason about this space and the trade-offs involved.

As an example of how it's useful to take a broader view here, using eps() would often not work well because of rounding errors. The Wikipedia section is good on this: https://en.wikipedia.org/wiki/Numerical_differentiation#Practical_considerations_using_floating-point_arithmetic

The rules we use are stated here: https://scicomp.stackexchange.com/a/14357 The link to Numerical Recipes will point you to a longer discussion of the topic.