Numerical Differentiation Goes Wrong

The mathematical definition of differentiation is simple. It is the slope of a graph, or, algebraicly, the limit as h tends to zero of

[f(x+h)-f(x)]/h

The numerical problem is how to choose h. Maths says that it should be as close to zero as possible, but computers calculate to a finite precision, so there will come a point when x and x+h are indistinguishable, and so f(x+h)-f(x) evaluates to zero.

This is not the only problem. The values of f(x+h) and f(x) will be very similar if h is small. Subtracting two small numbers leads to a loss of precision. If both are evaluated to ten significant figures, but the first six digits are identical, then their difference is accurate to just four digits.

A simple piece of python shows the problem. Here we try to find numerically the derivative of x*exp(-x) at x=0.5.

plot of x*exp(-x)

#!/usr/bin/python
from math import exp

def func(x):
    return(x*exp(-x))

h=0.1
x=0.5

for i in range(0,16):
    print("%e %.9f"%(h,(func(x+h)-func(x))/h))
    h=h/10

The output is

1.000000e-01 0.260216518
1.000000e-02 0.298741534
1.000000e-03 0.302810684
1.000000e-04 0.303219843
1.000000e-05 0.303260781
1.000000e-06 0.303264875
1.000000e-07 0.303265284
1.000000e-08 0.303265330
1.000000e-09 0.303265357
1.000000e-10 0.303265746
1.000000e-11 0.303268521
1.000000e-12 0.303257419
1.000000e-13 0.303645997
1.000000e-14 0.299760217
1.000000e-15 0.333066907
1.000000e-16 0.000000000

So for large h the answer is not very good as the step length along the function, which is gently curving, is too great. We need h to be less than about 1e-5 to get an answer correct to just five figures. But once h falls below 1e-12 then again the answer is no longer correct to five figures.

We can differentiate the expression analytically.

f(x)=x*exp(-x)

f'(x)=-x*exp(-x)+exp(-x)

f'(x)=(1-x)*exp(-x)

So at x=0.5 the correct answer is 0.303265329856...

The answer for h=1e-8 was correct to nine significant digits, but its two neighbours were not.

Move just a little further along the same function, to x=0.9999, and things get more difficult.

1.000000e-01 -0.017179201
1.000000e-02 -0.001790755
1.000000e-03 -0.000147062
1.000000e-04 0.000018395
1.000000e-05 0.000034952
1.000000e-06 0.000036608
1.000000e-07 0.000036774
1.000000e-08 0.000036793
1.000000e-09 0.000036804
1.000000e-10 0.000036637
1.000000e-11 0.000038858
1.000000e-12 0.000055511
1.000000e-13 0.000000000
1.000000e-14 0.000000000
1.000000e-15 0.055511151
1.000000e-16 0.555111512

The best answer is still at h=1e-8, and the correct answer starts 0.0000367916, so the best answer does not quite manage five significant digits of accuracy.

Where does the magic sequence 55511151... come from?

The value of the function x*exp(-x) is around 0.37 when x is close to one. So a computer, working in base two, will store it as m*2-1, so that m lies between 0.5 and 1. The value of m, around 0.74, is stored to 53 binary digits of precision. So for a given m, the next distinct number is m+2-53. If two numbers are as close as they can be without being regarded as identical, m will differ by 2-53. Here, with a binary exponent of -1, the numbers themselves will differ by 2-54. What is 2-54? It is 0.55511151231e-16, and, in the table above, appears divided by h when h is 1e-12, 1e-15 and 1e-16.

But sometimes it is possible to perform differentiation without subtraction.