hide long namesshow long names
hide short namesshow short names
Integer type:  int32  int64  nag_int  show int32  show int32  show int64  show int64  show nag_int  show nag_int

PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

NAG Toolbox: nag_numdiff_rcomm (d04ba)

Purpose

nag_numdiff_rcomm (d04ba) calculates a set of derivatives (up to order 1414) 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

[der, erest, ifail] = d04ba(xval, fval)
[der, erest, ifail] = nag_numdiff_rcomm(xval, fval)

Description

nag_numdiff_rcomm (d04ba) provides a set of approximations:
der(j) ,   j = 1,2,,14
derj ,   j=1,2,,14
to the derivatives:
f(j) (x0) ,   j = 1,2,,14
f(j) (x0) ,   j=1,2,,14
of a real valued function f(x)f(x) at a real abscissa x0x0, together with a set of error estimates:
erest(j) ,   j = 1,2,,14
erestj ,   j=1,2,,14
which hopefully satisfy:
|der(j)f(j)(x0)| < erest(j) ,   j = 1,2,,14 .
| derj - f (j) (x0) | < erestj ,   j= 1,2,,14 .
The results der(j)derj and erest(j)erestj are based on 2121 function values:
f(x0) , f (x0 ± (2i1)h) ,   i = 1,2,,10 .
f(x0) , f ( x0 ± (2i-1) h ) ,   i=1,2,,10 .
The abscissae xx and the corresponding function values f(x)f(x) should be passed into nag_numdiff_rcomm (d04ba) as the vectors xval and fval respectively. The step size hh is derived from the abscissae in xval. See Section [Further Comments] for a discussion of how the derived value of hh 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:
Tk,p,s ,   p = s,s + 1,,6 ,   k = 0,1,,9p
T k,p,s ,   p=s,s+1,,6 ,   k=0,1,,9-p
each of which is an approximation to f(2s + 1) (x0) / (2s + 1) ! f (2s+1) ( x0 ) / (2s+1) ! . A specific approximation Tk,p,sTk,p,s is of polynomial degree 2p + 22p+2 and is based on polynomial interpolation using function values f(x0 ± (2i1)h)f(x0±(2i-1)h), for k = k,,k + pk=k,,k+p. In the absence of round-off error, the better approximations would be associated with the larger values of pp and of kk. However, round-off error in function values has an increasingly contaminating effect for successively larger values of pp. This function proceeds to make a judicious choice between all the approximations in the following way.
For a specified value of ss, let:
Rp = Up Lp ,   p = s,s + 1,,6
Rp = Up - Lp ,   p=s,s+1,,6
where Up = maxk  (Tk,p,s) Up = maxk (Tk,p,s)  and Lp = mink  (Tk,p,s) Lp = mink (Tk,p,s) , for k = 0,1,,9pk=0,1,,9-p, and let pp- be such that Rp = minp  (Rp) Rp- = minp (Rp) , for p = s,,6p=s,,6.
This function returns:
der(2s + 1) = 1/(8p) ×
(9p )
T k, p, s UpLp
k = 0
(2s + 1) !
der2s+1 = 1 8-p- × { k=0 9-p- T k, p-, s - Up- - Lp- } (2s+1) !
and
erest(2s + 1) = Rp × (2s + 1) ! × K2s + 1
erest2s+1 = Rp- × (2s+1) ! × K 2s+1
where KjKj is a safety factor which has been assigned the values:
Kj = 1Kj=1, j9j9
Kj = 1.5Kj=1.5, j = 10,11j=10,11
Kj = 2Kj=2 j12j12,
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:     xval(2121) – double array
The abscissae at which the function has been evaluated, as described in Section [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 {x0, x0 ± (2j1) h } {x0, x0 ± (2j-1) h } , for j = 1,2,,10j=1,2,,10.
2:     fval(2121) – double array
fval(j)fvalj must contain the function value at xval(j)xvalj, for j = 1,2,,21j=1,2,,21.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     der(1414) – double array
The 14 derivative estimates.
2:     erest(1414) – double array
The 14 error estimates for the derivatives.
3:     ifail – int64int32nag_int scalar
ifail = 0ifail=0 unless the function detects an error (see [Error Indicators and Warnings]).

Error Indicators and Warnings

Errors or warnings detected by the function:
  ifail = 1ifail=1
On entry, the values of xval are not correctly spaced.
The derived hh is below tolerance.

Accuracy

The accuracy of the results is problem dependent. An estimate of the accuracy of each result der(j)derj is returned in erest(j)erestj (see Sections [Description], [Parameters] and [Further Comments]).
A basic feature of any floating point 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 der(14)der14 will be unusable. As an aid to this process, the sign of erest(j)erestj 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 der(j)derj 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 hh (see Section [Further Comments]). The only hard and fast rule is that for a given objective function and hh, the values of erest(j)erestj increase with increasing jj. The limit of 1414 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 hh. The overall accuracy is diminished as hh becomes small (because of the effect of round-off error) and as hh becomes large (because the discretization error also becomes large). If this function is used four or five times with different values of hh one can find a reasonably good value. A process in which the value of hh 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 x0x0 has a finite radius of convergence RR, the choices of h > R / 19h>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 x0x0 will be stored in xval(11)xval11.

Example

function nag_numdiff_rcomm_example
fprintf('\nFind the derivatives of the polygamma (psi) function\n');
fprintf('using function values generated by nag_specfun_psi_deriv_real.\n\n');
fprintf('Demonstrate the effect of successively reducing hbase.\n\n');
hbase = 0.0025;

% Set the target location and calculate the objective value
x_0 = 0.05;
[f_0, ifail] = nag_specfun_psi_deriv_real(x_0, int64(0));

% Compute the actual derivatives using nag_specfun_psi_deriv_real for comparison
actder = zeros(3, 1);
for j=1:3
  [actder(j), ifail] = nag_specfun_psi_deriv_real(x_0, int64(j));
end

der_comp = zeros(14*4, 3);
% Attempt 4 applications, reducing hbase by factor 0.1 each time
for j=1:4
  % Generate the abscissa xval using nag_numdiff_sample
  [xval] = nag_numdiff_sample(x_0, hbase);

  % Calculate the corresponding objective function values
  fval(11) = f_0;
  for k=1:10
    [fval(k),    ifail] = nag_specfun_psi_deriv_real(xval(k),    int64(0));
    [fval(k+11), ifail] = nag_specfun_psi_deriv_real(xval(k+11), int64(0));
  end

  % Call nag_numdiff_rcomm to calculate the derivative estimates
  [der, erest, ifail] = nag_numdiff_rcomm(xval, fval);
  if ifail ~= 0
    break;
  end

  % Store results in der_comp
  der_comp(j:4:14*4, 1) = hbase;
  der_comp(j:4:14*4, 2) = der;
  der_comp(j:4:14*4, 3) = erest;

  % Decrease hbase for next loop
  hbase = hbase*0.1;
end

% Display results for first 3 derivatives
if ifail == 0
  for i=1:3
    fprintf('\nderivative %d calculated using nag_specfun_psi_deriv_real: %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(4*(i-1)+j,:));
    end
  end
else
  fprintf('\nnag_numdiff_rcomm failed with ifail = %d\n', ifail);
end
 

Find the derivatives of the polygamma (psi) function
using function values generated by nag_specfun_psi_deriv_real.

Demonstrate the effect of successively reducing hbase.


derivative 1 calculated using nag_specfun_psi_deriv_real:  4.0153e+02
Derivative and error estimates for derivative 1
      hbase        der(1)      erest(1)
    2.5000e-03   4.0204e+02   1.3940e+02
    2.5000e-04   4.0153e+02   4.1666e-11
    2.5000e-05   4.0153e+02   2.1799e-10
    2.5000e-06   4.0153e+02   1.2417e-09

derivative 2 calculated using nag_specfun_psi_deriv_real: -1.6002e+04
Derivative and error estimates for derivative 2
      hbase        der(2)      erest(2)
    2.5000e-03  -1.6022e+04   5.5760e+03
    2.5000e-04  -1.6002e+04   1.3738e-07
    2.5000e-05  -1.6002e+04   6.0436e-06
    2.5000e-06  -1.6002e+04   9.5762e-04

derivative 3 calculated using nag_specfun_psi_deriv_real:  9.6001e+05
Derivative and error estimates for derivative 3
      hbase        der(3)      erest(3)
    2.5000e-03   9.1465e+05  -7.3750e+06
    2.5000e-04   9.6001e+05   2.1981e-04
    2.5000e-05   9.6001e+05   9.4922e-02
    2.5000e-06   9.6001e+05   5.9679e+01

function d04ba_example
fprintf('\nFind 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');
hbase = 0.0025;

% Set the target location and calculate the objective value
x_0 = 0.05;
[f_0, ifail] = s14ae(x_0, int64(0));

% Compute the actual derivatives using s14ae for comparison
actder = zeros(3, 1);
for j=1:3
  [actder(j), ifail] = s14ae(x_0, int64(j));
end

der_comp = zeros(14*4, 3);
% Attempt 4 applications, reducing hbase by factor 0.1 each time
for j=1:4
  % Generate the abscissa xval using d04bb
  [xval] = d04bb(x_0, hbase);

  % Calculate the corresponding objective function values
  fval(11) = f_0;
  for k=1:10
    [fval(k),    ifail] = s14ae(xval(k),    int64(0));
    [fval(k+11), ifail] = s14ae(xval(k+11), int64(0));
  end

  % Call d04ba to calculate the derivative estimates
  [der, erest, ifail] = d04ba(xval, fval);
  if ifail ~= 0
    break;
  end

  % Store results in der_comp
  der_comp(j:4:14*4, 1) = hbase;
  der_comp(j:4:14*4, 2) = der;
  der_comp(j:4:14*4, 3) = erest;

  % Decrease hbase for next loop
  hbase = hbase*0.1;
end

% Display results for first 3 derivatives
if ifail == 0
  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(4*(i-1)+j,:));
    end
  end
else
  fprintf('\nd04ba failed with ifail = %d\n', ifail);
end
 

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.5000e-03   4.0204e+02   1.3940e+02
    2.5000e-04   4.0153e+02   4.1666e-11
    2.5000e-05   4.0153e+02   2.1799e-10
    2.5000e-06   4.0153e+02   1.2417e-09

derivative 2 calculated using s14ae: -1.6002e+04
Derivative and error estimates for derivative 2
      hbase        der(2)      erest(2)
    2.5000e-03  -1.6022e+04   5.5760e+03
    2.5000e-04  -1.6002e+04   1.3738e-07
    2.5000e-05  -1.6002e+04   6.0436e-06
    2.5000e-06  -1.6002e+04   9.5762e-04

derivative 3 calculated using s14ae:  9.6001e+05
Derivative and error estimates for derivative 3
      hbase        der(3)      erest(3)
    2.5000e-03   9.1465e+05  -7.3750e+06
    2.5000e-04   9.6001e+05   2.1981e-04
    2.5000e-05   9.6001e+05   9.4922e-02
    2.5000e-06   9.6001e+05   5.9679e+01


PDF version (NAG web site, 64-bit version, 64-bit version)
Chapter Contents
Chapter Introduction
NAG Toolbox

© The Numerical Algorithms Group Ltd, Oxford, UK. 2009–2013