Curve fitting 2.

S. Arlinghaus

5.  Bounded curve fitting--using polynomial segments to link a finite set
of points.

     The following strategy enables one to link a finite set of points
with pieces
of a cubic curve (polynomial of degree 3).  Different cubics are spliced
together to
form a single smooth curve composed of segments of different cubics.  The
splice is
at the data point; at the splice, the slope of the tangent line to the
cubic on the left
should be the same as the slope of the tangent line to the cubic on the
right.  The smooth
 curve is often called an interpolating curve; estimates of unknown
y-values may be made
between known x-values; unlike the curve1.wk1 there is no logical way to
extrapolate, however.
     The technique for making all the equations work out rests on
elementary calculus,
to ensure continuous curvature (that the second derivatives be
continuous), and on linear 
algebra.  [Interested students might look at the handout in 2044.]
The following example shows how to perform the interpolation in black-box
mode.
The idea is to find five cubics which, when linked will form a single
smooth curve
linking the six points. The reason to use cubics is that from linear beam
theory, "the
fourth derivative of the displacement of a beam is zero along any interval
of the x-axis
that contains no external forces acting on the beam" (Applications of
Linear Algebra, p. 677).
  The idea is to find, simultaneously, values for a, b, c, d, in each
interval that will 
create a cubic piece which when linked to the next piece will be smooth;
generally, each
cubic segment is of the form ax^3 + bx^2 + cx + d.

CUBIC SPLINE INTERPOLATION
EXAMPLE--CHOOSE SIX POINTS
(x1, y1)=(1,1.25)
(x2, y2)=(2,1.75)
(x3, y3)=(3,3)
(x4, y4)=(4,2.5)
(x5, y5)=(5,2)
(x6, y6)=(6,1.75)       
        
Basic Theorem--see notes in 2044 for derivation--from textbook  
Given n points (x_1,y_1),...,(x_n,y_n) with x_(i+1)-x_i=h, i=1,2,...,n-1,       
the cubic spline        
        
        a_1(x-x_1)^3+b_1(x-x_1)^2+c_1(x-x_1)+d_1, x_1<= x <= x_2
        a_2(x-x_2)^3+b_2(x-x_2)^2+c_2(x-x_2)+d_2, x_2<= x <= x_3
S(x) =  ........................................................

a_(n-1)(x-x_(n-1))^3+b_(n-1)(x-x_(n-1))^2+c_(n-1)(x-x_(n-1))+d_(n-1),
x_(n-1) <= x <= x_n
        
that interpolates these points has coefficients given by        
        
a_i = (M_(i+1)-M_i)/6h  
b_i = (M_i)/2   
c_i = (y_(i+1)-y_i)/h-((M_(i+1)+2M_i)h/6)       
d_i = y_i                                                       
                                                        
for i = 1,2,...,n-1, where M_i=S''(x_i), i=1,2,...,n.                                                   
                                                        
Some substitutions and algebraic manipulations make it possible to                                                      
rewrite this system of equations as a single matrix equation.  When this                                                        
single matrix equation is coupled with the physical assumption that the                                                 
natural spline is allowed to extend freely beyond the interpolated
region--                                                        
that is, when the ends shoot off on a straight line so that S''(x)=0 and                                                        
M_1=M_n=0--the following matrix equation is the result.                                                 
                                                        
1  0  0  0  ...  0  0  0                                M_1
0 
1  4  1  0  ...  0  0  0                                M_2
y_1-2y_2+y_3
0  1  4  1  ...  0  0  0                                M_3
y_2-2y_3+y_4
                                                        
........................                           *    ...             =
............
                                                                                                        
0  0  0  0  ...  1  4  1                                M_(n-1)
y_(n-2)-2y_(n-1)+y_n                                            
0  0  0  0  ...  0  0  1                                M_n
0                                               
                                                                                                        

Use this latter equation to find coefficients for the cubic spline to fit
the sample points.                                                                                                      
Equation from Black Box:                                                                                                        
1       0       0       0       0       0               M1
0                       (x1, y1)=(1,1.25)
1       4       1       0       0       0               M2
0.75                    (x2, y2)=(2,1.75)
0       1       4       1       0       0         times M3      =       6
times   -1.75                   (x3, y3)=(3,3)
0       0       1       4       1       0               M4
0                       (x4, y4)=(4,2.5)
0       0       0       1       4       1               M5
0.25                    (x5, y5)=(5,2)
0       0       0       0       0       1               M6
0                       (x6, y6)=(6,1.75)


Solve the matrix equation for M1...M6:                                                                                                  

M1                      1       0       0       0       0       0
0 
M2                      -0.267943       0.267943        -0.072  0.019139
0       0               0.75 
M3      =       6 times 0.07177         -0.07177        0.2871  -0.076555
0.02    0       times   -1.75 
M4                      -0.019139       0.019139        -0.077  0.287081
-0.1    0.07            0 
M5                      0.004785        -0.00478        0.0191  -0.07177
0.27    -0.3            0.25 
M6                      0       0       0       0       0       1
0 

                        0               0                                       
                        0.325359                1.9522                                  
        =       6 times -0.551435             = -3.309                                  
                        0.130383                0.7823                                  
                        0.029904                0.1794                                  
                        0               0                                       

Matrix used:                                                                            
1       0       0       0       0       0       
1       4       1       0       0       0       
0       1       4       1       0       0       
0       0       1       4       1       0       
0       0       0       1       4       1       
0       0       0       0       0       1       

Inverse used (calculated in spreadsheet):                                               
        1       0       0       0       0       0 
        -0.2679         0.2679426       -0.07177        0.019139
-0.005  0.004785 
        0.0718  -0.07177        0.287081        -0.07656        0.0191
-0.019139 
        -0.0191         0.0191388       -0.076555       0.287081
-0.072  0.07177 
        0.0048  -0.004785       0.019139        -0.07177        0.2679
-0.267943 
        0       0       0       0       0       1 

Find coefficients, a's, b's, c's, and d's as in black box                                               

                a1      a2      a3      a4      a5                                                      
0       M1      0.3253589                                                                                       
1.952153        M2              -0.876794                                                                               
-3.308612       M3                      0.681818                                                                        
0.782297        M4                              -0.1                                                            
0.179426        M5                                      -0.029904                                                       
0       M6                                                                                              

                b1      b2      b3      b4      b5                                                      
0       M1      0                                                                                       
1.952153        M2              0.976077                                                                                
-3.308612       M3                      -1.65431                                                                        
0.782297        M4                              0.3911                                                          
0.179426        M5                                      0.089713                                                        
0       M6                                                                                              

                c1      c2      c3      c4      c5
1.25 
0       M1      0.1746411
1.75 
1.952153        M2              1.150718
3 
-3.308612       M3                      0.472488
2.5 
0.782297        M4                              -0.791
2 
0.179426        M5                                      -0.309809
1.75 
0       M6                                                                                              

                d1      d2      d3      d4      d5                                                      
0       M1      1.25                                                                                    
1.952153        M2              1.75                                                                            
-3.308612       M3                      3                                                                       
0.782297        M4                              2.5                                                             
0.179426        M5                                      2                                                       
0       M6                                                                                              

                a's     b's     c's     d's             
                0.325358        0       0.174641        1.25            
                -0.87679        0.976076        1.150717        1.75            
                0.681818        -1.6543         0.472488        3               
                -0.10047        0.391148        -0.79066        2.5             
                -0.0299         0.089712        -0.3098         2               

Equation of cubic spline to fit the six sample points.                                                  
Five equations to fit curve pieces between the six sample points.                                                       
        0.3254  (x-1)^3 +0      (x-1)^2 +0.174641       (x-1)   +1.25
        -0.8768         (x-2)^3 +0.97607        (x-2)^2 +1.15071
(x-2)   +1.75
S(x)=   0.6818  (x-3)^3 -1.6543 (x-3)^2 +0.47248        (x-3)   +3
        -0.1005         (x-4)^3 +0.39114        (x-4)^2 -0.79066
(x-4)   +2.5
        -0.0299         (x-5)^3 +0.08971        (x-5)^2 -0.3098 (x-5)   +2






Plot point by point; tenths of a unit between sample points.    
x value spline fit
1       1.25 
1.1     1.2678 
1.2     1.2875 
1.3     1.3112 
1.4     1.3407 
1.5     1.378 
1.6     1.4251 
1.7     1.4838 
1.8     1.5563 
1.9     1.6444  
2       1.75    same value using equations from above and below
2.1     1.874   
2.2     2.0122  
2.3     2.1594  
2.4     2.3103  
2.5     2.4598  
2.6     2.6024  
2.7     2.733   
2.8     2.8463  
2.9     2.9371  
3       3       same value using equations from above and below
3.1     3.0314  
3.2     3.0338  
3.3     3.0113  
3.4     2.9679  
3.5     2.9079  
3.6     2.8352  
3.7     2.754   
3.8     2.6683  
3.9     2.5823  
4       2.5     same value using equations from above and below
4.1     2.4247  
4.2     2.3567  
4.3     2.2953  
4.4     2.2399  
4.5     2.1899  
4.6     2.1447  
4.7     2.1037  
4.8     2.0664  
4.9     2.032   
5       2       same value using equations from above and below
5.1     1.9699 
5.2     1.9414 
5.3     1.9143 
5.4     1.8885 
5.5     1.8638 
5.6     1.84 
5.7     1.8168 
5.8     1.7943 
5.9     1.772 
6       1.75 




































1       0.325358*(C103-1)^3+0*(C103-1)^2+0.174641*(C103-1)+1.25
1.1     0.325358*(C104-1)^3+0*(C104-1)^2+0.174641*(C104-1)+1.25
1.2     0.325358*(C105-1)^3+0*(C105-1)^2+0.174641*(C105-1)+1.25
1.3     0.325358*(C106-1)^3+0*(C106-1)^2+0.174641*(C106-1)+1.25
1.4     0.325358*(C107-1)^3+0*(C107-1)^2+0.174641*(C107-1)+1.25
1.5     0.325358*(C108-1)^3+0*(C108-1)^2+0.174641*(C108-1)+1.25
1.6     0.325358*(C109-1)^3+0*(C109-1)^2+0.174641*(C109-1)+1.25
1.7     0.325358*(C110-1)^3+0*(C110-1)^2+0.174641*(C110-1)+1.25
1.8     0.325358*(C111-1)^3+0*(C111-1)^2+0.174641*(C111-1)+1.25
1.9     0.325358*(C112-1)^3+0*(C112-1)^2+0.174641*(C112-1)+1.25
2       -0.87679*(C113-2)^3+0.976076*(C113-2)^2+1.150717*(C113-2)+1.75
2.1     -0.87679*(C114-2)^3+0.976076*(C114-2)^2+1.150717*(C114-2)+1.75
2.2     -0.87679*(C115-2)^3+0.976076*(C115-2)^2+1.150717*(C115-2)+1.75
2.3     -0.87679*(C116-2)^3+0.976076*(C116-2)^2+1.150717*(C116-2)+1.75
2.4     -0.87679*(C117-2)^3+0.976076*(C117-2)^2+1.150717*(C117-2)+1.75
2.5     -0.87679*(C118-2)^3+0.976076*(C118-2)^2+1.150717*(C118-2)+1.75
2.6     -0.87679*(C119-2)^3+0.976076*(C119-2)^2+1.150717*(C119-2)+1.75
2.7     -0.87679*(C120-2)^3+0.976076*(C120-2)^2+1.150717*(C120-2)+1.75
2.8     -0.87679*(C121-2)^3+0.976076*(C121-2)^2+1.150717*(C121-2)+1.75
2.9     -0.87679*(C122-2)^3+0.976076*(C122-2)^2+1.150717*(C122-2)+1.75
3       0.681818*(C123-3)^3-1.6543*(C123-3)^2+0.472488*(C123-3)+3
3.1     0.681818*(C124-3)^3-1.6543*(C124-3)^2+0.472488*(C124-3)+3
3.2     0.681818*(C125-3)^3-1.6543*(C125-3)^2+0.472488*(C125-3)+3
3.3     0.681818*(C126-3)^3-1.6543*(C126-3)^2+0.472488*(C126-3)+3
3.4     0.681818*(C127-3)^3-1.6543*(C127-3)^2+0.472488*(C127-3)+3
3.5     0.681818*(C128-3)^3-1.6543*(C128-3)^2+0.472488*(C128-3)+3
3.6     0.681818*(C129-3)^3-1.6543*(C129-3)^2+0.472488*(C129-3)+3
3.7     0.681818*(C130-3)^3-1.6543*(C130-3)^2+0.472488*(C130-3)+3
3.8     0.681818*(C131-3)^3-1.6543*(C131-3)^2+0.472488*(C131-3)+3
3.9     0.681818*(C132-3)^3-1.6543*(C132-3)^2+0.472488*(C132-3)+3
4       -0.10044*(C133-4)^3+0.391148*(C133-4)^2-0.79066*(C133-4)+2.5
4.1     -0.10044*(C134-4)^3+0.391148*(C134-4)^2-0.79066*(C134-4)+2.5
4.2     -0.10044*(C135-4)^3+0.391148*(C135-4)^2-0.79066*(C135-4)+2.5
4.3     -0.10044*(C136-4)^3+0.391148*(C136-4)^2-0.79066*(C136-4)+2.5
4.4     -0.10044*(C137-4)^3+0.391148*(C137-4)^2-0.79066*(C137-4)+2.5
4.5     -0.10044*(C138-4)^3+0.391148*(C138-4)^2-0.79066*(C138-4)+2.5
4.6     -0.10044*(C139-4)^3+0.391148*(C139-4)^2-0.79066*(C139-4)+2.5
4.7     -0.10044*(C140-4)^3+0.391148*(C140-4)^2-0.79066*(C140-4)+2.5
4.8     -0.10044*(C141-4)^3+0.391148*(C141-4)^2-0.79066*(C141-4)+2.5
4.9     -0.10044*(C142-4)^3+0.391148*(C142-4)^2-0.79066*(C142-4)+2.5
5       -0.0299*(C143-5)^3+0.089712*(C143-5)^2-0.3098*(C143-5)+2
5.1     -0.0299*(C144-5)^3+0.089712*(C144-5)^2-0.3098*(C144-5)+2
5.2     -0.0299*(C145-5)^3+0.089712*(C145-5)^2-0.3098*(C145-5)+2
5.3     -0.0299*(C146-5)^3+0.089712*(C146-5)^2-0.3098*(C146-5)+2
5.4     -0.0299*(C147-5)^3+0.089712*(C147-5)^2-0.3098*(C147-5)+2
5.5     -0.0299*(C148-5)^3+0.089712*(C148-5)^2-0.3098*(C148-5)+2
5.6     -0.0299*(C149-5)^3+0.089712*(C149-5)^2-0.3098*(C149-5)+2
5.7     -0.0299*(C150-5)^3+0.089712*(C150-5)^2-0.3098*(C150-5)+2
5.8     -0.0299*(C151-5)^3+0.089712*(C151-5)^2-0.3098*(C151-5)+2
5.9     -0.0299*(C152-5)^3+0.089712*(C152-5)^2-0.3098*(C152-5)+2
6       -0.0299*(C153-5)^3+0.089712*(C153-5)^2-0.3098*(C153-5)+2