Cubic Spline Interpolation#

Just like linear interpolation, cubic spline interpolation joins the data in a piecewise fashion, but requires four data points for each piecewise cubic polynomial. For n points there are n-1 cubic functions, each requiring their own set of 4 coefficients.

\[S_i(x) = a_i \cdot x^3 + b_i \cdot x^2 + c_i \cdot x + d_i\]

Say we have 3 data points, this will require two piecewise (cubic spline) functions. Each pair of data points has their own function with four unknown coefficients, so we require 8 equations, (the middle point is used for both functions). So far the points only supply two equations for each function, giving 4 equations. The first and second derivative is the same for both functions, which gives two more equations. The last two equations are supplied by adding type. As more datapoints and their functions are added so the same argument shows that there is sufficient information to solve the equations.

Note

The cubic spline types are somewhat arbitrary and the selection is a value judgement based on appearance and data type

  • not-a-knot

    The first and second curve segments are the same polynomial - scipy default

  • periodic

    First and last y values are identical, f´(x0) = f´(xn); f´´(x0) = f´´(xn), with the first and second derivatives being equal at the first and last point.

  • clamped

    First derivative at curve ends is zero, f´(x0) = f´(xn) = 0

  • natural

    Second derivative at curve ends is zero, f´´(x0) = f´´(xn) = 0, this gives a straight line at both ends.

Each type affects the shape of the splines.

../_images/cubic_spline_interp.png

Cubic spline interpolation of 3 data points on a sine curve.#

All four types pass through all three data points but their paths differ slightly, knowledge of how the data is generated helps in finding the best fit.

The simplest way to implement is to use the scipy module Interpolate and import cubic_interpolate.

Note

Periodic Type

Setting this up may be a bit tricky if the original function is unknown. The first and second derivatives need to be equal, so probably best to use another type if unsure.

../_images/cubic_spline_interpolation.png

Using the cubic spline interpolation on a sine wave.#

Show/Hide Code cub_spline_interp.py

# https://gist.github.com/dankkom/01922dbf4f61d62cf78fa0f86560d14d

import numpy as np
from math import sqrt

def cubic_interp1d(x0, x, y):
    """
    Interpolate a 1-D function using cubic splines.
      x0 : a float or an 1d-array
      x : (N,) array_like
          A 1-D array of real/complex values.
      y : (N,) array_like
          A 1-D array of real values. The length of y along the
          interpolation axis must be equal to the length of x.

    Implement a trick to generate at first step the cholesky matrice L of
    the tridiagonal matrice A (thus L is a bidiagonal matrice that
    can be solved in two distinct loops).

    additional ref: www.math.uh.edu/~jingqiu/math4364/spline.pdf
    """
    x = np.asfarray(x)
    y = np.asfarray(y)

    # remove non finite values
    # indexes = np.isfinite(x)
    # x = x[indexes]
    # y = y[indexes]

    # check if sorted
    if np.any(np.diff(x) < 0):
        indexes = np.argsort(x)
        x = x[indexes]
        y = y[indexes]

    size = len(x)

    xdiff = np.diff(x)
    ydiff = np.diff(y)

    # allocate buffer matrices
    Li = np.empty(size)
    Li_1 = np.empty(size-1)
    z = np.empty(size)

    # fill diagonals Li and Li-1 and solve [L][y] = [B]
    Li[0] = sqrt(2*xdiff[0])
    Li_1[0] = 0.0
    B0 = 0.0 # natural boundary
    z[0] = B0 / Li[0]

    for i in range(1, size-1, 1):
        Li_1[i] = xdiff[i-1] / Li[i-1]
        Li[i] = sqrt(2*(xdiff[i-1]+xdiff[i]) - Li_1[i-1] * Li_1[i-1])
        Bi = 6*(ydiff[i]/xdiff[i] - ydiff[i-1]/xdiff[i-1])
        z[i] = (Bi - Li_1[i-1]*z[i-1])/Li[i]

    i = size - 1
    Li_1[i-1] = xdiff[-1] / Li[i-1]
    Li[i] = sqrt(2*xdiff[-1] - Li_1[i-1] * Li_1[i-1])
    Bi = 0.0 # natural boundary
    z[i] = (Bi - Li_1[i-1]*z[i-1])/Li[i]

    # solve [L.T][x] = [y]
    i = size-1
    z[i] = z[i] / Li[i]
    for i in range(size-2, -1, -1):
        z[i] = (z[i] - Li_1[i-1]*z[i+1])/Li[i]

    # find index
    #print(x0)
    index = x.searchsorted(x0)
    np.clip(index, 1, size-1, index)

    xi1, xi0 = x[index], x[index-1]
    yi1, yi0 = y[index], y[index-1]
    zi1, zi0 = z[index], z[index-1]
    hi1 = xi1 - xi0

    # calculate cubic
    f0 = zi0/(6*hi1)*(xi1-x0)**3 + \
         zi1/(6*hi1)*(x0-xi0)**3 + \
         (yi1/hi1 - zi1*hi1/6)*(x0-xi0) + \
         (yi0/hi1 - zi0*hi1/6)*(xi1-x0)
    return f0

if __name__ == '__main__':
    import matplotlib.pyplot as plt
    import seaborn as sns

    sns.set_theme()
    x = np.linspace(0, 10, 11)
    y = np.sin(x)
    plt.scatter(x, y)

    x_new = np.linspace(0, 10, 201)

    plt.plot(x_new, np.sin(x_new), 'y', label='sine')
    plt.plot(x_new, cubic_interp1d(x_new, x, y), 'k', label='cubic spline', linestyle=':')
    plt.xlabel(r'$x$')
    plt.ylabel(r'$f(x)$')
    plt.title('Cubic Spline Interpolation')
    plt.legend()
    plt.show()
    #plt.savefig('../../figures/cubic_spline_interpolation.png')

The cubic spline takes in its stride the data points that caused trouble with the Lagrange and Newton polynomials.

../_images/cubic_spline_knot_points2.png

Using the cubic spline interpolation on a the data set knot_points2.#

The result should be compared with the Lagrange polynomial.

This Page