Friday, April 08, 2005

Computing derivatives using Python

In a recent thread on Python edu-sig, Kirby Urner started a discussion about the use of python's decorators in mathematics. Some of the examples included numerically computing derivatives (and integrals) of functions. This generated a fair bit of discussion and the general conclusion was that this was not a good use of decorators. Some information contributed by yours truly required the reader to do a graph to properly visualise exactly what was going on. I thought it appropriate to present a summary that included the required graphical information here as a permanent record.

The mathematical definition of the derivative of a function f(x) at point x is to take a limit as "h" goes to zero of the following expression:

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

An example is given below: the red curve is the function f(x) and the green curve is the tangent line at x (with same slope as the derivative). The blue curve is a line going through x whose slope equals that of the above formula with a non-zero value of h.

As you can see, the slopes of the two straight lines are quite different. A better way of numerically computing the derivative for the same value of "h" is by taking a symmetrical interval around x as follows:

( f(x + h/2) - f(x - h/2) ) / h

This is illustrated in the next picture. As you can see, the two straight lines have much more similar slopes - hence, the value computed for the derivative is going to be more accurate.

The corresponding python code is as follows:

def derivative(f):
....Computes the numerical derivative of a function.
....def df(x, h=0.1e-5):
........return ( f(x+h/2) - f(x-h/2) )/h
....return df

And we use it as follows:

# sample function
def g(x): return x*x*x

# first derivative
dg = derivative(g)

# second derivative
d2g = derivative(dg) # == derivative(derivative(g))

# printing the value computed at a given point:

print dg(3)
print dg(3, 0.001)
print dg(3, 1e-10) # smaller h is not always more precise


Petri said...

Nice example. Care to add the name of the numerical method used?

André said...

I do not know if this method has a name. It is something I picked up more than 15 years ago. I tried to look it up in a few books but didn't find it anywhere.

Anonymous said...

Approximations of this form are called "finite differences", of which there are many variations. The particular difference formula your employing is often called a "second-order centered difference", while the one your contrasting against is a "first-order forward difference." The "second-order" and "first-order" refer to the accuracy. Second-order means that the error in the approximation decreases with decreasing h proportional h*h, while first-order decreases only proportional to h.

Anonymous said...

Directly from the time machine:

Matthew Strax-Haber said...

Check this out (possibly a bit clearer for debugging):

def D(f,h=1e-5):
'''Return derivative of function f'''
def df(x):
return (f(x+h)-f(x))/h
df.__name__ = f.__name__ + '_dx'
return df

Igor said...

You should be careful when using this formula since the limit of this formula might exist while the limit that defines the derivative might not. They are not the same thing although the "second-order centred difference" formula provides better numerical approximation.