For handling experimental data from Alcator C-Mod (e.g., density or temperature profiles), I commonly need to take a standard numerical derivative of the profile. However, since the abscissa of these data is commonly irregular (dictated by hardware constraints) the typical numerical solutions provided in Python are insufficient. The standard options I’ve encountered:
do not allow for arbitrary \(x\)- and \(y\)-axes.
I have implemented a simple solution to this, using Lagrange interpolating polynomials to generate the finite-difference matrix).
In general, \(n\) data points may be interpolated by a Lagrange polynomial of degree \(\le n-1\), with the polynomial passing through each point \((x_0,y_0)...(x_n,y_n)\). These polynomials are defined for the general case by
\[P(x) = \sum_{j=0}^{n} P_j(x)\]where
\[P_j(x) = y_j \prod_{m=0,m \ne j}^{n} \frac{x-x_m}{x_j-x_m}\]Intuitively, we may use this to generate a two-sided finite difference method - for each point, we use a \(n=2\) Lagrange polynomial for the points on either side to approximate a locally-parabolic curve through each point, and generate the difference matrix based on this polynomial. For a given point \(k\), this is
\[f_k(x) = y_{k-1} P_{k-1}(x) + y_k P_k(x) + y_{k+1} P_{k+1}(x)\] \[P_{k-1}(x) = \frac{(x-x_k)(x-x_{k+1})}{(x_{k-1}-x_{k})(x_{k-1}-x_{k+1})}\] \[P_k(x) = \frac{(x-x_{k-1})(x-x_{k+1})}{(x_k-x_{k-1})(x_k-x_{k+1})}\] \[P_{k+1}(x) = \frac{(x-x_{k-1})(x-x_{k})}{(x_{k+1}-x_{k-1})(x_{k+1}-x_{k})}\]The local derivative at \(x_k\), then, is
\[f'(x_k) = y_{k-1} P'_{k-1}(x_k) + y_k P'_k(x_k) + y_{k+1} P'_{k+1}(x_k)\]with each \(P'(x_k)\) factors providing the matrix elements for the difference matrix:
\[P'_{k-1}(x_k) = a_k = \frac{x_{k} - x_{k+1}}{(x_{k-1}-x_{k})(x_{k-1}-x_{k+1})}\] \[P'_k(x_k) = b_k = \frac{-x_{k+1} + 2x_k - x_{k-1}}{(x_k-x_{k-1})(x_k-x_{k+1})}\] \[P'_{k+1}(x_k) = c_k = \frac{x_k - x_{k-1}}{(x_{k+1}-x_{k-1})(x_{k+1}-x_{k})}\]Or, we may simplify this by defining the spacing between points, \(h_k = x_k - x_{k-1}\):
\[a_k = \frac{-h_{k+1}}{h_k(h_{k+1} + h_k)}\] \[b_k = \frac{h_{k+1} - h_k}{h_{k+1} h_k}\] \[c_k = \frac{h_k}{h_{k+1}(h_{k+1} + h_k)}\]It is evident that, in the special case of a uniform grid (all \(h_k\) equal) these reduce to the standard centered single difference matrix elements.
At the boundaries, we construct a single-sided difference using similar quadratic polynomials,
\[f_0(x) = y_0 P_0(x) + y_1 P_1(x) + y_2 P_2(x)\] \[f_n(x) = y_{n-2} P_{n-2} (x) + y_{n-1} P_{n-1}(x) + y_n P_n(x)\]resulting in matrix elements
\[a_0 = - \frac{2h_1 + h_2}{h_1(h_1 + h_2)}\] \[b_0 = \frac{h_1 + h_2}{h_1 h_2}\] \[c_0 = - \frac{h_1}{h_2(h_1 + h_2)}\]and
\[a_n = \frac{h_n}{h_{n-1}(h_{n-1} + h_n)}\] \[b_n = - \frac{h_{n-1} + h_n}{h_{n}h_{n-1}}\] \[c_n = \frac{2h_n + h_{n-1}}{h_n (h_{n-1} + h_n)}\]Using these elements to construct the difference matrix, I’ve implemented a finite-difference method in Python, compatible with typical numpy and scipy workflows, allowing numerical differentiation on an arbitrary grid.
This difference-matrix method is based on work developed on scientificpython.net (with some corrections and modifications). For finite-difference methods, see section 1.2 of Strang’s Computational Science and Engineering.