Spline Interpolation

Spline interpolation is useful for creating “continuous” curves out of discrete data points. These curves are can have C1 (first derivative) and C2 (second derivative) continuity, allowing for a smooth join between segments of the curve. An algorithm for computing natural cubic splines was listed on Wikipedia, and have reproduced it here so as to have the full math -> code translation in one place.

Definition of a Cubic Spline

Given a list of data values
Y = [y_0, y_1, \dots, y_n]

where y_n is an n-tuple defined as
y_n = (a_0, a_1, \dots, a_{n-1})

with corresponding location values X defined as
X = [x_0, x_1, \dots, x_n]

we wish to compute a value of Y for any value in the range [x_0, x_n]. We will do this for some x_i \leq x < x_j by evaluating a cubic spline S_i(x).

A cubic spline is of the form
{S}_{i}(x) = a_i + b_i (x-x_i) + c_i (x-x_i)^2 + d_i (x-x_i)^3.

where (x_i, a, b, c, d) is a 5-tuple describing the parameters of S_i(x).

Given our set of data Y and locations X, we wish to find n polynomials S_i (x) for i = 0, \dots , n-1

such that
S_i (x_i) = y_i = S_{i-1} (x_i),\quad i = 1 , \dots , n-1
S'_i (x_i) = S'_{i-1} (x_i),\quad i = 1 , \dots , n-1
S''_i (x_i) = S''_{i-1} (x_i),\quad i = 1 , \dots , n-1
S''_0 (x_0) = S''_{n-1} (x_n) = 0

We can compute the parameters of these n polynomials using the algorithm below.

Computation of Cubic Spline Parameters

Input: list of points Y \,, with \left | Y \right | =n+1, with locations X
Output: list spline parameters composed of n 5-tuples.

  1. Create new arrays a, b, d of size n
  2. Create new arrays c, h, l, u, z of size n+1
  3. Set l_0 = 1, u_0 = z_0 = 0
  4. Set h_0 = x_1 - x_0
  5. For i = 1, \dots, n-1
    1. Set h_i = x_{i+1} - x_i
    2. Set l_i = 2 \left ( x_{i+1} - x_{i-1} \right ) - h_{i-1}u_{i-1}
    3. Set u_i = \frac{h_i}{l_i}
    4. Set a_i = \frac{3}{h_i} \left ( y_{i+1} - y_i \right ) - \frac{3}{h_{i-1}} \left ( y_i - y_{i-1} \right )
    5. Set z_i = \frac{a_i - h_{i-1}z_{i-1}}{l_i}
  6. Set l_n = 1, z_n = c_n = 0
  7. For j = n-1, n-2, \dots, 0
    1. Set c_j = z_j - u_j c_{j+1}
    2. Set b_j = \frac{y_{j+1} - y_j}{h_j} - \frac{h_j \left( c_{j+1} + 2c_j \right)}{3}
    3. Set d_j = \frac{c_{j+1} - c_j}{3h_j}
  8. For i = 0, \dots, n-1
    1. Set S_{i,\left(x_i,a,b,c,d\right)} = \left( x_i, y_i, b_i, c_i, d_i \right)

For a C++ implementation of this algorithm, see Cubic Spline Interpolation