What's the accuracy difference between the two versions derivative calculation?

Question

The following is Python code aiming to calculate the derivative of a given function f.

1. Version one (Solution)

   x[ix] += h # increment by h
   fxh = f(x) # evalute f(x + h)
   x[ix] -= 2 * h 
   fxnh = f(x)
   x[ix] += h
   numgrad = (fxh - fxnh) / 2 / h
   

2. Version two (my version)

   fx = f(x) # evalute f(x)
   x[ix] += h 
   fxh = f(x) # evalute f(x+h)
   x[ix] -= h
   numgrad = (fxh - fx) / h
   

It has shown version one gives a better accuracy, could anyone explain why it is the case, what's the difference between the two calculations?

UPDATES I didn't realize it's a mathematical problem at the first place, I thought it was a problem relating to effects of floating accuracy. As suggested by MSeifert, I do agree that float-point noise matters, small magnitude result is more susceptible when exposing to the noise.

Answer 1

This is not a Python problem but a pure algorythmic one. Assuming that the function f has nice properties, you can look at its Taylor series development:

f(x+h) = f(x) + h f'(x) + h*h/2 f"(x) + h*h*h/6 f'''(x) + o(h3)

It comes that your first form gives for the error:

((f(x+h) - f(x)) /  h) - f'(x) = h/2 f"(x) + o(h)

that is an error in the magnitude order of h

If you use the second form, you get:

((f(x+h) - f(x-h)) / 2*h) - f'(x) = h*h/3 f'''(x) + o(h2)

the terms in h have fell and the error is in the magnitude order of h²

Of course it only makes sense if the required derivatives exist...

Answer 2

Your solution is "one-sided", you compare f(x+h) - f(x) and the general solution is "two-sided" f(x+h) - f(x-h).

It would be good to know what

It has shown version one gives a better accuracy

means. Because that's far too general.

But I think I have an example which might be appropriate here:

def double_sided_derivative(f, x, h):
    x_l, x_h = x - h, x + h
    return (f(x_h) - f(x_l)) / 2 / h

def one_sided_derivative(f, x, h):
    x_h = x + h
    return (f(x_h) - f(x)) / h

h = 1e-8

def f(x):
    return 1e-6 * x

# difference to real derivate:
double_sided_derivative(f, 10, h) - 1e-6, one_sided_derivative(f, 10, h) - 1e-6
# (6.715496481486314e-14, 1.5185825954029317e-13)

Note that the double-sided result is closer to the expected value. This might even lead to catastrophic cancellation. Then you may end up with a result that is mostly governed by floating-point-noise. This effect is further enhanced because the value is divided by a really small number.

By using both sides you increase (depending on your function!) the difference and thus the point at which the cancellation may occur. But in my opinion the biggest advantage is that you take the slope at both sides into account (averaging it somewhat). For example:

h = 1e-6

def f(x):
    return 4 + x + 5 * x**2

def fdx(x):
    return 1 + 10 * x

x = 10

double_sided_derivative(f, x, h) - fdx(x), one_sided_derivative(f, x, h) - fdx(x)
# (-2.7626811061054468e-08, 4.974594048690051e-06)

This is much closer to the true value (two orders of magnitude) than the one-sided approximation.