Back to November Calendar            ←Previous Entry                             Next Entry→

November 10, 2005

Spline Representations II

 

Yesterday we left off with this matrix equation involving the set of coefficients {ci} for the linear terms of the sequence of cubic functions xi(t) whose endpoints pass through  a given pair of points (ti, xi(ti)) and (ti+1, xi(ti+1)).  Together with the natural spline boundary constraints that the curvature (second derivative) at the extreme endpoints is zero, this can be written in matrix form as Ac = d:

   

 

where hk = tk+1-tk and wk is the slope of the segment connecting (ti, xi(ti)) and (ti+1, xi(ti+1)). 

 

An outline of a cubic spline algorithm with natural boundary constraints can be pieced together as follows:

 

  1. For k = 1 to n1 compute {
     
      
  2. Solve the mammoth tri-diagonal system whose matrix form, Ac = d, is shown above.  You could use Cholesky decomposition with an eye towards the fact that the coefficient matrix is nearly positive definite.
  3. For k = 1, 2, … , n-1, compute  
  4. For each k and , compute   where

              

A quick, simple example of implementing this algorithm is to fit three points, say (0,1), (1,4) and (2,3).   Following step 1 of the algorithm, we compute the step sizes h1 = h2 = 1 and then slopes .  This leads to (step 2) solving the tri-diagonal coefficient matrix system:

                                                           

 

The solution is consistent with the natural spline boundary conditions: c1 = c3 = 0, and it’s easy to solve for c2=−6.

 

Having the values for ci, we can proceed to step 3 and compute the coefficients of higher order terms:  , ,   and . 

 

The final step is to display (somehow) the cubic functions which fit together so smoothly.  There formulas are

                                        

and

                             

These two cubics and the points they fit are shown at right. 

 

The continuity conditions are met since,

  1. the control points are all hit,
  2. the natural spline boundary conditions are met :  and  and
  3. at the splicing of the splines,  and . 

 

Looks a bit like a parabola maybe, but, obviously, it’s not.

 

Alternatively, people sometimes start by finding the unique cubics fitting the first four and last four control points and then filling in cubics between to meet first order and second-order continuity conditions between successive control points. 

 

Searching JSTOR leads to an article, Cubic Splines from Simpson’s Rule, by Nishan Krikorian and Marc Ramras in the College Math Journal, Vol. 27, No. 2 (Mar., 1996), 124-126, which claims the more general result.

 

Where mk is the slope of the line tangent at control point k.  The paper points out that this is a rather tedious derivation in most standard texts, but there is a slick way of obtaining it.  Well, I can’t help but do it the dumb way first to help understand the slick way, maybe.

 

If we restrict the problem to the case with evenly spaced control points, , then the above system simplifies to

 

If you’ve been catching a waft of Simpson’s rule throughout this cubic spline business, it may have just gotten quite a bit stronger!  That 1-4-1 coefficient sequence sure looks like it! 

 

To see this, consider the piecewise cubic spline function C(x) on the double interval .  By the fundamental theorem of calculus, .

Now Simpson’s rule to compute the integral yields ,  which is, in fact, the second equation in the system described above.  All the other double intervals are similarly matched. 

 

It is well-known that Simpson’s rule is exact for cubics, so it seems like it would also be exact for the piecewise quadratics such as C’(x).  To be sure, if Q(x) is the quadratic function fitting (x0, s0) , (x1, s1) and (x2, s2), then, since C(x) is assumed to be twice differentiable (smooth), C’(x)  Q(x) is a differential piecewise quadratic interpolating the equally spaced data (x0, 0) , (x1, 0) and (x2, 0).  It’s at this point that I look back, look for what’s what and see, and say, “Oh yeah!  Why didn’t I think of this clever bit?” 

 

In fact, C’(x)  Q(x) is anti-symmetric about x1 (this is a good exercise.)

 

Maybe a quick example is called for.  Suppose we have the following table of points and  slopes to fit:

x

0

3

6

y

0

4

1

m

1

1/2

-1/3

Using the form

 

 

and applying the boundary conditions that

Pk(xk) = yk ,  Pk(xk+1) = yk+1, Pk’(xk) = mk ,  and Pk’ (sk+1) = yk+1

leads to the two systems, one for each spline:

 

 

Solving these systems for the coefficients, yields the two cubic polynomials

 

 

with derivatives

     

The quadratic fitting the points

x

0

3

6

m

1

1/2

-1/3

can be found by substituting x and m into x2a + xb + 1 = m to obtain

 

 

Thus the piecewise function

 

 

So, oh dear, it’s not anti-symmetric!  After re-reading the Math Mag article, I now detect that “the slopes are not given, but are arbitrary!  Ouch.  So…I did pick them out of a hat!