# Efficient Arbitrary Order Finite Differences in 1D

I am implementing on Matlab a high-order finite differences scheme to approximate the first derivative of $$f(x_i)$$ given $$x = [x(1), x(2),..., x(i),..., x(n)]$$ and $$f = [f(x(1)),..,f(x(n))]$$ with $$x$$ a non-uniform grid from $$x(1)$$ to $$x(n)$$ (in fact they are Legendre-Gauss-Lobatto nodes in $$[-1,1]$$). Right now I use the method detailed in chapter 2.1 at https://en.wikibooks.org/wiki/Using_High_Order_Finite_Differences/Definitions_and_Basics#approximation_of_a_first_derivative, which uses the Vandermonde matrix.

I decided to first set $$h = x(i+1)-x(i)$$ or $$x(i)-x(i-1)$$ depending on $$i$$, then I compute the Vandermonde matrix $$M$$ using $$\alpha(j) = (x(j)-x(i))/h$$ and I get the $$a_j$$'s coefficients with $$M\backslash [0,1,0,..,0]$$'. And finally $$f'(x_i) ~ (1/h)*(a_1*f(x_1) + ... + a_n*f(x_n))$$. It seems like this method returns the right result, unfortunately it makes the computations really slow...

Do you know another method to compute arbitrary order finite differences schemes ? Or a faster way to compute the $$a_j$$'s coefficients ?

I hope my explanations are clear enough, and already thank you for your answers.

• Since your grid is fixed you need to compute the coefficients of your finite-difference scheme only once, before the actual computation, so it would not slow down anything. Apr 18 '20 at 20:39
• Can you use MathJax for your equations? Apr 18 '20 at 23:06
• For many node distributions, the differentiation matrix can be computed easily. See Section 3.5 in the book Kopriva: Implementing Spectral Methods. Apr 20 '20 at 6:35

One of the key takeaway's in polynomial approximation is to avoid methods that involve Vandermonde matrices. This should also apply to the method you cited, which is basically doing a polynomial interpolation to the first derivative of a given function in a monomial basis.

One way to get around Vandermonde matrices is to use polynomial basis sets other than the monomial basis, and often the basis of Lagrange polynomials is the method of choice. A popular method in this regard has been derived by Fornberg, see e.g. here or here. I would recommend you to use this method instead.

Moreover, for Gauss-Lobatto nodes, the derivative matrix is known exactly. As shown in the second link above, however, this does not mean that the result has to be more accurate than using Fornberg's method (together with the diagonal correction also mentioned in the linked post).