Suppose we want to fit a smooth curve which actually goes through n + I given data points when n is quite large. Since an interpolating polynomial of correspondingly high degree n tends to be highly oscillatory and therefore is likely to give an unsatisfactory fit, at least in certain locations, an interpolation is often constructed by linking lower degree polynomials (piecewise polynomials) at certain or all of the given data points (called nodes or knots). This interpolation is smooth if we also insist that the piecewise polynomials have matching derivatives at the nodes, and this smoothness is enhanced by matching higher order derivatives.
Let the data points (x0, (f0), (x1(f1), . . ., (xn), (fn), be ordered according to their magnitude so that
We will seek a function S which is a polynomial of degree d on each subinterval [xj - 1,,xj], j = 1, 2, …,n, where
.
In order to achieve maximum smoothness at the nodes, we shall allow S to have up to d - 1 continuous derivatives. Such functions are referred to as splines. An example of a spline for the linear (d = 1) case is the polygon (Curve A) of Figure 11(b) in STEP 26. It is clear that this spline is continuous, but does not have a continuous first derivative. In practice, most popular are cubic splines, constructed from polynomials of degree three with continuous first and second derivatives at the nodes and discussed below in detail. Figure 13 below shows an example of a cubic spline S for n = 5. (The data points are taken from the table in Step 26) We see that S passes through all the data points. The function Sj on the subinterval [xj-1, xj] is a cubic. As has already been indicated, the first and second derivatives of Sj and Sj+1 match at (xj, fj), the point where they meet.
The term spline refers to the thin flexible rods, used in the past by draughtsmen to draw smooth curves, for example, in ship design. The graph of a cubic spline approximates the shape which forms when such a rod is forced to pass through given n + 1 nodes and corresponds, according to the theory of bending of thin rods, to minimum strain energy.
As has been indicated, a cubic spline S is constructed by fitting a cubic to each
subinterval [xj-1, xj]
for j = 1, 2, . . , n, whence it is
convenient to assume that S has values Sj(x)
for 
, where
FIGURE 13. Schematic example of a cubic spline over subintervals [x0, x1], [x1, x2], [x2, x3], [x3, x4], and [x4, x5], each function Sj on [xj-1, xj] being a cubic.
We now impose the condition S(xj) = fj, whence aj = fj for j = 1, 2, . . , n. For S to be continuous and to have continuous first and second derivatives at the given data points, we require
for j=1, 2, .. , n - 1.
Since we have a cubic with the four unknowns (aj, bj, cj, dj) on each of the n subintervals, and so a total of 4n unknowns, we need 4n equations to specify them. The requirements S(xj) = fj, j=0,1, 2, . . , n yield n + 1 equations, while 3(n -1) equations arise from the continuity requirement on S and its first two derivatives given above. This yields a total of n+1+3(n-1)=4n-2 equations, whence we need to impose two more conditions to specify S completely. The choice of these two extra conditions determines the type of the cubic spline obtained. Two common options are:
for
some given constants 
and 
If the values of f'(x0)
and f'''(xn) are known,
then and can be
set to these values.We shall not go here into the algebraic details; however, it turns out that if we write hj = xj - xj-1 and mj = S(xj), then the coefficients of Sj are given by
The spline is thus
determined by the values of 
which depend on whether we want a natural or a clamped
cubic spline.
For a natural cubic spline, we have m0 = mn = 0, and the equations:
for j = 1, 2, . . .
, n - 1. (Note that if all the values of hj
are the same, then the right-hand side of this
last equation is just 
. Setting 
, these
linear equations become the (n -1) ´ (n -
1) system:
,
where
,
Note that the coefficient matrix has non-zero entries only on the leading diagonal and the two sub-diagonals either side of it. Such a system is called a tri-diagonal system. Since most of the entries below the leading diagonal are zero, it is possible to modify Gauss elimination (cf. STEP 11) to produce a very efficient method for solving tri-diagonal systems.
For the clamped boundary conditions, the equations are:
,
and
,
It may be verified that these equations for m0, m1,. . . , mn can be written as an (n + 1 ) x (n + 1 ) tridiagonal system.
We will fit a natural cubic spline to data, taken from STEP 26:
,
Since the values of the xj are equally spaced, we find hj = 1 for j = 1, 2, . . , 5. Also, m0=m5 = 0 and the remaining values m1, m2, m3, and m4 satisfy the linear system:
.
Using Gauss elimination to solve this system to 5D, we find:
.
Calculating the coeficients, we then find that the spline S is given by:
,
where
,
The data points and the cubic spline have already been displayed in Figure 13.
The next example demonstrates graphically a comment made earlier, namely that the behaviour of interpolating polynomials of high degree tends to be very oscillatory. For this purpose, we consider the data of the function f(x) = 10/(1 + x2):
.
It is readily verified that the interpolating polynomial of degree 6 is given by
.
The function f is plotted as a solid line along with the interpolating polynomial as a dashed line in Figure 14 below. (Since both f and P6 are symmetric about the x-axis, only the (0,3) sectiont of the graph has been displayed.) The oscillatory behaviour of the interpolating polynomial of degree 6 is obvious.
Now fit a natural cubic spline to the data. In order to do so, we must solve the linear system
.
Gauss elimination yields (to 5D )
.
We now find that the natural spline S is symmetric about the y-axis. On [0, 3], it is given by
,
where
.
This spline is also plotted in Figure 14 as a dotted line. It is clear that it is a much better approximation to f than the interpolating polynomial.
FIGURE 14. The function f (x) =10/(1 +x2) (solid line) approximated by an interpolating polynomial (dashed line) and a natural cubic spline (dotted line).
Given the data points (0,1 ), ( 1, 4), (2,15), and (3, 40), find the natural cubic spline fitting these data. Use the spline to estimate the value of y at x = 2.3.