Polynomial interpolation

The formula

The setup: n+1 points (x₀,y₀), (x₁,y₁), …, (xₙ,yₙ) with all x-values distinct. There is exactly one polynomial p(x) of degree ≤ n through them — no more, no fewer. That uniqueness is the whole foundation of lookup-table engineering: whatever route you take to the polynomial, you land on the same one.

The Lagrange form builds it directly:

p(x) = Σ yₖ·Lₖ(x),   where  Lₖ(x) = Π_{j≠k} (x − xⱼ) / (xₖ − xⱼ)

Reading: each basis polynomial Lₖ equals 1 at its own node and 0 at every other node, so each data point contributes its y-value exactly where it should and nothing anywhere else.

The Newton form builds the same polynomial one point at a time:

p(x) = y₀ + [y₀,y₁]·(x−x₀) + [y₀,y₁,y₂]·(x−x₀)(x−x₁) + …

Reading: each new term is a correction that vanishes at all the earlier nodes, so adding a calibration point never disturbs the coefficients you already computed. The bracketed quantities are divided differences, built recursively from [yᵢ,yⱼ] = (yⱼ − yᵢ)/(xⱼ − xᵢ).

The degree-1 case is the one you use weekly, whether you call it by name or not:

p(x) = y₀ + (y₁ − y₀)·(x − x₀)/(x₁ − x₀)

Reading: linear interpolation between two table entries is polynomial interpolation with n = 1.

And the error term, the part that earns its keep:

f(x) − p(x) = f⁽ⁿ⁺¹⁾(ξ)/(n+1)! · (x−x₀)(x−x₁)···(x−xₙ)

Reading: the miss between the true function and your polynomial rides on the (n+1)-th derivative somewhere in the interval, times the product of your distances to every node ^[8]^[9]. Curvier function or wider node spacing, bigger error — and the formula says exactly how much bigger. For linear interpolation on a table with spacing h, it collapses to the bound every instrumentation engineer should have taped to the monitor:

|error| ≤ h²·max|f″| / 8

Where you meet it

The thermocouple table on the bench. The ITS-90 type K reference function gives 0.7981 mV at 20 °C and 1.2033 mV at 30 °C ^[10]. Your DMM reads 1.000 mV. Linear interpolation puts the junction at 24.98 °C; the reference function says 24.99 °C. That 0.01 °C miss is the h²·max|f″|/8 term with h = 10 °C on a gently curved function — which is why nobody carries a finer table for type K near room temperature, and why the same move gets dangerous on a curve that actually bends.

The five-point cal on a test stand. A pressure transducer comes back from the cal lab with five points and you need readings everywhere in between. Piecewise-linear between the points, a low-order fit, or — the classic mistake — a degree-4 polynomial threaded exactly through all five. The theory in this entry is how you decide, and the "How it works" section is why the degree-4 option loses.

Compressible-flow and atmosphere tables. Isentropic ratios, normal-shock tables, the standard atmosphere: all published at discrete Mach numbers or altitudes, all read in between by interpolation. Near a shock table's steep region, the second derivative is large and the linear-interpolation bound tells you the table spacing is doing real work.

Sizing a lookup table in embedded code. A sine LUT with 1° spacing and linear interpolation between entries has a worst-case error of (π/180)²/8 ≈ 3.8·10⁻⁵ — about 14–15 bits of accuracy from a 91-entry quarter-wave table. The error bound turns "how many entries do I need" from guesswork into arithmetic before you've written a line of firmware.

How it works

Everything downstream of the uniqueness theorem is a fight between two terms in the error formula: the derivative factor f⁽ⁿ⁺¹⁾(ξ)/(n+1)! and the node polynomial Π(x − xᵢ). You control the second one — where the nodes sit and how many there are. You are stuck with the first.

The intuition that fails people is "more points, higher degree, better fit." Carl Runge killed that idea in 1901 with a function that looks harmless: f(x) = 1/(1 + 25x²) on [−1, 1]. Interpolate it at evenly spaced points and the polynomial oscillates harder near the ends as the degree climbs. At degree 10 the worst-case error is about 1.9 — nearly twice the function's own peak value. At degree 20 it is about 60. The fit passes through every data point perfectly and is garbage everywhere between the outer ones. The derivative factor grows faster than the factorial kills it, and the equally spaced node polynomial is largest exactly where the trouble starts, near the interval ends.

Move the same 21 nodes to Chebyshev spacing — clustered toward the ends, xₖ = cos((2k+1)π/(2n+2)) — and the degree-20 error drops from 60 to about 0.015. Chebyshev nodes minimize the worst case of the node polynomial, which is why every serious approximation library uses them. But no fixed node layout saves you for every function: for any prescribed table of nodes there exists a continuous function whose interpolants diverge. Node choice is mitigation, not immunity.

The engineering conclusion is the one the numerical analysts reached decades ago: past degree 3 or 4, stop raising the degree and start splitting the interval. Piecewise-linear tables and cubic splines are polynomial interpolation used at the degree where the error formula behaves — which is why your DAQ software offers splines and does not offer a degree-19 fit.

Two implementation notes worth knowing. First, the obvious way to find the coefficients — solve the Vandermonde system — is numerically nasty; the Vandermonde matrix's condition number grows fast with degree, so small data noise becomes large coefficient error. The Newton form sidesteps the solve entirely, and the barycentric variant of the Lagrange form evaluates the polynomial in O(n) per point with excellent stability — it is the form production code actually uses. Second, interpolation and least-squares fitting are different contracts. Interpolation passes through every point and therefore reproduces every point's noise exactly. If your calibration data has scatter, you want a fit with fewer coefficients than points, not a polynomial that treats the noise as gospel.

Extrapolation deserves its own warning. Outside the node interval, the node polynomial Π(x − xᵢ) grows without bound, and the error formula stops being a bound on anything you'd sign your name to. An interpolating polynomial one span past the last calibration point is a random number generator with good posture.

History

The subject is older than the name on it. James Gregory had the finite-difference interpolation formula by 1670, in his correspondence with John Collins ^[1]. Newton, working independently, made it a system: his general formula for unevenly spaced data appears in the Principia (1687, Book III, Lemma V), where he needed it to thread a smooth path through discrete observations of a comet, and the full treatment — written in the mid-1670s — was published in 1711 as the Methodus Differentialis ^[2]^[3]. Divided differences are his machinery, and the equal-spacing special case still carries both names as the Gregory–Newton formula ^[1]^[2].

The formula every textbook calls Lagrange's was published in 1779 by Edward Waring, Lucasian Professor at Cambridge, as a cleaner way to write Newton's result without computing divided differences at all ^[2]^[4]. Euler gave a related form in 1783, and Lagrange published the same formula in 1795, apparently unaware of both ^[2]^[5]. The wrong name stuck.

The subject's cold shower came in 1901, when Carl Runge published the example above showing that equally spaced interpolation of a smooth, bounded function can diverge as the degree grows ^[6]. Runge's phenomenon reshaped the field: it is the reason node placement, Chebyshev spacing, and piecewise methods dominate modern practice. The last major practical development is recent — Berrut and Trefethen's 2004 SIAM Review paper made the case that the barycentric form of Lagrange interpolation is fast and stable enough to be the default, a century and a half after everyone had filed the Lagrange form under "elegant but impractical" ^[7].

Related tools

/tools/standard-atmosphere — atmosphere properties are published at discrete altitudes; reading between them is this entry
/tools/isentropic-flow — the calculator replaces the compressible-flow tables you'd otherwise interpolate by hand
/tools/loop-4-20ma — 4–20 mA scaling is degree-1 interpolation between the two calibration endpoints
/tools/psychrometric-dewpoint — computes directly what psychrometric charts and tables make you interpolate
/tools/wire-gauge-awg — a table-lookup calculator; AWG properties between listed gauges follow the same logic

Sources

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