Piecewise Interpolation: Cubic Spline Interpolation
Cubic Spline Interpolation
There are different schemes of piecewise cubic spline interpolation functions which vary according to the end conditions. The scheme presented here is sometimes referred to as “Not-a-knot” end condition in which the first cubic spline is defined over the interval and the last cubic spline is defined on the interval
(Figure 9). Given
data points with
intervals,
cubic polynomials are defined on the intervals:
.
each cubic polynomial is defined as:
. The
unknowns can be obtained using the following conditions(Figure 9):
- Each polynomial has to pass through the data points which are the bounds of the interval, therefore this results in
equations. In addition the first polynomial and the last (
) polynomial have to pass through the second and next-to-last data points respectively. In total, this produces
equations.
- On the
intermediate nodes between the polynomials, the first and second derivatives are equal to each other. This results in
equations
The above system of equations is guaranteed a unique solution as long as the data points are ordered: .

Figure 9. Scheme of piecewise cubic interpolation with continuity in the first and second derivatives.

As shown in Figure 10, this cubic interpolation scheme produces a very smooth interpolating function for the data points of the Runge function given in the previous section.
View Mathematica CodeData = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4,0.20}, {-0.2, 0.50}, {0, 1}, {0.2, 0.5}, {0.4, 0.2}, {0.6, 0.1}, {0.8, 0.058}, {1, 0.038}}; Clear[x, a, b, c, d] Spline3[Data_] := (k = Length[Data]; atable = Table[Subscript[a, i], {i, 1, k - 3}]; btable = Table[Subscript[b, i], {i, 1, k - 3}]; ctable = Table[Subscript[c, i], {i, 1, k - 3}]; dtable = Table[Subscript[d, i], {i, 1, k - 3}]; Poly = Table[atable[[i]] + btable[[i]]*x + ctable[[i]]*x^2 + dtable[[i]]*x^3, {i, 1, k - 3}]; intermediate = Table[Data[[i + 2]], {i, 1, k - 4}]; Equations = Table[0, {i, 1, 4 (k - 3)}]; Do[Equations[[i]] = (Poly[[i]] /. x -> Data[[i + 2, 1]]) == Data[[i + 2, 2]], {i, 1, k - 4}]; Do[Equations[[i + k - 4]] = (Poly[[i + 1]] /. x -> Data[[i + 2, 1]]) == Data[[i + 2, 2]], {i, 1, k - 4}]; Equations[[2 (k - 4) + 1]] = (Poly[[1]] /. x -> Data[[1, 1]]) == Data[[1, 2]]; Equations[[2 (k - 4) + 2]] = (Poly[[1]] /. x -> Data[[2, 1]]) == Data[[2, 2]]; Equations[[2 (k - 4) + 3]] = (Poly[[k - 3]] /. x -> Data[[k - 1, 1]]) == Data[[k - 1, 2]]; Equations[[2 (k - 4) + 4]] = (Poly[[k - 3]] /. x -> Data[[k, 1]]) == Data[[k, 2]]; Do[Equations[[2 (k - 3) + 2 + i]] = (D[Poly[[i]], x] /. x -> intermediate[[i, 1]]) == (D[Poly[[i + 1]], x] /. x -> intermediate[[i, 1]]), {i, 1, k - 4}]; Do[Equations[[2 (k - 3) + 2 + k - 4 + i]] = (D[Poly[[i]], {x, 2}] /. x -> intermediate[[i, 1]]) == (D[Poly[[i + 1]], {x, 2}] /.x -> intermediate[[i, 1]]), {i, 1, k - 4}]; un = Flatten[{atable, btable, ctable, dtable}]; f = Solve[Equations]; Poly = Poly /. f[[1]]; xs = intermediate; xs = Prepend[Append[xs, Data[[k]]], Data[[1]]]; pf = Table[{Poly[[i]], xs[[i, 1]] <= x <= xs[[i + 1, 1]]}, {i, 1, k - 3}]; Piecewise[pf]) y = Spline3[Data] yactual = 1/(1 + 25 x^2); a = Plot[y, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y3"}, PlotLabel -> "cubic spline"] a = Plot[{y, yactual}, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y3", "yactual"}, PlotLabel -> "cubic spline"]
import numpy as np import sympy as sp import matplotlib.pyplot as plt Data = [[-1, 0.038], [-0.8, 0.058], [-0.60, 0.10], [-0.4,0.20], [-0.2, 0.50], [0, 1], [0.2, 0.5], [0.4, 0.2], [0.6, 0.1], [0.8, 0.058], [1, 0.038]] def Spline3(x_val, Data): k = len(Data) a,b,c,d = sp.symbols('a:d', cls=sp.Function) x = sp.symbols('x') atable = [a(i) for i in range(k - 3)] btable = [b(i) for i in range(k - 3)] ctable = [c(i) for i in range(k - 3)] dtable = [d(i) for i in range(k - 3)] Poly = [atable[i] + btable[i]*x + ctable[i]*x**2 + dtable[i]*x**3 for i in range(k - 3)] intermediate = [Data[i + 2] for i in range(k - 4)] Equations = [0 for i in range(4*(k - 3))] for i in range(k - 4): Equations[i] = (Poly[i].subs(x, Data[i + 2][0])) - Data[i + 2][1] for i in range(k - 4): Equations[i + k - 4] = (Poly[i + 1].subs(x, Data[i + 2][0])) - Data[i + 2][1] Equations[2*(k - 4)] = (Poly[0].subs(x, Data[0][0])) - Data[0][1] Equations[2*(k - 4) + 1] = (Poly[0].subs(x, Data[1][0])) - Data[1][1] Equations[2*(k - 4) + 2] = (Poly[k - 4].subs(x, Data[k - 2][0])) - Data[k - 2][1] Equations[2*(k - 4) + 3] = (Poly[k - 4].subs(x, Data[k - 1][0])) - Data[k - 1][1] for i in range(k - 4): Equations[2*(k - 3) + i + 2] = (Poly[i].diff(x).subs(x, intermediate[i][0])) - (Poly[i + 1].diff(x).subs(x, intermediate[i][0])) for i in range(k - 4): Equations[2*(k - 2) + k + i - 4] = (Poly[i].diff(x, 2).subs(x, intermediate[i][0])) - (Poly[i + 1].diff(x, 2).subs(x, intermediate[i][0])) un = atable + btable + ctable + dtable f = list(sp.linsolve(Equations, un)) vars = {} for i in range(len(un)): vars[un[i]] = f[0][i]; Poly = [Poly[i].subs(vars) for i in range(len(Poly))] xs = np.vstack([Data[0], intermediate, Data[k - 1]]) display([['{:45}'.format(str(sp.N(Poly[i],4))), '{:15}'.format(str(round(xs[i][0],3))+' <=x<= '+str(round(xs[i + 1][0],3)))] for i in range(k - 3)]) return np.piecewise(x_val, [(x_val >= xs[i][0])&(x_val <= xs[i + 1][0]) for i in range(len(Poly))], [sp.lambdify(x, Poly[i]) for i in range(len(Poly))]) x = np.arange(-1, 1, 0.01) y = Spline3(x, Data) yactual = 1/(1 + 25*x**2) plt.plot(x, yactual, label="yactual") plt.plot(x, y, label="y3") plt.scatter([point[0] for point in Data],[point[1] for point in Data], c='k') plt.title("Cubic spline") plt.xlabel("x"); plt.ylabel("y3") plt.legend(); plt.grid(); plt.show()
The following MATLAB code fits a cubic spline interpolation function to data extracted using the Runge function:

Figure 10. Behaviour of the cubic spline interpolation scheme when applied to the Runge function data points
Example
Use the cubic spline interpolation scheme presented above to find the interpolating function given the following 5 data points (-1, 0.038) (-0.8, 0.058), (-0.60, 0.10), (-0.4, 0.20), (-0.2, 0.5)
Solution
There are data points, therefore, there are
cubic polynomials required for interpolation. These are:
To find the associated unknown coefficients, the following system of 8 linear equations is solved:
Solving the above system, yields the following interpolation function:
Figure 11 shows the resulting interpolation function along with the five data points. The Mathematica code is provided below the figure.

Figure 11. Cubic spline interpolation applied to five data points.
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4, 0.20}, {-0.2, 0.50}}; Clear[x, a, b, c, d] Spline3[Data_] := (k = Length[Data]; atable = Table[Subscript[a, i], {i, 1, k - 3}]; btable = Table[Subscript[b, i], {i, 1, k - 3}]; ctable = Table[Subscript[c, i], {i, 1, k - 3}]; dtable = Table[Subscript[d, i], {i, 1, k - 3}]; Poly = Table[atable[[i]] + btable[[i]]*x + ctable[[i]]*x^2 + dtable[[i]]*x^3, {i, 1, k - 3}]; intermediate = Table[Data[[i + 2]], {i, 1, k - 4}]; Equations = Table[0, {i, 1, 4 (k - 3)}]; Do[Equations[[i]] = (Poly[[i]] /. x -> Data[[i + 2, 1]]) == Data[[i + 2, 2]], {i, 1, k - 4}]; Do[Equations[[i + k - 4]] = (Poly[[i + 1]] /. x -> Data[[i + 2, 1]]) == Data[[i + 2, 2]], {i, 1, k - 4}]; Equations[[2 (k - 4) + 1]] = (Poly[[1]] /. x -> Data[[1, 1]]) == Data[[1, 2]]; Equations[[2 (k - 4) + 2]] = (Poly[[1]] /. x -> Data[[2, 1]]) == Data[[2, 2]]; Equations[[2 (k - 4) + 3]] = (Poly[[k - 3]] /. x -> Data[[k - 1, 1]]) == Data[[k - 1, 2]]; Equations[[2 (k - 4) + 4]] = (Poly[[k - 3]] /. x -> Data[[k, 1]]) == Data[[k, 2]]; Do[Equations[[2 (k - 3) + 2 + i]] = (D[Poly[[i]], x] /. x -> intermediate[[i, 1]]) == (D[Poly[[i + 1]], x] /. x -> intermediate[[i, 1]]), {i, 1, k - 4}]; Do[Equations[[2 (k - 3) + 2 + k - 4 + i]] = (D[Poly[[i]], {x, 2}] /. x -> intermediate[[i, 1]]) == (D[Poly[[i + 1]], {x, 2}] /. x -> intermediate[[i, 1]]), {i, 1, k - 4}]; f = Solve[Equations]; Poly = Poly /. f[[1]]; xs = intermediate; xs = Prepend[Append[xs, Data[[k]]], Data[[1]]]; pf = Table[{Poly[[i]], xs[[i, 1]] <= x <= xs[[i + 1, 1]]}, {i, 1, k - 3}]; Piecewise[pf]) Equations // MatrixForm y = Spline3[Data] a = Plot[y, {x, -1, 0}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y3"}, PlotLabel -> "cubic spline"]
import numpy as np import sympy as sp import matplotlib.pyplot as plt Data = [[-1, 0.038], [-0.8, 0.058], [-0.60, 0.10], [-0.4, 0.20], [-0.2, 0.50]] def Spline3(x_val, Data): k = len(Data) a,b,c,d = sp.symbols('a:d', cls=sp.Function) x = sp.symbols('x') atable = [a(i) for i in range(k - 3)] btable = [b(i) for i in range(k - 3)] ctable = [c(i) for i in range(k - 3)] dtable = [d(i) for i in range(k - 3)] Poly = [atable[i] + btable[i]*x + ctable[i]*x**2 + dtable[i]*x**3 for i in range(k - 3)] intermediate = [Data[i + 2] for i in range(k - 4)] Equations = [0 for i in range(4*(k - 3))] for i in range(k - 4): Equations[i] = (Poly[i].subs(x, Data[i + 2][0])) - Data[i + 2][1] for i in range(k - 4): Equations[i + k - 4] = (Poly[i + 1].subs(x, Data[i + 2][0])) - Data[i + 2][1] Equations[2*(k - 4)] = (Poly[0].subs(x, Data[0][0])) - Data[0][1] Equations[2*(k - 4) + 1] = (Poly[0].subs(x, Data[1][0])) - Data[1][1] Equations[2*(k - 4) + 2] = (Poly[k - 4].subs(x, Data[k - 2][0])) - Data[k - 2][1] Equations[2*(k - 4) + 3] = (Poly[k - 4].subs(x, Data[k - 1][0])) - Data[k - 1][1] for i in range(k - 4): Equations[2*(k - 3) + i + 2] = (Poly[i].diff(x).subs(x, intermediate[i][0])) - (Poly[i + 1].diff(x).subs(x, intermediate[i][0])) for i in range(k - 4): Equations[2*(k - 2) + k + i - 4] = (Poly[i].diff(x, 2).subs(x, intermediate[i][0])) - (Poly[i + 1].diff(x, 2).subs(x, intermediate[i][0])) un = atable + btable + ctable + dtable f = list(sp.linsolve(Equations, un)) vars = {} for i in range(len(un)): vars[un[i]] = f[0][i]; Poly = [Poly[i].subs(vars) for i in range(len(Poly))] xs = np.vstack([Data[0], intermediate, Data[k - 1]]) display([['{:45}'.format(str(sp.N(Poly[i],4))), '{:15}'.format(str(round(xs[i][0],3))+' <=x<= '+str(round(xs[i + 1][0],3)))] for i in range(k - 3)]) return np.piecewise(x_val, [(x_val >= xs[i][0])&(x_val <= xs[i + 1][0]) for i in range(len(Poly))], [sp.lambdify(x, Poly[i]) for i in range(len(Poly))]) x = np.arange(-1, -0.2, 0.001) y = Spline3(x, Data) yactual = 1/(1 + 25*x**2) plt.plot(x, yactual, label="yactual") plt.plot(x, y, label="y3") plt.scatter([point[0] for point in Data],[point[1] for point in Data], c='k') plt.title("Cubic spline") plt.xlabel("x"); plt.ylabel("y3") plt.legend(); plt.grid(); plt.show()
The following MATLAB code fits a cubic spline interpolation function to the data from this example: