Numerical derivative and finite difference coefficients: any update of the Fornberg method?

Question

When one want to compute numerical derivatives, the method presented by Bengt Fornberg here (and reported here) is very convenient (both precise and simple to implement). As the original paper date from 1988, I would like to know whether there is a better alternative today (as (or almost as) simple and more precise) ?

Answer 1

Overview

Good question. There is a paper entitled "Improving the accuracy of the matrix differentiation method for arbitrary collocation points" by R. Baltensperger. It's no big deal in my opinion, but it has a point (that already was known before the appearance in 2000): it stresses the importance of an accurate representation of the fact that the derivative of the constant function $f(x)=1$ should be zero (this holds exactly in the mathematical sense, but not necessarily in the numerical representation).

It is simple to see that this requires the row sums of the n-th derivative matrices $D^{(n)}$ to be zero. It is common to enforce this constraint by adjusting the diagonal entry, i.e. by setting $$ \tag 1 D^{(n)}_{jj} := -\sum_{\substack{i=1\\i \neq j }}^N D_{ij}\,. $$ It's clear that this feature does not hold exactly when working on a computer due to roundoff errors in floating point calculations. What is more surprising is that these errors are even more severe when using the analytical formulas for the derivative matrix (which are available for many classical collocation points, e.g. Gauss-Lobatto).

Now, the paper (and references therein) states that the error of the derivative is in the order of the deviation of the row sums from zero. The goal is therefore to make these numerically as small as possible.

Numerical tests

The good point is that the Fornberg procedure seems to be quite good in this regard. In the picture below I've compared the behaviour of the exact, i.e. analytical, first derivative matrix and the one derived by the Fornberg algorithm, for varying number of Chebyshev-Lobatto collocation points.

Again, believing the statement in the cited paper, this implies that the Fornberg algorithm will yield more accurate results for the derivative.

In order to prove that, I'll use the same function as in the paper, $$\tag 2 f(x) = \frac{1}{1+x^2}\,.$$ and evaluate the error $$\tag 3 E_n = \max_{i\in \{0,\ldots,n\}} \bigg| f^\prime(x_i) - \sum_{j=1}^n D_{ij} f(x_j) \bigg|\,.$$ This is done for (i) the analytically obtained derivative matrix, (ii) the Fornberg derivative matrix, and (iii) an adjusted version of the Fornberg matrix, where the above Eq. (1) is enforced quite straightforwardly via $$\tag 4 \tilde D_{jj} = D_{jj} - \left(\sum_{i=1}^n D_{ji} \right), \qquad \text{for all}\ j\,. $$ Here is what I get (again for the example of Gauss-Lobatto abscissas):

Conclusion

In conclusion, Fornberg's method seems to be quite accurate, in the case $N=512$ even by about 3 orders of magnitude more accurate than results from the analytical formulas. This should be accurate enough for most applications. Moreover, this is remarkable because Fornberg does not seem to explicitly include this fact in his method (at least there is no mention in the two Fornberg papers).

Another order of magnitude can be gained for this example through a straightforward inclusion of Eq.(4). As this is a quite simple approach and applied only once for each derivative, I see no reason in not using it.

The method from the Baltensperger paper -- which uses a more sophisticated approach for evaluating the sum in Eq.(1) in order to reduce roundoff errors -- yields about the same order of magnitude for the error. So, at least for this example, it's roughly equivalent to the "Adjusted Fornberg" method above.

Answer 2

Assuming you are trying to differentiate a numerical implementation of a continuous function, there are a large number of methods:

1) Automatic differentiation. The most accurate and general method. Painful to code, requiring operator overloading and argument dependent lookup. Puts a burden on the user to understand these concepts. Also struggles with removable singularities, such as differentiating sinc at $x= 0$.

2) A Chebyshev transform. Project your function onto a span of Chebyshev polynomials and differentiate the three term recurrence. Super fast, very accurate. But requires that you have a compactly supported domain of interest; outside the selected domain $[a,b]$, the three term recurrence is unstable.

3) Finite differencing. Underrated in 1D; see Nick Higham's Tips and Tricks in Numerical Computing. The idea is that if you balance the truncation error and the roundoff error then you don't need to select a stepsize; it can be chosen automatically. In Boost, this idea is used to recover (by default) 6/7th of the correct digits for the type. (Higham only shows the idea for the simpler case of 1/2 the correct digits, but the idea is easily extended.) The coefficients are from Fornberg's equispaced table, but the stepsize is chosen under the assumption that the function can be evaluated to 1ULP accuracy. The disadvantage is that it requires 2 function evaluations to recover half the digits of the type, 4 to recover 3/4th the digits, so on. In 1D, not a bad deal. In higher dimensions, it's catastrophic.

4) The complex step derivative. Use $f'(x) \approx \Im(f(x+ih))$. Take $h$ to be the unit roundoff and this will recover almost every bit correct. However, it's kinda cheating, because it's generally harder to implement a function in the complex plane than to hand code its real derivative. Still a cool idea and useful in certain circumstances.

Answer 3

I'm not aware of anyone having improved Fornberg's algorithm (see also his slightly more recent paper). As an aside, it seems to me that looking at his algorithm as a way to compute numerical derivatives is not right. All he's done is derive an efficient algorithm to compute weights for finite-difference methods. The advantage of his method is that it gives you the weights for all derivatives up to the desired derivative in one go.

Answer 4

A simpler scheme

In addition to my other answer which is more about an extension of the Fornberg method, I'll address here the question for more simple alternatives.

For this I sketch an alternative scheme which produces the derivative coefficients the of Lagrangian interpolation more directly. Its implementation requires only a few lines of code, works for arbitrary grids, and according to my first experiments, is as accurate as Fornberg's.

Basis of the implementation is the imaginary-step derivative $$ f'(x) \ = \ \frac{1}{h}\,\text{Im}\big(f(x+ih)\big)\,, $$ where $h$ is a variable in the order of machine-precision. The imaginary-step derivative is known to stably produce derivative values and does not suffer the numerical instability of a Finite difference implementation with $h\rightarrow 0$.

Second ingredient is the Lagrange interpolation polynomial $L_i(x)$ on the grid $\{x_1,\ldots, x_N\}$ evaluated with one of the barycentric forms, e.g. $$ L_i(z) \ = \ \begin{cases}1 & \text{if }z = x_i\\\ \frac{\frac {\mu_i}{z-x_k}}{ \sum_k\frac {\mu_k}{z-x_k} } & \text{otherwise}\,. \end{cases}$$ where $\mu_i = \frac{1}{\prod_{k\neq i} (x_i - x_k)}$ In order to use the complex-step derivative, one has to make sure these formulas also work for complex arguments $z$. Moreover, for a given function f(x) and vector of coefficients $f_i = f(x_i)$, we denote the interpolation polynomial through $(x_1,f_i)$ by $$ P(x;f) = \sum_{i=1}^N f_i L_i(x)\,. $$

---

Algorithm

The algorithm is sketched in the following. It has the same input and output parameters as the one of Fornberg, but is much more a no-brainer.

Input:

• $x$ : a grid with N distinct gridpoints
• $ord$ : derivative order
• z : a point at which the derivative is to be evaluated
• Possibly: a function $f(x)$ or its funciton values $f_i$ at the gridpoints (only required for output variant 2.)

Initialization

• Initialize the Lagrange polynomials $L$ via barycentric interpolation.
• Initialize an array of $N\times N$-matrices $D^{(o)}$, for $o=0, \ldots, ord$
• Set $D^{(0)} := I_N$, i.e. the unit matrix of dimension $N\times N$
• Set $o := 0$

Algorithm

While $o < ord$:

• Calculate $D^{(o+1)}_{ik} = CSD \bigg( P\big(x_k; D^{(o)}_{i}\big) \bigg)$ by the complex-step derivative $CSD$ for all $i$ and $k$.
  Here, $D^{(o)}_{i}$ denotes the $i$-th row of $D^{(o)}$.

• Set o = o+1;

Decide what to Output:

1. The vector $d^{(ord)}$ of Finite Difference coefficients at point $z$, where $d_i = P(z;D^{(ord)}_i)$. This is what Fornberg does.

2. The interpolation function $p^{(ord)}(x) = \sum_{i=1}^N f(x_i) P(x;D^{(ord)}_i)$ to the derivate $f^{(ord)}(x)$ of order $ord$. For this, you have to input a function resp. the function values $f_i$ at $x_i$ to the algorithm.

3. A meta-function $\text{diff}(\text{int } ord, \text{function } g)$ which returns the interpolation function of variant 2., but for an arbitrary function $g$ that is to be interpolated at the gridpoints.

Personally, I like variant 3. the most.

---

Analysis of the algorithm

As Fornberg's, this algorithm is $\mathcal O(ord \cdot N^2)$. I'll post more empirical results regarding accuracy, stability, etc. if I find the time.

Answer 5

To increase the precision of numerical differentiation do the following:

1) Chose your favorite high-precision "standard" method based on some step size h.

2) Compute the value of the derivative with the method chosen in 1) many times with different but reasonable step sizes h. Each time you may pick h as a random number from the interval (0.5*H/10, 1.5*H/10) where H is an appropriate step size for the method you use.

3) Average the results.

Your result may gain 2-3 orders of magnitude in the absolute error wrt. the non-averaged result.

https://arxiv.org/abs/1706.10219