Integer, Intent (In)	::	nder
Integer, Intent (Inout)	::	ifail
Real (Kind=nag_wp), External	::	f
Real (Kind=nag_wp), Intent (In)	::	xval, hbase
Real (Kind=nag_wp), Intent (Out)	::	der(14), erest(14)

C Header Interface

#include <nag.h>

void	d04aaf_ (const double xval, const Integer nder, const double hbase, double der[], double erest[], double (NAG_CALL f)(const double x), Integer ifail)

The routine may be called by the names d04aaf or nagf_numdiff_fwd.

3 Description

d04aaf provides a set of approximations:

der (j), j = 1, 2, \dots, n

to the derivatives:

f^{(j)} (x_{0}), j = 1, 2, \dots, n

of a real valued function

f (x)

at a real abscissa

x_{0}

, together with a set of error estimates:

erest (j), j = 1, 2, \dots, n

which hopefully satisfy:

|der (j) - f^{(j)} (x_{0})| < erest (j), j = 1, 2, \dots, n .

You must provide the value of

x_{0}

, a value of

n

(which is reduced to

14

should it exceed

14

), a subroutine which evaluates

f (x)

for all real

x

, and a step length

h

. The results

der (j)

and

erest (j)

are based on

21

function values:

f (x_{0}), f (x_{0} \pm (2 i - 1) h), i = 1, 2, \dots, 10 .

Internally d04aaf 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 routine 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, d04aaf calculates a set of numbers:

T_{k, p, s}, p = s, s + 1, \dots, 6, k = 0, 1, \dots, 9 - p

each of which is an approximation to

f^{(2 s + 1)} (x_{0}) / (2 s + 1)!

. A specific approximation

T_{k, p, s}

is of polynomial degree

2 p + 2

and is based on polynomial interpolation using function values

f (x_{0} \pm (2 i - 1) h)

, for

k = k, \dots, 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 routine proceeds to make a judicious choice between all the approximations in the following way.

For a specified value of

s

, let:

R_{p} = U_{p} - L_{p}, p = s, s + 1, \dots, 6

where

U_{p} = \max_{k} (T_{k, p, s})

and

L_{p} = \min_{k} (T_{k, p, s})

, for

k = 0, 1, \dots, 9 - p

, and let

\bar{p}

be such that

R_{\bar{p}} = \min_{p} (R_{p})

, for

p = s, \dots, 6

The routine returns:

der (2 s + 1) = \frac{1}{8 - \bar{p}} \times \{\sum_{k = 0}^{9 - \bar{p}} T_{k, \bar{p}, s} - U_{\bar{p}} - L_{\bar{p}}\} (2 s + 1)!

and

erest (2 s + 1) = R_{\bar{p}} \times (2 s + 1)! \times K_{2 s + 1}

where

K_{j}

is a safety factor which has been assigned the values:

$K_{j} = 1$ ,	$j \leq 9$
$K_{j} = 1.5$ ,	$j = 10, 11$
$K_{j} = 2$	$j \geq 12$ ,

on the basis of performance statistics.

The even order derivatives are calculated in a precisely analogous manner.

4 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

5 Arguments

1: $xval$ – Real (Kind=nag_wp) Input

On entry: the point at which the derivatives are required,

x_{0}

2: $nder$ – Integer Input

On entry: must be set so that its absolute value is the highest order derivative required.

$nder > 0$: All derivatives up to order $\min (nder, 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$ – Real (Kind=nag_wp) Input

On entry: the initial step length which may be positive or negative. For advice on the choice of hbase see Section 9.

Constraint:

hbase \neq 0.0

4: $der (14)$ – Real (Kind=nag_wp) array Output

On exit:

der (j)

contains an approximation to the

j

th derivative of

f (x)

x = xval

, so long as the

j

th derivative is one of those requested by you when specifying nder. For other values of

j

der (j)

is unused.

5: $erest (14)$ – Real (Kind=nag_wp) array Output

On exit: an estimate of the absolute error in the corresponding result

der (j)

so long as the

j

th derivative is one of those requested by you when specifying nder. The sign of

erest (j)

is positive unless the result

der (j)

is questionable. It is set negative when

|der (j)| < |erest (j)|

or when for some other reason there is doubt about the validity of the result

der (j)

(see Section 6). For other values of

j

erest (j)

is unused.

6: $f$ – real (Kind=nag_wp) Function, supplied by the user. External Procedure

f must evaluate the function

f (x)

at a specified point.

The specification of f is:

Fortran Interface

Function f (

Real (Kind=nag_wp)	::	f
Real (Kind=nag_wp), Intent (In)	::	x

C Header Interface

double

f_ (const double *x)

1: $x$ – Real (Kind=nag_wp) Input

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 d04aaf the only values of $x$ for which $f (x)$ will be required are $x = xval$ and $x = xval \pm (2 j - 1) hbase$ , for $j = 1, 2, \dots, 10$ ; and
(ii) $f (x_{0})$ is always computed, but it is disregarded when only odd order derivatives are required.

f must either be a module subprogram USEd by, or declared as EXTERNAL in, the (sub)program from which d04aaf is called. Arguments denoted as Input must not be changed by this procedure.

Note: f should not return floating-point NaN (Not a Number) or infinity values, since these are not handled by d04aaf. If your code inadvertently does return any NaNs or infinities, d04aaf is likely to produce unexpected results.

7: $ifail$ – Integer Input/Output

On entry: ifail must be set to

0

- 1 or 1

. If you are unfamiliar with this argument you should refer to Section 4 in the Introduction to the NAG Library FL Interface for details.

For environments where it might be inappropriate to halt program execution when an error is detected, the value

- 1 or 1

is recommended. If the output of error messages is undesirable, then the value

1

is recommended. Otherwise, if you are not familiar with this argument, the recommended value is

0

. When the value $- 1 or 1$ is used it is essential to test the value of ifail on exit.

On exit:

ifail = 0

unless the routine detects an error or a warning has been flagged (see Section 6).

6 Error Indicators and Warnings

If on entry

ifail = 0

- 1

, explanatory error messages are output on the current error message unit (as defined by x04aaf).

Errors or warnings detected by the routine:

$ifail = 1$: On entry, $hbase = 0.0$ .
Constraint: $hbase \neq 0.0$ .

On entry, $nder = 0$ .
Constraint: $nder \neq 0$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 7 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 8 in the Introduction to the NAG Library FL Interface for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 9 in the Introduction to the NAG Library FL Interface for further information.

If ifail has a value zero on exit then d04aaf has terminated successfully, but before any use is made of a derivative

der (j)

the value of

erest (j)

must be checked.

7 Accuracy

The accuracy of the results is problem dependent. An estimate of the accuracy of each result

der (j)

is returned in

erest (j)

(see Sections 3, 5 and 9).

A basic feature of any floating-point routine 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)

will be unusable. As an aid to this process, the sign of

erest (j)

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 routine indicates that the corresponding value of

der (j)

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

f (xval)

and hbase, the values of

erest (j)

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.

8 Parallelism and Performance

d04aaf is not threaded in any implementation.

9 Further Comments

The time taken by d04aaf 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 routine 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

f (x)

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.

10 Example

This example evaluates the odd-order derivatives of the function:

f (x) = \frac{1}{2} e^{2 x - 1}

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.