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?
$\begingroup$
$\endgroup$
7
-
9$\begingroup$ You can always construct a (piecewise) polynomial interpolant passing through your points, and then differentiate that. $\endgroup$– J. M.Dec 19, 2011 at 9:00
-
$\begingroup$ 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. $\endgroup$– rcollyerDec 19, 2011 at 11:11
-
$\begingroup$ What kind of function is that? $\endgroup$– mbqDec 19, 2011 at 12:33
-
$\begingroup$ 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. $\endgroup$– David ZDec 20, 2011 at 20:46
-
1$\begingroup$ 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. $\endgroup$– David ZDec 21, 2011 at 0:11
|
Show 2 more comments
6 Answers
$\begingroup$
$\endgroup$
10
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.
answered Dec 19, 2011 at 13:56
-
5$\begingroup$ 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... $\endgroup$– J. M.Dec 19, 2011 at 14:00
-
1$\begingroup$ I wonder if Koev's work and Fornberg's technique could be related? $\endgroup$ Dec 19, 2011 at 16:48
-
1$\begingroup$ 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. $\endgroup$– J. M.Dec 20, 2011 at 6:15
-
3$\begingroup$ 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. $\endgroup$ Dec 20, 2011 at 14:56
-
1$\begingroup$ Nice answer. Just for information, this paper gives an alternative approach to Fornberg's. It follows the same principle, but gives a different algorithm. $\endgroup$ Apr 20, 2017 at 1:08
$\begingroup$
$\endgroup$
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.
$\begingroup$
$\endgroup$
6
The above answers are great in terms of giving you a code to use, but aren't as good in terms of theory. If you want to delve deeper into interpolating polynomials, take a look at this theoretical treatment with a few concrete examples:
Singh, Ashok K., and B. S. Bhadauria. "Finite difference formulae for unequal sub-intervals using lagrange’s interpolation formula." International Journal of Mathematics and Analysis 3.17 (2009): 815-827. (Link to PDF)
The authors use Lagrangian Interpolation (see the Wikipedia article) to calculate 3-point, 4-point and 5-point interpolating polynomials, as well as their first, second and third derivatives. They have expressions for the truncation error as well, which is important to consider when using any finite difference scheme. They also have the generic formula for calculating interpolating polynomials using N points.
Lagrangian interpolating polynomials are useful because they and their derivatives can be very accurate in the domain you are interpolating, and they do not assume an even grid spacing. Due to the nature of Lagrangian interpolating polynomials, you can never have more orders of derivatives than you have grid points.
I think this answers your question well because the paper I cited has formulae for arbitrarily high-order finite difference schemes, which by nature are for uneven grids and are limited only by the number of grid points you include in your stencil. The paper also has a generic formula for the truncation error, which will help you evaluate the Lagrangian interpolating polynomial scheme against other schemes you might be considering. The author's paper should give the same results as Fornberg's method. Their contribution is really just tallying a few examples and giving an estimate of the error, which you may find useful.
I found both the paper I cited and Fornberg's work to be useful for my own research.
-
3$\begingroup$ sorry that I have to state it, but your cited reference looks weird -- they use horrible formulas and solve only a few special cases. In contrast, Fornberg has solved the general problem by giving a simple algorithm, and that already in the 80's. See here $\endgroup$ Mar 27, 2015 at 0:35
-
-
3$\begingroup$ and a last comment to disrespect this paper. In "an excellent theoretical treatment", you can't have 9 references, where 7 refer to you own work and one to a general numeric analysis book. At least not if you didn't invent the topic by yourself, which those authors have not. $\endgroup$ Mar 27, 2015 at 0:48
-
$\begingroup$ You're absolutely right. I wouldn't say the formulas are horrible, though they could be improved. The special cases is actually rather nice as tests/comparisons, and they do give a general formula, which must be the same as Fornberg's. $\endgroup$– jvriesemMar 30, 2015 at 19:23
-
2$\begingroup$ @jvriesem Please note that the cited paper has the wrong sign in the third term in Eqn. (15b) $\endgroup$– TarekSep 6, 2016 at 7:58
$\begingroup$
$\endgroup$
I once worked on this for a bit. It makes a big difference where your data comes from. Typically it will contain errors from experimental data, or is quoted only with finite precision. In such cases do not fit a curve through the points; this will just pick up the errors. It is MUCH better to take a range of data around the point of interest, fit a smooth curve to it, and differentiate that. The details (how big an interval, what order of curve-fitting) really depend on what sort of data you have, and which derivative you want. If this is important to you, I recommend experimenting with different methods and choosing the one you like.
A small illustrative problem. Suppose you have 100 points nonuniformly distributed in [0,1] and you want to know the second derivative at x=1/2. Should you fit exactly to the three nearest points and reconcile yourself to their errors? Or should you fit a parabola to every point you have, including the ones a long way from x=1/2? You will soon see that this is a hard problem. The question is what works best IN YOUR CASE
$\begingroup$
$\endgroup$
I found this paper on finite difference formulas with unequal sub-intervals. I'm going to use this instead of interpolation. Once I type all the formulas out, I'll post them here.
$\begingroup$
$\endgroup$
3
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
-
1$\begingroup$ 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. $\endgroup$– DanDec 21, 2011 at 6:54
-
8
-
$\begingroup$ And also because it just does not work in practice because there is no limit involved. If you take h too big you get truncation error and if you take h too small you get rounding error. $\endgroup$ Mar 15, 2020 at 15:35