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