Numerical Differentiation

Finite-difference formulas estimate derivatives from function samples, balancing truncation error that falls with step size against rounding error that rises.

Approximating slopes from samples

A derivative is a limit of difference quotients, and numerical differentiation keeps the step h finite. The forward difference (f(x+h)-f(x))/h has error proportional to h. The central difference (f(x+h)-f(x-h))/(2h) has error proportional to h^2 and is usually preferred for its higher accuracy at the same cost.

Higher-order formulas

Combining more sample points cancels more Taylor terms, raising the order of accuracy. The second derivative is commonly estimated by (f(x+h) - 2f(x) + f(x-h))/h^2, accurate to order h^2. These stencils are the building blocks of finite-difference methods for differential equations.

The competing errors

Numerical differentiation faces a fundamental tension. Truncation error shrinks as h decreases, but rounding error grows because subtracting nearly equal function values loses significance and is then divided by a tiny h. The total error is minimized at an intermediate h, roughly the square root of machine epsilon for central differences.

import math
def central_diff(f, x, h):
    return (f(x+h) - f(x-h)) / (2*h)
f = math.sin
for h in (1e-1, 1e-4, 1e-8, 1e-12):
    print(h, abs(central_diff(f, 1.0, h) - math.cos(1.0)))

Better alternatives

Where accuracy is critical, alternatives avoid the subtraction problem. Richardson extrapolation combines estimates at two step sizes to cancel leading error. The complex-step derivative Imag(f(x+ih))/h has no subtractive cancellation and reaches machine precision. Automatic differentiation computes exact derivatives from code directly.

Reliable derivatives feed gradient-based optimization and Jacobians for implicit solvers throughout physics simulation, including the sensitivity studies conducted for the breeder Hyperion.