NAG Library Chapter Introduction

D04 – Numerical Differentiation

This chapter is concerned with calculating approximations to derivatives of a function $f$.

The problem of numerical differentiation does not receive very much attention nowadays. Although the Taylor series plays a key role in much of classical analysis, the poor reputation enjoyed by numerical differentiation has led numerical analysts to construct techniques for most problems which avoid the explicit use of numerical differentiation.

One may briefly and roughly define the term numerical differentiation as any process in which numerical values of derivatives $f^{\left(s\right)}\left(x_0\right)$ are obtained from evaluations $f\left(x_i\right)$ of the function $f\left(x\right)$ at several abscissae $x_i$ near $x_0$. This problem can be stable, well-conditioned, and accurate when complex-valued abscissae are used. This was first pointed out by Lyness and Moler (1967). An item of numerical software for this appears in Lyness and Ande (1971). However, in many applications the use of complex-valued abscissae is either prohibitive or prohibited. The main difficulty in using real abscissae is that amplification of round-off error can completely obliterate meaningful results. In the days when one relied on hand calculating machines and tabular data, the frustration caused by this effect led to the abandonment of serious use of the process.

There are several reasons for believing that, in the present-day computing environment, numerical differentiation might find a useful role. The first is that, by present standards, it is rather a small-scale calculation, so its cost may well be negligible compared with any overall saving in cost that might result from its use. Secondly, the assignment of a step length $h$ is now generally open. One does not have to rely on tabular data. Thirdly, although the amplification of round-off error is an integral part of the calculation, its effect can be measured reliably and automatically by the routine at the time of the calculation.

Thus you do not have to gauge the round-off level (or noise level) of the function values or assess the effect of this on the accuracy of the results. A routine does this automatically, returning with each result an error estimate which has already taken round-off error amplification into account.

We now illustrate, by means of a very simple example, the importance of the round-off error effect. A very simple approximation of $f^{\prime }\left(0\right)$, based on the identity is If there were no precision problem, this formula would be the only one needed in the theory of first-order numerical differentiation. We could simply take $h={10}^{-40}$ (or $h={10}^{-1000}$) to obtain an excellent approximation based on two function values. It is only when we consider in detail how a machine with finite precision comes to calculate $\left(f\left(h\right)-f\left(-h\right)\right)/2h$ that the necessity for a sophisticated theory becomes apparent.

To simplify the subsequent description we shall assume that the quantities involved are neither so close to zero that the machine underflow characteristics need be considered nor so large that machine overflow occurs. The approximation mentioned above involves the function values $f\left(h\right)$ and $f\left(-h\right)$. In general no computer has available these numbers exactly. Instead it uses approximations $\hat{f}\left(h\right)$ and $\hat{f}\left(-h\right)$ whose relative accuracy is less than some tolerance ${\epsilon }_{\mathrm{f}}$. If the function $f\left(x\right)$ is a library function, for example $\sin x$, ${\epsilon }_f$ may coincide with the machine accuracy parameter ${\epsilon }_{\mathrm{m}}$. More generally the function $f\left(x\right)$ is calculated in a user-supplied routine and ${\epsilon }_f$ is larger than ${\epsilon }_m$ by a small factor $5$ or $6$ if the calculation is short or by some larger factor in an extended calculation. This factor is not usually known by you.

$\hat{f}\left(h\right)$ and $\hat{f}\left(-h\right)$ are related to $f\left(h\right)$ and $f\left(-h\right)$ by Thus even if the rest of the calculation were carried out exactly, it is trivial to show that The difference between the quantity actually calculated and the quantity which one attempts to calculate may be as large as ${\epsilon }_{f}f\left(\xi \right)/h$; for example using $h={10}^{-40}$ and ${\epsilon }_m={10}^{-7}$ this gives a ‘conditioning error’ of ${10}^{33}f\left(\xi \right)$.

In practice much more sophisticated formulae than (1) above are used, and for these and for the corresponding higher-derivative formulae the error analysis is different and more complicated in detail. But invariably the theory contains the same overall feature. In a finite length calculation, the error is composed of two main parts: a discretization error which increases as $h$ becomes large and is zero for $h=0$; and a ‘conditioning’ error due to amplification of round-off error in function evaluation, which increases as $h$ becomes small and is infinite for $h=0$.

The routines in this chapter have to take into account internally both these sources of error in the results. Thus they return pairs of results, $\mathbf{DER}\left(j\right)$ and $\mathbf{EREST}\left(j\right)$ where $\mathbf{DER}\left(j\right)$ is an approximation to $f^{\left(j\right)}\left(x_0\right)$ and $\mathbf{EREST}\left(j\right)$ is an estimate of the error in the approximation $\mathbf{DER}\left(j\right)$. If the routine has not been misled, $\mathbf{DER}\left(j\right)$ and $\mathbf{EREST}\left(j\right)$ satisfy In this respect, numerical differentiation routines are fundamentally different from other routines. You do not specify an error criterion. Instead the routine provides the error estimate and this may be very large.

We mention here a terminological distinction. A fully automatic routine is one in which you do not need to provide any information other than that required to specify the problem. Such a routine (at a cost in computing time) decides an appropriate internal parameter such as the step length $h$ by itself. On the other hand a routine which requires you to provide a step length $h$, but automatically chooses from several different formulae each based on the specified step length, is termed a semi-automatic routine.

The situation described above must have seemed rather depressing when hand machines were in common use. For example in the simple illustration one does not know the values of the quantities $f^{\prime \prime \prime }\left(\xi \right)$ or ${\epsilon }_f$ involved in the error estimates, and the idea of altering the value of $h$ and starting again must have seemed appalling. However, by present-day standards, it is a relatively simple matter to write a program which carries out all the previously considered time-consuming calculations which may seem necessary. None of the routines in this chapter are particularly revolutionary in concept. They simply utilize the computer for the sort of task for which it was originally designed. It carries out simple tedious calculations which are necessary to estimate the accuracy of the results of other even simpler tedious calculations.

(a)	At the present, there is only one numerical differentiation algorithm available in this chapter accessible using direct communication in D04AAF, and using reverse communication in D04BAF (see Section 3.2.3 in the Essential Introduction for a description of the difference between these two conventions). These are semi-automatic routines for obtaining approximations to the first fourteen derivatives of a real valued function $f\left(x\right)$ at a specified point $x_0$. For D04AAF, you must provide a FUNCTION representing $f\left(x\right)$, the value of $x_0$, an upper limit $n\le 14$ on the order of the derivatives required and a step length $h$. For D04BAF, you must supply the value of $x_0$, $20$ other abscissae and the function values at those abscissae. Both routines return a set of approximations $\mathbf{DER}\left(j\right)$ and corresponding error estimates $\mathbf{EREST}\left(j\right)$ which hopefully satisfy $$\left|\mathbf{DER}\left(j\right)-f^{\left(j\right)}\left(x_0\right)\right|\le \mathbf{EREST}\left(j\right)\text{, \hspace{1em}}j=1,2,\dots ,n\le 14\text{.}$$   It is important that the error estimate $\mathbf{EREST}\left(j\right)$ be taken into consideration before any use of $\mathbf{DER}\left(j\right)$ is made. The actual size of $\mathbf{EREST}\left(j\right)$ depends on the analytic structure of the function, on the computational precision used and on the step size $h$, and is difficult to predict. Thus you have to run the routine to find out how accurate the results are. Usually you will find the higher-order derivatives are successively more inaccurate and that past a certain order, say $10$ or $11$, the size of $\mathbf{EREST}\left(j\right)$ actually exceeds $\mathbf{DER}\left(j\right)$. Clearly when this happens the approximation $\mathbf{DER}\left(j\right)$ is unusable.
(b)	There is available in the algorithm section of CACM (see Lyness and Ande (1971)) a semi-automatic Fortran routine for numerical differentiation of an analytical function $f\left(z\right)$ at a possibly complex abscissa $z_0$. This is a stable problem. It can be used for any problem for which D04AAF might be used and produces more accurate results, and derivatives of arbitrary high order. However, even if $z_0$ is real and $f\left(z\right)$ is real for $z$, this routine requires a user-supplied FUNCTION which evaluates $f\left(z\right)$ for complex values of $z$ and it makes use of complex arithmetic.
(c)	Routines are available in Chapter E02 to differentiate functions which are polynomials (in Chebyshev series representation) or cubic splines (in B-spline representation).

None.

None.

Lyness J N and Ande G (1971) Algorithm 413, ENTCAF and ENTCRE: evaluation of normalised Taylor coefficients of an analytic function Comm. ACM 14(10) 669–675

Lyness J N and Moler C B (1967) Numerical differentiation of analytic functions SIAM J. Numer. Anal. 4(2) 202–210