Finite Differences

Derivatives without differentiating: slope from neighbouring values. Central beats one-sided, and shrinking the step eventually backfires.

On x cubed at 2 the forward secant overshoots at 12.61 while the symmetric central secant lands at 12.01, nearly on the true tangent of 12

Why Does This Exist?

Autodiff is unavailable more often than textbooks admit: black-box simulators, legacy Fortran, hardware in the loop, quick sanity checks of hand-derived gradients. Finite differences estimate $f'(x)$ from neighbouring values of $f$ alone, no symbolic work. Forward, backward and central formulas trade one extra evaluation for an order of accuracy, and the central difference's $O(h^2)$ error makes it the default checker behind every gradient verification routine.

The method is also the seed of numerical differential equations: replace derivatives with differences and a differential equation becomes algebra on a grid. Derivatives and partial derivatives define the target; Taylor series supply the error analysis; gradients shows the checking workflow this enables.

The one idea you need first

Taylor's expansion $f(x+h) = f(x) + hf'(x) + h^2 f''(x)/2 + \cdots$ rearranges into slope estimates: each formula below is Taylor solved for $f'(x)$ with the remainder kept as the error term. The order of the first surviving remainder term names the method's accuracy.

Think of It Like This

Speed from two photographs

A car's speed at noon cannot be photographed, but two timestamped photos give average speed over the gap. Both photos after noon (forward difference) measures where the car is heading, including its acceleration over the gap: biased when speed changes. One photo before and one after (central difference) straddles noon, and the acceleration contaminates both sides nearly equally, cancelling most of the bias.

That cancellation is the whole accuracy gap: symmetric sampling kills the leading error term for free. The price is one extra "past" evaluation, and the limit is photo resolution: timestamps closer than the clock's precision subtract into noise.

How It Actually Works

The three formulas with their errors:

Forward: $(f(x+h) - f(x))/h = f'(x) + O(h)$ . One-sided, first-order. The only option at domain boundaries.
Backward: $(f(x) - f(x-h))/h = f'(x) + O(h)$ . Mirror image, same accuracy.
Central: $(f(x+h) - f(x-h))/(2h) = f'(x) + O(h^2)$ . Symmetric, second-order: the $h$ terms cancel, leaving $h^2$ .

Halving $h$ halves forward error but quarters central error. Total error is U-shaped: truncation falls as $h$ shrinks while rounding (catastrophic cancellation in the subtraction) rises, so floats bottom out near $h \approx 10^{-5}$ for forward and $10^{-4}$ -ish for central in doubles. Smaller is not better.

Worked example

Differentiate $f(x) = x^3$ at $x = 2$ (truth $f'(2) = 12$ ) with $h = 0.1$ . Values: $f(2.1) = 9.261$ , $f(1.9) = 6.859$ .

Forward: $(9.261 - 8)/0.1 = 12.61$ . Error $0.61$ .
Central: $(9.261 - 6.859)/0.2 = 2.402/0.2 = 12.01$ . Error $0.01$ .

Same step, sixty times less error, one extra evaluation. The $O(h)$ versus $O(h^2)$ labels predict exactly this ratio: $h = 0.1$ against $h^2 = 0.01$ .

Code

The point is the check: confirm the accuracy gap and the U-shaped failure at tiny $h$ .

1def diffs(f: list[float], h: float) -> tuple[float, float]:2    fwd = (f[1] - f[0]) / h  # forward needs f(x), f(x+h)3    cen = (f[2] - f[0]) / (2 * h)  # central needs f(x-h), f(x+h)4    return (fwd, cen)5
6
7vals = [8.0, 9.261, 6.859]  # f(2), f(2.1), f(1.9) for x cubed8print([round(v, 2) for v in diffs(vals, 0.1)])  # -> [12.61, 12.01]

Watch Out For

Shrinking h past the rounding floor

Below $h \approx 10^{-8}$ in doubles, $f(x+h)$ and $f(x-h)$ share their first eight digits and the subtraction keeps mostly rounding noise, amplified by division with $2h$ . The check then fails on correct gradients and misleads the debugging that follows. Park forward checks at $h = 10^{-5}$ , central at $10^{-4}$ , and compare relatively.

Checking in float32

Single precision carries 7 digits, so the rounding floor arrives a thousand times sooner and no $h$ is simultaneously small enough for truncation and large enough for rounding. Gradient checks belong in float64, period. Cast the function up for the check even when production runs in float32 or below.

The Quick Version

Forward and backward differences err at $O(h)$ ; central differences at $O(h^2)$ for one extra evaluation.
On $x^3$ at 2 with $h = 0.1$ : forward 12.61, central 12.01, truth 12.
Total error is U-shaped in $h$ : truncation versus cancellation, bottoming near $10^{-5}$ .
Never gradient-check in float32; cast up to float64 for the check.
Symmetric sampling cancels the leading error: straddle the point whenever values exist both sides.