Back to November Calendar            ←Previous Entry                             Next Entry→

November 13, 2005

Spline Representations V

=========================  Based on http://www.mathworks.com/moler/interp.pdf

Moler on Splines and Bezier

Many of the most effective interpolation techniques are based on piecewise cubic polynomials. Let h_k denote the length of the k^th subinterval:

h_k = x_k+1  x_k

Then the first divided difference, μ_k, is given by

Let m_k denote the slope of the line tangent to the interpolant at x_k:

m_k = P′(x_k)

For the piecewise linear interpolant m_k = μ_k or μ_k-1, but this is not necessarily true for higher order interpolants. 

Consider the following function (pulled magically from a jaws of hat-hidden rabbit) on the interval x_k ≤  x  ≤  x_k+1, expressed in terms of local variables 

s = x  x_k and h = h_k:

This is a cubic polynomial in s, and hence in x, that satisfies four interpolation conditions, two on function values and two on the possibly unknown derivative values:

Functions that satisfy interpolation conditions on derivatives are known as Hermite or osculatory interpolants, because of the higher order contact at the interpolation sites. (“Osculari” means “to kiss” in Latin.)

If we happen to know both function values and first derivative values at a set of data points, then piecewise cubic Hermite interpolation can reproduce those data. But if we are not given the derivative values, we need to define the slopes d_k somehow. Of the many possible ways to do this, we will describe two, which Matlab calls pchip and spline. 

SPLINE

Our other piecewise cubic interpolating function is a cubic spline. The term “spline" refers to an instrument used in drafting; a thin, flexible wooden or plastic tool that is passed through given data points and defines a smooth curve in between.  The physical spline minimizes potential energy subject to the interpolation constraints.  The corresponding mathematical spline must have a continuous second derivative and satisfy the same interpolation constraints. The breakpoints of a spline are also referred to as its knots.

The world of splines extends far beyond the basic one-dimensional, cubic, interpolatory spline we are describing here. There are multidimensional, high-order, variable knot, and approximating splines. A valuable expository and reference text for both the mathematics and the software is A Practical Guide to Splines by Carl de Boor [1]. De Boor is also the author of the spline function and the Spline Toolbox for Matlab.

The first derivative P′(x) of our piecewise cubic function is defined by different formulas on either side of a knot x_k.  Both formulas yield the same value d_k at the knots, so P′ (x) is continuous.

2^nd Order Continuity Too

On the kth subinterval, the second derivative is a linear function of s = x  x_k and (show this, it’s easy, but worth it):

If x = x_k, then s = 0 and so, approaching the knot from the right:

And approaching x_k+1 from the left:

The formula for the second derivative, approaching x from the left can then be obtained by a shift of indices:

Requiring P′′(x) to be continuous at x = x_k means we can write the following equivalent equations:

Aha!  That’s the derivation I’ve been looking for!  The parenthetical factor on the RHS is

Summing up, we pulled a crazy function out thin air, showed it satisfied zero and first order continuity conditions and then, imposing the second order condition, showed that the slopes m_k for 1 ≤ k ≤ n1 of the tangent lines at the control points must satisfy the equation boxed above.  If we add end conditions for the slopes, we’ve verified the matrix equation given without justification in the CMJ article.  As some of the examples on this page illustrate, if you don’t specify either end condition, then the system is underdetermined, but if you specify both, it may be over determined.

The slopes m_k of a spline are closely related to the differences μ_k, they are a kind of running average of the μ_k 's.  So we’ve got n2 equations involving the n unknowns, d_k.  There are a variety of approaches used at the ends of the interval.  One strategy, known

as “not-a-knot,”  uses a single cubic on the first two subintervals, x₁ ≤ x ≤ x₃, and on the last two subintervals, x_n−2 ≤ x ≤ x_n. In effect, x₂ and x_n−1 are not knots.

If the knots are equally spaced with, say, all h_k = 1, the interior formula becomes

While at the ends (using the not knot scheme) the equations are

and

Of course, If there are only n = 3 or 4 points, then a single cubic can be made to satisfy all necessary conditions.  These methods make best sense when n is at least 5, more like 500.  

In the general case, the system can be written in matrix form this way:

The two values r₁ and r_n are associated with the end conditions.

If the knots are equally spaced with all h_k = 1, the coefficient matrix is quite simple:

In the Moler textbook function, splinetx, the linear equations defining the slopes are solved with the tridisolve function introduced in chapter 2 of Moler, Linear Equations.  In the spline functions distributed with Matlab and the Spline Toolbox, the slopes are computed by the MatLab backslash operator d = A\r. 

In the second screen shot, you can see the result of doing the rref (reduced row-echelon form) function on the augmented matrix, thus producing the slopes we want as the right hand side.

At this stage, we’ve got the slopes to fit the not knot end conditions, but these are only a means to finding the coefficients we need to determine the cubic splines.

Consider the general case:

Solving the 2x2 system that arises from the last two equations yields  and   

Thus, we can write a TI-89 program like this

:coeffs()

:Prgm

:{0,2,1,3,1}→d 

:1→h

:rref(tri)→rreftri

:For i,1,5,1

:  rreftri[i,6]→c[i]

:EndFor

:For i,1,4,1

:  (c[i]+c[i+1])/h^2-2*(d[i+1]-d[i])/h^3→a[i]

:  3*(d[i+1]-d[i])/h^2-(2*c[i]+c[i+1])/h→b[i]

:EndFor

:Graph seq(a[i]*(x-(i-1)*h)^3+b[i]*(x-(i-1)*h)^2+c[i]*(x-(i-1)*h)+d[i],i,1,4)

:EndPrgm

This produces the coefficients below and the graph shown to the right.

Are there 4 different cubics?  From the table, you might think so, but, a careful examination of the graph shows that there are only 2 different cubics drawn. So there’s not knots at (1,2) at (3,3).  In fact the first two and the last two coefficient arrays represent the same cubics (here I’ve written the coefficients in exact rational form):

Hmm. Below is a larger graph with the spline and knots shown in black and the cubics drawn in red and blue dash/dots.