# 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.

locationordinatefirst differencesecond differenceyδyδ-2 y^{2}y_{-2}y_{-1}- y_{-2}-1 y_{-1}δy_{-1}- δy_{-2}y_{0}- y_{-1}0 y_{0}δy_{0}- δy_{-1}y_{1}- y_{0}1 y_{1}δy_{1}- δy_{0}y_{2}- y_{1}2 y_{2}δy_{2}- δy_{1}y_{3}- y_{2}3 y_{3}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.

abcissaordinatefirst differencesecond differencexyδyδ0.05 0.005005 0.005018 0.010 0.010023 0.000047 0.005065 0.015 0.015088 0.000096 0.005161 0.020 0.020249 0.000273 0.005434 0.025 0.025683 0.000380 0.005814 0.030 0.031697^{2}y

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) ≈ y_{0}+ (y_{1}- y_{0}) p

It is easy to verify that substituting values `p=0.0`

or
`p=1.0`

produces results that match the tabulated
`y`

and _{0}`y`

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. _{1}

y(0.0 + p) ≈ y_{0}+ (y_{1}- y_{0}) 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) ≈ y_{0}+ (-y_{1}+ y_{0}) + 2.0 [coeff]

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

y_{1}- 2.0 * y_{0}+ y_{-1}[coeff] = -------------------------- 2.0

Using this new coefficient, complete the extended interpolation formula.

(y_{1}- 2 y_{0}+ y_{-1}) y(0.0 + p) ≈ y_{0}+ (y_{1}- y_{0}) (p) + ------------------ (p^{2}- p) 2.0

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

(δ_{0}+ δ_{-1}) (δ_{0}- δ_{-1}) y(0.0 + p) ≈ y_{0}+ ------------- (p) + ----------- (p^{2}- p) 2.0 2.0

Or collecting terms,

δ_{0}^{2}y(0.0 + p) ≈ y_{0}+ δ_{0}p + -------- p^{2}2.0

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 y_{k-2}- 2/3 y_{k-1}+ 2/3 y_{k+1}- 1/12 y_{k+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.

Footnotes:

[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.