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 (d04aa)

Purpose

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

[der, erest, ifail] = d04aa(xval, nder, hbase, fun)
[der, erest, ifail] = nag_numdiff(xval, nder, hbase, fun)

Description

nag_numdiff (d04aa) provides a set of approximations:
der(j),  j = 1,2,,n
derj,  j=1,2,,n
to the derivatives:
f(j)(x0),   j = 1,2,,n
f (j) (x0),   j= 1,2,,n
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,,n
erestj,  j=1,2,,n
which hopefully satisfy:
|der(j)f(j)(x0)| < erest(j),   j = 1,2,,n.
|derj-f (j) (x0)|<erestj,   j= 1,2,,n.
You must provide the value of x0x0, a value of nn (which is reduced to 1414 should it exceed 1414), a function which evaluates f(x)f(x) for all real xx, and a step length hh. 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.
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 2020 function values, nag_numdiff (d04aa) calculates a set of numbers:
Tk,p,s,   p = s,s + 1,,6,   k = 0,1,,9p
Tk,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.
The 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 (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:     xval – double scalar
The point at which the derivatives are required, x0x0.
2:     nder – int64int32nag_int scalar
Must be set so that its absolute value is the highest order derivative required.
nder > 0nder>0
All derivatives up to order min (nder,14)min(nder,14) are calculated.
nder < 0nder<0 and nder is even
Only even order derivatives up to order min (nder,14)min(-nder,14) are calculated.
nder < 0nder<0 and nder is odd
Only odd order derivatives up to order min (nder,13)min(-nder,13) are calculated.
3:     hbase – double scalar
The initial step length which may be positive or negative. For advice on the choice of hbase see Section [Further Comments].
Constraint: hbase0.0hbase0.0.
4:     fun – function handle or string containing name of m-file
fun must evaluate the function f(x)f(x) at a specified point.
[result] = fun(x)

Input Parameters

1:     x – double scalar
The value of the argument xx.
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 xx for which f(x)f(x) will be required are x = xvalx=xval and x = xval ± (2j1)hbasex=xval±(2j-1)hbase, for j = 1,2,,10j=1,2,,10; and
(ii) f(x0)f(x0) is always computed, but it is disregarded when only odd order derivatives are required.

Output Parameters

1:     result – double scalar
The result of the function.

Optional Input Parameters

None.

Input Parameters Omitted from the MATLAB Interface

None.

Output Parameters

1:     der(1414) – double array
der(j)derj contains an approximation to the jjth derivative of f(x)f(x) at x = xvalx=xval, so long as the jjth derivative is one of those requested by you when specifying nder. For other values of jj, der(j)derj is unused.
2:     erest(1414) – double array
An estimate of the absolute error in the corresponding result der(j)derj so long as the jjth derivative is one of those requested by you when specifying nder. The sign of erest(j)erestj is positive unless the result der(j)derj is questionable. It is set negative when |der(j)| < |erest(j)||derj|<|erestj| or when for some other reason there is doubt about the validity of the result der(j)derj (see Section [Error Indicators and Warnings]). For other values of jj, erest(j)erestj is unused.
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,nder = 0nder=0,
orhbase = 0.0hbase=0.0.
If ifail has a value zero on exit then nag_numdiff (d04aa) has terminated successfully, but before any use is made of a derivative der(j)derj the value of erest(j)erestj must be checked.

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 the function 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 user-supplied step length hbase (see Section [Further Comments]). The only hard and fast rule is that for a given fun(xval)fun(xval) and hbase, 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 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 2121 times and calculate a finite difference table having about 200200 entries in total.
The results depend very critically on the choice of the user-supplied step length hbase. The overall accuracy is diminished as hbase becomes small (because of the effect of round-off 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 fun(x)fun(x) about xval has a finite radius of convergence RR, the choices of hbase > R / 19hbase>R/19 are not likely to lead to good results. In this case some function values lie outside the circle of convergence.

Example

function nag_numdiff_example
xval = 0.5;
nder = int64(-7);
hbase = 0.5;
fun = @(x) 0.5*exp(2.0*x-1.0);
[der, erest, ifail] = nag_numdiff(xval, nder, hbase, fun)
 

der =

   1.0e+04 *

    0.1392
         0
   -0.3139
         0
    0.8762
         0
   -2.4753
         0
         0
         0
         0
         0
         0
         0


erest =

   1.0e+05 *

   -1.0734
         0
   -1.4378
         0
   -2.4790
         0
   -4.4838
         0
         0
         0
         0
         0
         0
         0


ifail =

                    0


function d04aa_example
xval = 0.5;
nder = int64(-7);
hbase = 0.5;
fun = @(x) 0.5*exp(2.0*x-1.0);
[der, erest, ifail] = d04aa(xval, nder, hbase, fun)
 

der =

   1.0e+04 *

    0.1392
         0
   -0.3139
         0
    0.8762
         0
   -2.4753
         0
         0
         0
         0
         0
         0
         0


erest =

   1.0e+05 *

   -1.0734
         0
   -1.4378
         0
   -2.4790
         0
   -4.4838
         0
         0
         0
         0
         0
         0
         0


ifail =

                    0



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