Numerical differentiation

Today I needed to the derivative of the zeta function. SciPy implements the zeta function, but not its derivative, so I needed to write my own version.

The most obvious way to approximate a derivative would be to simply stick a small step size into the definition of derivative:

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

However, we could do much better using

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

To see why, expand f(x) in a power series:

f(x + h) = f(x) + h f'(x) + h²f"(x)/2 + O(h³)

A little rearrangement shows that the error in the one-sided difference, the first approximation above, is O(h). Now if you replace h with -h and do a little algebra you can also show that the two-sided difference is O(h²). When h is small, h² is very small, so the two-sided version will be more accurate for sufficiently small h.

So how small should h be? The smaller the better, in theory. In computer arithmetic, you lose precision whenever you subtract two nearly equal numbers. The more bits two numbers share, the more bits of precision you may lose in the subtraction. In my application, h = 10^-5 works well: the precision after the subtraction in the numerator is comparable to the precision of the (two-sided) finite difference approximation. The following code was adequate for my purposes.

 from scipy.special import zeta def zeta_prime(x): h = 1e-5 return (zeta(x+h,1) - zeta(x-h,1))/(2*h)

The zeta function in SciPy is Hurwitz zeta function, a generalization of the Riemann zeta function. Setting the second argument to 1 gives the Riemann zeta function.

There's a variation on the method above that works for real-valued functions that extend to a complex analytic function. In that case you can use the complex step differentiation trick to use

Im( f(x+ih)/h )

to approximate the derivative. It amounts to the two-sided finite difference above, except you don't need to have a computer carry out the subtraction, and so you save some precision. Why's that? When x is real, x + ih and x - ih are complex conjugates, and f(x - ih) is the conjugate of f(x + ih), i.e. conjugation and function application commute in this setting. So (f(x+ih) - f(x-ih)) is twice the imaginary part of f(x + ih).

SciPy implements complex versions many special functions, but unfortunately not the zeta function.