Prev
Unit 8 · Numerical Differentiation and Integration
Finite Differences and Richardson Extrapolation
55m100 XP
Next
∫Numerical Differentiation and Integration

Finite Differences and Richardson Extrapolation

Difference formulas from Taylor series, the step-size trade-off, noisy data, Richardson extrapolation, complex step

A static fire, one sample per millisecond
A solid rocket motor is bolted horizontally to a test stand. When it fires, a load cell records thrust and a transducer records chamber pressure, a thousand times a second, for three seconds. The engineers want more than the curves. The burning rate of the propellant depends on how fast pressure rises; ignition transients are judged by dp/dtdp/dt; nozzle erosion shows up as a slow drift in the ratio of thrust to pressure. All of these are derivatives — of numbers sampled every millisecond, each carrying a little electrical noise. The first attempt is always the same: subtract consecutive samples and divide by the time step. On a pressure signal with noise of 0.2% of full scale, the result is a derivative whose noise is ten times larger than the derivative itself. Nothing was programmed wrongly. Differentiation of data is simply one of the most ill-conditioned operations in this course.
Every difference formula comes from Taylor's theorem (lesson 1-3). Expanding f(x±h)f(x \pm h) about xx:
f(x+h)=f+hf′+h22f′′+h36f′′′+…,f(x−h)=f−hf′+h22f′′−h36f′′′+…f(x+h) = f + hf' + \tfrac{h^2}{2}f'' + \tfrac{h^3}{6}f''' + \dots, \qquad f(x-h) = f - hf' + \tfrac{h^2}{2}f'' - \tfrac{h^3}{6}f''' + \dots
Subtracting cancels the even terms. Combining more points cancels more terms:
f(x+h)−f(x)h=f′+O(h),f(x+h)−f(x−h)2h=f′+O(h2),−f2+8f1−8f−1+f−212h=f′+O(h4)\frac{f(x+h)-f(x)}{h} = f' + O(h),\qquad \frac{f(x+h)-f(x-h)}{2h} = f' + O(h^2),\qquad \frac{-f_{2}+8f_{1}-8f_{-1}+f_{-2}}{12h} = f' + O(h^4)

Forward, central and fourth-order central differences (fₖ = f(x + kh)). Second derivative: (f(x+h) − 2f(x) + f(x−h))/h² = f″ + O(h²) — the stencil of every heat and beam equation in Units 3 and 4.

💡The step-size trade-off
Truncation error shrinks as h→0h \to 0, but f(x+h)−f(x)f(x+h) - f(x) is a difference of nearly equal numbers, and each carries a rounding error of about ε∣f∣\varepsilon|f|. Divided by hh, that becomes ε∣f∣/h\varepsilon|f|/h — growing as hh shrinks. For the forward difference the total error is about h2∣f′′∣+ε∣f∣/h\tfrac{h}{2}|f''| + \varepsilon|f|/h, smallest at h≈ε≈10−8h \approx \sqrt{\varepsilon} \approx 10^{-8} with about 8 correct digits. Central differences do better: h≈ε1/3≈6×10−6h \approx \varepsilon^{1/3} \approx 6 \times 10^{-6} and about 11 digits; the fourth-order formula, h≈10−3h \approx 10^{-3} and about 13. No finite-difference formula reaches full precision, and pushing hh smaller makes things worse.
Noisy data is the same problem, worse. Measurement noise σ\sigma plays the role of ε∣f∣\varepsilon|f|, and it is millions of times larger. A central difference turns noise σ\sigma into a derivative error of about σ/(2 h)\sigma/(\sqrt2\,h): with σ=0.014\sigma = 0.014 MPa and h=1h = 1 ms, about 10 MPa/s — while the true dp/dtdp/dt during the burn is under 1 MPa/s. The cure is to smooth while differentiating: fit a low-degree polynomial to a window of samples by least squares (Unit 5) and differentiate the fit. That is the Savitzky–Golay filter (1964); with a 51-sample cubic window the derivative error drops thirtyfold, to 0.33 MPa/s. Wider windows remove more noise but blur fast events — the step-size trade-off again, now with σ\sigma instead of ε\varepsilon.
Derivation
How much a difference amplifies noise
Let each sample carry an independent random error of standard deviation σ\sigma: f~k=fk+εk\tilde f_k = f_k + \varepsilon_k. The central difference adds the error
ε+1−ε−12h.\frac{\varepsilon_{+1} - \varepsilon_{-1}}{2h}.
Variances of independent errors add (lesson 1-2's random walk), and a minus sign does not change a variance, so ε+1−ε−1\varepsilon_{+1} - \varepsilon_{-1} has variance 2σ22\sigma^2 and standard deviation 2 σ\sqrt2\,\sigma. Dividing by 2h2h:
noise in the derivative≈2 σ2h=σ2 h.\text{noise in the derivative} \approx \frac{\sqrt2\,\sigma}{2h} = \frac{\sigma}{\sqrt2\,h}.
With σ=0.014\sigma = 0.014 MPa and h=0.001h = 0.001 s that is 9.99.9 MPa/s. Halving the step doubles it — the opposite of what truncation error wants.
💡Richardson extrapolation
If an approximation D(h)D(h) has error c2h2+c4h4+…c_2h^2 + c_4h^4 + \dots, compute it twice: 4D(h/2)−D(h)3=f′+O(h4).\dfrac{4D(h/2) - D(h)}{3} = f' + O(h^4). The h2h^2 terms cancel without our ever knowing c2c_2. Lewis Fry Richardson called it 'the deferred approach to the limit' (1911, 1927). For Squire and Trapp's test function at x=1.5x = 1.5, central differences with h=0.05h = 0.05 are off by 1.3×10−21.3 \times 10^{-2}; one Richardson step brings that to 8×10−68 \times 10^{-6}. The same idea powers Romberg integration (lesson 8-3) and error estimation in ODE solvers (Unit 9).
Derivation
Richardson's combination, and the fourth-order stencil it produces
The central difference has an error expansion with only even powers (lesson 1-3 showed the odd terms cancel):
D(h)=f′+c2h2+c4h4+…,D(h/2)=f′+c2h24+c4h416+…D(h) = f' + c_2h^2 + c_4h^4 + \dots, \qquad D(h/2) = f' + \frac{c_2h^2}{4} + \frac{c_4h^4}{16} + \dots
Multiply the second by 4 and subtract the first: the c2c_2 terms cancel, and 4f′−f′=3f′4f' - f' = 3f' remains, so
4D(h/2)−D(h)3=f′−c4h44+…\frac{4D(h/2) - D(h)}{3} = f' - \frac{c_4h^4}{4} + \dots
Nothing about c2c_2 was needed except that the error goes like h2h^2. In general, for an error chpch^p, the combination is (2pD(h/2)−D(h))/(2p−1)\big(2^pD(h/2) - D(h)\big)/(2^p - 1); for a first-order error (p=1p = 1) it is 2D(h/2)−D(h)2D(h/2) - D(h), the extrapolation used for seepage in lesson 4-4.
The fourth-order formula. Apply the same combination to central differences with steps hh and 2h2h, written out with fk=f(x+kh)f_k = f(x + kh):
13[4 f1−f−12h−f2−f−24h]=−f2+8f1−8f−1+f−212h.\frac13\Big[4\,\frac{f_1 - f_{-1}}{2h} - \frac{f_2 - f_{-2}}{4h}\Big] = \frac{-f_2 + 8f_1 - 8f_{-1} + f_{-2}}{12h}.
That is the five-point stencil above: every higher-order difference formula is a Richardson extrapolation of a lower-order one.
They are neither finite Quantities nor Quantities infinitely small, nor yet nothing. May we not call them the Ghosts of departed Quantities?
– George Berkeley, The Analyst (1734)
🕰The complex-step trick
For an analytic function implemented with complex arithmetic, Taylor's theorem in the imaginary direction gives f(x+ih)=f(x)+ihf′(x)−h22f′′(x)−…f(x + ih) = f(x) + ihf'(x) - \tfrac{h^2}{2}f''(x) - \dots, so f′(x)≈Im⁡f(x+ih)/hf'(x) \approx \operatorname{Im} f(x + ih)/h with error O(h2)O(h^2) and no subtraction at all. There is nothing to cancel, so hh can be 10−2010^{-20} or 10−10010^{-100} and the derivative is correct to machine precision. James Lyness and Cleve Moler described it in 1967; William Squire and George Trapp popularised it in 1998, and it became a standard tool for computing sensitivities in aerospace design codes.
💡Where you will use this
Rocket Propulsion: burn rates and ignition transients from test data — always smoothed. Structures: strains from displacement measurements (digital image correlation, extensometers). Optimal Guidance and Control: finite-difference and complex-step Jacobians for Newton and gradient methods (Units 2 and 11). Every PDE course: the difference stencils of this lesson are the finite-difference method.
📚Further reading
  • Bengt Fornberg, 'Generation of finite difference formulas on arbitrarily spaced grids', Mathematics of Computation 51(184), 699–706 (1988) — a short algorithm that produces every stencil in this lesson, and more.
  • Abraham Savitzky and Marcel Golay, 'Smoothing and differentiation of data by simplified least squares procedures', Analytical Chemistry 36(8), 1627–1639 (1964).
  • Joaquim Martins, Peter Sturdza and Juan Alonso, 'The complex-step derivative approximation', ACM Transactions on Mathematical Software 29(3), 245–262 (2003) — the complex step in engineering design codes.
↓But there is a problem
Differentiation divides small differences by small numbers and amplifies every error. Its inverse does the opposite: integration adds many small contributions, and random errors partly cancel. The total impulse from the same noisy thrust record is accurate to a few hundredths of a percent. The oldest integration rule still in daily use was written down because of a wedding. In 1613 Johannes Kepler, astronomer to the Emperor, bought barrels of wine for his marriage in Linz — and did not like the way the merchant measured how much was in them.