Document Sample

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” Deﬁnition: A cubic spline S(x) is a piecewise-deﬁned function that satisﬁes 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 coefﬁcients. October 16, 2008 c Steven Rauch and John Stockie Cubic Splines MACM 316 2/14 Derivation Using the deﬁnition, let’s determine equations relating the coef- ﬁcients 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 simpliﬁed considerably by deﬁning ′′ 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 ﬂat 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 speciﬁed): ′ 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 ﬁrst/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 ﬁrst/last RHS entries are modiﬁed 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 ﬁrst/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 coefﬁcients 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 coefﬁcients: 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 deﬁning 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 satisﬁes 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 ﬂat re- gion – non-smooth data cause some problems for cubic splines. October 16, 2008 c Steven Rauch and John Stockie

DOCUMENT INFO

Shared By:

Categories:

Tags:
cubic splines, cubic spline, natural cubic spline, control points, data points, natural spline, spline interpolation, second derivative, Cubic Spline Interpolation, first derivative

Stats:

views: | 104 |

posted: | 1/22/2011 |

language: | English |

pages: | 14 |

OTHER DOCS BY dfsdf224s

How are you planning on using Docstoc?
BUSINESS
PERSONAL

By registering with docstoc.com you agree to our
privacy policy and
terms of service, and to receive content and offer notifications.

Docstoc is the premier online destination to start and grow small businesses. It hosts the best quality and widest selection of professional documents (over 20 million) and resources including expert videos, articles and productivity tools to make every small business better.

Search or Browse for any specific document or resource you need for your business. Or explore our curated resources for Starting a Business, Growing a Business or for Professional Development.

Feel free to Contact Us with any questions you might have.