Why cubic splines?

Interpolation is the problem of estimating Y values for X positions that lie between known data points. For a set of points like (0, 1), (1, 3), (2, 2), a single high-degree polynomial can fit all points exactly, but such polynomials tend to overfit and suffer from numerical instability. An alternative is to fit piecewise polynomials between adjacent points, and the most common variant is linear interpolation: fast and reasonable for large datasets, but crude when only a few points are available.

Cubic splines offer a better tradeoff. By using third-degree polynomials between each pair of adjacent points, we can enforce continuity of both the first and second derivatives at the interior points where polynomials meet. This yields a curve that looks natural and smooth, without the oscillation problems of high-degree polynomials. For N points, we need N-1 cubic polynomials, each of the form:

Pi(x) = aix3 + bix2 + cix + di

Each polynomial contributes four unknown coefficients, giving 4N-4 unknowns total. Solving for these coefficients requires the same number of independent constraints.

Setting up the constraint equations

Using the three-point example above, we have two cubic polynomials and eight coefficients to determine. The constraints come from four sources, applied generically regardless of N:

  • Interpolation at endpoints: each polynomial must pass through the two original points it connects, yielding 2N-2 equations.
  • First-derivative continuity at interior points: where polynomials meet, their slopes must match, giving N-2 equations. The first derivative of a cubic is P′(x) = 3ax² + 2bx + c.
  • Second-derivative continuity at interior points: the curvature must also match, yielding another N-2 equations. The second derivative is P″(x) = 6ax + 2b.
  • Boundary conditions: for a natural spline, the second derivative is set to zero at the first and last original points, providing the final 2 equations.

For the concrete example, these eight equations can be written as a linear system. Solving it yields the coefficient vector (-0.75, 0, 2.75, 1, 0.75, -4.5, 7.25, -0.5), which defines the two interpolating polynomials. Once these are known, any X between consecutive original points can be evaluated by the appropriate polynomial.

Solving with Gaussian elimination

In general, spline interpolation with more than a few points quickly produces a large system. The matrix formulation Ax = b lets us apply standard linear algebra techniques. The implementation used in the accompanying code performs Gaussian elimination to row-echelon form, using partial pivoting for numerical stability, then applies Gauss-Jordan reduction to extract the solution.

The resulting algorithm is straightforward and dependency-free, implemented in readable, commented JavaScript. A complete code sample is available on GitHub, including functions to construct the equation system from a set of points, solve it, and render the results via an SVG-based demo page. A related line-plotting demo shows the same spline interpolation applied to arbitrary mathematical functions when sampling is sparse.