Why do different step sizes produce large numerical error while approximating derivatives?

user_14974457 0
Viewed 119

I'm attempting a bit of code that would help me approximate the derivatives of an arbitrary function. I saw there were four options on another post:

  1. Finite Differences
  2. Automatic Derivatives
  3. Symbolic Differentiation
  4. Compute derivatives by hand

I saw that my approach falls best in the first option, which had the note, "prone to numerical error". So I'm aware that this method isn't expected to be exact, which is fine.

That being said, I did some research into what size numbers can be stored by different data types, and found in this post that it can be quite small, (on the order of 10–308) and that "In the normal range, results of elementary operations will be accurate within the normal precision of the format".

That being said, I seem to be getting extremely bad outcomes for this following code snippet where I explore different sized intervals; the smallest difference shouldn't be much smaller than 10–27 (10–9, cubed), which is much larger than the limiting value. I would appreciate maybe a more specific response?

epsilon = 0.01 # is "small" w.r.t. to 3

def approx_derivative(func): # rough derivative factory function
  return lambda x : (func(x + epsilon) - func(x)) / epsilon

while epsilon > 10**-9:
  nth_deriv = lambda x : x ** 3 # 0th derivative
  for i in range(5): # should read about 27, 27, 18, 6, 0
    print(nth_deriv(3), end=', ')
    nth_deriv = approx_derivative(nth_deriv) # take derivative
  print('\n')
  epsilon *= 0.1

The output is:

27, 27.090099999999495, 18.0599999999842, 6.000000002615025, -3.552713678800501e-07, 

27, 27.009000999996147, 18.00600000123609, 6.000000496442226, -0.007105427357601002, 

27, 27.00090001006572, 18.000599766310188, 6.004086117172847, -71.05427357601002, 

27, 27.000090000228735, 18.000072543600254, 3.5527136788005005, 355271.36788005003, 

27, 27.000009005462285, 17.998047496803334, 0.0, 3552713678.8005, 

27, 27.000000848431675, 18.11883976188255, 0.0, -35527136788004.99, 

27, 27.0000001023618, 0.0, 0.0, 3.552713678800497e+17, 

27, 27.000002233990003, 0.0, 0.0, 0.0,

As we can see in the first couple of examples, the results aren't exact but are pretty good. For certain interval sizes, though, some values are blown up; others go to 0; and some are just plain wrong, like giving half the value, despite the intuition that they should become more accurate for smaller epsilons. What main things can I attribute to this error? What should I be looking out for/be cautious of? Are there errors I should be worried about catching with a block (like division by 0)?

Is there a value for epsilon that is generally considered "best" for doing computations with floats? Or is there a "rule-of-thumb" for choosing a good-sized epsilon based on your input? Is there a preferred definition of the derivative to use over the one I implemented?

derivative numeric ieee-754 precision python
1 Answers

First, the minimum value representable in the floating-point type is irrelevant here. The precision is a concern. The Python specification is not specific about what floating-point format a Python implementation uses, but many use the IEEE-754 binary64 format, which has 53 bits in its significand, meaning the smallest representable change in value is about 2−52 times the magnitude of the number represented. We will look at how that affects your calculations.

Computing the approximation of the derivative, f', computes (x+e)3−x3. The results is 3ex2 + 3e2x + e3, approximately 3ex2, which is then divided by e to yield approximately 3x2. Note that two terms around x3 in magnitude were involved, and they differed by about 3ex2. Since you are taking the derivative at x = 3, 3ex2 happens to equal ex3, so this is lower than x3 by a factor of e. As long as e is modest, two binary64 numbers near x3 can readily differ by ex3 without losing accuracy.

Computing the approximation of the second derivative, f'', computes ((x+2e)3−(x+e)3)−((x+e)3−x3) = (x+2e)3−2(x+e)3+x3 = 6e2x + 6e3. Along the way, this gets divided by e2, but we can ignore that for error analysis for now. The original terms had magnitude around x3, but the result of our subtractions is around 6e2x. So this relies on the numbers being able to differ by about e2. Since 2−52 is about 2•10−16, this fails considerably when e is around 10−8, with lesser errors before that. When e is around this size, the binary64 format simply cannot record the differences between the numbers precisely enough to retain accuracy when they are later subtracted.

For the third derivative, you need to represent differences around e3 relative to the primary magnitudes, so it fails considerably around e = 10−5, since then e3 is getting near the limit of 2•10−16.

user_2982250
Related