Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

problems with finite difference #132

Open
JianghuiDu opened this issue Apr 19, 2019 · 1 comment
Open

problems with finite difference #132

JianghuiDu opened this issue Apr 19, 2019 · 1 comment

Comments

@JianghuiDu
Copy link

JianghuiDu commented Apr 19, 2019

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
Copy link
Collaborator

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

2 participants