How can I numerically differentiate an unevenly sampled function?

Standard finite difference formulas are usable to numerically compute a derivative under the expectation that you have function values $f(x_k)$ at evenly spaced points, so that $h \equiv x_{k+1} - x_k$ is a constant. What if I have unevenly spaced points, so that $h$ now varies from one pair of adjacent points to the next? Obviously I can still compute a first derivative as $f'(x) \approx \frac{1}{h_k}[f(x_{k+1}) - f(x_k)]$, but are there numerical differentiation formulas at higher orders and accuracies that can adapt to variation in the grid size?

-
You can always construct a (piecewise) polynomial interpolant passing through your points, and then differentiate that. –  Ｊ. Ｍ. Dec 19 '11 at 9:00
Or, you can reconstruct the finite difference formulas without the simplification $h = x_{k+1} - x_k$. Often times this must be done for integration, but it is likely that J.M.'s suggestion is more stable. –  rcollyer Dec 19 '11 at 11:11
What kind of function is that? –  mbq Dec 19 '11 at 12:33
The example that prompted this question is a function sampled at logarithmically spaced values $x_k = x_0 \delta^k$, but calculating the second derivative of the log-transformed data gives funny results and I wanted a check on it. Plus I figured I'd ask as general a question as possible. –  David Z Dec 20 '11 at 20:46
As far as I'm concerned, something that works for only first and second derivatives would be a perfectly fine answer to the question. I wrote the question as I did to allow for a general answer if someone had one, but of course in practice it's the first and second derivs that are most useful. –  David Z Dec 21 '11 at 0:11

J.M's comment is right: you can find an interpolating polynomial and differentiate it. There are other ways of deriving such formulas; typically, they all lead to solving a van der Monde system for the coefficients. This approach is problematic when the finite difference stencil includes a large number of points, because the Vandermonde matrices become ill-conditioned. A more numerically stable approach was devised by Fornberg, and is explained more clearly and generally in a second paper of his.

Here is a simple MATLAB script that implements Fornberg's method to compute the coefficients of a finite difference approximation for any order derivative with any set of points. For a nice explanation, see Chapter 1 of LeVeque's text on finite difference methods.

A bit more on FD formulas: Suppose you have a 1D grid. If you use the whole set of grid points to determine a set of FD formulas, the resulting method is equivalent to finding an interpolating polynomial through the whole grid and differentiating that. This approach is referred to as spectral collocation. Alternatively, for each grid point you could determine a FD formula using just a few neighboring points. This is what is done in traditional finite difference methods.

As mentioned in the comments below, using finite differences of very high order can lead to oscillations (the Runge phenomenon) if the points are not chosen carefully.

-
On the other hand, when you use interpolating polynomials, one must always remember things like Runge's phenomenon possibly happening with your data, if your data is perversely configured enough. I'd say piecewise polynomials might be less susceptible to this... –  Ｊ. Ｍ. Dec 19 '11 at 14:00
Recent work by Plamen Koev (building on work by Björck and Pereyra), however, suggests that Vandermonde systems can be solved with good accuracy even if the systems can be ill-conditioned. –  Ｊ. Ｍ. Dec 19 '11 at 14:28
I wonder if Koev's work and Fornberg's technique could be related? –  David Ketcheson Dec 19 '11 at 16:48
Interestingly enough, there seems to be a resemblance between Fornberg's formulae, and earlier formulae developed by Lyness and Moler based on the classical Neville method for generating the interpolating polynomial. They might in fact be the same formulae in different notation, but I haven't checked thoroughly. –  Ｊ. Ｍ. Dec 20 '11 at 6:15
Polynomial interpolation with many points requires special point distributions to be well-conditioned. In general, for non uniform point distributions it is not recommended to do interpolation and then differentiation the interpolation polynomial because it can be highly oscillatory (think "Runge phenomenon" as mentioned by J. M.). Depending on your needs, it might be a better idea to just use cubic splines which for many practical purposes can provide good answers to the approximation problem of approximating derivatives. –  Allan P. Engsig-Karup Dec 20 '11 at 14:56

http://mathformeremortals.wordpress.com/2013/01/12/a-numerical-second-derivative-from-three-points/

This addresses your question and shows the formula you are looking for, for the second derivative. Higher-Order derivatives follow the same pattern.

-

The simplest method is to use finite difference approximations.

A simple two-point estimation is to compute the slope of a nearby secant line through the points (x,f(x)) and (x+h,f(x+h)).[1] Choosing a small number h, h represents a small change in x, and it can be either positive or negative. The slope of this line is

$$f(x+h)-f(x)\over h$$

This expression is Newton's difference quotient.

The slope of this secant line differs from the slope of the tangent line by an amount that is approximately proportional to h. As h approaches zero, the slope of the secant line approaches the slope of the tangent line. Therefore, the true derivative of f at x is the limit of the value of the difference quotient as the secant lines get closer and closer to being a tangent line

-
I think you're getting downvoted because David Zaslavsky specifically mentioned the difference quotient formula, and the question is asking if there are any better approximations. –  Dan Dec 21 '11 at 6:54
Also because it's a direct copy-and-paste from Wikipedia, except for the spam link that was originally part of the answer. –  David Z Dec 21 '11 at 22:37