Introduction to Numerical Analysis: Piecewise Interpolation
Introduction
In the previous section, it was shown that when the order of the interpolating polynomial increases — which is natural when there is a large number of data points — the interpolating polynomial function can highly oscillate or fluctuate between the data points. With a large number of data points, it is better to keep the order of interpolation low. This can be done using piecewise interpolation using splines. Splines are polynomials on subintervals that are connected in a smooth manner. The order of the polynomials is not associated with the number of data points in this case, so, we can have data points, and use a spline of degree
where
. Assuming
data points
such that
, a spline function
of degree
satisfies the following:
- On each subinterval
,
is a polynomial of degree
- The derivatives of
up to and including the
derivative are continuous on the domain
.
The data points, when chosen such that the spline function passes through them, are called knots. For example, a spline of degree 0 is a piecewise constant function constructed by simply assuming to have a constant value on any subinterval
(Figure 1). A spline of degree 1 is a piecewise linear (piecewise affine) function constructed by connecting a straight line between every two data points (Figure 1). A spline of degree 2 is a piecewise quadratic function on each subinterval. The first derivative exists and is continuous, but it is not smooth. As shown in Figure 2, the first derivative is a broken line at the knots and thus the derivatives at the knots do not exist. A spline of degree 3 is a piecewise cubic function on each subinterval. The first and second derivatives are continuous functions. The second derivative is not smooth. Figure 1 shows the resulting 0 order, linear, quadratic, and cubic spline interpolating functions for the data points generated by the Runge function described in the previous section. Figure 2 on the other hand shows the generated derivatives. The 0 order spline function gives a zero derivative every where. The linear spline function has derivatives that are constant on each subinterval. The quadratic spline gives derivatives that are not smooth at the data points. The cubic spline provides smooth derivatives that are very close to the actual derivative of the Runge function that was used to generate the data!
The chosen order of the spline function depends on the accuracy sought. Linear splines can provide acceptable interpolation between the knots. However, their derivatives are constant. Linear splines are very useful when the data points are very close to each other, and there is no interest in their derivatives. Quadratic and cubic splines provide very smooth interpolation of the data points. However, the derivatives of the cubic splines are “smoother” than the quadratic splines. The choice of the spline order depends on the oscillatory nature of the data points to be interpolated.
The “Interpolation” built-in function generates an interpolation function through the data. To choose splines, the option: Method->”Spline” has to be specified as well as the interpolation order. The following code generates the 0 order, linear, quadratic, and cubic spline functions through the data points of the Runge function shown in the previous section. The code also compares the first derivatives of the spline functions with the actual first derivative of the Runge function.
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}}; y0 = Interpolation[Data, Method -> "Spline", InterpolationOrder -> 0]; y1 = Interpolation[Data, Method -> "Spline", InterpolationOrder -> 1]; y2 = Interpolation[Data, Method -> "Spline", InterpolationOrder -> 2]; y3 = Interpolation[Data, Method -> "Spline", InterpolationOrder -> 3]; yactual = 1/(1 + 25 x^2); a = Plot[{y0[x], yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y0"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y0", "yactual"}, PlotLabel -> "0 order splines"]; b = Plot[{y1[x], yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y1"}, ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y1", "yactual"}, PlotLabel -> "Linear splines"]; c = Plot[{y2[x], yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y2"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y2", "yactual"}, PlotLabel -> "Quadratic splines"]; d = Plot[{y3[x], yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y3", "yactual"}, PlotLabel -> "Cubic splines"]; Grid[{{a, b}, {c, d}}] y0prime = D[y0[x], x]; y1prime = D[y1[x], x]; y2prime = D[y2[x], x]; y3prime = D[y3[x], x]; yactualprime = D[yactual, x]; a = Plot[{y0prime, yactualprime}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y0"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y0'", "y'actual"}, PlotLabel -> "Derivative of 0 order splines"]; b = Plot[{y1prime, yactualprime}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y1"}, ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y1'", "y'actual"}, PlotLabel -> "Linear splines"]; c = Plot[{y2prime, yactualprime}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]},AxesLabel -> {"x", "y2"}, ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y2'", "y'actual"}, PlotLabel -> "Quadratic splines"]; d = Plot[{y3prime, yactualprime}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y3"}, ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y3'", "y'actual"}, PlotLabel -> "Cubic splines"]; Grid[{{a, b}, {c, d}}]

Figure 1. Piecewise interpolation using: 0 order, linear, quadratic, and cubic splines

Figure 2. The derivatives of piecewise interpolation using 0 order, linear, quadratic, and cubic splines. The derivative of the 0 order splines is the zero function.
0 Degree Spline Interpolation
Given a set of data points
, a spline of degree 0 can be defined as:
In the previous spline definition, the interpolated value in the interval is equal to the data value
. Alternatively, the interpolated value in the interval
can be set equal to the data value
. In this case,
can be defined as follows:
Piecewise Linear (Piecewise Affine) Spline Interpolation
Given a set of data points
, a piecewise linear (piecewise affine) spline can be defined as:
The data points have
intervals. The linear function for each interval is defined using two coefficients, and therefore, we need to find
coefficients
. The coefficients
and
can be found by solving the two equations:
Therefore:
Example
Consider the following 11 data points: (-1.00,0.038),(-0.80,0.058),(-0.60,0.100),(-0.4,0.200),(-0.20,0.500),(0.00,1.000),(0.20,0.500),(0.40,0.200),(0.60,0.100),(0.80,0.0580),(1.00,0.038) which were generated using the Runge function. Find the explicit representation of the linear spline interpolating function for these data points.
Solution
There are 11 data points surrounding 10 intervals. For each interval , a set of coefficients
and
are required for the linear interpolation. For the first interval, the coefficients are:
Similarly:
Repeating for the other coefficients, the required explicit equation has the form:
Notice that for the procedure to work, the components of the data points have to satisfy:
. The following Mathematica code utilizes a user defined procedure “Spline1” that creates the required piecewise linear function.
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}}; Spline1[Data_] := ( k1 = Length[Data] - 1; atable = Table[(Data[[i, 2]] - Data[[i + 1, 2]])/(Data[[i, 1]] -Data[[i + 1, 1]]), {i, 1, k1}]; btable = Table[(Data[[i + 1, 2]]*Data[[i, 1]] -Data[[i, 2]]*Data[[i + 1, 1]])/(Data[[i, 1]] - Data[[i + 1, 1]]), {i, 1, k1}]; pf = Table[{atable[[i]] x + btable[[i]],Data[[i, 1]] <= x <= Data[[i + 1, 1]]}, {i, 1, k1}]; Piecewise[pf] ) y = Spline1[Data] a = Plot[{y, yactual}, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y1"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y1", "yactual"}, PlotLabel -> "Linear splines"]
Quadratic Spline Interpolation
For quadratic spline interpolation, we present two possible quadratic interpolation schemes.
Scheme 1:
In the first scheme, the intervals between the data points are used as intervals on which a quadratic function is defined. On the intermediate nodes, the slope of the parabola on the left is equated to the slope of the parabola on the right (Figure 3). Consider a set of data points
, with
intervals. We seek to find a quadratic polynomial
on the
interval
. In total, there are
coefficients. On each interval, the two equations
and
ensure that the spline passes through the data points which also ensure the continuity of the splines. These continuity equations will provide
equations. In addition, the connection of the quadratic polynomials at intermediate nodes will be assumed to be smooth, and so at every intermediate point the first derivative is assumed to be continuous, that is
. These first-derivative continuity equations provide another
equations. Finally, the second derivative at the starting point can be set to 0 to provide one additional condition so that the number of equations would be equal to
which is equal to the number of unknowns. The following are the resulting equations:

Figure 3. Scheme 1 of piecewise quadratic interpolation with continuity in the interpolating function and its first derivative
The following Mathematica code implements this procedure for a set of data (the procedure will only work if ).
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}}; Spline2[Data_] := (k1 = Length[Data] - 1; y = Table[If[i <= 2 k1,If[EvenQ[i], Data[[i/2 + 1, 2]], Data[[(i + 1)/2, 2]]], 0], {i,1, 3 k1}]; M = Table[0, {i, 1, 3 k1}, {j, 1, 3 k1}]; For[i = 1, i <= k1, i = i + 1, M[[2 i - 1, 3 (i - 1) + 1]] = M[[2 i, 3 (i - 1) + 1]] = 1; M[[2 i - 1, 3 (i - 1) + 2]] = Data[[i, 1]]; M[[2 i - 1, 3 (i - 1) + 1 + 2]] = Data[[i, 1]]^2; M[[2 i, 3 (i - 1) + 1 + 1]] = Data[[i + 1, 1]]; M[[2 i, 3 (i - 1) + 1 + 2]] = Data[[i + 1, 1]]^2]; For[i = 1, i <= k1 - 1, i = i + 1, M[[2 k1 + i, 3 (i - 1) + 1 + 1]] = 1; M[[2 k1 + i, 3 (i - 1) + 1 + 4]] = -1; M[[2 k1 + i, 3 (i - 1) + 1 + 2]] = 2 Data[[i + 1, 1]]; M[[2 k1 + i, 3 (i - 1) + 1 + 5]] = -2 Data[[i + 1, 1]]]; M[[3 k1, 3]] = 2; Coef = Inverse[M].y; pf = Table[{Coef[[3 (i - 1) + 1]] + Coef[[3 (i - 1) + 2]]*x + Coef[[3 (i - 1) + 3]]*x^2, Data[[i, 1]] <= x <= Data[[i + 1, 1]]}, {i, 1, k1}]; Piecewise[pf]) y = Spline2[Data] M // MatrixForm yactual=1/(1+25x^2); a = Plot[{y, yactual}, {x, -1, 1},Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y1"},ImageSize -> Medium, PlotRange -> All,PlotLegends -> {"y2", "yactual"}, PlotLabel ->"Quadratic spline"]
Example
Use scheme 1 of quadratic interpolation to find the interpolating function given the following 4 data points (-1, 0.038) (-0.8, 0.058), (-0.60, 0.10), (-0.4, 0.20)
Solution
The four data points () with three intervals (
), therefore, there are
unknowns. The following are the nine equations according to scheme 1 to be solved to find the unknowns:
Solving the above equations produces the following interpolation function:
Figure 4 shows the resulting interpolation function along with the four data points. The Mathematica code is provided below the figure.

Figure 4. Scheme 1 quadratic interpolation applied to four data points
Data = {{-1, 0.038}, {-0.8, 0.058}, {-0.60, 0.10}, {-0.4, 0.20}}; Spline2[Data_] := (k1 = Length[Data] - 1; y = Table[ If[i <= 2 k1, If[EvenQ[i], Data[[i/2 + 1, 2]], Data[[(i + 1)/2, 2]]], 0], {i, 1, 3 k1}]; M = Table[0, {i, 1, 3 k1}, {j, 1, 3 k1}]; For[i = 1, i <= k1, i = i + 1, M[[2 i - 1, 3 (i - 1) + 1]] = M[[2 i, 3 (i - 1) + 1]] = 1; M[[2 i - 1, 3 (i - 1) + 2]] = Data[[i, 1]]; M[[2 i - 1, 3 (i - 1) + 1 + 2]] = Data[[i, 1]]^2; M[[2 i, 3 (i - 1) + 1 + 1]] = Data[[i + 1, 1]]; M[[2 i, 3 (i - 1) + 1 + 2]] = Data[[i + 1, 1]]^2]; For[i = 1, i <= k1 - 1, i = i + 1, M[[2 k1 + i, 3 (i - 1) + 1 + 1]] = 1; M[[2 k1 + i, 3 (i - 1) + 1 + 4]] = -1; M[[2 k1 + i, 3 (i - 1) + 1 + 2]] = 2 Data[[i + 1, 1]]; M[[2 k1 + i, 3 (i - 1) + 1 + 5]] = -2 Data[[i + 1, 1]]]; M[[3 k1, 3]] = 2; Coef = Inverse[M].y; pf = Table[{Coef[[3 (i - 1) + 1]] + Coef[[3 (i - 1) + 2]]*x + Coef[[3 (i - 1) + 3]]*x^2, Data[[i, 1]] <= x <= Data[[i + 1, 1]]}, {i, 1, k1}]; Piecewise[pf]) y = Spline2[Data] M // MatrixForm a = Plot[{y}, {x, -1, -0.4}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y2"},ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y2"}, PlotLabel -> "Quadratic spline"]
Disadvantages
In certain cases, such as fluctuating data, this scheme might lead to some oscillations within the intervals between the data points. As an example, when this scheme is used to fit through the 11 data points of the Runge function given in the previous section, a nice smooth interpolating function is produced for the first few intervals (Figure 5). The condition of having zero second derivative on the initial interval causes the interpolating in the first interval to follow the line connecting the first two data points. However, in the last five intervals, the interpolating function oscillates and deviates away from the expected behaviour.

Figure 5. Behaviour of scheme 1 for quadratic interpolation of the Runge function data points.
Scheme 2
For the second scheme, the intervals where the quadratic interpolation functions are defined, are chosen different from the intervals between the data points in a manner to ensure symmetry of the quadratic interpolation. The new intervals where the piecewise quadratic functions are defined are bounded by the points (Figure 6). Consider a set of
data points
, with
intervals. On each of the
intervals bounded by
a quadratic polynomial
is defined. In total, there are
coefficients. First,the quadratic polynomials have to pass through the
data poitns resulting in
equations. These quadratic polynomials have to be continuous and differentiable at the
intermediate points that are the bounds of the intervals resulting in
equations. Therefore, in total, there are
equations which is equal to the number of unknowns. The following are the resulting equations:

Figure 6. Scheme 2 of piecewise quadratic interpolation with continuity in the interpolating function and its first derivative
The following Mathematica code implements this procedure for a set of data (the procedure will only work if ). As shown in Figure 7, the oscillations produced by scheme 1 disappear if this scheme is used for interpolating the data points of the Runge function given in the previous section.
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}}; Spline2[Data_] := (k = Length[Data]; NewData = Table[0, {i, 1, k - 3}]; Do[NewData[[i]] = (Data[[i + 1]] + Data[[i + 2]])/2, {i, 1, k - 3}]; NNewData = Prepend[Append[NewData, Data[[k]]], Data[[1]]]; y = Table[If[i <= k - 1 + 1, Data[[i, 2]], 0], {i, 1, 3 k - 6}]; M = Table[0, {i, 1, 3 k - 6}, {j, 1, 3 k - 6}]; M[[1, 1]] = M[[2, 1]] = 1; M[[1, 2]] = Data[[1, 1]]; M[[1, 3]] = Data[[1, 1]]^2; M[[2, 2]] = Data[[2, 1]]; M[[2, 3]] = Data[[2, 1]]^2; M[[k - 1, 3 k - 8]] = M[[k, 3 k - 8]] = 1; M[[k - 1, 3 k - 7]] = Data[[k - 1, 1]]; M[[k - 1, 3 k - 6]] = Data[[k - 1, 1]]^2; M[[k, 3 k - 7]] = Data[[k, 1]]; M[[k, 3 k - 6]] = Data[[k, 1]]^2; For[i = 3, i <= k - 2, i = i + 1, M[[i, 3 (i - 2) + 1]] = 1; M[[i, 3 (i - 2) + 2]] = Data[[i, 1]]; M[[i, 3 (i - 2) + 3]] = Data[[i, 1]]^2]; For[i = 1, i <= k - 3, i = i + 1, M[[k + i, 3 (i - 1) + 1]] = 1; M[[k + i, 3 (i - 1) + 4]] = -1; M[[k + i, 3 (i - 1) + 2]] = NewData[[i, 1]]; M[[k + i, 3 (i - 1) + 5]] = -NewData[[i, 1]]; M[[k + i, 3 (i - 1) + 3]] = NewData[[i, 1]]^2; M[[k + i, 3 (i - 1) + 6]] = -1*NewData[[i, 1]]^2; M[[k + i + k - 3, 3 (i - 1) + 2]] = 1; M[[k + i + k - 3, 3 (i - 1) + 5]] = -1; M[[k + i + k - 3, 3 (i - 1) + 3]] = 2*NewData[[i, 1]]; M[[k + i + k - 3, 3 (i - 1) + 6]] = -2*NewData[[i, 1]]]; Coef = Inverse[M].y; pf = Table[{Coef[[3 (i - 1) + 1]] + Coef[[3 (i - 1) + 2]]*x + Coef[[3 (i - 1) + 3]]*x^2,NNewData[[i, 1]] <= x <= NNewData[[i + 1, 1]]}, {i, 1, k - 2}]; Piecewise[pf]) y = Spline2[Data] yactual = 1/(1 + 25 x^2); M//MatrixForm a = Plot[{y, yactual}, {x, -1, 1}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y2"},ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y2", "yactual"}, PlotLabel -> "Quadratic spline"]

Figure 7. Behaviour of scheme 2 of quadratic interpolation when applied to the Runge function data points
Example
Use scheme 2 of quadratic interpolation 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
The five data points () with four intervals (
).
intervals for the quadratic interpolation functions are defined between the points
. Therefore, there are
unknowns. The following are the nine equations according to scheme 2 to be solved to find the unknowns:
Solving the above equations produces the following interpolation function:
Figure 8 shows the resulting interpolation function along with the five data points. The Mathematica code is provided below the figure.

Figure 8. Scheme 1 quadratic 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}}; Spline2[Data_] := (k = Length[Data]; NewData = Table[0, {i, 1, k - 3}]; Do[NewData[[i]] = (Data[[i + 1]] + Data[[i + 2]])/2, {i, 1, k - 3}]; NNewData = Prepend[Append[NewData, Data[[k]]], Data[[1]]]; y = Table[If[i <= k - 1 + 1, Data[[i, 2]], 0], {i, 1, 3 k - 6}]; M = Table[0, {i, 1, 3 k - 6}, {j, 1, 3 k - 6}]; M[[1, 1]] = M[[2, 1]] = 1; M[[1, 2]] = Data[[1, 1]]; M[[1, 3]] = Data[[1, 1]]^2; M[[2, 2]] = Data[[2, 1]]; M[[2, 3]] = Data[[2, 1]]^2; M[[k - 1, 3 k - 8]] = M[[k, 3 k - 8]] = 1; M[[k - 1, 3 k - 7]] = Data[[k - 1, 1]]; M[[k - 1, 3 k - 6]] = Data[[k - 1, 1]]^2; M[[k, 3 k - 7]] = Data[[k, 1]]; M[[k, 3 k - 6]] = Data[[k, 1]]^2; For[i = 3, i <= k - 2, i = i + 1, M[[i, 3 (i - 2) + 1]] = 1; M[[i, 3 (i - 2) + 2]] = Data[[i, 1]]; M[[i, 3 (i - 2) + 3]] = Data[[i, 1]]^2]; For[i = 1, i <= k - 3, i = i + 1, M[[k + i, 3 (i - 1) + 1]] = 1; M[[k + i, 3 (i - 1) + 4]] = -1; M[[k + i, 3 (i - 1) + 2]] = NewData[[i, 1]]; M[[k + i, 3 (i - 1) + 5]] = -NewData[[i, 1]]; M[[k + i, 3 (i - 1) + 3]] = NewData[[i, 1]]^2; M[[k + i, 3 (i - 1) + 6]] = -1*NewData[[i, 1]]^2; M[[k + i + k - 3, 3 (i - 1) + 2]] = 1; M[[k + i + k - 3, 3 (i - 1) + 5]] = -1; M[[k + i + k - 3, 3 (i - 1) + 3]] = 2*NewData[[i, 1]]; M[[k + i + k - 3, 3 (i - 1) + 6]] = -2*NewData[[i, 1]]]; Coef = Inverse[M].y; pf = Table[{Coef[[3 (i - 1) + 1]] + Coef[[3 (i - 1) + 2]]*x + Coef[[3 (i - 1) + 3]]*x^2, NNewData[[i, 1]] <= x <= NNewData[[i + 1, 1]]}, {i, 1, k - 2}]; Piecewise[pf]) y = Spline2[Data] M // MatrixForm a = Plot[{y}, {x, -1, -0.2}, Epilog -> {PointSize[Large], Point[Data]}, AxesLabel -> {"x", "y2"}, ImageSize -> Medium, PlotRange -> All, PlotLegends -> {"y2"}, PlotLabel -> "Quadratic spline"]
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.
The following Mathematica code implements this procedure for a set of data (the procedure will only work if ).
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.
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}}; 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"]

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"]