Math 104A Homework #4 Instructor: Lihui Chai 1. Write a code to compute a natural spline S(x) which interpolates a collection of given points (x0 , y0 ), (x1 , y1 ), ..., (xn , yn ) where x0 < x1 < x2 < ... < xn (do not assume they are equidistributed). Your code should have a triadiagonal solver for the resulting linear system of equations (youre not allowed to use Matlabs \\ operator to solve the linear system). 2. One important application of spline interpolation is the construction of smooth curves that are not necessarily the graph of a function but that have a parametric representation x = x(t) and y = y(t) for t [a, b]. Hence one needs to determine two splines interpolating (tj , xj ) and (tj , yj ) (j = 0, 1, ...n). The arc length of the curve is a natural choice for the parameter t. However, this is not known a priori and instead the tj s are usually chosen as the distances of consecutive points: q t0 = 0, tj = tj1 + (xj xj1 )2 + (yj yj1 )2 , j = 1, 2, ...n. Use the values in Table 1 to construct a smooth parametric representation of a curve passing through the points (xj , yj ), j = 0, 1, ..., 8 by finding the two natural cubic splines interpolating (tj , xj ) and (tj , yj ), j = 0, 1, ...8, respectively. Tabulate the coefficients of the splines and plot the resulting (parametric) curve. All course materials (class lectures and discussions, handouts, homework assignments, examinations, web materials) and the intellectual content of the course itself are protected by United States Federal Copyright Law, the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lectures notes and all other course materials without the prior written permission of Prof. Hector D. Ceniceros. 1 Table 1 j 0 1 2 3 4 5 6 7 8 tj 0 0.618 0.935 1.255 1.636 1.905 2.317 2.827 3.330 xj 1.50 0.90 0.60 0.35 0.20 0.10 0.50 1.00 1.50 2 yj 0.75 0.90 1.00 0.80 0.45 0.20 0.10 0.20 0.25 Piece-wise Linear Interpolation and Splines Hector D. Ceniceros 1 Piece-wise Linear Interpolation One way to reduce the error in linear interpolation is to divide [a, b] into small subintervals [x0 , x1 ], ..., [xn1 , xn ]. In each of the subintervals [xj , xj+1 ] we approximate f by (1) f (xj+1 ) f (xj ) (x xj ), xj+1 xj x [xj , xj+1 ]. 1 f (x) p(x) = f 00 ()(x xj )(x xj+1 ), 2 x [xj , xj+1 ] p(x) = f (xj ) + We know that (2) where is some point between xj and xj+1 . Suppose that |f 00 (x)| M2 , x [a, b] then (3) 1 |f (x) p(x)| M2 max |(x xj )(x xj+1 )|. xj xxj+1 2 Now the max at the right hand side is attained at the midpoint (xj + xj+1 )/2 and \u0012 \u00132 xj+1 xj 1 (4) max |(x xj )(x xj+1 )| = = h2j , xj xxj+1 2 4 These are lecture notes for Math 104 A. These notes and all course materials are protected by United States Federal Copyright Law, the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lectures notes and all other course materials without the prior written permission of the instructor. 1 where hj = xj+1 xj . Therefore (5) 1 |f (x) p(x)| M2 h2j . xj xxj+1 8 max If we want this error to be smaller than a prescribed tolerance we can take sufficiently small subintervals. Namely, we can pick hj such that 18 M2 h2j which implies that r 8 hj (6) . M2 2 Cubic Splines Several applications require a smoother curve than that provided by a piecewise linear approximation. Continuity of the first and second derivatives provide that required smoothness. One of the most frequently used such approximations is a cubic spline, which is is a piecewise cubic function, s(x), which interpolates a set points (x0 , f0 ), (x1 , f1 ), . . . (xn , fn ), and has two continuous derivatives. In each subinterval [xj , xj+1 ], s(x) is a cubic polynomial, which we may represent as (7) sj (x) = aj (x xj )3 + bj (x xj )2 + cj (x xj ) + dj . Let (8) hj = xj+1 xj . The spline s(x) interpolates the given data: (9) (10) sj (xj ) = fj = dj , sj (xj+1 ) = aj h3j + bj h2j + cj hj + dj = fj+1 . Now s0j (x) = 3aj (x xj )2 + 2bj (x xj ) + cj and s00j (x) = 6aj (x xj ) + 2bj . Therefore (11) (12) (13) (14) s0j (xj ) = cj , s0j (xj+1 ) = 3aj h2j + 2bj hj + cj , s00j (xj ) = 2bj , s00j (xj+1 ) = 6aj hj + 2bj . 2 We are going to write the spline coefficients aj , bj , cj , and dj in terms of fj and fj+1 and the unknown values zj = s00j (xj ) and zj+1 = s00j (xj+1 ). We have dj = fj , 1 bj = zj , 2 6aj hj + 2bj = zj+1 aj = 1 (zj+1 zj ) 6hj and substituting these values in (10) we get cj = 1 1 (fj+1 fj ) hj (zj+1 + 2zj ). hj 6 Let us collect all our formulas for the spline coefficients: (15) (16) (17) (18) 1 (zj+1 zj ), 6hj 1 bj = zj , 2 1 1 cj = (fj+1 fj ) hj (zj+1 + 2zj ), hj 6 dj = f j . aj = Note that the second derivative of s is continuous, s00j (xj+1 ) = zj+1 = s00j+1 (xj+1 ), and by construction s interpolates the given data. We are now going to use the condition of continuity of the first derivative of s to determine equations for the values of the second derivative. s0j (xj+1 ) = 3aj h2j + 2bj hj + cj 1 1 1 1 =3 (zj+1 zj )h2j + 2 zj hj + (fj+1 fj ) hj (zj+1 + 2zj ) 6hj 2 hj 6 1 1 = (fj+1 fj ) + hj (2zj+1 + zj ), hj 6 Decreasing the index by 1 we get (19) s0j1 (xj ) = 1 (fj fj1 ) + hj1 (2zj + zj1 ) hj1 6 1 3 Continuity of the first derivative at an interior node means s0j1 (xj ) = s0j (xj ) for j = 1, 2, ..., n 1. Therefore 1 1 1 1 (fj fj1 ) + hj1 (2zj + zj1 ) = cj = (fj+1 fj ) hj (zj+1 + 2zj ) hj1 6 hj 6 which can be written as (20) hj1 zj1 + 2(hj1 + hj )zj + hj zj+1 = 6 hj1 (fj fj1 ) + 6 (fj+1 fj ) hj for j = 1, 2, . . . , n 1. This is a linear system of n 1 equations for the n 1 unknowns z1 , z2 , . . . , zn1 . In matrix form 2(h0 + h1 ) h1 0 d1 z1 . z2 0 h1 2(h1 + h2 ) h2 . . d2 (21) , = . . .. .. .. .. .. . . hn2 . ... dn1 0 h 2(h + h ) zn1 n2 n2 n1 where (22) d1 d2 .. . h60 (f1 f0 ) + h61 (f2 f1 ) h0 f000 h61 (f2 f1 ) + h62 (f3 f2 ) .. . . = 6 6 dn2 hn3 (fn2 fn3 ) + hn2 (fn1 fn2 ) 6 6 dn1 hn2 (fn1 fn2 ) + hn1 (fn fn1 ) hn1 fn00 Note that z0 = f000 and zn = fn00 are unspecified. The values z0 = zn = 0 defined what is called a Natural Spline. The matrix of the linear system (21) is diagonally dominant and hence nonsingular so there is a unique solution for z1 , z2 , . . . , zn1 corresponding the second derivative values of the spline at the interior nodes. If the nodes are equidistributed, xj = x0 + jh for j = 0, 1, . . . , n then the linear system ((21), after dividing by h simplifies to 6 (f 2f + f ) z z1 0 1 2 0 4 1 0 h 6 z2 (f 2f2 + f3 ) .. h 1 1 4 1 . 0 . .. . . (23) = . . . ... .. 1 1 6 h (fn3 2fn2 + fn1 ) z n2 . 6 0 .. 1 4 zn1 (f 2fn1 + fn ) zn h n2 4 Once the z1 , z2 , . . . , zn1 are found the spline coefficients can be computed from (15)-(18). Example 1. Find the natural cubic spline that interpolates (0, 0), (1, 1), (2, 16). We know z0 = 0 and z2 = 0. We only need to find z1 (only 1 interior node). The system (23) degenerates to just one equation and h = 1, thus z0 + 4z1 + z2 = 6[f0 2f1 + f2 ] z1 = 21 In [0, 1] we have 1 1 7 (z1 z0 ) = 21 = , 6h 6 2 1 b0 = z0 = 0 2 1 1 5 1 c0 = (f1 f0 ) h(z1 + 2z0 ) = 1 21 = , h 6 6 2 d0 = f0 = 0. a0 = Thus, s0 (x) = a0 (x 0)3 + b0 (x 0)3 + c0 (x 0) + d0 = 27 x3 25 x. Now in [1, 2] 1 1 7 (z2 z1 ) = (21) = , 6h 6 2 21 1 b1 = z1 = , 2 2 1 1 1 c1 = (f2 f1 ) h(z2 + 2z1 ) = 16 1 (2 21) = 8, h 6 6 d1 = f1 = 1. a1 = and s1 (x) = 72 (x 1)3 + 21 (x 1)2 + 8(x 1) + 1. Therefore the spline is 2 given by \u001a 7 3 x 25 x x [0, 1], 2 s(x) = 7 21 3 2 (x 1) + 2 (x 1)2 + 8(x 1) + 1 x [1, 2]. 5 2.1 Solving the Tridiagonal System The matrix of coefficients of the linear system (21) has the tridiagonal form (24) a1 b 1 c 1 a2 b 2 . . c2 . . . . ... ... ... A= ... ... .. . ... .. . cN 1 bN 1 aN where for the natural splines (n 1) (n 1) system (21), the non-zero tridiagonal entries are (25) (26) (27) j = 1, 2, . . . , n 1 j = 1, 2, . . . , n 2 j = 1, 2, . . . , n 2. aj = 2(hj1 + hj ), bj = hj , cj = hj1 N.B. Do not confuse these matrix coefficients with the spline coefficients aj , bj , cj , and dj given by (15)-(18). We can solve the corresponding linear system of equations Ax = d by factoring this matrix A into the product of a lower triangular matrix L and an upper triangular matrix U . To illustrate the idea let us take N = 5. A 5 5 tridiagonal linear system has the form a1 b 1 0 0 0 x1 d1 c1 a2 b2 0 0 x2 d2 0 c2 a3 b3 0 x3 = d3 (28) 0 0 c3 a4 b4 x4 d4 x5 d5 0 0 0 c 4 a5 and we a1 c1 0 0 0 seek a factorization of the form b1 0 0 0 1 0 0 l1 1 0 a2 b 2 0 0 c 2 a3 b 3 0 = 0 l2 1 0 c3 a4 b4 0 0 l3 0 0 0 0 0 c 4 a5 6 0 0 0 1 l4 0 m1 u1 0 0 0 0 0 0 m2 u2 0 0 0 0 m u 0 3 3 0 0 0 0 m4 u4 1 0 0 0 0 m5 Note the the first matrix on the right hand side is lower triangular and the second one is upper triangular. Performing the product of the matrices and comparing with the corresponding entries of the left hand side matrix we have 1st row: 2nd row: 3rd row: 4th row: 5th row: a1 c1 c2 c3 c4 = m1 , b1 = u1 , = m1 l1 , a2 = l1 u1 + m2 , b2 = m2 l2 , a3 = l2 u2 + m3 , b3 = m3 l3 , a4 = l3 u3 + m4 , b4 = m4 l4 , a5 = l4 u4 + m5 , b5 = u2 , = u3 , = u4 , = u5 . So we can determine the unknowns in the following order m1 ; u1 , l1 , m2 ; u2 , l2 , m3 ; . . . , u4 , l4 , m5 . Since uj = bj for all j we can write down the algorithm for general N as % Determine the factorization coefficients m1 = a1 for j = 1 : N 1 lj = cj /mj mj+1 = aj+1 lj bj end % Forward substitution on Ly = d y 1 = d1 for j = 2 : N yj = dj lj1 yj1 end % Backward substitution to solve U x = y xN = yN /mN for j = N 1 : 1 xj = (yj bj xj+1 )/mj end 7 2.2 Complete Splines Sometimes it is more appropriate to specify the first derivative at the end points than the second derivative. This is called a complete or clamped spline. In this case z0 = f000 and zn = fn00 become unknowns together with z1 , z2 , . . . , zn1 . We need to more equations to have a complete system for all the n + 1 unknown values of the second derivative in a complete spline. Recall that sj (x) = aj (x xj )3 + bj (x xj )2 + cj (x xj ) + dj and so sj (x) = 3aj (x xj )2 + 2bj (x xj ) + cj . Therefore (29) (30) s00 (x0 ) = c0 = f00 s0n1 (xn ) = 3an1 h2n1 + 2bn1 hn1 + cn1 = fn0 . Substituting c0 , an1 , bn1 , and cn1 from (15)-(17) we get (31) (32) 6 (f1 f0 ) 6f00 , h0 6 hn1 zn1 + 2hn1 zn = (fn fn1 ) + 6fn0 . hn1 2h0 z0 + h0 z1 = These two equations together with (21) uniquely determine the second derivative values at all the nodes. The resulting (n + 1) (n + 1) is also tridiagonal and diagonally dominant (hence nonsingular). Once the values z0 , z1 , . . . , zn are found the splines coefficients are obtained from (15)-(18). 2.3 Parametric Curves In graphics and computer animation applications it is often required to construct smooth curves that are not necessarily the graph of a function but that have a parametric representation x = x(t) and y = y(t) for t [a, b]. Hence one needs to determine two splines interpolating (tj , xj ) and (tj , yj ) (j = 0, 1, . . . n), respectively. 8 The arc length of the curve is a natural choice for the parameter t. However, this is not known a priori and instead the nodes tj 's are usually chosen as the distances of consecutive, judiciously chosen points: q (33) t0 = 0, tj = tj1 + (xj xj1 )2 + (yj yj1 )2 , j = 1, 2, . . . n. 9