Skip to the content of the web site.

5.8 Cubic Spline Interpolation

Introduction Theory HOWTO Error Analysis Examples Questions Applications in Engineering Matlab Maple


A piecewise linear interpolating function does not suffer the same problems which an interpolating polynomial may exhibit, for example, polynomial wiggle and the Runge phenomenon; however, the interpolating function is only continuous and not differentiable everywhere.

In the previous section on matching derivatives, we saw that we could define not only the values but also the values of higher-order derivatives at various values of x. A spline is peicewise interpolating function and a cubic spline is a piecewise interpolating function where on each segment the interpolating function is a cubic polynomial. By placing constraints on the derivatives and concavity at the interpolating points, we get a function which is twice differentiable.


The interpolating polynomials which have been seen to this point have been defined on for all the n points (x1, y1), (x2, y2), ..., (xn, yn). An alternative approach is to define a different interpolating polynomial on each sub-interval under the assumption that the x values are given in order.

The simplest means is to take each pair of adjacent points and find an interpolating polynomial between the points which using Newton polynomials is

This can be expanded to reduce the number of required operations by reducing it to a form ax + b which can be computed immediately. The reader may note that if the value x = xk + 1 is substituted into the above equation that the value is yk + 1.

A significant issue with piecewise linear interpolation is that the interpolant is not differentiable or smooth. A non-differentiable function can introduce new issues in a system almost as easily as a non-continuous function.


Error Analysis



Applications to Engineering

Very applicable....


The following code finds the

% return two (n-1)-dimensional vectors so that a(k)*x + b(k) is the
% linear interpolating polynomial between points x(k) and x(k + 1)
function [a, b] = piecewise( x, y )
	n = length( x );
	a = (y( 2:n ) - y( 1:(n - 1) )) ./ (x( 2:n ) - x( 1:(n - 1) ));
	b = (y( 1:(n - 1) ).*x( 2:n ) - y( 2:n ).*x( 1:(n - 1) )) ./ (x( 2:n ) - x( 1:(n - 1) ));


Figure 1 shows a Maple plot of the function sin(4x)e-x/2u(x).

> plot( 2*Heaviside(x)*exp(-x/2)*cos(4*x), x = -2..6 );

Figure 2. A Maple plot.

The plot is generated by linearly interpolating a number of points where first fifty points are sampled and then additional subsamples where the points do not appear to be sufficiently smooth according to certain heuristics. The interpolated points are shown in Figure 3.

> plot( 2*Heaviside(x)*exp(-x/2)*cos(4*x), x = -2..6, style = point, symbol = circle );

Figure 3. The points interpolated in making the maple plot.

It is possible to turn off the adaptive subsampling to reveal the first fifty samples as shown in Figure 4.

> plot( 2*Heaviside(x)*exp(-x/2)*cos(4*x), x = -2..6, adaptive = false );

Figure 4. The initial 50 points without adaptive subsampling.

Interpolating an insufficient number of points will result in a less pleasing result shown in Figure 5.

> plot( 2*Heaviside(x)*exp(-x/2)*cos(4*x), x = -2..6, style = point, symbol = circle, adaptive = false );

Figure 5. The linear interpolation of 50 points not using subsampling.

Maple also allows the user to state the number of initially sampled points and require the plotting function to first determine discontinuties.

> plot( 2*Heaviside(x)*exp(-x/2)*cos(4*x), x = -2..6, symbol = circle, style = point, numpoints = 1000, discont = true );

Figure 6. Starting with a larger number of initial points and allowing Maple to determine discontinuities.