The Classic 'Central Differences' Method

Numerical Methods for Derivative Estimation

If you try to research the topic "numerical estimation of derivatives," what you are likely to find — maybe the only thing that you can find — is the method of Central Differences. Because of its widespread presence, it is important to understand in general terms how this method works. And then, you should have a good basis for understanding why its results are so lousy so much of the time.


What are Central Differences?

Central Differences (often confused with various formulas that use them) are a mathematical obscurity once considered essential for calculating advanced mathematical functions. That was back in the day when there were no electronic calculators or digital computers. Back then, calculating highly precise values of functions as commonplace as an exponential or cosine function could require hours of tedious effort. To reduce this labor to a manageable level, tables of mathematical functions were typically developed and printed by experts. This worked great when the values that you wanted happened to be the ones provided directly as table entries. But what could you do for the values that lie between the table entries?

Instead of giving up on the table and performing direct calculations, it was much more efficient to interpolate the in-between values, based on the values present in the table. There are some popular and effective methods for this; and there are many, many other equivalent formulas. The interpolated results can be very satisfactory as long as the approximation errors are acceptably small.

The function tables will have abscissa values (the function's inputs) at small equally-spaced steps, so that the ordinate values (the function's outputs) change smoothly and only a small amount from line to line. The interpolation formulas can use as few as two tabulated points, but more accurate methods will use more points. The best methods for interpolation within the body of the table will select multiple tabulated points on each side of the point to be evaluated. Formulas can be expressed in terms of tabulated function values directly, or in terms of difference values obtained by subtracting pairs of tabulated terms. Central differences are a type of differences that cluster around the evaluation point in a balanced manner.

Some notations to help define the central differences:

    x === the input value

    f === the function value to be calculated

    y === f(x), the value returned by the function 

    h === the "step size" for the input variable x

    p === a partial step to reach the desired evaluation position

    δ === a calculated central difference

The table is often normalized so that the step size h is considered "one unit". If this is done, the results can be scaled later to restore the original tabulation step size. Input values can then be identified by an integer position index. You will very often see derivative estimator formulas showing a 1/h multiplier in front of it, as a reminder to apply this scaling after applying the interpolation formula.

Differences in the table are simply the result of subtracting consecutive tabulated values from the column to the immediate left. In particular, central differences are calculated as follows.

     location   ordinate   first difference  second difference
                   y            δy              δ2y
        -2        y-2     
                            y-1 - y-2
        -1        y-1                       δy-1 - δy-2
                            y0 - y-1
         0        y0                        δy0 - δy-1            
                            y1 - y0
         1        y1                        δy1 - δy0 
                            y2 - y1
         2        y2                        δy2 - δy1     
                            y3 - y2
         3        y3                                     etc.

For example, suppose we have input values x and corresponding function values y. As shown below, we can easily calculate the difference terms to populate the table.

     abcissa   ordinate   first difference  second difference
       x          y            δy              δ2y
     0.05     0.005005     
     0.010    0.010023                    0.000047
     0.015    0.015088                    0.000096             
     0.020    0.020249                    0.000273
     0.025    0.025683                    0.000380      
     0.030    0.031697

These calculations are messy but easy. Notice that higher order differences tend to have fewer significant digits, which makes hand calculations much simpler — though this is not nearly as important for calculations by automatic computers.


Polynomial Approximations

The basic idea of the polynomial interpolation process is to identify a polynomial of low order that passes through selected tabulated points exactly. Polynomials have the highly desirable features that they are very simple to calculate, and they have only a few nonzero derivatives, so they are very smooth. Then, this constructed polynomial is evaluated at the particular position of interest between the tabulated ordinates to yield the interpolated function value. The coefficients of the constructed interpolation polynomials can be expressed directly in terms of tabulated values, or in terms of differences. The results are the same; it is all a matter of notation convenience.

A simple linear interpolation assumes that polynomial is first order, and the interpolated points lie along a straight line between tabulated points. When the function shows very little curvature, this crude approximation is often good enough. For example, here is a linear interpolation to estimate a value from the table above.

    y(0.0156)  ≈ 0.015088  +  0.005161 * (0.001/0.005)   

The main weakness is that the linear formula has no way to represent any sort of curvature , which any sort of interesting function will have. We could get a much more accurate fit using a polynomial that represents at least some of the curvature features.

For example, the interpolation could be based on a second-order polynomial instead of a straight line. In the tabulation, a quadratic interpolation function would pass through three values in a neighborhood of the point of evaluation. The well known Newton's Interpolation Formula[1] can be used to generate suitable polynomials of any order. The principal disadvantage of this formula is that it is one sided it does not use tabulated information in a balanced manner. So we will not discuss it further, but instead, illustrate how constructing interpolation polynomials through a balanced set of points is straightforward.


Constructing balanced polynomial approximations

Start with a linear interpolation formula. The linear polynomial curve passes through two ordinate points. We want to build upon this and extend the fit to a higher-order polynomial by including information from an additional table line. To simplify this presentation, assume that the input variables are scaled and shifted so that the table abscissa step is 1.0 unit, and the location of the evaluation is near the point x = 0.0.

    y(0.0 + p) ≈ y0  +  (y1 - y0) p

It is easy to verify that substituting values p=0.0 or p=1.0 produces results that match the tabulated y0 and y1 values. We want to preserve this property when extending this formula with a new quadratic term. We can see that this is accomplished for any nonzero value of the new coefficient in the following formula.

    y(0.0 + p) ≈ y0  +  (y1 - y0) p   +   [coeff] · p (p-1)

Now we need to find a value for [coeff] that also forces the approximation to match the tabulated value at a new position. Let's make the formula match at the point where p = -1. Insert p = -1 into the above formula.

    y(0.0 + p) ≈ y0  +  (-y1 + y0) +  2.0 [coeff]

It is now possible to solve directly for the new coefficient.

                  y1  -  2.0 * y0  +  y-1
    [coeff]  =   --------------------------

Using this new coefficient, complete the extended interpolation formula.

                                               (y1 - 2 y0 + y-1)
    y(0.0 + p)  ≈  y0  + (y1 - y0) (p)   +    ------------------  (p2 - p)

The coefficient terms can now be expressed using the notation of central differences.

0 + δ-1)               (δ0 - δ-1)
    y(0.0 + p)  ≈  y0  +  -------------  (p)   +   -----------  (p2 - p)
                             2.0                       2.0

Or collecting terms,

    y(0.0 + p)  ≈  y0  +   δ0 p    +   --------  p2 

This process can be applied repeatedly, adding a new term with new unknown coefficient, and evaluating the new coefficient value, thus extending the polynomial curve to as high an order as desired.

The details of all of this are not nearly as important as the understanding that there is nothing particularly difficult about it. Having constructed the fit, expressing the polynomial coefficients in terms of central differences is just a matter of notation. The Central Differences formulas happen to provide a particularly elegant representation with a balanced set of points surrounding the point where an interpolated value is needed.


Approximating derivatives

Once you have a polynomial that matches the function values, the derivative of this polynomial should be a good estimate of the derivative of the underlying function producing the tabulated data. The formula for differentiating a polynomial function is simple and very well known.

We are going to concentrate on the important case that the derivative value is to be evaluated at a point that corresponds to a tabulated function value. Define the indexing for the interpolation polynomial so that "index point zero" of the table aligns with the evaluation point of interest. Since the polynomial argument p is zero at this point, this results in a drastic simplification when the polynomial is evaluated.

Though theoretically simple, it is rather laborious to construct all of the derivative estimator formulas you might need in this manner. Fortunately, all of this has been done for you in the past, and you do not need to repeat all of that work. Here is a table[2] of coefficients for the derivative estimators that result from using central differences notations.

Derivative Fit order -4 -3 -2 -1 0 1 2 3 4
1 2   _ _ _ -1/2 0 1/2 _ _ _
1 4   _ _ 1/12 -2/3 0 2/3 -1/12 _ _
1 6   _ -1/60 3/20 -3/4 0 3/4 -3/20 1/60 _
1 8   1/280 -4/105 1/5 -4/5 0 4/5 -1/5 4/105 1/280

For example, the second estimator formula from the table above is implemented as:

    y'k  =  (1/h) ( 1/12 yk-2 - 2/3 yk-1 + 2/3 yk+1 - 1/12 yk+2 )

where the 1/h term accounts for the actual size of step in the independent (input) variable.

You can think of this formula as a vector of numerical estimator coefficients that are lined up with a sequence of consecutive function values. When the estimator is evaluated, corresponding numbers are pairwise multiplied, and then these intermediate products are summed to produce the desired estimate. To evaluate a succession of points, the vector of coefficients "slides" along the function data from one position to the next. This kind of mathematical operation is known as a convolution.[3]


What could possibly go wrong?

Though you will find many many sources telling you about the central differences estimators, what they are not likely to tell you is that results you get from them are generally pretty awful[4]. In fact, the only special case that produces high quality results is the case when the numerical values in the table are accurate to a very high degree of precision — which is the purpose for which Central Differences were devised. For data obtained any other way, for example by digitization, measurements, or approximate observation, you can expect some seriously disappointing results. The next section will discuss why.



[1] See the Wikipedia article Newton Polynomial for information about the Newton Interpolation Formula and related topics.

[2] See the Wikipedia article Finite difference coefficient for a relatively complete compilation.

[3] See the Wikipedia article Convolution for a more complete and formal discussion of convolution operations.

[4] Not everybody is fooled. See for example the notes by Stephen Roberts, ENGINEERING COMPUTATION Lecture 6, Computing derivatives and integrals.


Contents of the "Numerical Estimation of Derivatives" section of the website, including this page, are licensed under a Creative Commons Attribution-ShareAlike 4.0 International License unless otherwise specifically noted. For complete information about the terms of this license, see The license allows usage and derivative works for individual or commercial purposes under certain restrictions. For permissions beyond the scope of this license, see the contact information page for this site.