This implementation could certainly use some love, but finite difference differentiation is always unstable, and the situation gets worse as you take higher order derivatives.
You might want to consider using the differentiate method to take your derivatives symbolically (if this works for the functions you're using), or have a look at the differentiation example in PowerSeries.jl. To do a symbolic second derivative, you can do, e.g. julia> using Calculus julia> @eval d2(x) = $(differentiate(differentiate(:(sin(x))))) julia> d2(1.0) -0.8414709848078965 On Monday, January 20, 2014 10:40:28 AM UTC-8, Hans W Borchers wrote: > > I looked into the *Calculus* package and its derivative functions. First, > I got errors when running examples from the README file: > > julia> second_derivative(x -> sin(x), pi) > ERROR: no method eps(DataType,) > in finite_difference at > /Users/HwB/.julia/Calculus/src/finite_difference.jl:27 > in second_derivative at > /Users/HwB/.julia/Calculus/src/derivative.jl:67 > > Then I was a bit astonished to see not too accurate results such as > > julia> abs(second_derivative(sin, 1.0) + sin(1.0)) > 6.647716624952338e-7 > > while, when applying the standard central formula for second derivatives, > (f(x+h) - 2*f(x) + f(x-h)) / h^2 with the (by theory) suggested step > length eps^0.25 (for second derivatives) will result in a much better > value: > > julia> h = eps()^0.25; > > julia> f = sin; x = 1.0; > > julia> df = (sin(x+h) - 2*sin(x) + sin(x-h)) / h^2 > -0.8414709866046906 > > julia> abs(df + sin(1.0)) > 1.7967940468821553e-9 > > The functions for numerical differentiation in *Calculus* look quite > involved, maybe it would be preferable to apply known approaches derived > from Taylor series. Even the fourth order derivative will in this case lead > to an absolute error below 1e-05! >
