Numerical Differentiation

Numerical differentiation estimates $f^{\prime }(x)$ from evaluations of $f$. It is useful when an analytic derivative is unavailable, but it is more sensitive to floating-point rounding than evaluating $f$ itself.

Choosing a method

Use an analytic or automatic derivative when available; use complex-step for complex-safe analytic functions; use a central difference for a simple robust estimate; and use Richardson extrapolation when higher accuracy justifies more evaluations.

Forward difference

The derivative is defined by

$$f^{\prime }(x)=\underset{h\to 0}{\lim }\frac{f(x+h)-f(x)}{h}.$$

For a small nonzero step $h$, use the forward-difference approximation

$$D_{+}(h)=\frac{f(x+h)-f(x)}{h}.$$

Taylor expansion gives

$$f(x+h)=f(x)+f^{\prime }(x)h+\frac{1}{2}{f}^{\prime \prime }(x)h^2+O(h^3),$$

therefore

$$D_{+}(h)-f^{\prime }(x)=\frac{1}{2}{f}^{\prime \prime }(x)h+O(h^2).$$

Its truncation error is $O(h)$: halving $h$ approximately halves this error.

Smaller $h$ is not always better

When $h$ becomes very small, $f(x+h)$ and $f(x)$ are nearly equal. Their subtraction loses significant digits, and division by $h$ amplifies rounding error.

import numpy as np
def forward_difference(f, x, h):
    return (f(x + h) - f(x)) / h

Complex-step differentiation

If $f$ has an analytic extension to complex arguments, evaluate it at a purely imaginary displacement:

$$f(x+ih)=f(x)+ih{f}^{\prime }(x)-\frac{1}{2}{h}^2{f}^{\prime \prime }(x)-\frac{i}{6}{h}^3{f}^{\prime \prime \prime }(x)+O(h^4).$$

Taking the imaginary part yields

$$\frac{\mathrm{Im}(f(x+ih))}{h}=f^{\prime }(x)-\frac{1}{6}{f}^{\prime \prime \prime }(x)h^2+O(h^4),$$

so $f^{\prime }(x)\approx \mathrm{Im}(f(x+ih))/h$. There is no subtraction of nearly equal numbers, so $h$ can often be extremely small (for example $10^{-20}$ in double precision), giving near machine-precision derivatives.

Applicability

Complex-step differentiation requires a complex-safe, analytic implementation of $f$. It can fail for abs, comparisons or branching based on the input, max, min, and code that discards imaginary parts.

def complex_step(f, x, h=1e-20):
    return np.imag(f(x + 1j * h)) / h

Richardson extrapolation

Central differences have an error expansion containing only even powers of $h$. Write the unknown derivative as $L=f^{\prime }(x)$ and the base approximation as $D(h)=D_0(h)$:

$$D(h)=\frac{f(x+h)-f(x-h)}{2h}=L+a_1{h}^2+a_2{h}^4+a_3{h}^6+\cdots .$$

Using half the step changes each error term predictably:

$$\begin{array}{r}D(h/2)=L+a_1\frac{h^2}{4}+a_2\frac{h^4}{16}+a_3\frac{h^6}{64}+\cdots \\ 4D(h/2)-D(h)=3L-\frac{3}{4}{a}_2{h}^4-\frac{15}{16}{a}_3{h}^6+\cdots \\ R_{1,1}=\frac{4D(h/2)-D(h)}{3}=L-\frac{1}{4}{a}_2{h}^4+O(h^6)\end{array}$$

The leading error has improved from $O(h^2)$ to $O(h^4)$. This is Richardson extrapolation: combine two estimates with known error scaling to eliminate the leading error term.

For the next step, start with $R_{1,1}(h)$ and $R_{1,1}(h/2)$. Their leading errors scale as $h^4$ and $(h/2{)}^4=h^4/16$, so the same cancellation uses $16=4^2$:

$$R_{2,2}=\frac{16R_{2,1}-R_{1,1}}{15},$$

where $R_{2,1}$ is the first-order extrapolation built from $D(h/2)$ and $D(h/4)$. The result has error $O(h^6)$.

In general, define the first column by

$$R_{\ell ,0}=D\left(\frac{h}{2^{\ell }}\right),\ell =0,1,2,\dots$$

and use the recurrence

$$R_{\ell ,k}=\frac{4^k{R}_{\ell ,k-1}-R_{\ell -1,k-1}}{4^k-1},1\le k\le \ell .$$

Here $k$ is the extrapolation level. Since $R_{\ell ,k-1}$ has leading error proportional to $(h/2^{\ell }{)}^{2k}$, its error differs by a factor $4^k$ from the corresponding entry in the preceding row. The recurrence removes that term, so $R_{\ell ,k}$ has error $O((h/2^{\ell }{)}^{2k+2})$.

The first entries form a triangular table:

$\ell$	$k=0$	$k=1$	$k=2$
$0$	$R_{0,0}=D(h)$
$1$	$R_{1,0}=D(h/2)$	$R_{1,1}=\frac{4R_{1,0}-R_{0,0}}{3}$
$2$	$R_{2,0}=D(h/4)$	$R_{2,1}=\frac{4R_{2,0}-R_{1,0}}{3}$	$R_{2,2}=\frac{16R_{2,1}-R_{1,1}}{15}$

Rounding-error limit

Extrapolation removes truncation error, but increasingly small steps eventually amplify floating-point error.

def richardson_derivative(f, x, h0=1.0, rtol=1e-12, atol=1e-12, max_levels=30):
    """Central differences accelerated by Richardson extrapolation."""
    estimates = []
    for level in range(max_levels):
        h = h0 / (2 ** level)
        row = [central_difference(f, x, h)]
        for k in range(1, level + 1):
            factor = 4 ** k
            row.append((factor * row[k - 1] - estimates[level - 1][k - 1]) / (factor - 1))
        estimates.append(row)
        if level and abs(row[-1] - estimates[level - 1][-1]) < max(rtol * abs(row[-1]), atol):
            return row[-1], h
    return estimates[-1][-1], h0 / (2 ** (max_levels - 1))