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