Docstoc

Cubic Splines

Document Sample
Cubic Splines Powered By Docstoc
					Cubic Splines                      MACM 316                                        1/14



                            Cubic Splines
Given the following list of points:
          x :      a = x0 < x1 < · · · < xn = b

          y :          y0    y1     ···       yn

  ® Basic idea: “piecewise polynomial interpolation”
           • Use lower order polynomials that interpolate on each subin-
             terval [xi, xi+1].
           • Force the polynomials to join up as smoothly as possible.
  ® Simplest example: a linear spline just “connects the dots”

Definition: A cubic spline S(x) is a piecewise-defined function
that satisfies the following conditions:
  1. S(x) = Si (x) is a cubic polynomial on each subinterval [xi, xi+1]
     for i = 0, 1, . . . , n − 1.
  2. S(xi) = yi for i = 0, 1, . . . , n (S interpolates all the points).
  3. S(x), S ′(x), and S ′′ (x) are continuous on [a, b] (S is smooth).


So we write the n cubic polynomial pieces as
        Si (x) = ai + bi (x − xi ) + ci(x − xi)2 + di (x − xi)3,
for i = 0, 1, . . . , n − 1, where ai, bi , ci, and di represent 4n un-
known coefficients.



October 16, 2008                                           c Steven Rauch and John Stockie
Cubic Splines                              MACM 316                                           2/14



                                     Derivation
Using the definition, let’s determine equations relating the coef-
ficients and keep a count of the number of equations:
                                                                                # eqns.
   Interpolation and continuity:
     Si (xi) = yi           for i = 0, 1, . . . , n − 1                                n
       Si (xi+1) = yi+1                 for i = 0, 1, . . . , n − 1                    n
       (based on both of these, S(x) is continuous)

   Derivative continuity:
     ′           ′
    Si (xi+1) = Si+1 (xi+1) for i = 0, 1, . . . , n − 2                           n−1
        ′′          ′′
       Si (xi+1) = Si+1 (xi+1) for i = 0, 1, . . . , n − 2                        n−1

                                        Total # of equations:                   4n − 2

There are still 2 equations missing! (later)




October 16, 2008                                                      c Steven Rauch and John Stockie
Cubic Splines                           MACM 316                                            3/14



                              Derivation (cont’d)
We need expressions for the derivatives of Si :
          Si(x) = ai + bi(x − xi) + ci(x − xi)2 + di (x − xi)3

          Si (x) = bi + 2ci(x − xi) + 3di(x − xi)2
           ′


           ′′
          Si (x) = 2ci + 6di(x − xi)
It is very helpful to introduce the hi = xi+1 − xi. Then the spline
conditions can be written as follows:
 • Si(xi) = yi for i = 0, 1, . . . , n − 1:
                   ai = y i

 • Si(xi+1) = yi+1 for i = 0, 1, . . . , n − 1:

                   ai + hi bi + h2ci + h3di = yi+1
                                 i      i

    ′           ′
 • Si (xi+1) = Si+1 (xi+1) for i = 0, 1, . . . , n − 2:

                   bi + 2hici + 3h2di − bi+1 = 0
                                  i

    ′′          ′′
 • Si (xi+1) = Si+1 (xi+1) for i = 0, 1, . . . , n − 2:

                   2ci + 6hidi − 2ci+1 = 0

The boxed equations above can be written as a large linear sys-
tem for the 4n unknowns
        [a0, b0, c0, d0, a1, b1, c1, d1 , . . . , an−1, bn−1, cn−1, dn−1 ]T




October 16, 2008                                                    c Steven Rauch and John Stockie
Cubic Splines                                     MACM 316                                              4/14



                             Alternate formulation
This linear system can be simplified considerably by defining
                            ′′                                             mi
                      mi = Si (xi) = 2ci                  or        ci =   2

and thinking of the mi as unknowns instead. Then:

    ′′          ′′
 • Si (xi+1) = Si+1 (xi+1) for i = 0, 1, . . . , n − 2:
     =⇒            2ci + 6hidi − 2ci+1 = 0
     =⇒            mi + 6hidi − mi+1 = 0

                           mi+1 −mi
     =⇒             di =     6hi


 • Si(xi) = yi and Si (xi+1) = yi+1 for i = 0, 1, . . . , n − 1:
     =⇒            yi + hibi + h2ci + h3di = yi+1
                                i      i

     Substitute ci and di from above:

                           yi+1 −yi       hi          hi
     =⇒             bi =      hi
                                      −   2
                                             mi   −   6
                                                         (mi+1   − mi)

    ′           ′
 • Si (xi+1) = Si+1 (xi+1) for i = 0, 1, . . . , n − 2:
     =⇒            bi + 2hici + 3h2di = bi+1
                                  i

     Substitute bi, ci and di from above and simplify:

                                                                         yi+2 −yi+1            yi+1 −yi
         himi + 2(hi + hi+1 )mi+1 + hi+1 mi+2 = 6                           hi+1
                                                                                         −        hi


Notice: These are n − 1 linear equations for n + 1 unknowns
        [m0, m1, m2, . . . , mn]T                  where               ′′
                                                                 mi = gi (xi).
October 16, 2008                                                                c Steven Rauch and John Stockie
Cubic Splines                            MACM 316                                            5/14



              Endpoint Conditions: Natural Spline
  ® Now, we’ll deal with the two missing equations . . .
  ® Note: Derivative matching conditions are applied only at in-
    terior points, which suggests applying some constraints on
    the derivative(s) at x0 and xn.
  ® There is no unique choice, but several common ones are
          1. Natural Spline (zero curvature):
                      m0 = 0     and       mn = 0
                   for taken together gives an (n + 1) × (n + 1) system:
                                                                             
             1       0           0                        ···            0                        
            h0 2(h0 + h1 )     h1           0                                 m0
                                                                              m                 
            0      h1      2(h1 + h2 )      h2            0                   1                 
                                                                         .   m2
                                                                                                
        
             0       0          h2      2(h2 + h3 )       h3             . 
                                                                         .
                                                                                                   
                                                                                m3
                                                                                                  
             .
              .                         ...          ...         ...           .
                                                                                                   
             .                                                                .
                                                                                                   
                                                                                   .              
                                0          hn−2    2(hn−2 + hn−1 ) hn−1 
                                                                                   mn
             0     ···                       0             0             1
                                                                                   
                                                                       0
                                                                y2 −y1
                                                          
                                                                 h1
                                                                       − y1h0 0
                                                                             −y     
                                                                                    
                                                                y3 −y2    y2 −y1
                                                                       − h1
                                                                                   
                                                                 h2                
                                                               y4 −y3    y3 −y2    
                                                       = 6       h3
                                                                       − h2         
                                                                       .
                                                                                   
                                                          
                                                                      .
                                                                       .            
                                                                                    
                                                           yn −yn−1     yn−1 −yn−2 
                                                           hn−1 − hn−2             
                                                                       0


. . . or you could just drop the equations for m0 and mn and write
it as an (n − 1) × (n − 1) system.


October 16, 2008                                                     c Steven Rauch and John Stockie
Cubic Splines                    MACM 316                                      6/14



                   Interpretation of a Natural Spline




    • The word “spline” comes from the draftsman’s spline, which
      is shown here constrained by 3 pegs.
    • There is no force on either end to bend the spline beyond the
      last pegs, which results in a flat shape (zero curvature).
    • This is why the S ′′ = 0 end conditions are called natural.




October 16, 2008                                       c Steven Rauch and John Stockie
Cubic Splines                            MACM 316                                            7/14



           Endpoint Conditions: Clamped Spline
          2. Clamped Spline, sometimes called a “complete spline” (deriva-
             tive is specified):
                    ′
                   S0 (x0) = A     =⇒    b0 = A
                                              y1 − y0       h0          h0
                                   =⇒    A=             −        m0 −        (m1 − m0)
                                                  h0        2           6
                                                                   y1 − y0
                                   =⇒    2h0m0 + h0m1 = 6                  −A
                                                                      h0
                   and
                    ′
                   Sn−1 (xn) = B    =⇒     bn−1 = B
                                                                    yn − yn−1
                      =⇒     hn−1mn−1 + 2hn−1mn = 6 B −
                                                                      hn−1

                   These two equations are placed in the first/last rows
                                                                                          
                        2h0    h0          0      ···       ···         0
                                                                        .
                                                                        .
                       h0 2(h0 + h1)     h1       0                    .
                                                                                          
                                                                                           
                      
                       0                                               .
                                                                        .
                                                                                           
                              h1     2(h1 + h2) h2          0          .                  
                                                                                           
                       .                 ...     ...       ...
                       ..      0                                       0
                                                                                           
                                                                                           
                       0                  0     hn−2 2(hn−2 + hn−1 ) hn−1
                                                                                          
                              ···                                                          
                         0    ···        ···       0       hn−1       2hn−1

                   and the first/last RHS entries are modified accordingly.




October 16, 2008                                                     c Steven Rauch and John Stockie
Cubic Splines                               MACM 316                                           8/14



                    Endpoint Conditions: Not-A-Knot
          3. Not-A-Knot Spline (third derivative matching):
                       ′′′       ′′′                      ′′′          ′′′
                      S0 (x1) = S1 (x1)         and      Sn−2(xn−1) = Sn−1 (xn−1)

                          ′′′                          mi+1 −mi
                   Using Si (x) = 6di and di =            6
                                                                ,   these conditions
                   become

                      m1 − m0 = m2 − m1

                   and

                      mn−1 − mn−2 = mn − mn−1
                   The matrix in this case is
                                                                                         
                          −1      2         −1        ···      ···         0
                                                                            .
                                                                            .
                          h0 2(h0 + h1 )     h1        0                    .
                                                                                         
                                                                                         
                                                                           .
                                                                            .
                                                                                          
                      
                         0      h1      2(h1 + h2 ) h2          0          .             
                           .
                                                                                          
                      
                          .
                           .      0          ...      ...       ...        0
                                                                                          
                                                                                          
                          0                   0      hn−2 2(hn−2 + hn−1 ) hn−1
                                                                                         
                               ···                                                       
                          0     ···         ···       −1         2         −1

                   and the first/last RHS entries are zero.

                   Note: Matlab’s spline function implements the not-a-knot
                   condition (by default) as well as the clamped spline, but
                   not the natural spline. Why not? (see Homework #4).




October 16, 2008                                                       c Steven Rauch and John Stockie
Cubic Splines                                   MACM 316                                          9/14



                          Summary of the Algorithm
     Starting with a set of n + 1 data points
        (x0, y0), (x1, y1), (x2, y2), . . . , (xn, yn)
  1. Calculate the values hi = xi+1 − xi for i = 0, 1, 2, . . . , n − 1.
  2. Set up the matrix A and right hand side vector r as on page 6
     for the natural spline. If clamped or not-a-knot conditions are
     used, then modify A and r appropriately.
  3. Solve the (n + 1) × (n + 1) linear system Am = r for the
     second derivative values mi.
  4. Calculate the spline coefficients for i = 0, 1, 2, . . . , n − 1 us-
     ing:
                   ai = y i
                          yi+1 − yi       hi          hi
                   bi =               −        mi −        (mi+1 − mi )
                               hi         2           6
                          mi
                   ci =
                        2
                        mi+1 − mi
                   di =
                           6hi
  5. On each subinterval xi ≤ x ≤ xi+1, construct the function
                   gi (x) = ai + bi(x − xi) + ci(x − xi)2 + di (x − xi)3




October 16, 2008                                                          c Steven Rauch and John Stockie
Cubic Splines                                MACM 316                                           10/14



                             Natural Spline Example
                             x 4.00 4.35 4.57 4.76 5.26 5.88
                             y 4.19 5.77 6.57 6.23 4.90 4.77
                             h 0.35 0.22 0.19 0.50 0.62
Solve Am = r when
                         1    0        0    0    0
                                                       
                                                                                           
                                                                                0
                      0.35 1.14 0.22 0     0    0      
                                                                       −5.2674           
                        0 0.22 0.82 0.19 0      0                                        
        A=                                                and   r=    −32.5548           
                        0    0 0.19 1.38 0.50 0                                          
                                                                       −5.2230           
                        0    0    0 0.50 2.24 0.62     
                                                                         14.70180
                         0    0    0        0    1

        =⇒             m = [0, 3.1762, −40.4021, −0.6531, 6.7092, 0]T


Calculate the spline coefficients:
                 xi y i = a i               bi     ci         di      Interval
          S0(x) 4.00 4.19              4.3290           0    1.5125 [4.00, 4.35]

          S1(x) 4.35          5.77     4.8848     1.5881 −33.0139 [4.35, 4.57]
          S2(x) 4.57          6.57     0.7900 −20.2010      34.8675 [4.57, 4.76]

          S3(x) 4.76          6.23   −3.1102     −0.3266     2.4541 [4.76, 5.26]
          S4(x) 5.26          4.90   −1.5962      3.3546    −1.8035 [5.26, 5.88]

Then the spline functions are easy to read off, for example:
        S2 (x) = 6.57 + 0.7900(x − 4.57) − 20.2010(x − 4.57)2 + 34.8675(x − 4.57)3

        S3 (x) = 6.23 − 3.1102(x − 4.76) − 0.3266(x − 4.76)2 + 2.4541(x − 4.76)3


Exercise: Verify that S2 and S3 satisfy the conditions defining a
cubic spline.
October 16, 2008                                                         c Steven Rauch and John Stockie
Cubic Splines                                  MACM 316                                             11/14



                             Splines vs. Interpolation
                        7

                       6.5

                        6

                       5.5
                   y




                        5

                       4.5
                                                                    Polynomial
                                                                    Natural Spline
                        4
                         4            4.5                 5           5.5
                                                      x

     The Newton divided difference table for the data points is
                    x    y       a1           a2              a3     a4         a5
                   4.00 4.19
                               4.5143
                   4.35 5.77                -1.5402
                               3.6364                 -15.3862
                   4.57 6.57             -13.2337                  22.6527
                               -1.7895                13.1562                -15.7077
                   4.76 6.23                -1.2616                -6.8778
                               -2.6600                2.6331
                   5.26 4.90                2.1878
                               -0.2097
                   5.88 4.77

=⇒           P5 (x) = 4.19 + 4.5143(x − 4) − 1.5402(x − 4)(x − 4.35) − . . .

October 16, 2008                                                              c Steven Rauch and John Stockie
Cubic Splines                             MACM 316                                            12/14



                           Clamped Spline Example
                   x 4.00 4.35 4.57 4.76 5.26 5.88
                   y 4.19 5.77 6.57 6.23 4.90 4.77                 (same data)
                   h 0.35 0.22 0.19 0.50 0.62

Assume S ′ (4.00) = −1.0 and S ′(5.88) = −2.0 (clamped).

Solve Am = r where
                       0.70 0.35 0     0    0    0                      33.0858
                                                                                       
          
                      0.35 1.14 0.22 0     0    0    
                                                      
                                                                  
                                                                       −5.2674           
                                                                                          
                        0 0.22 0.82 0.19 0      0                   −32.5548           
        A=                                              and   r=                       
          
                        0    0 0.19 1.38 0.50 0      
                                                      
                                                                  
                                                                       −5.2230           
                                                                                          
                        0    0    0 0.50 2.24 0.62                   14.7018           
                         0    0    0   0 0.62 1.24                     −10.7418

=⇒           m = [54.5664, −14.6022, −35.0875, −3.0043, 11.1788, −14.2522]T


                 xi y i = a i bi                ci       di      Interval
          S0(x) 4.00 4.19 −1.0000             27.2832 −32.9375 [4.00, 4.35]
          S1(x) 4.35         5.77    5.9937   −7.3011 −15.5191 [4.35, 4.57]

          S2(x) 4.57         6.57    0.5279 −17.5437      28.1431 [4.57, 4.76]
          S3(x) 4.76         6.23   −3.0908   −1.5021      4.7277 [4.76, 5.26]

          S4(x) 5.26         4.90   −1.0472    5.5894     −6.8363 [5.26, 5.88]


Exercise: Verify that the spline satisfies the clamped end condi-
tions:
         ′                    ′
        S0 (4.00) = −1.0 and S4(5.88) = −2.0.



October 16, 2008                                                       c Steven Rauch and John Stockie
Cubic Splines                      MACM 316                                  13/14



                             Comparison
Below is a comparison of the splines with various endpoint con-
ditions:
                        7
                                                Natural
                                                Clamped
                       6.5                      Not−A−Knot

                        6

                       5.5
                   y




                        5

                       4.5

                        4
                         4   4.5            5    5.5
                                        x




October 16, 2008                                       c Steven Rauch and John Stockie
Cubic Splines                           MACM 316                                         14/14



                                  Second Example
Consider the following data representing                   the thrust of a model
rocket versus time.
             t    T                      t                   T
           0.00 0.00                   0.60                16.00
           0.05 1.00
                                       0.70                16.00
           0.10 5.00
           0.15 15.00                  0.80                16.00
           0.20 33.50                  0.85                16.00
           0.30 33.00                  0.90                 6.00
           0.40 16.50                  0.95                 2.00
           0.50 16.00
                                       1.00                 0.00
Below is a plot of the natural spline:
                            40

                            35

                            30

                            25
                   Thrust




                            20

                            15

                            10

                             5

                             0
                              0   0.2   0.4          0.6     0.8         1
                                              time

Notice how wiggles are introduced near the ends of the flat re-
gion – non-smooth data cause some problems for cubic splines.

October 16, 2008                                                   c Steven Rauch and John Stockie

				
DOCUMENT INFO