What numerical differentiation is for
You differentiate numerically when you cannot differentiate symbolically. That happens in three common situations: the function only exists as a table of measured values, the function is the output of a simulation you cannot see inside, or the expression is differentiable in principle but so ugly that a finite difference is faster to code and easier to trust than a hand-derived formula. Every optimiser, every implicit ODE solver and every sensitivity analysis in engineering practice contains a finite difference somewhere.
The idea is the definition of the derivative with the limit removed. Instead of lim(h→0) [f(x+h) − f(x)]/h, you stop at a specific small h and accept the resulting error. What separates a numerical-methods answer from a guess is that you can say how large that error is and how it behaves as you change h.
Numerical differentiation is the awkward member of the numerical-methods family, and it is worth knowing why. Integration is a smoothing operation, so errors in the data average out; differentiation is a roughening operation, so errors in the data get amplified by 1/h. Halving the step improves the truncation error but doubles the amplification of noise. That tension is the whole subject, and the step-size sweep in the table above shows it happening.
Where the four formulas come from
Every finite-difference formula is a linear combination of Taylor expansions chosen so that the terms you do not want cancel. Expand f(x+h) = f(x) + h f'(x) + h²f''(x)/2 + h³f'''(x)/6 + …. Rearranging for f'(x) gives the forward difference with a leading error term of −h f''(x)/2: first order, because halving h halves the error.
Now also expand f(x−h) = f(x) − h f'(x) + h²f''(x)/2 − h³f'''(x)/6 + … and subtract the two. The f(x) terms cancel, and so does the whole h²f'' term, because it appears with the same sign in both. Dividing by 2h leaves the central difference with a leading error of −h²f'''(x)/6: second order, so halving h divides the error by four. That extra order costs nothing — it uses the same two function evaluations as a forward difference plus one — which is why central differencing is the default everywhere it is available.
The five-point stencil repeats the trick with a second pair of nodes at x ± 2h. Choosing the weights −1, 8, −8, 1 over 12h cancels the h² and h³ terms as well, leaving a leading error of h⁴f⁽⁵⁾(x)/30. Because the error term involves the fifth derivative, the formula is exact for any polynomial of degree four or lower — which is the sharpest test you can run on an implementation, and one of the test vectors on this page.
Backward differencing is the mirror of forward. It exists because at the last row of a data table there is no f(x+h) to use, exactly as forward differencing exists for the first row. Whenever you have neighbours on both sides, prefer central.
Worked example: f(x) = x³ at x = 2 with h = 0.1
The exact derivative is f'(x) = 3x², so f'(2) = 12. Work each scheme by hand.
- Sample the function. f(1.8) = 5.832, f(1.9) = 6.859, f(2.0) = 8, f(2.1) = 9.261, f(2.2) = 10.648.
- Forward. (9.261 − 8)/0.1 = 1.261/0.1 = 12.61. The error is +0.61. Compare the predicted leading term h f''(x)/2 = 0.1 × 12 / 2 = 0.6 — the prediction accounts for almost all of it.
- Backward. (8 − 6.859)/0.1 = 1.141/0.1 = 11.41. The error is −0.59, the same size and the opposite sign, which is exactly why averaging the two helps.
- Central. (9.261 − 6.859)/(2 × 0.1) = 2.402/0.2 = 12.01. The error has dropped from 0.61 to 0.01, a factor of 61. The predicted term is h²f'''(x)/6 = 0.01 × 6 / 6 = 0.01, matching exactly.
- Five-point. [−10.648 + 8(9.261) − 8(6.859) + 5.832] / (12 × 0.1). The bracket is −10.648 + 74.088 − 54.872 + 5.832 = 14.4, and 14.4/1.2 = 12 exactly, because the error term uses the fifth derivative of a cubic, which is zero.
Now halve the step to h = 0.05 and redo the central difference: f(2.05) = 8.615125, f(1.95) = 7.414875, and (8.615125 − 7.414875)/0.1 = 1.20025/0.1 = 12.0025. The error fell from 0.01 to 0.0025, a factor of four, which is the second-order behaviour the formula promises.
Choosing h, and reading the sweep table
The total error of a finite difference is the sum of two competing pieces. Truncation error falls like hᵖ. Round-off error grows like ε|f|/h, where ε ≈ 2.2 × 10⁻¹⁶ is double-precision machine epsilon, because subtracting two nearly equal function values destroys significant digits and then you divide by a small number. Minimising the sum gives an optimal step of roughly ε^(1/2) ≈ 10⁻⁸ for a first-order scheme, ε^(1/3) ≈ 6 × 10⁻⁶ for a central difference, and ε^(1/5) ≈ 7 × 10⁻⁴ for the five-point stencil, each scaled by the size of x.
That is the single most useful piece of practical knowledge here, and it surprises people: for a central difference, h = 10⁻¹² is far worse than h = 10⁻⁵. Smaller is not better past the optimum. The sweep table makes the crossover visible — errors fall by the predicted factor per row while truncation dominates, then stall or climb when cancellation takes over.
The best achievable accuracy also differs by scheme. At its optimum a forward difference gets about half the available digits, a central difference about two thirds, and the five-point stencil about four fifths. If you need a derivative to ten significant figures from a smooth function, only a high-order stencil or Richardson extrapolation will reach it — which is exactly what the reference value on this page uses.
With measured data the arithmetic is the same but the constants are not. Your noise floor is set by the instrument, not by machine epsilon, so if readings are good to ±0.001 and the spacing is 0.01, a central difference inherits an uncertainty near 0.001/0.01 = 0.1 in the slope regardless of how smooth the underlying physics is. In that regime, fitting a low-order polynomial over several points and differentiating the fit beats differencing adjacent points.
Reference: the four schemes at a glance
| Scheme | Nodes used | Weights ÷ denominator | Leading error term | Order | Evaluations |
|---|---|---|---|---|---|
| Forward | x, x+h | (−1, 1) ÷ h | −h f''(x)/2 | 1 | 2 |
| Backward | x−h, x | (−1, 1) ÷ h | +h f''(x)/2 | 1 | 2 |
| Central | x−h, x+h | (−1, 1) ÷ 2h | −h² f'''(x)/6 | 2 | 2 |
| Five-point | x±h, x±2h | (1, −8, 8, −1) ÷ 12h | +h⁴ f⁽⁵⁾(x)/30 | 4 | 4 |
The forward and backward leading terms differ only in sign, which is why their average is the central difference and why the second-order term disappears from it.
Pitfalls and things this calculator does not do
- Chasing accuracy by shrinking h. Below the optimum the error rises. If your answer changes wildly as you shrink the step, you are looking at cancellation, not convergence.
- Differencing across a kink. Taylor's theorem needs the derivatives that appear in the error term to exist. At x = 0 the function |x| has no second derivative, so no order claim holds and a central difference returns 0 for a function whose one-sided slopes are −1 and +1.
- Using a central difference at the edge of a data table. There is no node beyond the end; use the forward or backward form, and accept that the endpoint slope is the least accurate value in the whole table.
- Assuming the step is yours to choose. With tabulated data h is fixed by the sampling interval. Then the only lever left is a higher-order stencil or a smoothing fit.
- Ignoring the scale of x. A step of 0.001 is small near x = 1 and enormous near x = 10⁻⁶. All the optimal-step rules above scale with the magnitude of x.
- Reading the reference as exact. It is itself a numerical estimate, good to roughly ten significant figures for a smooth function. Differences at that level are the reference's own error.
Richardson extrapolation in one line
If D(h) is a second-order estimate, then [4D(h/2) − D(h)]/3 cancels the h² term and gives a fourth-order estimate for the price of two extra evaluations. Applying it repeatedly is Romberg's method, and it is how the reference value here is produced. The same idea applied to the trapezoid rule is the fastest route to a high-accuracy integral — see the Riemann sum calculator for where that family of rules starts.
Related methods and where to go next
Finite differences underpin more of numerical analysis than any other single idea. A root-finder that uses a numerical derivative instead of an analytic one is the secant method, a cousin of Newton's; when even that is unreliable, the fallback is the unconditionally convergent bisection method, which needs no derivative at all. Replacing the derivatives in a differential equation by differences is the finite-difference method for PDEs, and the same Taylor cancellation argument sets its order of accuracy.
In several variables the same formulas are applied one coordinate at a time. Doing that for a scalar field gives the gradient; combining the partial derivatives of a vector field gives the divergence and the curl. A gradient computed by finite differences costs n + 1 function evaluations in n dimensions, which is why large-scale optimisation reaches for automatic differentiation or an adjoint method instead.
If what you want is not a slope at a point but the location of a slope-zero, use the critical points calculator, which differentiates polynomials symbolically and so avoids the noise issue entirely. And when a limit is genuinely indeterminate rather than merely awkward, differentiating numerator and denominator is a rule rather than an approximation — that is the L'Hôpital's rule calculator.
One assumption is worth restating: everything here presumes f is smooth near x to at least the order of the error term you are claiming. When it is not — a discontinuity, a corner, a jump in a data set — no finite difference is meaningful, and the honest answer is that the derivative does not exist there.