Back to November Calendar            slog-11-18-05.htm                             slog-11-20-05.htm

November 19, 2005

Spline Representations XI

A Spline with Irregular Spacing.

Let’s compare the natural cubic spline with this parametric curve.  First, we seek the spline fitting A, B, C and D.  Recall that with the natural spline, the boundary conditions are that the second derivatives are zero there. 

Finding the spline is significantly complicated by the fact that the x-values are not evenly spaced in this example.  I resort to crude matrix methods…no lovely tridiagonalization… 

We want to fit a cubic to each of the three segments determined by consecutive points
A(6, 7), B(0,0), C(12, 9) and D(14,2).   Each cubic has 4 degrees of freedom (four parameters are needed to define the cubic) so there are 3*4 = 12 degrees of freedom to play with.  So let’s count their usages.  The 2 position conditions for each of the 3 cubics account for 2*3 = 6 degrees of freedom.  Then we have both slope and curvature (1^st and 2^nd derivative) conditions at each of the two knots for 2*2 = 4 more degrees of freedom.  This leaves 2 degrees for the end conditions, y″ = 0.

Let the three cubics be written in this form:

The position conditions show that, for i = 1,2,3

and also

The slope conditions are

and the curvature conditions are

Finally, the end conditions (for a natural spline) are

Clearly, b₁ = 0.  To find the remaining 8 unknown parameters, we can solve the following system of equation (in matrix form):

In MatLab, I enter the augmented coefficient matrix like so:

>> A=[0 0 0 0 0 48 2 0 0;   216 6 0 0 0 0 0 0 7;  0 0 1728 144 12 0 0 0 -9;
    108 1 0 0 -1 0 0 0 0; 0 0 432 24 1 0 0 -1 0; 36 0 0 -2 0 0 0 0 0;
     0 0 72 2 0 0 -2 0 0; 0 0 0 0 0 8 4 2 11] 

A =

     0     0     0     0     0    48     2     0     0

   216     6     0     0     0     0     0     0     7

     0     0  1728   144    12     0     0     0    -9

   108     1     0     0    -1     0     0     0     0

     0     0   432    24     1     0     0    -1     0

    36     0     0    -2     0     0     0     0     0

     0     0    72     2     0     0    -2     0     0

     0     0     0     0     0     8     4     2    11

The system is then solved by

>> B=rref(A);

yielding the solutions to system in the 9^th column:

>> B(:,9)

  -0.02388059701493

   2.02641509433962

   0.03445378151261

  -0.42987804878049

  -0.55284552845528

  -0.03378378378378

   0.81045751633987

   4.01415094339623

At this stage, we have the coefficients for functions

The trouble is, to plot these using MatLab’s polyval function, you need to expand and collect to get the standard coefficients.  This is perhaps clumsy, but my first goal is to just get the job done somehow, so here’s what I do.  First, to simplify further commands, I copy the 9^th column of B to a single vector I call c9:

>> c9 = B(:,9);

Then, since I need to expand two binomials to simplify f₁ and f3 I first create the binomials:

>> binomial_1=[1 6];

>> binomial_3=[1 -12];

and then set about squaring and cubing these like so:

>> bi_square_1 = conv(binomial_1, binomial_1);

>> bi_square_3 = conv(binomial_3, binomial_3);

>> bi_cube_1   = conv(binomial_1, bi_square_1);

>> bi_cube_3   = conv(binomial_3, bi_square_3);

Then I pad the various vectors with zeros so they’re all compatible:

>> binomial_1 = [0 0 binomial_1]

binomial_1 =

     0     0     1     6

>> bi_square_1 = [0 bi_square_1]

biquare =

     0     1    12    36

And take the linear combination of these to form the standard coefficient array.  Note that the coefficient of the square term is zero, so that one isn’t included, and the constant term needs to be added in.  Really the constant term should be zero, but I keep the round off error, just as an indication of the background noise that is present in all the terms.

>> spline1=c9(1)*bi_cube_1+c9(2)*binomial_1

spline1 =

   -0.0239   -0.4299   -0.5527    7.0003

>> spline1(4)=spline1(4)-7

spline1 =

   -0.0239   -0.4299   -0.5527    0.0003

Similarly, the third piece is computed:

>> spline3 = c9(6)*bi_cube_3+c9(7)*bi_square_3+c9(8)*binomial_3

spline3 =

   -0.0338    2.0267  -30.0314  126.9144

>> spline3(4) = spline3(4)-9

spline3 =

   -0.0338    2.0267  -30.0314  117.9144

The second (middle) piece is somewhat simpler, since no shifting/expanding is involved:

>> spline2 = [c9(3) c9(4) c9(5) 0]

spline2 =

    0.0345   -0.4299   -0.5528         0

Armed with these, I partition the various x-intervals like so:

>> x1=[-6:.1:0];

>> x2=[0:.1:12];

>> x3=[12:.1:14];

and then plot splines and the control points:

>>a=[-6 0 12 14];

>>b = [-7 0 -9 2];

>>plot(x1,polyval(p1,x1),
       x2,polyval(p2,x2),
       x3,polyval(p3,x3),'LineWidth',2)

>>hold on

>>plot(a,b,'or')

This yields the irregularly spaced natural spline plot below: