5.3: Cubic Spline Interpolation
( \newcommand{\kernel}{\mathrm{null}\,}\)
The n+1 points to be interpolated are again
(x0,y0),(x1,y1),…(xn,yn)
Here, we use n piecewise cubic polynomials for interpolation,
gi(x)=ai(x−xi)3+bi(x−xi)2+ci(x−xi)+di,i=0,1,…,n−1,
with the global interpolation function written as
g(x)=gi(x), for xi≤x≤xi+1.
To achieve a smooth interpolation we impose that g(x) and its first and second derivatives are continuous. The requirement that g(x) is continuous (and goes through all n+1 points) results in the two constraints
gi(xi)=yi,i=0 to n−1gi(xi+1)=yi+1,i=0 to n−1
The requirement that g′(x) is continuous results in
g′i(xi+1)=g′i+1(xi+1),i=0 to n−2
And the requirement that g′′(x) is continuous results in
g′′i(xi+1)=g′′i+1(xi+1),i=0 to n−2
There are n cubic polynomials gi(x) and each cubic polynomial has four free coefficients; there are therefore a total of 4n unknown coefficients. The number of constraining equations from (5.1)-(5.4) is 2n+2(n−1)=4n−2. With 4n−2 constraints and 4n unknowns, two more conditions are required for a unique solution. These are usually chosen to be extra conditions on the first g0(x) and last gn−1(x) polynomials. We will discuss these extra conditions later.
We now proceed to determine equations for the unknown coefficients of the cubic polynomials. The polynomials and their first two derivatives are given by
gi(x)=ai(x−xi)3+bi(x−xi)2+ci(x−xi)+dig′i(x)=3ai(x−xi)2+2bi(x−xi)+cig′′i(x)=6ai(x−xi)+2bi
We will consider the four conditions (5.1)-(5.4) in turn. From (5.1) and (5.5), we have
di=yi,i=0 to n−1,
which directly solves for all of the d-coefficients.
To satisfy (5.2), we first define
hi=xi+1−xi
and
fi=yi+1−yi
Now, from (5.2) and (5.5), using (5.8), we obtain the n equations
aih3i+bih2i+cihi=fi,i=0 to n−1
From (5.3) and (5.6) we obtain the n−1 equations
3aih2i+2bihi+ci=ci+1,i=0 to n−2
From (5.4) and (5.7) we obtain the n−1 equations
3aihi+bi=bi+1i=0 to n−2
It is will be useful to include a definition of the coefficient bn, which is now missing. (The index of the cubic polynomial coefficients only go up to n−1. ) We simply extend (5.11) up to i=n−1 and so write
3an−1hn−1+bn−1=bn
which can be viewed as a definition of bn.
We now proceed to eliminate the sets of a - and c-coefficients (with the d-coefficients already eliminated in (5.8)) to find a system of linear equations for the b-coefficients. From (5.11) and (5.12), we can solve for the na-coefficients to find
ai=13hi(bi+1−bi),i=0 to n−1
From (5.9), we can solve for the n c-coefficients as follows:
ci=1hi(fi−aih3i−bih2i)=1hi(fi−13hi(bi+1−bi)h3i−bih2i)=fihi−13hi(bi+1+2bi),i=0 to n−1.
We can now find an equation for the b-coefficients by substituting (5.8), (5.13) and (5.14) into (5.10) :
3(13hi(bi+1−bi))h2i+2bihi+(fihi−13hi(bi+1+2bi))=(fi+1hi+1−13hi+1(bi+2+2bi+1))
which simplifies to
13hibi+23(hi+hi+1)bi+1+13hi+1bi+2=fi+1hi+1−fihi
an equation that is valid for i=0 to n−2. Therefore, (5.15) represent n−1 equations for the n+1 unknown b-coefficients. Accordingly, we write the matrix equa- tion for the b-coefficients, leaving the first and last row absent, as
(⋯⋯⋯⋯ missing ⋯⋯13h023(h0+h1)13h1⋯000⋮⋮⋮⋱⋮⋮⋮000⋯13hn−223(hn−2+hn−1)13hn−1⋯⋯⋯⋯ missing ⋯⋯)(b0b1⋮bn−1bn)(missingf1h1−f0h0⋮fn−1hn−1−fn−2hn−2missing)
Once the missing first and last equations are specified, the matrix equation for the b-coefficients can be solved by Gaussian elimination. And once the b-coefficients are determined, the a - and c-coefficients can also be determined from (5.13) and (5.14), with the d-coefficients already known from (5.8). The piecewise cubic polynomials, then, are known and g(x) can be used for interpolation to any value x satisfying x0≤x≤xn
The missing first and last equations can be specified in several ways, and here we show the two ways that are allowed by the MATLAB function spline.m. The first way should be used when the derivative g′(x) is known at the endpoints x0 and xn; that is, suppose we know the values of α and β such that
g′0(x0)=α,g′n−1(xn)=β.
From the known value of α, and using (5.6) and (5.14), we have
α=c0=f0h0−13h0(b1+2b0)
Therefore, the missing first equation is determined to be
23h0b0+13h0b1=f0h0−α.
From the known value of β, and using (5.6),(5.13), and (5.14), we have
β=3an−1h2n−1+2bn−1hn−1+cn−1=3(13hn−1(bn−bn−1))h2n−1+2bn−1hn−1+(fn−1hn−1−13hn−1(bn+2bn−1)),
which simplifies to
13hn−1bn−1+23hn−1bn=β−fn−1hn−1
to be used as the missing last equation.
The second way of specifying the missing first and last equations is called the not- a-knot condition, which assumes that
g0(x)=g1(x),gn−2(x)=gn−1(x)
Considering the first of these equations, from (5.5) we have
a0(x−x0)3+b0(x−x0)2+c0(x−x0)+d0=a1(x−x1)3+b1(x−x1)2+c1(x−x1)+d1
Now two cubic polynomials can be proven to be identical if at some value of x, the polynomials and their first three derivatives are identical. Our conditions of continuity at x=x1 already require that at this value of x these two polynomials and their first two derivatives are identical. The polynomials themselves will be identical, then, if their third derivatives are also identical at x=x1, or if
a0=a1
From (5.13), we have
13h0(b1−b0)=13h1(b2−b1)
or after simplification
h1b0−(h0+h1)b1+h0b2=0
which provides us our missing first equation. A similar argument at x=xn−1 also provides us with our last equation,
hn−1bn−2−(hn−2+hn−1)bn−1+hn−2bn=0.
The MATLAB subroutines spline.m and ppval.m can be used for cubic spline interpolation (see also interp1.m). I will illustrate these routines in class and post sample code on the course web site.