PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_numdiff (d04aa)
Purpose
nag_numdiff (d04aa) calculates a set of derivatives (up to order $14$) of a function of one real variable at a point, together with a corresponding set of error estimates, using an extension of the Neville algorithm.
Syntax
Description
nag_numdiff (d04aa) provides a set of approximations:
to the derivatives:
of a real valued function
$f\left(x\right)$ at a real abscissa
${x}_{0}$, together with a set of error estimates:
which hopefully satisfy:
You must provide the value of
${x}_{0}$, a value of
$n$ (which is reduced to
$14$ should it exceed
$14$), a function which evaluates
$f\left(x\right)$ for all real
$x$, and a step length
$h$. The results
${\mathbf{der}}\left(j\right)$ and
${\mathbf{erest}}\left(j\right)$ are based on
$21$ function values:
Internally
nag_numdiff (d04aa) calculates the odd order derivatives and the even order derivatives separately. There is an option you can use for restricting the calculation to only odd (or even) order derivatives. For each derivative the function employs an extension of the Neville Algorithm (see
Lyness and Moler (1969)) to obtain a selection of approximations.
For example, for odd derivatives, based on
$20$ function values,
nag_numdiff (d04aa) calculates a set of numbers:
each of which is an approximation to
${f}^{\left(2s+1\right)}\left({x}_{0}\right)/\left(2s+1\right)!$. A specific approximation
${T}_{\mathit{k},p,s}$ is of polynomial degree
$2p+2$ and is based on polynomial interpolation using function values
$f\left({x}_{0}\pm \left(2i1\right)h\right)$, for
$\mathit{k}=\mathit{k},\dots ,\mathit{k}+p$. In the absence of roundoff error, the better approximations would be associated with the larger values of
$p$ and of
$k$. However, roundoff error in function values has an increasingly contaminating effect for successively larger values of
$p$. This function proceeds to make a judicious choice between all the approximations in the following way.
For a specified value of
$s$, let:
where
${U}_{p}={\displaystyle \underset{\mathit{k}}{\mathrm{max}}}\phantom{\rule{0.25em}{0ex}}\left({T}_{\mathit{k},p,s}\right)$ and
${L}_{p}={\displaystyle \underset{\mathit{k}}{\mathrm{min}}}\phantom{\rule{0.25em}{0ex}}\left({T}_{\mathit{k},p,s}\right)$, for
$\mathit{k}=0,1,\dots ,9p$, and let
$\stackrel{}{\mathit{p}}$ be such that
${R}_{\stackrel{}{\mathit{p}}}={\displaystyle \underset{\mathit{p}}{\mathrm{min}}}\phantom{\rule{0.25em}{0ex}}\left({R}_{\mathit{p}}\right)$, for
$\mathit{p}=s,\dots ,6$.
The function returns:
and
where
${K}_{j}$ is a safety factor which has been assigned the values:
${K}_{j}=1$, 
$j\le 9$ 
${K}_{j}=1.5$, 
$j=10,11$ 
${K}_{j}=2$ 
$j\ge 12$, 
on the basis of performance statistics.
The even order derivatives are calculated in a precisely analogous manner.
References
Lyness J N and Moler C B (1966) van der Monde systems and numerical differentiation Numer. Math. 8 458–464
Lyness J N and Moler C B (1969) Generalised Romberg methods for integrals of derivatives Numer. Math. 14 1–14
Parameters
Compulsory Input Parameters
 1:
$\mathrm{xval}$ – double scalar

The point at which the derivatives are required, ${x}_{0}$.
 2:
$\mathrm{nder}$ – int64int32nag_int scalar

Must be set so that its absolute value is the highest order derivative required.
 ${\mathbf{nder}}>0$
 All derivatives up to order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{nder}},14\right)$ are calculated.
 ${\mathbf{nder}}<0$ and nder is even
 Only even order derivatives up to order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{nder}},14\right)$ are calculated.
 ${\mathbf{nder}}<0$ and nder is odd
 Only odd order derivatives up to order $\mathrm{min}\phantom{\rule{0.125em}{0ex}}\left({\mathbf{nder}},13\right)$ are calculated.
 3:
$\mathrm{hbase}$ – double scalar

The initial step length which may be positive or negative. For advice on the choice of
hbase see
Further Comments.
Constraint:
${\mathbf{hbase}}\ne 0.0$.
 4:
$\mathrm{fun}$ – function handle or string containing name of mfile

fun must evaluate the function
$f\left(x\right)$ at a specified point.
[result] = fun(x)
Input Parameters
 1:
$\mathrm{x}$ – double scalar

The value of the argument
$x$.
If you have equally spaced tabular data, the following information may be useful:
(i) 
in any call of nag_numdiff (d04aa) the only values of $x$ for which $f\left(x\right)$ will be required are $x={\mathbf{xval}}$ and
$x={\mathbf{xval}}\pm \left(2\mathit{j}1\right){\mathbf{hbase}}$, for $\mathit{j}=1,2,\dots ,10$; and 
(ii) 
$f\left({x}_{0}\right)$ is always computed, but it is disregarded when only odd order derivatives are required. 
Output Parameters
 1:
$\mathrm{result}$ – double scalar

The value of $f\left(x\right)$ at the specified point.
Optional Input Parameters
None.
Output Parameters
 1:
$\mathrm{der}\left(14\right)$ – double array

${\mathbf{der}}\left(j\right)$ contains an approximation to the
$j$th derivative of
$f\left(x\right)$ at
$x={\mathbf{xval}}$, so long as the
$j$th derivative is one of those requested by you when specifying
nder. For other values of
$j$,
${\mathbf{der}}\left(j\right)$ is unused.
 2:
$\mathrm{erest}\left(14\right)$ – double array

An estimate of the absolute error in the corresponding result
${\mathbf{der}}\left(j\right)$ so long as the
$j$th derivative is one of those requested by you when specifying
nder. The sign of
${\mathbf{erest}}\left(j\right)$ is positive unless the result
${\mathbf{der}}\left(j\right)$ is questionable. It is set negative when
$\left{\mathbf{der}}\left(j\right)\right<\left{\mathbf{erest}}\left(j\right)\right$ or when for some other reason there is doubt about the validity of the result
${\mathbf{der}}\left(j\right)$ (see
Error Indicators and Warnings). For other values of
$j$,
${\mathbf{erest}}\left(j\right)$ is unused.
 3:
$\mathrm{ifail}$ – int64int32nag_int scalar
${\mathbf{ifail}}={\mathbf{0}}$ unless the function detects an error (see
Error Indicators and Warnings).
Error Indicators and Warnings
Errors or warnings detected by the function:
 ${\mathbf{ifail}}=1$

On entry,  ${\mathbf{nder}}=0$, 
or  ${\mathbf{hbase}}=0.0$. 
 ${\mathbf{ifail}}=99$
An unexpected error has been triggered by this routine. Please
contact
NAG.
 ${\mathbf{ifail}}=399$
Your licence key may have expired or may not have been installed correctly.
 ${\mathbf{ifail}}=999$
Dynamic memory allocation failed.
If
ifail has a value zero on exit then
nag_numdiff (d04aa) has terminated successfully, but before any use is made of a derivative
${\mathbf{der}}\left(j\right)$ the value of
${\mathbf{erest}}\left(j\right)$ must be checked.
Accuracy
The accuracy of the results is problem dependent. An estimate of the accuracy of each result
${\mathbf{der}}\left(j\right)$ is returned in
${\mathbf{erest}}\left(j\right)$ (see
Description,
Arguments and
Further Comments).
A basic feature of any floatingpoint function for numerical differentiation based on real function values on the real axis is that successively higher order derivative approximations are successively less accurate. It is expected that in most cases ${\mathbf{der}}\left(14\right)$ will be unusable. As an aid to this process, the sign of ${\mathbf{erest}}\left(j\right)$ is set negative when the estimated absolute error is greater than the approximate derivative itself, i.e., when the approximate derivative may be so inaccurate that it may even have the wrong sign. It is also set negative in some other cases when information available to the function indicates that the corresponding value of ${\mathbf{der}}\left(j\right)$ is questionable.
The actual values in
erest depend on the accuracy of the function values, the properties of the machine arithmetic, the analytic properties of the function being differentiated and the usersupplied step length
hbase (see
Further Comments). The only hard and fast rule is that for a given
${\mathbf{fun}}\left({\mathbf{xval}}\right)$ and
hbase, the values of
${\mathbf{erest}}\left(j\right)$ increase with increasing
$j$. The limit of
$14$ is dictated by experience. Only very rarely can one obtain meaningful approximations for higher order derivatives on conventional machines.
Further Comments
The time taken by nag_numdiff (d04aa) depends on the time spent for function evaluations. Otherwise the time is roughly equivalent to that required to evaluate the function $21$ times and calculate a finite difference table having about $200$ entries in total.
The results depend very critically on the choice of the usersupplied step length
hbase. The overall accuracy is diminished as
hbase becomes small (because of the effect of roundoff error) and as
hbase becomes large (because the discretization error also becomes large). If the function is used four or five times with different values of
hbase one can find a reasonably good value. A process in which the value of
hbase is successively halved (or doubled) is usually quite effective. Experience has shown that in cases in which the Taylor series for
${\mathbf{fun}}\left({\mathbf{x}}\right)$ about
xval has a finite radius of convergence
$R$, the choices of
${\mathbf{hbase}}>R/19$ are not likely to lead to good results. In this case some function values lie outside the circle of convergence.
Example
This example evaluates the oddorder derivatives of the function:
up to order
$7$ at the point
$x=\frac{1}{2}$. Several different values of
hbase are used, to illustrate that:
(i) 
extreme choices of hbase, either too large or too small, yield poor results; 
(ii) 
the quality of these results is adequately indicated by the values of erest. 
Open in the MATLAB editor:
d04aa_example
function d04aa_example
fprintf('d04aa example results\n\n');
fprintf(['Four separate runs to calculate the first four odd order ',...
'derivatives of\n fun(x) = exp(2*x1)/2 at x = 0.5.\n']);
fprintf('The exact results are 1, 4, 16 and 64\n\n');
xval = 0.5;
nder = int64(7);
hbase = 0.5;
fun = @(x) exp(2*x1)/2;
for k = 1:4
hbase = 5*10^(k);
[der, erest, ifail] = d04aa( ...
xval, nder, hbase, fun);
fprintf('\nwith step length %9.4f the results are:\n',hbase);
fprintf('Order Derivative Error estimate\n');
for i = 1:2:abs(nder)
fprintf('%d%21.3e%21.4e\n',i, der(i), erest(i));
end
end
d04aa example results
Four separate runs to calculate the first four odd order derivatives of
fun(x) = exp(2*x1)/2 at x = 0.5.
The exact results are 1, 4, 16 and 64
with step length 0.5000 the results are:
Order Derivative Error estimate
1 1.392e+03 1.0734e+05
3 3.139e+03 1.4378e+05
5 8.762e+03 2.4790e+05
7 2.475e+04 4.4838e+05
with step length 0.0500 the results are:
Order Derivative Error estimate
1 1.000e+00 1.5294e11
3 4.000e+00 2.1127e09
5 1.600e+01 3.8162e07
7 6.400e+01 7.3906e05
with step length 0.0050 the results are:
Order Derivative Error estimate
1 1.000e+00 3.5527e14
3 4.000e+00 4.9597e10
5 1.600e+01 1.4335e05
7 6.404e+01 2.8513e01
with step length 0.0005 the results are:
Order Derivative Error estimate
1 1.000e+00 1.4289e13
3 4.000e+00 3.0894e07
5 1.599e+01 6.3314e01
7 3.826e+04 1.9644e+06
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015