JuliaMath / JuliaMath/Calculus.jl

problems with finite difference

Open
#132 1 comment 0 reactions 0 assignees View on GitHub
Dominant language
Julia
Stars
285
Forks
78
PR merge metrics
No merged PRs in 30d

Description

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?

Contributor guide

No contributing guide indexed for this repository

Assessment

This issue has not been assessed yet.

Get new issues in your inbox

A short digest of beginner-friendly GitHub issues.