How can I obtain a one dimensional finite difference formula for $U_{xx}$ with unevenly spaced nodes?

Question

I know that if I had evenly spaced points, I can use $U_{xx}\approx \frac{U_{i-1}-2U_{i}+U_{i+1}}{dx^2}$.

But if my gridpoints are unevenly spaced, I assume that I can obtain the finite difference formula in the following way:

$D^-[D^+(U)]$, where $D^-$ and $D^+$ are the backward and forward difference formulas. If I define $h_i=x_{i}-x_{i-1}$, then I think I would obtain:

$D^-[D^+(U)] = D^-[\frac{U_{i+1}-U_{i}}{h_{i+1}}] = \frac{\frac{U_{i+1}-U_{i}}{h_{i+1}} - \frac{U_{i}-U_{i-1}}{h_{i+1}}}{h_i}=\frac{U_{i-1}-2U_{i}+U_{i+1}}{h_ih_{i+1}}$.

However, I have some doubts about this formula. Particularly, the term $\frac{U_{i}-U_{i-1}}{h_{i+1}}$ might be wrong. Any feedback on this would be greatly appreciated :)

Answer 1

Since the finite difference approximation to the gradient is $$ u'(x) \approx \frac{\Delta u}{\Delta x} $$ it is clear that in your original formula, it needs to be $$ u''(x) \approx \frac{\frac{U_{i+1}-U_i}{x_{i+1}-x_{i}} -\frac{U_{i}-U_{i-1}}{x_{i}-x_{i-1}}} {\frac 12 (x_{i+1}+x_i) - \frac 12 (x_{i}+x_{i-1})} $$ which in your notation is then equal to $$ u''(x) \approx \frac{\frac{U_{i+1}-U_i}{h_{i+1}} -\frac{U_{i}-U_{i-1}}{h_{i}}} {\frac 12 (h_{i+1}+h_i)} = \frac{(U_{i+1}-U_i)h_i -(U_{i}-U_{i-1}){h_{i+1}}} {\frac 12 (h_{i+1}+h_i)h_ih_{i+1}}. $$ Conveniently, and not by accident, this reduces to the well-known formula for uniform mesh sizes.

Answer 2

This paper, Generation of Finite Difference Formulas on Arbitrarily Spaced Grids, by Bengt Fornberg, provides Pseudo-Code for the... well... the title is a bit of a give-away.

If you pack it into a loop, you can create Differentiation matrices of arbitrary degrees over arbitrary point sets. If you really need it, I've got a Matlab implementation thereof that creates said matrices, which I could provide on request.

Update

Alexander has just pointed-out fdcoeffF(k,xbar,x), an implementation of Fornberg's algorithm which computes the weights of the stencil for the kth derivative at a point xbar over a set of data at points x.