Splines
The previous page ended on a warning. Thread a single high-degree polynomial through many
equally-spaced points and it can erupt into wild oscillations near the ends —
Runge's phenomenon.
Adding more points makes it worse. So the great idea of this page is to stop asking one
stiff, global polynomial to do everything, and instead lay down many short, low-degree
pieces — one per gap between data points — and glue them together so smoothly that the joins
are invisible. That patchwork of polynomials is a spline.
The name is not a metaphor. Before computers, a draughtsman drawing the curved hull of a ship pinned a
thin, flexible strip of wood or steel — a spline — through the marked points and let it
bend into a fair, smooth curve. The strip settled into the shape that a thin elastic beam
naturally takes, and that shape turns out to be, almost exactly, a chain of cubics. This page teaches you
the star of the family, the natural cubic spline: what it is, the smoothness conditions
that pin it down, the neat tridiagonal system you solve to build it, and the beautiful fact that it is
literally the smoothest curve — the one of least bending — that passes through your points.
Warm-up: the linear spline connects the dots
The crudest spline just joins each point to the next with a straight line. On the interval
[x_i, x_{i+1}] it is the degree-1 interpolant
S_i(x) = y_i + \frac{y_{i+1} - y_i}{x_{i+1} - x_i}\,(x - x_i),
a simple weighted average of the two endpoint heights. It is honest and it never overshoots — the value
on each interval stays between its two data points, so a linear spline can never produce the
crazy humps of a high-degree polynomial. Runge's disaster simply cannot happen. Let's evaluate one:
// Linear spline through (x_i, y_i): find the bracketing interval, then blend.
const xs: number[] = [0, 1, 2, 3, 4];
const ys: number[] = [0, 2, 1, 3, 1];
function linearSpline(X: number): number {
// locate the interval [xs[i], xs[i+1]] containing X
let i = 0;
while (i < xs.length - 2 && X > xs[i + 1]) i++;
const t = (X - xs[i]) / (xs[i + 1] - xs[i]); // fraction across the gap, 0..1
return (1 - t) * ys[i] + t * ys[i + 1]; // weighted average of endpoints
}
for (const X of [0, 0.5, 1.5, 2.25, 3.75, 4]) {
console.log(`S(${X}) = ${linearSpline(X).toFixed(3)}`);
}
It reproduces the data at the nodes and interpolates linearly between them. But run your eye along the
curve and you feel the problem instantly: at every interior point the slope jumps. The
curve is continuous but kinked — a corner at every data point. A car following this road would
jerk the wheel at each node. The value matches across the join, but the first derivative does not. To
get a curve that looks and behaves smoothly, we need the pieces to agree not just in height but in
slope, and even in curvature. That is the leap to cubics.
A quadratic per interval can match the value and the first derivative at each join, which already
kills the kinks. But a piece-wise quadratic has an awkward asymmetry: fixing the slope at the left
end of each interval forces the slope at the right end, so information marches one-way across the data
and the curve tends to develop a lopsided wobble. Cubics are the sweet spot — the lowest degree that
can match value, slope and curvature at both ends of every interval at once, giving a curve
with no built-in direction and no visible seams. That is why "spline" almost always means
cubic spline.
The cubic spline: smooth to the second derivative
A cubic spline puts a separate cubic polynomial on each subinterval
[x_i, x_{i+1}],
S_i(x) = a_i + b_i(x - x_i) + c_i(x - x_i)^2 + d_i(x - x_i)^3,
and then demands that neighbouring pieces meet as gracefully as possible. Two adjacent cubics don't just
touch at the shared knot — we force their heights, their slopes, and their
curvatures to agree there. Matching the second derivative is the secret ingredient: it
is what removes the last visible trace of a seam, because the eye reads a discontinuity in curvature as a
subtle crease. A curve whose value, first and second derivatives are all continuous is called
C^2, and a cubic spline is the standard way to get a C^2
interpolant.
-
Piecewise cubic. On each of the n subintervals of
n+1 knots x_0 < x_1 < \cdots < x_n,
S is a cubic — that is 4n unknown coefficients
in all.
-
Interpolation. Each piece matches the data at both its ends,
S_i(x_i) = y_i and S_i(x_{i+1}) = y_{i+1}.
-
C^2 smoothness. At every one of the
n-1 interior knots the value, the first derivative
S', and the second derivative S'' all match
across the join.
-
Two spare degrees of freedom. Those conditions supply
4n - 2 equations for 4n unknowns, leaving
exactly two free — fixed by a boundary condition at each end:
- Natural: S''(x_0) = S''(x_n) = 0 (the ends
straighten out).
- Clamped: prescribe the end slopes S'(x_0) and
S'(x_n).
- Not-a-knot: force the third derivative to match at the first and last
interior knots, so the two end pieces are a single cubic.
Let's count the bookkeeping once, slowly, because it is the heart of why the spline is uniquely
determined. Interpolation at both ends of every piece is 2n equations.
Continuity of S' and of S'' at each of the
n-1 interior knots adds 2(n-1) more. That is
2n + 2(n-1) = 4n - 2 equations against 4n unknowns
— a deficit of two, closed by the two boundary conditions. With them, the system has exactly one
solution.
Building it: a tridiagonal system
Solving 4n equations sounds heavy, but a change of variables shrinks it
dramatically. Take the unknowns to be the second derivatives at the knots,
M_i = S''(x_i). Because each piece is a cubic, its second derivative is
linear, so on [x_i, x_{i+1}] it slides straight from
M_i to M_{i+1}. Integrate that twice, pin the two
constants using S(x_i)=y_i and S(x_{i+1})=y_{i+1},
and impose that the slopes agree at each interior knot. With
h_i = x_{i+1} - x_i everything collapses to one tidy relation at each
interior knot i = 1, \ldots, n-1:
h_{i-1}M_{i-1} + 2(h_{i-1}+h_i)M_i + h_i M_{i+1} = 6\!\left(\frac{y_{i+1}-y_i}{h_i} - \frac{y_i-y_{i-1}}{h_{i-1}}\right).
Look at what each equation touches: only M_{i-1},
M_i and M_{i+1} — three neighbours. Stack them up
and the coefficient matrix is tridiagonal: non-zero only on the diagonal and the two
bands hugging it. For the natural spline we set the ends
M_0 = M_n = 0, leaving n-1 equations in
n-1 unknowns. A tridiagonal system is a joy to solve: the
Thomas algorithm — Gaussian elimination that only ever touches the three bands — cracks
it in O(n) time, a sweep down and a sweep back up, no more work than reading
the data once. Here is the whole pipeline: assemble the system, solve it, and evaluate the spline at a
query point.
// Natural cubic spline: build the tridiagonal system for M_i = S''(x_i),
// solve it with the Thomas algorithm, then evaluate at a query point.
const xs: number[] = [0, 1, 2, 3, 4];
const ys: number[] = [0, 2, 1, 3, 1];
const n: number = xs.length - 1; // number of intervals
const h: number[] = [];
for (let i = 0; i < n; i++) h.push(xs[i + 1] - xs[i]);
// Tridiagonal system for the interior M_1..M_{n-1}; natural ends M_0 = M_n = 0.
// Bands: lo (sub-diagonal), di (diagonal), up (super-diagonal); rhs = right-hand side.
const lo: number[] = [], di: number[] = [], up: number[] = [], rhs: number[] = [];
for (let i = 1; i < n; i++) {
lo.push(h[i - 1]);
di.push(2 * (h[i - 1] + h[i]));
up.push(h[i]);
rhs.push(6 * ((ys[i + 1] - ys[i]) / h[i] - (ys[i] - ys[i - 1]) / h[i - 1]));
}
// Thomas algorithm: forward-eliminate the sub-diagonal, then back-substitute.
const m: number = di.length;
for (let k = 1; k < m; k++) {
const w = lo[k] / di[k - 1];
di[k] -= w * up[k - 1];
rhs[k] -= w * rhs[k - 1];
}
const Min: number[] = new Array(m);
Min[m - 1] = rhs[m - 1] / di[m - 1];
for (let k = m - 2; k >= 0; k--) Min[k] = (rhs[k] - up[k] * Min[k + 1]) / di[k];
// Reassemble the full second-derivative vector with the natural zero ends.
const M: number[] = [0, ...Min, 0];
console.log(`second derivatives M_i = [${M.map((v) => v.toFixed(3)).join(", ")}]`);
// Evaluate the spline on the piece containing X.
function spline(X: number): number {
let i = 0;
while (i < n - 1 && X > xs[i + 1]) i++;
const A = xs[i + 1] - X, B = X - xs[i], hi = h[i];
return (
(M[i] * A * A * A + M[i + 1] * B * B * B) / (6 * hi) +
(ys[i] / hi - (M[i] * hi) / 6) * A +
(ys[i + 1] / hi - (M[i + 1] * hi) / 6) * B
);
}
for (const X of [0, 1, 1.5, 2.5, 3.5, 4]) {
console.log(`S(${X}) = ${spline(X).toFixed(4)}`);
}
The spline nails the data at the knots (check S(0), S(1), S(4)) and threads
smoothly between them. Notice M_0 and M_n print as
exactly 0 — that is the natural boundary condition, visible in the output.
Proving it's smooth: the join is seamless
A claim as strong as "value, slope and curvature all match at every knot" deserves a direct
test. Let's take the spline we just built, walk up to an interior knot from the piece on its
left and from the piece on its right, and compare the height, the first derivative,
and the second derivative from each side. If the cubic spline is genuinely C^2,
all three must agree to full machine precision:
// Verify C^2: at an interior knot, value / S' / S'' agree from both pieces.
const xs: number[] = [0, 1, 2, 3, 4];
const ys: number[] = [0, 2, 1, 3, 1];
const M: number[] = [0, -6.6, 6.9, -6.6, 0]; // natural M_i from the previous solve (rounded)
const h = (i: number): number => xs[i + 1] - xs[i];
// Piece i as a function of x, plus its first and second derivatives.
function piece(i: number, X: number): [number, number, number] {
const A = xs[i + 1] - X, B = X - xs[i], hi = h(i);
const S =
(M[i] * A * A * A + M[i + 1] * B * B * B) / (6 * hi) +
(ys[i] / hi - (M[i] * hi) / 6) * A +
(ys[i + 1] / hi - (M[i + 1] * hi) / 6) * B;
const S1 =
(-M[i] * A * A + M[i + 1] * B * B) / (2 * hi) -
(ys[i] / hi - (M[i] * hi) / 6) +
(ys[i + 1] / hi - (M[i + 1] * hi) / 6);
const S2 = (M[i] * A + M[i + 1] * B) / hi; // linear in x, as promised
return [S, S1, S2];
}
const knot = 2; // interior knot x = 2, shared by pieces 1 and 2
const left = piece(1, xs[knot]); // approach from the left piece
const right = piece(2, xs[knot]); // approach from the right piece
const labels = ["value S ", "slope S'", "curv S''"];
for (let k = 0; k < 3; k++) {
console.log(`${labels[k]}: left = ${left[k].toFixed(4)} right = ${right[k].toFixed(4)}`);
}
All three lines match from both sides — the join really is seamless in value, slope and curvature. (The
tiny differences you may see are only because M_i was rounded when pasted in;
with the exact solve they vanish entirely.) That C^2 continuity is precisely
what your eye reads as "smooth". A linear spline would pass the value test and fail the slope test at
every knot.
Linear kinks versus cubic smoothness, side by side
Here are the same five data points joined two ways: by the connect-the-dots linear spline
and by the natural cubic spline. The linear version is a run of straight segments with a
corner at every knot; the cubic version glides through the very same points with no visible seam, curving
to keep its slope and curvature continuous. Both hit every data point exactly — they differ only in
how they travel between.
Look near the peak at x = 3: the cubic swings a touch above the straight
segments as it eases in and out. That gentle bulge is the price of smoothness — and, occasionally, a
trap.
Smooth does not mean well-behaved. To keep its curvature continuous, a cubic spline is free to bulge
past your data between the knots. If your five points are all rising, the spline can still dip
below the lowest of two neighbours or rear above the highest — an overshoot. For most
drawing that is invisible and harmless, but if the quantity is physically bounded — a concentration
that must stay \ge 0, a probability that must stay
\le 1 — the spline can hand you an impossible value between data points.
Two things to remember. First, a spline is only as good as its boundary conditions:
the natural choice S''=0 is convenient but assumes the curve straightens at
the ends, which may be flatly wrong for your data — a clamped spline with the true end slopes is often
far better. Second, when monotonicity matters, reach for a monotone spline (such as a
PCHIP / Fritsch–Carlson Hermite spline) that trades a little of the C^2
smoothness for a guarantee never to overshoot.
Why cubics? The draughtsman's beam
The cubic spline is not an arbitrary choice of degree — it is the answer to a physics problem. A thin
elastic beam bent to pass through fixed points settles into the shape that minimises its stored
bending energy, which to good approximation is
E[S] = \int_{x_0}^{x_n} \big(S''(x)\big)^2 \, dx,
the total squared curvature along the curve. Among all smooth functions that interpolate the
given points, the one that makes this integral smallest is exactly the natural cubic
spline. That is why the draughtsman's flexible strip traced a cubic spline without anyone
solving a single equation: nature minimised the energy for them, and the boundary term
S''=0 at the free ends is precisely the "natural" condition that a beam left
loose at its tips adopts. The cubic spline is, in a rigorous sense, the smoothest possible curve
through your data.
-
Among all C^2 functions g that interpolate the
data g(x_i) = y_i, the natural cubic spline S
uniquely minimises the bending energy
\displaystyle \int_{x_0}^{x_n} \big(g''\big)^2\,dx.
-
Equivalently, S is the interpolant of least total squared curvature — the
mathematical idealisation of the flexible drafting spline that gave the method its name.
Setting S'' = 0 at the two ends says the curve stops bending right at the
boundary — it leaves the interval as a straight line. That is exactly how a physical beam behaves when
its ends are free rather than pinned to a slope, and it is what makes the natural spline the true
energy minimiser. The catch is that real data rarely straightens out at the edges; if you happen to
know the slopes there, a clamped spline honours them and usually tracks the underlying
function more faithfully near the ends, where — just as with polynomial interpolation — the worst
errors always lurk.