D401 Numerical Differentiation

Subroutine subprograms RDERIV and DDERIV compute an approximate numerical value of the derivative $f\text{'}(x)$ of a function $f(x)$ at a specified point x.

On computers other than CDC and Cray, only the double-precision version DDERIV is available. On CDC and Cray computers, only the single-precision version RDERIV is available.

Structure:

SUBROUTINE subprograms
User Entry Names : RDERIV, DDERIV
Obsolete User Entry Names: DERIV $\equiv$ RDERIV
Files Referenced : Unit 6
External References: MTLMTR (N002), ABEND (Z035), user-supplied FUNCTION subprogram

Usage:

For $t=R$ (type REAL), $t=D$ (type DOUBLE PRECISION),

    CALL tDERIV(F,X,DELTA,DFDX,RERR)

F

(type according to t) Name of a user-supplied FUNCTION subprogram, declared EXTERNAL in the calling program. This subprogram must set $F(X)=\; f(X)$ .

X

(type according to t) The specified point x for which the derivative is to be calculated.

DELTA

(type according to t) On entry, DELTA must contain a scaling factor $\delta$ , which can usually be set equal to 1. On exit, it contains the last value of this factor (see Method).

DFDX

(type according to t) On exit, DFDX contains an approximation to $f\text{'}(X)$ .

RERR

(type according to t) On exit, RERR contains an estimate of the relative error of DFDX.

Method:

The method is based on an extension to numerical differentiation of Romberg's principle of sequence extrapolation, originally developed for numerical integration. The subroutine starts by computing the 10 numbers

$T$₀^(k)= [f(x+h)-f(x-h)]/(2h), &quad;(k = 0,1,...,9),

with

$h\; =\; 0.0256*2-k/2$	$(k\; even)$
$h\; =\; 0.0192*2-(k-1)/2$	$(k\; odd),$

where the scaling factor $\delta$ is initially set to DELTA. This procedure is repeated up to 9 times, with $\delta$ replaced by $\delta /10$ each time, until the sequence $T$₀^(k) is found to be monotone. A Romberg-like triangular table

$T$_m^(k)= u_mT_m-1^(k+1)- w_mT_m-1^(k).

with appropriate weights $u$_m,w_m is then computed for $m\; =\; 1,2,...,9;\; k\; =\; 0,1,...,9-m,$ and DFDX is set equal to $T$₉⁽⁰⁾ .

Restrictions:

The function $f(x)$ must be differentiable at $x=X$ and in a neighbourhood of X. Misleading results will be obtained if this is not true.

Error handling:

Error D401.1: If the function $f(x)$ is such that, after 9 successive reductions of $\delta$ by a factor 1/10, the sequence $T$₀^(k) is not monotone, an error message is written Unit 6, unless subroutine MTLSET (N002) has been called. DFDX is set equal to zero, RERR is set equal to one in this case.

References:

1. H. Rutishauser, Ausdehnung des Rombergschen Prinzips, Numer. Math. 5 (1963) 48--54.

D501