NAG C Library Function Document

nag_numdiff_1d_real (d04aac)


nag_numdiff_1d_real (d04aac) 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.


#include <nag.h>
#include <nagd04.h>
void  nag_numdiff_1d_real (double xval, Integer nder, double hbase, double der[], double erest[],
double (*fun)(double x, Nag_Comm *comm),
Nag_Comm *comm, NagError *fail)


nag_numdiff_1d_real (d04aac) provides a set of approximations:
der[j-1],  j=1,2,,n  
to the derivatives:
f j x0,   j= 1,2,,n  
of a real valued function fx at a real abscissa x0, together with a set of error estimates:
erest[j-1],  j=1,2,,n  
which hopefully satisfy:
der[j-1]-f j x0<erest[j-1],   j= 1,2,,n.  
You must provide the value of x0, a value of n (which is reduced to 14 should it exceed 14), a function which evaluates fx for all real x, and a step length h. The results der[j-1] and erest[j-1] are based on 21 function values:
fx0,fx0±2i-1h,  i=1,2,,10.  
Internally nag_numdiff_1d_real (d04aac) 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_1d_real (d04aac) calculates a set of numbers:
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!. A specific approximation Tk,p,s is of polynomial degree 2p+2 and is based on polynomial interpolation using function values fx0±2i-1h, for k=k,,k+p. In the absence of round-off error, the better approximations would be associated with the larger values of p and of k. However, round-off 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:
Rp = Up - Lp ,   p=s,s+1,,6  
where Up = maxk Tk,p,s  and Lp = mink Tk,p,s , for k=0,1,,9-p, and let p- be such that Rp- = minp Rp , for p=s,,6.
The function returns:
der[2s] = 1 8-p- × k=0 9-p- T k, p-, s - Up- - Lp- 2s+1 !  
erest[2s] = Rp- × 2s+1 ! × K 2s+1  
where Kj is a safety factor which has been assigned the values:
Kj=1, j9
Kj=1.5, j=10,11
Kj=2 j12,
on the basis of performance statistics.
The even order derivatives are calculated in a precisely analogous manner.


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


1:     xval doubleInput
On entry: the point at which the derivatives are required, x0.
2:     nder IntegerInput
On entry: must be set so that its absolute value is the highest order derivative required.
All derivatives up to order minnder,14 are calculated.
nder<0 and nder is even
Only even order derivatives up to order min-nder,14 are calculated.
nder<0 and nder is odd
Only odd order derivatives up to order min-nder,13 are calculated.
3:     hbase doubleInput
On entry: the initial step length which may be positive or negative. For advice on the choice of hbase see Section 9.
Constraint: hbase0.0.
4:     der[14] doubleOutput
On exit: der[j-1] contains an approximation to the jth derivative of fx at x=xval, so long as the jth derivative is one of those requested by you when specifying nder. For other values of j, der[j-1] is unused.
5:     erest[14] doubleOutput
On exit: an estimate of the absolute error in the corresponding result der[j-1] so long as the jth derivative is one of those requested by you when specifying nder. The sign of erest[j-1] is positive unless the result der[j-1] is questionable. It is set negative when der[j-1]<erest[j-1] or when for some other reason there is doubt about the validity of the result der[j-1] (see Section 6). For other values of j, erest[j-1] is unused.
6:     fun function, supplied by the userExternal Function
fun must evaluate the function fx at a specified point.
The specification of fun is:
double  fun (double x, Nag_Comm *comm)
1:     x doubleInput
On entry: 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_1d_real (d04aac) the only values of x for which fx will be required are x=xval and x=xval±2j-1hbase, for j=1,2,,10; and
(ii) fx0 is always computed, but it is disregarded when only odd order derivatives are required.
2:     comm Nag_Comm *
Pointer to structure of type Nag_Comm; the following members are relevant to fun.
userdouble *
iuserInteger *
The type Pointer will be void *. Before calling nag_numdiff_1d_real (d04aac) you may allocate memory and initialize these pointers with various quantities for use by fun when called from nag_numdiff_1d_real (d04aac) (see Section in How to Use the NAG Library and its Documentation).
Note: fun should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by nag_numdiff_1d_real (d04aac). If your code inadvertently does return any NaNs or infinities, nag_numdiff_1d_real (d04aac) is likely to produce unexpected results.
7:     comm Nag_Comm *
The NAG communication argument (see Section in How to Use the NAG Library and its Documentation).
8:     fail NagError *Input/Output
The NAG error argument (see Section 3.7 in How to Use the NAG Library and its Documentation).

Error Indicators and Warnings

Dynamic memory allocation failed.
See Section in How to Use the NAG Library and its Documentation for further information.
On entry, argument value had an illegal value.
On entry, nder=0.
Constraint: nder0.
An internal error has occurred in this function. Check the function call and any array sizes. If the call is correct then please contact NAG for assistance.
See Section 2.7.6 in How to Use the NAG Library and its Documentation for further information.
Your licence key may have expired or may not have been installed correctly.
See Section 2.7.5 in How to Use the NAG Library and its Documentation for further information.
On entry, hbase=0.0.
Constraint: hbase0.0.


The accuracy of the results is problem dependent. An estimate of the accuracy of each result der[j-1] is returned in erest[j-1] (see Sections 3, 5 and 9).
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[13] will be unusable. As an aid to this process, the sign of erest[j-1] 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-1] 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 9). The only hard and fast rule is that for a given funxval and hbase, the values of erest[j-1] 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.

Parallelism and Performance

nag_numdiff_1d_real (d04aac) is not threaded in any implementation.

Further Comments

The time taken by nag_numdiff_1d_real (d04aac) 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 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 funx about xval has a finite radius of convergence R, the choices of hbase>R/19 are not likely to lead to good results. In this case some function values lie outside the circle of convergence.


This example evaluates the odd-order derivatives of the function:
fx = 12 e 2x-1  
up to order 7 at the point x=12 . 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.

Program Text

Program Text (d04aace.c)

Program Data


Program Results

Program Results (d04aace.r)