) 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 routine differs from d04aaf, in that the abscissae and corresponding function values must be calculated before this routine is called. The abscissae may be generated using d04bbf.

2

Specification

Fortran Interface

Subroutine d04baf (

xval, fval, der, erest, ifail)

Integer, Intent (Inout)	::	ifail
Real (Kind=nag_wp), Intent (In)	::	xval(21), fval(21)
Real (Kind=nag_wp), Intent (Out)	::	der(14), erest(14)

C Header Interface

#include nagmk26.h

void	d04baf_ ( const double xval[], const double fval[], double der[], double erest[], Integer *ifail)

3

Description

d04baf provides a set of approximations:

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

to the derivatives:

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

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, 14

which hopefully satisfy:

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

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 .

The abscissae

x

and the corresponding function values

f (x)

should be passed into d04baf as the vectors xval and fval respectively. The step size

h

is derived from the abscissae in xval. See Section 9 for a discussion of how the derived value of

h

may affect the results of d04baf. 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 d04bbf if desired. For each derivative d04baf employs an extension of the Neville Algorithm (see Lyness and Moler (1969)) to obtain a selection of approximations.

For example, for odd derivatives, this routine 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

This 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 (1969) Generalised Romberg methods for integrals of derivatives Numer. Math. 14 1–14

5

Arguments

1: $xval (21)$ – Real (Kind=nag_wp) arrayInput: On entry: the abscissae at which the function has been evaluated, as described in Section 3. These can be generated by calling d04bbf. The order of the abscissae is irrelevant.

Constraint: the values in xval must span the set $\{x_{0}, x_{0} \pm (2 j - 1) h\}$ , for $j = 1, 2, \dots, 10$ .
2: $fval (21)$ – Real (Kind=nag_wp) arrayInput: On entry: $fval (j)$ must contain the function value at $xval (j)$ , for $j = 1, 2, \dots, 21$ .
3: $der (14)$ – Real (Kind=nag_wp) arrayOutput: On exit: the 14 derivative estimates.
4: $erest (14)$ – Real (Kind=nag_wp) arrayOutput: On exit: the 14 error estimates for the derivatives.
5: $ifail$ – IntegerInput/Output: On entry: ifail must be set to $0$ , $- 1 or 1$ . If you are unfamiliar with this argument you should refer to Section 3.4 in How to Use the NAG Library and its Documentation 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, the values of xval are not correctly spaced.
Derived $h = 〈value〉$ .

The derived $h$ is below tolerance.
Derived $h > 〈value〉$ is required. Derived $h = 〈value〉$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.
See Section 3.9 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.
See Section 3.8 in How to Use the NAG Library and its Documentation for further information.

$ifail = - 999$: Dynamic memory allocation failed.
See Section 3.7 in How to Use the NAG Library and its Documentation for further information.

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 d04baf 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 step length

h

(see Section 9). The only hard and fast rule is that for a given objective function and

h

, 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

d04baf is not threaded in any implementation.

9

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 round-off error) and as

h

becomes large (because the discretization error also becomes large). If this routine 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 d04bbf, then they will be in ascending order. In particular, the target abscissa

x_{0}

will be stored in

xval (11)

10

Example

This example evaluates the derivatives of the polygamma function, calculated using s14aef, and compares the first

3

derivatives calculated to those found using s14aef.

10.1

Program Text

Program Text (d04bafe.f90)

10.2

Program Data

None.

10.3

Program Results

Program Results (d04bafe.r)

d04baf (PDF version)

D04 (numdiff) Chapter Contents

D04 (numdiff) Chapter Introduction

NAG Library Manual

NAG Library Routine Document

d04baf (rcomm)

▸▿ Contents

1 Purpose

2 Specification

3 Description

4 References

5 Arguments

6 Error Indicators and Warnings

7 Accuracy

8 Parallelism and Performance

9 Further Comments

10 Example

10.1 Program Text

10.2 Program Data

10.3 Program Results

1

Purpose

2

Specification

3

Description

4

References

5

Arguments

6

Error Indicators and Warnings

7

Accuracy

8

Parallelism and Performance

9

Further Comments

10

Example

10.1

Program Text

10.2

Program Data

10.3

Program Results