Cubic Spline Interpolation

This documentation was generated with the assistance of AI. Please report any inaccuracies.

Local, low-order polynomial interpolation gives a function whose first and higher derivatives are discontinuous wherever the interpolation switches which tabulated points it uses. A cubic spline avoids this by fitting a cubic polynomial between every pair of tabulated points whose coefficients are chosen — using information slightly beyond just that interval — so that the resulting function is continuous through its second derivative everywhere, including at the tabulated points themselves. It is described in Press et al., 2007, Section 3.3 (page 120), and implemented as CubicSplineInterpolator in the com.irurueta.numerical.interpolation package.

From linear interpolation to a spline

Linear interpolation within an interval is:

This has zero second derivative inside each interval and an undefined one at the tabulated points. Suppose, in addition to the , the second derivatives were also known at every tabulated point. Adding a cubic term to the formula above — built to vanish at both endpoints of the interval, and whose second derivative varies linearly from to — gives a continuous second derivative without disturbing the values at the endpoints:

Solving for the second derivatives

The are not actually known in advance — they are exactly what the spline needs to solve for. Requiring the first derivative (obtained by differentiating the formula above) to also be continuous across every interior tabulated point gives, for :

This is linear equations in the unknown , so two more conditions — typically boundary conditions at and — are needed for a unique solution. The two common choices, both supported by this library’s implementation, are:

  • Set and/or to zero, giving the natural cubic spline.

  • Set them to whatever value makes the first derivative of the interpolating function match a known, specified value at that boundary.

Crucially, each is coupled only to its immediate neighbors, so this system is tridiagonal and can be solved in operations, rather than the a generic linear solver would need.

Algorithm

flowchart TD A["Boundary conditions: natural (y''=0) or\nspecified first derivative at each end"] --> B["Tridiagonal decomposition:\nforward elimination over N-2 interior equations"] B --> C["Back-substitution: solve for all y_i''"] C --> D["Store y'' array once, at construction time"] D --> E["interpolate(x): locate interval [x_j, x_j+1],\nevaluate y = A*y_j + B*y_j+1 + C*y_j'' + D*y_j+1''"]

The tridiagonal solve happens once, at construction time, and its result — the array of second derivatives — is stored for the lifetime of the interpolator. Each subsequent call to interpolate(x) is then cheap: locate the bracketing interval (via the same locate/hunt search described in Polynomial Interpolation) and evaluate the closed-form cubic formula above.

Library implementation

CubicSplineInterpolator extends BaseInterpolator with a fixed order of M = 2 (it always brackets between two adjacent tabulated points, unlike the higher-order polynomial interpolators):

  • CubicSplineInterpolator(double[] x, double[] y) — natural spline (zero second derivative at both ends).

  • CubicSplineInterpolator(double[] x, double[] y, double yp1, double ypn) — clamped spline, with the first derivative at the lower and upper endpoints specified explicitly.

Unlike PolynomialInterpolator, this class has no error estimate — but it guarantees the interpolated curve passes exactly through every tabulated point, and is smoother across interval boundaries than piecewise low-order polynomial interpolation.

A natural spline is the right choice when nothing is known about the true first derivative at the endpoints. If the derivative is known there — analytically, or accurately enough by numerical differencing — supplying it via yp1/ypn produces a more accurate spline near the boundaries.

Example

// natural cubic spline
final var interpolator = new CubicSplineInterpolator(x, y);
final var value = interpolator.interpolate(target);