PDF version (NAG web site
, 64bit version, 64bit version)
NAG Toolbox: nag_numdiff_rcomm (d04ba)
Purpose
nag_numdiff_rcomm (d04ba) calculates a set of derivatives (up to order
$14$) of a function at a point with respect to a single variable. A corresponding set of error estimates is also returned. Derivatives are calculated using an extension of the Neville algorithm. This function differs from
nag_numdiff (d04aa), in that the abscissae and corresponding function values must be calculated before this function is called. The abscissae may be generated using
nag_numdiff_sample (d04bb).
Syntax
Description
nag_numdiff_rcomm (d04ba) 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:
The results
${\mathbf{der}}\left(j\right)$ and
${\mathbf{erest}}\left(j\right)$ are based on
$21$ function values:
The abscissae
$x$ and the corresponding function values
$f\left(x\right)$ should be passed into
nag_numdiff_rcomm (d04ba) as the vectors
xval and
fval respectively. The step size
$h$ is derived from the abscissae in
xval. See
Further Comments for a discussion of how the derived value of
$h$ may affect the results of
nag_numdiff_rcomm (d04ba). The order in which the abscissae and function values are stored in
xval and
fval is irrelevant, provided that the function value at any given index corresponds to the value of the abscissa at the same index. Abscissae may be automatically generated using
nag_numdiff_sample (d04bb) if desired. For each derivative
nag_numdiff_rcomm (d04ba) employs an extension of the Neville Algorithm (see
Lyness and Moler (1969)) to obtain a selection of approximations.
For example, for odd derivatives, this function 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{}{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$.
This 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 (1969) Generalised Romberg methods for integrals of derivatives Numer. Math. 14 1–14
Parameters
Compulsory Input Parameters
 1:
$\mathrm{xval}\left(21\right)$ – double array

The abscissae at which the function has been evaluated, as described in
Description. These can be generated by calling
nag_numdiff_sample (d04bb). The order of the abscissae is irrelevant.
Constraint:
the values in
xval must span the set
$\left\{{x}_{0},{x}_{0}\pm \left(2\mathit{j}1\right)h\right\}$, for
$\mathit{j}=1,2,\dots ,10$.
 2:
$\mathrm{fval}\left(21\right)$ – double array

${\mathbf{fval}}\left(\mathit{j}\right)$ must contain the function value at ${\mathbf{xval}}\left(\mathit{j}\right)$, for $\mathit{j}=1,2,\dots ,21$.
Optional Input Parameters
None.
Output Parameters
 1:
$\mathrm{der}\left(14\right)$ – double array

The 14 derivative estimates.
 2:
$\mathrm{erest}\left(14\right)$ – double array

The 14 error estimates for the derivatives.
 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, the values of
xval are not correctly spaced.
The derived $h$ is below tolerance.
 ${\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.
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 nag_numdiff_rcomm (d04ba) 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 step length
$h$ (see
Further Comments). The only hard and fast rule is that for a given objective function and
$h$, 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 results depend very critically on the choice of the step length $h$. The overall accuracy is diminished as $h$ becomes small (because of the effect of roundoff error) and as $h$ becomes large (because the discretization error also becomes large). If this function is used four or five times with different values of $h$ one can find a reasonably good value. A process in which the value of $h$ is successively halved (or doubled) is usually quite effective. Experience has shown that in cases in which the Taylor series for the objective function about ${x}_{0}$ has a finite radius of convergence $R$, the choices of $h>R/19$ are not likely to lead to good results. In this case some function values lie outside the circle of convergence.
As mentioned, the order of the abscissae in
xval does not matter, provided the corresponding values of
fval are ordered identically. If the abscissae are generated by
nag_numdiff_sample (d04bb), then they will be in ascending order. In particular, the target abscissa
${x}_{0}$ will be stored in
${\mathbf{xval}}\left(11\right)$.
Example
This example evaluates the derivatives of the polygamma function, calculated using
nag_specfun_psi_deriv_real (s14ae), and compares the first
$3$ derivatives calculated to those found using
nag_specfun_psi_deriv_real (s14ae).
Open in the MATLAB editor:
d04ba_example
function d04ba_example
fprintf('d04ba example results\n\n');
fprintf('Find the derivatives of the polygamma (psi) function\n');
fprintf('using function values generated by s14ae.\n\n');
fprintf('Demonstrate the effect of successively reducing hbase.\n\n');
x_0 = 0.05;
[f_0, ifail] = s14ae(x_0, int64(0));
actder = zeros(3, 1);
for j=1:3
[actder(j), ifail] = s14ae(x_0, int64(j));
end
der_comp = zeros(14, 3, 4);
for j=1:4
hbase = 0.025*10^(j);
[xval] = d04bb(x_0, hbase);
fval(11) = f_0;
for k=[1:10,12:21]
[fval(k), ifail] = s14ae(xval(k), int64(0));
end
[der, erest, ifail] = d04ba(xval, fval);
der_comp(:, 1, j) = hbase;
der_comp(:, 2, j) = der;
der_comp(:, 3, j) = erest;
end
for i=1:3
fprintf('\nderivative %d calculated using s14ae: %11.4e\n', i, actder(i));
fprintf('Derivative and error estimates for derivative %d\n', i);
fprintf(' hbase der(%d) erest(%d)\n', i, i);
for j=1:4
fprintf(' %12.4e %12.4e %12.4e\n', der_comp(i,:,j));
end
end
d04ba example results
Find the derivatives of the polygamma (psi) function
using function values generated by s14ae.
Demonstrate the effect of successively reducing hbase.
derivative 1 calculated using s14ae: 4.0153e+02
Derivative and error estimates for derivative 1
hbase der(1) erest(1)
2.5000e03 4.0204e+02 1.3940e+02
2.5000e04 4.0153e+02 4.9170e11
2.5000e05 4.0153e+02 2.1799e10
2.5000e06 4.0153e+02 1.1826e09
derivative 2 calculated using s14ae: 1.6002e+04
Derivative and error estimates for derivative 2
hbase der(2) erest(2)
2.5000e03 1.6022e+04 5.5760e+03
2.5000e04 1.6002e+04 1.2831e07
2.5000e05 1.6002e+04 6.0543e06
2.5000e06 1.6002e+04 9.5762e04
derivative 3 calculated using s14ae: 9.6001e+05
Derivative and error estimates for derivative 3
hbase der(3) erest(3)
2.5000e03 9.1465e+05 7.3750e+06
2.5000e04 9.6001e+05 2.3718e04
2.5000e05 9.6001e+05 4.2253e02
2.5000e06 9.6001e+05 5.9679e+01
PDF version (NAG web site
, 64bit version, 64bit version)
© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2015