Spline interpolation

The formula

On each span between knots, the spline is just a cubic:

Sᵢ(x) = aᵢ + bᵢ·(x − xᵢ) + cᵢ·(x − xᵢ)² + dᵢ·(x − xᵢ)³,   xᵢ ≤ x ≤ xᵢ₊₁

Reading: nothing exotic per piece — the intelligence is entirely in how the pieces are stitched. Requiring S, S′, and S″ to be continuous at every interior knot gives 4n − 2 conditions on 4n coefficients; the two left over are yours to spend at the boundaries, and that choice matters more than most people realize.

The stitching reduces to one equation per interior knot for the second derivatives Mᵢ = S″(xᵢ), with hᵢ = xᵢ₊₁ − xᵢ:

hᵢ₋₁·Mᵢ₋₁ + 2·(hᵢ₋₁ + hᵢ)·Mᵢ + hᵢ·Mᵢ₊₁ = 6·[ (yᵢ₊₁ − yᵢ)/hᵢ − (yᵢ − yᵢ₋₁)/hᵢ₋₁ ]

Reading: each knot's curvature is set by a balance with its two neighbors and the local change in slope of the data. Neighbors-only coupling means a tridiagonal system — solved in O(n) by the Thomas algorithm, which is why a spline through ten thousand points is instant.

The property that earned the natural spline its keep (Holladay's theorem):

∫ [S″(x)]² dx  ≤  ∫ [F″(x)]² dx    for every C² interpolant F through the same points

Reading: among all twice-differentiable curves through your data, the natural cubic spline carries the least total bending energy — it is the mathematical version of the draftsman's flexible wooden strip ^[7]^[8].

The error bound, for a clamped spline with correct end slopes on a function with a bounded fourth derivative:

max |f(x) − S(x)| ≤ (5/384)·h⁴·max |f⁗(x)|

Reading: fourth-order convergence — halve the knot spacing and the worst-case error drops by a factor of 16 ^[9]^[10]. A clamped spline of sin(x) on [0, π] with 21 knots lands within 1.6·10⁻⁶ everywhere; the bound says 7.9·10⁻⁶.

The Hermite form, when you know slopes m₀, m₁ as well as values on a span of width h, with t = (x − x₀)/h:

S(x) = (2t³ − 3t² + 1)·y₀ + (t³ − 2t² + t)·h·m₀ + (−2t³ + 3t²)·y₁ + (t³ − t²)·h·m₁

Reading: the unique cubic honoring all four numbers. Every cubic spline is secretly this formula applied span by span — the schemes differ only in where the slopes come from.

Where you meet it

The cam profile on its way to the grinder. A valve or indexing cam is specified as a displacement table, but the follower feels derivatives: contact force rides on acceleration, which is the second derivative of displacement. A curve that is merely continuous in slope puts a step in acceleration at every knot — an impact, cycle after cycle. The cubic spline's C² continuity is not a nicety here; it is the difference between a cam that runs quiet and one that hammers its follower into early retirement. Third-derivative (jerk) jumps at the knots are the residual compromise, and high-speed cam designers reach for higher-order splines precisely to chase them out.

Airfoil coordinates into the CAD model. Published airfoil sections arrive as 60-some coordinate pairs, and the CFD mesh or the CNC toolpath needs a smooth curve through them. Two traps: the leading edge has a vertical tangent, so you cannot fit y(x) — you fit x(s) and y(s) parametrically against arc length. And curvature continuity is exactly what the flow cares about; a slope-continuous-only fit puts a curvature kink at each data point that shows up as a spurious pressure blip in the solution.

The calibration curve from the cal lab. An RTD or load cell comes back with a dozen certified points and the DAQ needs readings in between. Spline it and the interior error is superb — but read the "How it works" section before you accept the default boundary condition, and before you spline anything with a plateau in it.

Resampling test data. Two instruments logged at different rates and you need both channels on one time base before you can cross-plot them. When the DAQ package aligns nonuniform timestamps onto a common base, it is interpolating under the hood — often with a spline, if you ask for one — and its behavior at the record's ends is the boundary-condition question again. (Uniform-ratio rate converters are a different animal: those run anti-aliasing filters, and their edge artifacts are filter startup transients, not spline boundaries.)

How it works

Splines exist because high-degree polynomial interpolation fails. Runge's function 1/(1 + 25x²) sampled at 21 evenly spaced points on [−1, 1] gives a degree-20 interpolating polynomial with a worst-case error near 60 — the fit oscillates wildly between the outer nodes. A cubic spline through the same 21 points has a worst-case error of about 0.003. Four orders of magnitude, same data. The mechanism: a polynomial's stiffness is global, so forcing it through a point on the left makes it swing on the right; a spline's pieces are local, coupled only through the tridiagonal chain, and a disturbance at one knot decays fast as it propagates down the line.

The boundary condition is where practitioners get hurt. "Natural" (S″ = 0 at both ends) sounds like the safe default and is anything but: unless your underlying function genuinely has zero curvature at the ends, you've imposed a wrong value there, and the penalty is real. Splining eˣ on [0, 1] with knots every 0.1: the natural spline misses by 1.3·10⁻³ near the ends, while the clamped spline — fed the true end slopes — misses by 7·10⁻⁷. Three-plus orders of magnitude, and the natural spline's convergence degrades from h⁴ to h² near the boundaries. If you know the end slopes, clamp. If you don't, use not-a-knot (the default in MATLAB and a bc_type option in SciPy), which recovers full fourth-order accuracy without pretending to know anything. Reserve "natural" for when zero end curvature is physically true — which, for the drafting strip that named the method, it was.

Second gotcha: splines overshoot. Run a cubic spline through data that steps from a run of 0s to a run of 1s and it undershoots to −0.11 and overshoots to 1.11 — an 11% excursion the data never made. On a calibration curve this manufactures non-monotonic behavior: two different inputs reading back as the same output. The fix is monotone piecewise-cubic interpolation — Fritsch and Carlson's 1980 scheme, shipped in every major package as PCHIP — which sacrifices C² continuity to guarantee no overshoot between points ^[11]^[12]. The engineering rule: smooth physics, use the spline; plateaus, saturation, or steps, use PCHIP.

Third: interpolation is a contract to pass through every point, which means reproducing every point's noise exactly. A spline through scattered data is a curve that treats scatter as signal, and its derivative is worse — differentiation amplifies exactly what interpolation preserved. Noisy data wants a smoothing spline or a least-squares fit with fewer coefficients than points, not an interpolant.

And the standing warning that applies to all interpolation: outside the last knot, a spline is just the final cubic continuing on its own opinion. Extrapolate one span past your calibration range and you are quoting a polynomial that no data ever constrained.

History

The word came off the shop floor. A spline was the thin flexible strip — wood or metal — that ship loftsmen bent through fixed points to draw hull curves at full scale, pinned in place by hooked lead weights the trade called ducks ^[4]^[5]^[6]. Elasticity theory says a thin strip bent through fixed points settles into the shape of least bending energy — which, for small deflections, is mathematically a piecewise cubic with continuous second derivative. The tool was solving the variational problem for a century before anyone wrote it down.

The writing down happened at the Army's Ballistic Research Laboratory in Aberdeen, Maryland. Isaac Jacob Schoenberg, released from the University of Pennsylvania for war work from 1943 to 1945, developed the theory while working on the approximation of equally spaced ballistics data, and published it in a two-part 1946 paper in the Quarterly of Applied Mathematics — the paper that gave the mathematical object the loftsman's name ^[1]^[2]^[3]. In 1957 John Holladay proved the property that connects tool to theorem: the natural cubic spline minimizes the integrated squared second derivative among all smooth interpolants — it really is the mathematical drafting strip ^[7]^[8].

Industry got there in parallel, three times. In the late 1950s and early 1960s, Paul de Casteljau at Citroën, Pierre Bézier at Renault, and Garrett Birkhoff with Carl de Boor at General Motors Research each built piecewise-polynomial machinery for describing car body surfaces numerically — the work that grew into B-splines, NURBS, and the CAD kernels every mechanical engineer now draws in ^[4]^[5]. The sharp error constants — the 5/384 — were pinned down by Hall and Meyer in 1976 ^[9]^[10], and Fritsch and Carlson closed the overshoot loophole with monotone piecewise cubics in 1980 ^[11]^[12]. From a bent stick in a ship loft to the default interpolator in every analysis package, in about one working lifetime.

Related tools

/tools/standard-atmosphere — atmosphere tables are the classic spline customer: smooth physics published at discrete altitudes
/tools/isentropic-flow — computes what you'd otherwise spline out of compressible-flow tables
/tools/psychrometric-dewpoint — psychrometric charts are hand-drawn spline territory; the calculator goes straight to the equations
/tools/loop-4-20ma — the two-point linear scaling that a multi-point spline calibration generalizes
/tools/wire-gauge-awg — table lookup between listed gauges is the same read-between-the-points problem

Sources

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