📈Interpolation and Approximation

Splines and Curves

Piecewise interpolation, cubic splines and their tridiagonal system, shape preservation, Bézier curves

Paris, 1959

For centuries, the lines of a ship's hull were laid out full size on the floor of a loft. The draughtsman pinned a long, thin, flexible strip of wood — a spline — at the points the design required and held it down with heavy lead weights called ducks. The strip settled into the smoothest shape it could, the one that stores the least bending energy, and he traced along it. Car bodies were designed the same way, from plaster and clay models measured by hand. In 1959 the carmaker Citroën hired a young mathematician, Paul de Casteljau, to find a way to describe such curves to the numerically controlled machines that were starting to cut metal. His algorithm — repeated linear interpolation between control points — was kept as an internal company secret. A few years later at rival Renault, the engineer Pierre Bézier developed the same curves independently for the UNISURF design system, and published them. The curves bear Bézier's name; the algorithm that evaluates them bears de Casteljau's. Every font on your screen, every vector drawing and every CAD surface is built from them.

Instead of one high-degree polynomial through all the data, use a low-degree polynomial on each interval and join the pieces. The simplest is piecewise linear — continuous, but with a kink at every data point. Piecewise cubic Hermite interpolation also matches slopes at the nodes, giving a continuous first derivative. A cubic spline goes one step further: on each interval it is a cubic, and at every interior node the value, slope and curvature are continuous. That is exactly the shape of the draughtsman's wooden strip, which by beam theory has piecewise-cubic deflection between the pegs.

💡The cubic spline system

Let

M_i = s''(x_i)

be the unknown curvatures at the nodes, and

h_i = x_{i+1} - x_i

. Requiring continuous slopes at each interior node gives, for

i = 1, \dots, n-1

h_{i-1}M_{i-1} + 2(h_{i-1} + h_i)M_i + h_iM_{i+1} = 6\left(\dfrac{y_{i+1} - y_i}{h_i} - \dfrac{y_i - y_{i-1}}{h_{i-1}}\right).

A tridiagonal, diagonally dominant system — solved in

O(n)

by the Thomas algorithm of lesson 3-3. Two extra conditions close it: natural (

M_0 = M_n = 0

, a free strip), clamped (given end slopes), or not-a-knot (the default in SciPy and MATLAB). Among all twice-differentiable interpolants, the natural spline minimises

\int (s'')^2\,dx

— the bending energy of the strip.

Derivation

Deriving the spline equations

On one interval. On

[x_i, x_{i+1}]

the spline is a cubic, so its second derivative is a straight line, running from

M_i

M_{i+1}

s''(x) = M_i\,\frac{x_{i+1} - x}{h_i} + M_{i+1}\,\frac{x - x_i}{h_i}.

Integrate twice. The two integration constants are easiest to write as multiples of

(x_{i+1} - x)

and

(x - x_i)

s(x) = M_i\frac{(x_{i+1} - x)^3}{6h_i} + M_{i+1}\frac{(x - x_i)^3}{6h_i} + C_i\,(x_{i+1} - x) + D_i\,(x - x_i).

Requiring

s(x_i) = y_i

and

s(x_{i+1}) = y_{i+1}

fixes them:

C_i = y_i/h_i - M_ih_i/6

and

D_i = y_{i+1}/h_i - M_{i+1}h_i/6

. Values and curvatures are now continuous by construction.

Matching slopes. Differentiate and evaluate at the left end of interval

i

, then at the right end of interval

i - 1

s'(x_i^+) = \frac{y_{i+1} - y_i}{h_i} - \frac{h_i}{3}M_i - \frac{h_i}{6}M_{i+1}, \qquad s'(x_i^-) = \frac{y_i - y_{i-1}}{h_{i-1}} + \frac{h_{i-1}}{3}M_i + \frac{h_{i-1}}{6}M_{i-1}.

Set them equal and multiply by 6:

h_{i-1}M_{i-1} + 2(h_{i-1} + h_i)M_i + h_iM_{i+1} = 6\Big(\frac{y_{i+1} - y_i}{h_i} - \frac{y_i - y_{i-1}}{h_{i-1}}\Big).

One equation per interior node, each linking a curvature to its two neighbours: tridiagonal. The diagonal entry

2(h_{i-1} + h_i)

is larger than the sum of the off-diagonal ones,

h_{i-1} + h_i

— strictly diagonally dominant, so the Thomas algorithm needs no pivoting (lesson 3-3).

📐A cam profile

Cam lift measured every 30° from 0 to 180°: 0, 1.2, 4.6, 8.0, 9.1, 6.3, 2.0 mm. The hand-built natural spline — tridiagonal system, Thomas algorithm — gives 2.70875 mm at 45°, 8.99875 mm at 105° and 4.21125 mm at 165°, identical to scipy.interpolate.CubicSpline(x, y, bc_type='natural') at every one of 361 test points. The curvature is continuous, which matters: a cam follower's acceleration is proportional to it, and a jump in curvature is a jump in force.

⚠Splines overshoot; shape-preserving cubics do not

Smoothness has a price: a spline through a step must overshoot it. Interpolating a solid-motor thrust curve (rise to 20 kN, flat burn, tail-off), the cubic spline peaks at 23.4 kN and dips below zero — thrust the motor never produced, and a 1.8% error in total impulse. PCHIP (Fritsch and Carlson, 1980) chooses the slopes at the nodes so that each piece stays monotone wherever the data are: maximum exactly 20.0 kN, minimum exactly 0, at the cost of a discontinuous second derivative. For physical quantities that must not go negative or exceed measured bounds — concentrations, thrust, probabilities — use a shape-preserving interpolant.

💡Bézier curves and de Casteljau's algorithm

A cubic Bézier curve is defined by four control points

P_0, \dots, P_3

B(t) = (1-t)^3P_0 + 3(1-t)^2tP_1 + 3(1-t)t^2P_2 + t^3P_3, \qquad 0 \le t \le 1.

It starts at

P_0

, ends at

P_3

, leaves along

P_0P_1

and arrives along

P_2P_3

; the inner points pull the curve without lying on it, and the curve stays inside their convex hull. De Casteljau's algorithm evaluates it by three rounds of linear interpolation — between consecutive control points, then between those points, then once more. It is perfectly stable (only convex combinations), and the intermediate points split the curve into two Bézier curves that together trace the original exactly — the basis of adaptive drawing and intersection algorithms.

Derivation

End tangents, and why repeated interpolation gives the same curve

Tangents. Differentiate

B(t)

term by term and collect:

B'(t) = 3\Big[(1-t)^2(P_1 - P_0) + 2(1-t)t\,(P_2 - P_1) + t^2(P_3 - P_2)\Big].

t = 0

only the first term survives,

B'(0) = 3(P_1 - P_0)

; at

t = 1

only the last,

B'(1) = 3(P_3 - P_2)

. The curve leaves along the first leg of the control polygon and arrives along the last — which is how a designer sets a tangent by dragging one point.

De Casteljau. For three control points, one round of interpolation gives

P_{01} = (1 - t)P_0 + tP_1

and

P_{12} = (1 - t)P_1 + tP_2

, and the second round gives

(1 - t)P_{01} + tP_{12} = (1 - t)^2P_0 + 2(1 - t)t\,P_1 + t^2P_2,

the quadratic Bézier curve. A cubic needs one more round, and expanding it produces the coefficients

(1 - t)^3,\ 3(1 - t)^2t,\ 3(1 - t)t^2,\ t^3

— the binomial pattern of

\big((1 - t) + t\big)^3

. Each round uses only weights

1 - t

and

t

between 0 and 1, so nothing is ever subtracted: no cancellation.

As long as algebra and geometry proceeded along separate paths, their advance was slow and their applications limited. But when these sciences joined company, they drew from each other fresh vitality and thenceforward marched on at a rapid pace toward perfection.

– Joseph-Louis Lagrange, Leçons élémentaires sur les mathématiques (1795)

💡Where you will use this

Mechanical design and CAD/CAM: B-splines and NURBS (the generalisation of Bézier curves) define every surface in a modern CAD model and every toolpath. Aircraft Engines and Aerodynamics: blade and airfoil shapes are parameterised by splines for optimisation. Rocket Propulsion: thrust and pressure curves interpolated with shape-preserving cubics. Optimal Guidance and Robotics: smooth trajectories through waypoints are splines, with continuous acceleration for the actuators.

📚Further reading

Carl de Boor, A Practical Guide to Splines (Springer, 1978; revised ed. 2001) — the standard reference, by the author of much of the theory.
Gerald Farin, Curves and Surfaces for CAGD: A Practical Guide (5th ed., Morgan Kaufmann, 2002) — Bézier curves, de Casteljau's algorithm and B-splines for design.
F. N. Fritsch and R. E. Carlson, 'Monotone piecewise cubic interpolation', SIAM Journal on Numerical Analysis 17(2), 238–246 (1980) — the PCHIP paper.

↓But there is a problem

Splines are local: moving one data point changes the curve only nearby. That is their strength for shapes. For signals that repeat — vibration of a gearbox, the alternating current in a grid, the Earth's daily rotation in a GNSS clock — the natural description is the opposite: a few global sine waves, each spread over the whole record. In the early 1960s the United States and the Soviet Union were negotiating a ban on nuclear tests, and the question was whether an underground explosion could be told apart from an earthquake on a distant seismometer. The answer depended on taking the frequency content of thousands of long records — and the computation was hopelessly slow.