All previously discussed methods of polynomial interpolation fit
a set of given points
by an
nth degree polynomial, and a higher degree polynomial is needed to
fit a larger set of data points. A major drawback of such methods
is *overfitting*, as domonstrated by the following example.

**Example:**

Vased on equally spaced points from to
with increment of 1, a function
can be approximated
by any of the interpolation methods discussed above by polynomial
of degree , as shown in the figure below. We note that the
approximation is very poor towards to the two ends where the error
is disappointingly high. This is known as
*Runge's phenomenon*, indicating the fact that higher degree
polynomial interpolation does not necessarily always produce more
accurate result, as the degree of the interpolating polynomial may
become unnecessarily high and the polynomial may become oscillatory.

This Runge's phenominon is a typical example of overfitting, due to an excessively complex model with too many parameters relative to the observed data, here specifically a polynomial of a degree too high (requiring too many coefficients) to model the given data points.

Now we consider a different method of spline interpolation, which
fits the given points by a piecewise polynomial function ,
known as the *spline*, a composite function formed by
low-degree polynomials each fitting in the interval
between and
:

(61) |

- For to be continuous, two consecutive polynomials
and
must join at :
(62) (63) - For to be smooth, they need to have the same
derivatives at the point they joint, i.e.,
(64)

In the following we consider approximating between any two consecutive points and by a linear, quadratic, and cubic polynomial (of first, second, and third degree).

**Linear spline**: with two parameters and can only satisfy the following two equations required for to be continuous:(65) (66)

(67)

which is represented in the form of by the first expression, or a linear interpolation of the two end points and in the second expression.As , the linear spline is continuous at . But as in general

(68) **Quadratic spline:**with three parameters and can satisfy the following three equations required for to be smooth () as well as continuous:(69) (70) (71) (72) (73) (74) (75) (76) **Cubic spline:**with four parameters , and can satisfy the following four equations required for to be continuous and smooth ():and (77) (78) (79) (80) (81) (82)

(83)

To find , we take derivative of and rearrange terms to get

(84)

which, when evaluated at and , becomes:

(85)

Replacing by in the second equation, we also get(86) (87) (88) (89) *divided differences*. We can rewrite the equation above as(90) (91) - Assume the first order derivatives at both ends
and are known. Specially
is called
*clamped boundary condition*. At the front end, we set(92) (93) (94) (95) (96) - Alternatively, we can also assume and
are known. Specially
is called
*natural boundary condition*. Now we can simply get and and solve the following system for the unknowns :(97)

- Assume the first order derivatives at both ends
and are known. Specially
is called

**Example:**

A function is sampled at the following points:

The interpolation results based on linear, quadratic and cubic splines are shown in the figure below, together with the original function , and the interpolating polynomials , used as the ith segment of between and .

For the quadratic interpolation, based on we get . For the cubic interpolation, we solve the following equation

The errors of these three methods are , , and , respectively. Obviously the higher the degree of the interpolating polynomial, the higher the accuracy. The error of the cubic spline method is significantly smaller than of the polynomial interpolation.

The Matlab code that implements the cubic spline method is listed below.

function [S C]=Spline3(u,x,y,dya,dyb) % vectors x and y contain n+1 points and the corresponding function values % vector u contains all discrete samples of the continuous argument of f(x) % dya and dyb are the derivatives f'(x_0) and f'(x_n), respectively n=length(x); % number of interpolating points k=length(u); % number of discrete sample points C=zeros(n,k); % the n-1 cubic interpolating polynomials A=2*eye(n); % coefficient matrix on left-hand side A(1,2)=1; A(n,n-1)=1; d=zeros(n,1); % vector on right-hand side d(1)=((y(2)-y(1))/(x(2)-x(1))-dya)/h0; % first element of d for i=2:n-1 h0=x(i)-x(i-1); h1=x(i+1)-x(i); h2=x(i+1)-x(i-1); A(i,i-1)=h0/h2; A(i,i+1)=h1/h2; d(i)=((y(i+1)-y(i))/h1-(y(i)-y(i-1))/h0)/h2; % 2nd divided difference end d(n)=(dyb-(y(n)-y(n-1))/h1)/h1; % last element of d M=6*inv(A)*d; % solving linear equation system for M's for i=2:n h=x(i)-x(i-1); x0=u-x(i-1); x1=x(i)-u; C(i-1,:)=(x1.^3*M(i-1)+x0.^3*M(i))/6/h... % the ith cubic polynomial -(M(i-1)*x1+M(i)*x0)*h/6+(y(i-1)*x1+y(i)*x0)/h; idx=find(u>x(i-1) & u<=x(i)); % indices between x(i-1) and x(i) S(idx)=C(i-1,idx); % constructing spline by cubic polynomials end end

**Example:**

The function used before is now approximated by both the Newton's method and the cubic spline method, with very different results as shown below. The Runge's phenomenon suffered by Newton's method is certainly avoided by the cubic spline method.