Numerical differentiation

The formula

Forward difference:

f'(x) ≈ ( f(x+h) − f(x) ) / h          error ≈ (h/2)·f''(x)

Reading: slope of the line through this sample and the next one. The error shrinks in proportion to the step h — halve the step, halve the error.

Central difference:

f'(x) ≈ ( f(x+h) − f(x−h) ) / (2h)     error ≈ (h²/6)·f'''(x)

Reading: slope of the line through the neighbors on either side, skipping the point itself. The error shrinks with h² — halve the step, quarter the error. Same two function evaluations as forward difference, one order better, which is why central is the default whenever you have data on both sides.

Second derivative, central:

f''(x) ≈ ( f(x+h) − 2·f(x) + f(x−h) ) / h²     error ≈ (h²/12)·f⁗(x)

Reading: how much the middle point sags below the average of its neighbors, scaled by the step squared. This three-point stencil is the beating heart of nearly every finite-difference PDE solver ever written.

The one that governs everything else — the total error of a computed difference:

E(h) ≈ C·h^p  +  ε/h

Reading: truncation error (first term) falls as the step shrinks; noise and roundoff error (second term) grow as the step shrinks, because you are dividing an increasingly uncertain difference by an increasingly small number. There is a floor, and it is not at h = 0.

Where you meet it

The test stand, computing a rise rate. Chamber pressure sampled at 1 kHz, and the customer wants dP/dt during startup. Take a central difference — subtract the neighbors on either side of each sample and divide by the two milliseconds between them — and 1 mV of channel noise becomes roughly 0.7 V/s of derivative noise; difference each sample from the last instead and it is twice that. Either way, a clean-looking trace turns into grass. Every data-reduction engineer learns this the first time they plot a raw finite-difference rate.

Velocity and acceleration from position. Radar tracks, encoder counts, photogrammetry, GPS — position is what gets measured, and rates are what the analysis needs. Differencing position samples is the obvious move; differencing them twice for acceleration is where the noise catch stops being theoretical and starts being a review-board finding.

Gradients for an optimizer. A trade study wraps an optimizer around a CFD or FEA run that has no analytic derivatives. The optimizer perturbs each design variable by h and finite-differences the objective. Pick h too large and the gradient is biased; too small and solver noise swamps it, and the optimizer wanders. Half of "the optimizer won't converge" tickets are really step-size tickets.

The D in PID. A derivative term in a digital controller is a finite difference of the error signal, executed every sample period. This is why real controllers ship the D term with a low-pass filter bolted on, and why an unfiltered derivative gain on a noisy process variable makes an actuator chatter itself to death.

How it works

Everything comes from Taylor's series. Expand f(x+h) about x, solve for f'(x), and the terms you drop are the truncation error. The forward difference drops a term proportional to h·f''; the central difference is a forward and a backward expansion subtracted, so the even-order terms cancel and the leading error rides on h²·f'''. The improvement is not cosmetic. Differentiating eˣ at x = 1 with h = 0.01: forward difference misses by 1.4·10⁻², central by 4.5·10⁻⁵ — three hundred times better for the same arithmetic. Wider stencils continue the pattern; the five-point central rule buys h⁴ accuracy for two more evaluations.

So shrink h to zero and win? No — and this is the trap. The difference f(x+h) − f(x−h) subtracts two nearly equal numbers, so whatever uncertainty each carries — roundoff at the sixteenth digit, or sensor noise at the third — survives the subtraction and then gets divided by 2h. In double precision the total error curve is a V: for the central difference it bottoms out near h ≈ ε^(1/3) — around 10⁻⁵ for machine epsilon ε ≈ 2.2·10⁻¹⁶ — and the best error you can buy is on the order of ε^(2/3), about 10⁻¹¹, no matter how much smaller you make the step ^[1]. Run the experiment and you can watch it happen: the error falls cleanly to about 3·10⁻¹¹ near h = 10⁻⁵·⁵, then climbs back up and hits order-one garbage by h = 10⁻¹⁶. Eleven good digits is the ceiling. There is no such thing as "let h go to zero" on a computer.

With measured data the same V exists, but the floor is far higher because ε is now sensor noise, not roundoff. Independent noise of standard deviation σ per sample passes through a central difference as σ/(√2·h) of derivative noise — the 1 kHz pressure example above, worked exactly. A sample-to-sample first difference is worse, √2·σ/h, because it divides by h instead of 2h; that factor of two is the "twice that" in the test-stand example. The frequency-domain view says the same thing in one line: differentiation has gain ω, so it amplifies a 100 Hz noise component a hundred times more than a 1 Hz signal component. Differentiation is a high-boost filter applied to data whose junk lives at high frequency.

The mistake people make follows directly: sampling faster to get a "better" derivative. Faster sampling shrinks h, and with fixed per-sample noise that makes the raw difference worse, not better. The fixes all amount to trading bandwidth for noise, deliberately. Low-pass filter before differencing, and difference across a wider stride than one sample. Or fit a low-order polynomial through a sliding window and differentiate the fit — the Savitzky–Golay approach, which is what your scope's "smoothed derivative" button is doing. Or, offline, fit a spline and differentiate that. For gradients on smooth code rather than noisy data, two sharper tools exist: Richardson extrapolation, which combines the answers at h and h/2 to cancel the leading error term, and the complex-step method — f'(x) ≈ Im( f(x + i·h) )/h — which involves no subtraction at all, so h can be 10⁻²⁰ and the derivative comes out at machine precision ^[1]. Complex step only works on code you can run in complex arithmetic, which is exactly the case in optimization wrappers, and it is why modern aircraft MDO frameworks lean on it.

Limits of validity worth keeping straight: the error terms assume the derivative they ride on exists and is bounded — difference across a discontinuity, a table lookup boundary, or a solver's own convergence noise and the formulas mean nothing. And a central difference needs both neighbors, so at the edges of a record you are forced to one-sided formulas with their bigger error; the standard three-point one-sided rule (−3f(x) + 4f(x+h) − f(x+2h))/(2h) recovers h² accuracy at a boundary.

History

The machinery is old — older than most of the calculus notation we use to describe it. Isaac Newton and James Gregory each built interpolation formulas on tables of forward differences in the late 1600s, and Newton put a compressed account into the Principia in 1687; the Gregory–Newton forward difference formula still carries both names ^[5]^[6]. Brook Taylor's 1715 Methodus Incrementorum Directa et Inversa — the same book that stated Taylor's theorem — organized the subject into what is now called the calculus of finite differences ^[3]^[4].

The twentieth century turned it from a hand-computation trick into a numerical method with an error theory. Lewis Fry Richardson used finite differences to attack physical problems — stresses in a masonry dam — in a 1911 paper, and in 1927, with J. A. Gaunt, published "the deferred approach to the limit": compute at two step sizes and extrapolate the error away, the scheme now called Richardson extrapolation ^[7]^[8]. In 1964 Abraham Savitzky and Marcel Golay, working on noisy spectrometer data, published sliding least-squares smoothing and differentiation in Analytical Chemistry; it became one of the most-cited papers the journal ever ran, because every instrument maker had the same noise-amplification problem ^[9]^[10]. The complex-step idea traces to Lyness and Moler in 1967 and was revived as a practical recipe by Squire and Trapp in a two-page 1998 SIAM Review note ^[11]^[12] — proof that in this field a good idea can wait thirty years for the hardware and the audience to catch up.

Related tools

/tools/rc-filter
/tools/adc-resolution
/tools/snr-enob
/tools/strain-gauge-bridge

Sources

Written by HE in our own words from the cited sources — engineering judgment included, your stamp still required. All entries →