Timo Denk's Blog

Cubic Spline Interpolation

· Timo Denk

Cubic spline interpolation is a mathematical method commonly used to construct new points within the boundaries of a set of known points. These new points are function values of an interpolation function (referred to as spline), which itself consists of multiple cubic piecewise polynomials. This article explains how the computation works mathematically.

After an introduction, it defines the properties of a cubic spline, then it lists different boundary conditions (including visualizations), and provides a sample calculation. Furthermore, it acts as a reference for the mathematical background of the cubic spline interpolation tool on tools.timodenk.com which is introduced at the end of the article.

Introduction

Cubic spline interpolation is the process of constructing a spline f:[x1,xn+1]R which consists of n polynomials of degree three, referred to as f1 to fn. A spline is a function defined by piecewise polynomials. Opposed to regression, the interpolation function traverses all n+1 pre-defined points of a data set D. The resulting function has the following structure:f(x)={a1x3+b1x2+c1x+d1if x[x1,x2]a2x3+b2x2+c2x+d2if x(x2,x3]anx3+bnx2+cnx+dnif x(xn,xn+1].Note that all polynomials are just valid within an interval; they compose the interpolation function. While extrapolation predicts a development outside the range of the data, interpolation works just within the data boundaries [x1,xn+1].

With properly chosen coefficients ai, bi, ci, and di for the polynomials, the resulting function traverses the points smoothly. For determining the coefficients, several equations are formulated which all together compose a uniquely solvable system of equations.

Plot of the data set D which will be used throughout this article except for the example section, defined as D={(1.5,1.2),(.2,0),(1,.5),(5,1),(10,1.2),(15,2),(20,1)}

Table of Notation

SymbolDescription
Dn+1 data points, a set of ordered pairs (x,y)
(xi,yi)ith element in D under the condition xi<xi+1 (leftmost point: (x1,y1))
f:[x1,xn+1]RInterpolation function (spline) which predicts new points and consists of n polynomials
fi:[xi,xi+1]Rith polynomial of third degree: aix3+bix2+cix+di

Properties of a Cubic Spline

In order to determine the 4n coefficients of all polynomials, several equations must be stated. First of all, every polynomial is known to pass through exactly two points. Therefore, the 2n equations f1(x1)=y1f1(x2)=y2f2(x2)=y2f2(x3)=y3fn(xn)=ynfn(xn+1)=yn+1,written-out asa1x13+b1x12+c1x1+d1=y1a1x23+b1x22+c1x2+d1=y2a2x23+b2x22+c2x2+d2=y2a2x33+b2x32+c2x3+d2=y3anxn3+bnxn2+cnxn+dn=ynanxn+13+bnxn+12+cnxn+1+dn=yn+1,are applicable. They express that at x=x1, the value of the first polynomial is equal to y1, and at x=x2, the value is y2. For the point where the second polynomial begins (x=x2), which is exactly where the first polynomial has ended, the second polynomial’s value is y2. At x=x3, it is y3, and so forth.

Furthermore, the first and second derivative of all polynomials are identical in the points where they _touch _their adjacent polynomial. The derivatives of polynomials of degree three are ddxfi(x)=3aix2+2bix+ci and d2dx2fi(x)=6aix+2bi. f1 touches the second polynomial f2 at x2, the second the third at x3, and so forth until fn1 touches fn at xn:ddxf1(x)=ddxf2(x)|x=x2ddxf2(x)=ddxf3(x)|x=x3ddxfn1(x)=ddxfn(x)|x=xn. Evaluated at the corresponding x values, that gives
3a1x22+2b1x2+c1=3a2x22+2b2x2+c23a2x32+2b2x3+c2=3a3x32+2b3x3+c33an1xn2+2bn1xn+cn1=3anxn2+2bnxn+cn. Equivalently, for the second derivative, the same is done by statingd2dx2f1(x)=d2dx2f2(x)|x=x2d2dx2f2(x)=d2dx2f3(x)|x=x3d2dx2fn1(x)=d2dx2fn(x)|x=xn, written-out as
6a1x2+2b1=6a2x2+2b26a2x3+2b2=6a3x3+2b36an1xn+2bn1=6anxn+2bn.
This adds another 2(n1) equations.

Boundary Conditions

In order to be able to solve the system of equations, two more pieces of information are required. Arbitrary constraints like setting the third derivative in the (say) fourth point to zero may be used. However, the selection of a boundary condition, consisting of a pair of equations, is the commonly used method. The four conditions “natural spline”, “not-a-knot spline”, “periodic spline”, and “quadratic spline”, are described in detail below.

Natural Spline

The natural spline is defined as setting the second derivative of the first and the last polynomial equal to zero in the interpolation function’s boundary points:
6a1x1+2b1=06anxn+1+2bn=0.The visual interpretation is that the function’s change of steepness approaches zero in the first and last point.

Boundary condition “natural”. Intuitively, the result looks best compared to the other methods.

Not-a-Knot Spline

For the not-a-knot boundary condition, the first two polynomials’ third derivatives are equal in the point where they touch (x2) each other. The same applies to the last two polynomials which touch at xn, giving d3dx3f1(x)=d3dx3f2(x)|x=x2d3dx3fn1(x)=d3dx3fn(x)|x=xn.
After calculating the derivatives this can be broken down to setting a1=a2an1=an;that is making the a-coefficients equal to each other. The resulting plot can be seen in the following figure:

Boundary condition “not-a-knot”. The consistent change of steepness results in a slightly shifted maximum next to the second rightmost point, compared to the natural spline.

Periodic Spline

For a periodic interpolation, the last polynomial’s first and second derivative are set to be equal to the first polynomial’s first and second derivative. Graphically, this means that the function is tileable. This can be particularly helpful if f(x1)=f(xn+1), but works just as well with a function where the y values of the first and last point differ. The figure below shows that f1 can be copied and attached to the last point and the last polynomial fn to the first point, without making the function look odd.

The corresponding equations are
3a1x12+2b1x1+c1=3anxn+12+2bnxn+1+cn6a1x1+2b1=6anxn+1+2bn.

Boundary condition “periodic”. The copied polynomials f1 and f6, colored in gray, fit the function’s shape smoothly after being moved to the the opposite side.

Quadratic Spline

The quadratic boundary condition defines the first and the last polynomial to be quadratic which makes this method the simplest one. The two parabolas are highlighted in red in the plot below. Mathematically spoken, the first coefficient of both polynomials is equal to zero: a1=an=0.

Boundary condition “quadratic”. The first and the last polynomial are not cubic but quadratic (parabola pieces; colored in red). For the given set of data, this does not result in a counter-intuitive shape of the interpolation function.

Sample Calculation

This section performs a sample calculation with a real set of data and the “natural” boundary condition. Its purpose is to demonstrate how the equations are determined and actually turned into an interpolation function.

Given the set of points D={(0,21),(1,24),(2,24),(3,18),(4,16)}, the interpolation function will consist of |D|1=4 polynomials, namely f1, f2, f3, and f4. Each is of degree three: fi=aix3+bix2+cix+di. The first polynomial f1 is known to pass through the first and the second point, therefore f1(0)=21 and f1(1)=24. Equivalently, the second polynomial traverses the second and third point and so forth, leading to the equations
f2(1)=24f2(2)=24f3(2)=24f3(3)=18f4(3)=18f4(4)=16.

Furthermore, the first and second derivatives of all adjacent polynomials are equal in the point where they touch. The f1 touches f2 at x=1, leading to ddxf1(x)=ddxf2(x)|x=1, which is, after taking the derivative, equal to 3a1+2b1+c1=3a2+2b2+c2. Again, the same can be done for the other polynomials, giving 12a2+8b2+c2=12a3+8b3+c327a3+18b3+c3=27a4+18b4+c4. The second derivative is handled in the same way so the following three equations apply:
6a1+2=6a2+212a2+2=12a3+218a3+2=18a4+2

The final two pieces of information are contained in the chosen boundary condition which in this example is the natural spline. The corresponding equations are
2b1=024a4+2b4=0.

All equations are forming a linear system of equations which can be written in matrix notation.

  • Row 1 to 8: Equations expressing that the polynomials traverse the points in D.
  • Row 9 to 11: The first derivative of two adjacent polynomials is equal in the point where they touch.
  • Row 12 to 14: The second derivative of two adjacent polynomials is equal in the point where they touch.
  • Row 15 and 16: The boundary condition (natural spline).

[00010000000000001111000000000000000011110000000000008421000000000000000084210000000000002793100000000000000002793100000000000064164132103210000000000000124101241000000000000027610276106200620000000000000012200122000000000000001820018200020000000000000000000000000024200][a1b1c1d1a2b2c2d2a3b3c3d3a4b4c4d4]=[212424242418181600000000]

This system of equations can be solved (e.g. with Gaussian elimination) so that the function coefficients are calculated and the interpolation function
f(x)={0.30357x3+3.3036x+21if x[0,1]1.4821x3+3.5357x20.23214x+22.179if x(1,2]3.2321x324.750x2+56.339x15.536if x(2,3]1.4464x3+17.357x269.982x+110.79if x(3,4].can be written down. f(x) is plotted in the figure below.

The interpolation function f(x) with boundary condition “natural”.
The polynomials (from top left to bottom right) f1, f2, f3, and f4 compose the interpolation function, within their intervals (colored in black). Outside of these bounds the polynomial function values do not match the points.

Online Tool

The cubic spline interpolation tool on tools.timodenk.com performs a cubic spline interpolation and visualizes the resulting function derived from a data set provided by the user. It offers selection of the boundary conditions and of the axes scaling. If needed, the corresponding LATEX code for rendering the interpolation function in a document may be exported. The latter was also used to create the figures of this article. Through tools.timodenk.com, there is also access to linearquadratic, and polynomial interpolation tools.

As a side note for developers, it is worth mentioning that a user of the tool reported a reproducible, weird interpolation result for a specific set of data. As it turned out, this was caused by the floating point precision of JavaScript which was too low for computing the coefficients from the system of equations without significant loss of precision. The issue was resolved by using the library Math.js. It features among other things a BigNumber class for very high precision numbers, which are now in use.