Back to November Calendar            ←Previous Entry                             Next Entry→

November 17, 2005 Spline Representations IX   Simpon’s Rule and Splines.   A week ago I tried to demonstrate the connection between Simpson’s Rule and Cubic Splines by following the CMJ article, Cubic Splines from Simpson’s Rule, by Nishan Krikorian and Marc Ramras, Vol. 27, No. 2 (Mar., 1996), 124-126, first mentioned in these pages 11/10.    To see the spline/Simpson connection, start with a definition:   Defintion:  C(x) is a piecewise cubic spline function satisfying the 2nd order continuity condition 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 of equations for the 2nd order cubic spline.  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.)   I found that I had an error in my formula for the coefficient bk of the second order term in the splines coefficients.  Also the example failed...so I’m going to try a new example, this time with 7 control points, 5 knots:     x y slope RHS -6 0 -1 -1.5 -4 2 0 0 -2 0 1 1.5 0 3 0 0 2 0 -1 -1.5 4 2 0 0 6 0 1 1.5

 

In MatLab, I implement the complete splines solution as follows:

 

>> A = [1 0 0 0 0 0 0 ; 1 4 1 0 0 0 0 ; 0 1 4 1 0 0 0 ; 0 0 1 4 1 0 0 ; 0 0 0 1 4 1 0 ; 0 0 0 0 1 4 1 ; 0 0 0 0 0 0 1]

A =

     1     0     0     0     0     0     0

     1     4     1     0     0     0     0

     0     1     4     1     0     0     0

     0     0     1     4     1     0     0

     0     0     0     1     4     1     0

     0     0     0     0     1     4     1

     0     0     0     0     0     0     1

>> y = [0 2 0 3 0 2 0];

>> RHS = 1.5.*[-1 0 1 0 -1 0 1]

>> RHS = RHS';

>> m=A\RHS;

>> xk

xk =

  Columns 1 through 7

   -6.0000   -5.9000   -5.8000   -5.7000   -5.6000   -5.5000   -5.4000

   -4.0000   -3.9000   -3.8000   -3.7000   -3.6000   -3.5000   -3.4000

   -2.0000   -1.9000   -1.8000   -1.7000   -1.6000   -1.5000   -1.4000

         0    0.1000    0.2000    0.3000    0.4000    0.5000    0.6000

    2.0000    2.1000    2.2000    2.3000    2.4000    2.5000    2.6000

    4.0000    4.1000    4.2000    4.3000    4.4000    4.5000    4.6000

  Columns 8 through 14

   -5.3000   -5.2000   -5.1000   -5.0000   -4.9000   -4.8000   -4.7000

   -3.3000   -3.2000   -3.1000   -3.0000   -2.9000   -2.8000   -2.7000

   -1.3000   -1.2000   -1.1000   -1.0000   -0.9000   -0.8000   -0.7000

    0.7000    0.8000    0.9000    1.0000    1.1000    1.2000    1.3000

    2.7000    2.8000    2.9000    3.0000    3.1000    3.2000    3.3000

    4.7000    4.8000    4.9000    5.0000    5.1000    5.2000    5.3000

  Columns 15 through 20

   -4.6000   -4.5000   -4.4000   -4.3000   -4.2000   -4.1000

   -2.6000   -2.5000   -2.4000   -2.3000   -2.2000   -2.1000

   -0.6000   -0.5000   -0.4000   -0.3000   -0.2000   -0.1000

    1.4000    1.5000    1.6000    1.7000    1.8000    1.9000

    3.4000    3.5000    3.6000    3.7000    3.8000    3.9000

    5.4000    5.5000    5.6000    5.7000    5.8000    5.9000

>> hold on

>> for i=1:6 plot(xk(i,:),splines04(xi(i,:),ak,bk,ck,dk,i,h)); end;

>> for k = 1:6 ak(k) = -2*(y(k+1)-y(k))/h^3+(m(k)+m(k+1))/h^2; end;

>> for k = 1:6 bk(k) =  3*(y(k+1)-y(k))/h^2-(2*m(k)+m(k+1))/h; end;

>> ck = m;

>> dk = y;

>> hold on

>> for i=1:6 plot(xk(i,:),splines04(xk(i,:),ak,bk,ck,dk,i,h,-6)); end;

>> xp=[-6:2:6]

>> hold on

>> plot(xp,dk,'o')

 

This MatLab code produces the graph shown below:

The coefficients of the cubics   are listed here:

 

a_k = -0.8000    0.6500   -0.6750    0.6750   -0.6500    0.8000

b_k =  2.8500   -1.9500    1.9500   -2.1000    1.9500   -1.9500

c_k = -1.5000    0.3000    0.3000    0.0000   -0.3000   -0.3000

d_k =  0         2         0         3         0         2

 

So, to check the anti-symmetry claim,  We need to find the quadratics fitting the (x_k, m_k) in groups of 3.  The quadratic fitting {( 6,-1.5),( 4,0.3),( 2,0.3)} and {(2,0.3),(0,0),(2, 0.3)} are

and {(2, 0.3),(4, 0.3),(6, 1.5)}, which should be a simple translation of the first, by symmetry.

 

I used the TI89 to solve the find the first parabola’s equation, as outlined in the following three screenshots showing the general form of the Vandermonde coefficient matrix, it’s augmented form after entering values for the parameters x₁, x₂ and x₃. and the result of solving the system by obtaining the reduced row-echelon form (rref) of the augmented matrix:

 

Thus  and on first double interval then we have

To visualize what’s going on here, I draw for you pictures.  The graph on the left includes three parabolas: the derivative functions for the cubic splines from -6 to -4 and from -4 to -2 and the parabola fitting the point they have in common and the two points at the extremes.  The red part of the second graph is the difference between the blue and amber parts of the first graph: C’(x)  Q(x):

To nail this down algebraically, it might help to write the two pieces of C’(x)  Q(x) in factored and vertex forms:  

So just by anti-symmetry of this function, it’s clear that there’s as much area below the graph above [-6,-4] as there is above the graph below [-4,-2] so  

 

This was demonstrated in this particular instance, but it is true in general, and since Simpson’s rule is exact for quadratics, it is also exact for piecewise quadratics and thus