Interpolation

Chapter 3 Interpolation Interpolation is the problem of fitting a smooth curve through a given set of points, generally as the graph of a function. It is useful at least in data analysis (interpolation is a form of regression), industrial design, signal processing (digital-to-analog conversion) and in numerical analysis. It is one of those important recurring concepts in applied mathematics. In this chapter, we will immediately put interpolation to use to formulate high-order quadrature and differentiation rules. 3.1 Polynomial interpolation Given N + 1 points xj ∈ R, 0 ≤ j ≤ N , and sample values yj = f (xj ) of a function at these points, the polynomial interpolation problem consists in finding a polynomial pN (x) of degree N which reproduces those values: yj = pN (xj ), j = 0, . . . , N. In other words the graph of the polynomial should pass through PN the points (xj , yj ). A degree-N polynomial can be written as pN (x) = n=0 an xn for some coefficients a0 , . . . , aN . For interpolation, the number of degrees of freedom (N + 1 coefficients) in the polynomial matches the number of points where the function should be fit. If the degree of the polynomial is strictly less than N , we cannot in general pass it through the points (xj , yj ). We can still try to pass a polynomial (e.g., a line) in the “best approximate manner”, but this is a problem in approximation rather than interpolation; we will return to it later in the chapter on least-squares. 1 CHAPTER 3. INTERPOLATION Let us first see how the interpolation problem can be solved numerically in a direct way. Use the expression of pN into the interpolating equations yj = pN (xj ): N X an xnj = yj , j = 0, . . . , N. n=0 In these N + 1 equations indexed by j, the unknowns are the coefficients a0 , . . . , aN . We are in presence of a linear system V a = y, N X ⇔ Vjn an = yj , n=0 with V the so-called Vandermonde matrix,  1 x0 · · ·  1 x1 · · ·  V =  .. .. . . 1 xN · · · Vjn = xnj , i.e.,  xN 0  xN 1  ..  . .  xN N We can then use a numerical software like Matlab to construct the vector of abscissas xj , the right-hand-side of values yj , the V matrix, and numerically solve the system with an instruction like a = V \ y (in Matlab). This gives us the coefficients of the desired polynomial. The polynomial can now be plotted in between the grid points xj (on a finer grid), in order to display the interpolant. Historically, mathematicians such as Lagrange and Newton did not have access to computers to display interpolants, so they found explicit (and elegant) formulas for the coefficients of the interpolation polynomial. It not only simplified computations for them, but also allowed them to understand the error of polynomial interpolation, i.e., the difference f (x) − pN (x). Let us spend a bit of time retracing their steps. (They were concerned with applications such as fitting curves to celestial trajectories.) We’ll define the interpolation error from the uniform (L∞ ) norm of the difference f − pN : kf − pN k∞ := max |f (x) − pN (x)|, x where the maximum is taken over the interval [x0 , xN ]. 2 3.1. POLYNOMIAL INTERPOLATION Call PN the space of real-valued degree-N polynomials: PN = { N X an xn : an ∈ R}. n=0 Lagrange’s solution to the problem of polynomial interpolation is based on the following construction. Lemma 1. (Lagrange elementary polynomials) Let {xj , j = 0, . . . , N } be a collection of disjoint numbers. For each k = 0, . . . , N , there exists a unique degree-N polynomial Lk (x) such that  1 if j = k; Lk (xj ) = δjk = 0 if j 6= k. Proof. Fix k. Lk has roots at xj for j 6= k, so Lk must be of the form1 Y Lk (x) = C (x − xj ). j=k 6 Evaluating this expression at x = xk , we get 1=C Y (xk − xj ) ⇒ C=Q j=k 6 1 . j=k 6 (xk − xj ) Hence the only possible expression for Lk is Q j=k 6 (x − xj ) Lk (x) = Q . j=k 6 (xk − xj ) These elementary polynomials form a basis (in the sense of linear algebra) for expanding any polynomial interpolant pN . 1 That’s because, if we fix j, we can divide Lk (x) by (x − xj ), j 6= k. We obtain Lk (x) = (x − xj )q(x) + r(x), where r(x) is a remainder of lower order than x − xj , i.e., a constant. Since Lk (xj ) = 0 we must have r(x) = 0. Hence (x − xj ) must be a factor of Lk (x). The same is true of any (x − xj ) for j 6= k. There are N such factors, which exhausts the degree N . The only remaining degree of freedom in Lk is the multiplicative constant. 3 CHAPTER 3. INTERPOLATION Theorem 4. (Lagrange interpolation theorem) Let {xj , j = 0, . . . , N } be a collection of disjoint real numbers. Let {yj , j = 0, . . . , N } be a collection of real numbers. Then there exists a unique pN ∈ PN such that pN (xj ) = yj , j = 0, . . . , N. Its expression is pN (x) = N X yk Lk (x), (3.1) k=0 where Lk (x) are the Lagrange elementary polynomials. Proof. The justification that (3.1) interpolates is obvious: pN (xj ) = N X yk Lk (xj ) = k=0 N X yk Lk (xj ) = yj . k=0 It remains to see that pN is the unique interpolating polynomial. For this purpose, assume that both pN and qN take on the value yj at xj . Then rN = pN − qN is a polynomial of degree N that has a root at each of the N + 1 points x0 , . . . , xN . The fundamental theorem of algebra, however, says that a nonzero polynomial of degree N can only have N (complex) roots. Therefore, the only way for rN to have N + 1 roots is that it is the zero polynomial. So pN = qN . By definition, pN (x) = N X f (xk )Lk (x) k=0 is called the Lagrange interpolation polynomial of f at xj . Example 7. Linear interpolation through (x1 , y1 ) and (x2 , y2 ): x − x2 x − x1 L1 (x) = , L2 (x) = , x1 − x2 x2 − x1 p1 (x) = y1 L1 (x) + y2 L2 (x) y2 − y1 y1 x2 − y2 x1 = x+ x2 − x1 x2 − x1 y2 − y1 = y1 + (x − x1 ). x2 − x1 4 3.1. POLYNOMIAL INTERPOLATION Example 8. (Example (6.1) in Suli-Mayers) Consider f (x) = ex , and interpolate it by a parabola (N = 2) from three samples at x0 = −1, x1 = 0, x2 = 1. We build (x − x1 )(x − x2 ) 1 L0 (x) = = x(x − 1) (x0 − x1 )(x0 − x2 ) 2 Similarly, 1 L2 (x) = x(x + 1). 2 L1 (x) = 1 − x2 , So the quadratic interpolant is p2 (x) = e−1 L0 (x) + e0 L1 (x) + e1 L2 (x), = 1 + sinh(1) x + (cosh(1) − 1) x2 , ' 1 + 1.1752 x + 0.5431 x2 . Another polynomial that approximates ex reasonably well on [−1, 1] is the Taylor expansion about x = 0: x2 t2 (x) = 1 + x + . 2 Manifestly, p2 is not very different from t2 . (Insert picture here) Let us now move on to the main result concerning the interpolation error of smooth functions. Theorem 5. Let f ∈ C N +1 [a, b] for some N > 0, and let {xj : j = 0, . . . , N } be a collection of disjoint reals in [a, b]. Consider pN the Lagrange interpolation polynomial of f at xj . Then for every x ∈ [a, b] there exists ξ(x) ∈ [a, b] such that f (N +1) (ξ(x)) f (x) − pN (x) = πN +1 (x), (N + 1)! where πN +1 (x) = N +1 Y (x − xj ). j =1 An estimate on the interpolation error follows directly from this theorem. Set MN +1 = max |f (N +1) (x)| x∈[a,b] 5 CHAPTER 3. INTERPOLATION (which is well defined since f (N +1) is continuous by assumption, hence reaches its lower and upper bounds.) Then |f (x) − pN (x)| ≤ MN +1 |πN +1 (x)| (N + 1)! In particular, we see that the interpolation error is zero when x = xj for some j, as it should be. Let us now prove the theorem. Proof. (Can be found in Suli-Mayers, Chapter 6) In conclusion, the interpolation error: • depends on the smoothness of f via the high-order derivative f (N +1) ; • has a factor 1/(N + 1)! that decays fast as the order N → ∞; • and is directly proportional to the value of πN +1 (x), indicating that the interpolant may be better in some places than others. The natural follow-up question is that of convergence: can we always expect convergence of the polynomial interpolant as N → ∞? In other words, does the factor 1/(N + 1)! always win over the other two factors? Unfortunately, the answer is no in general. There are examples of very smooth (analytic) functions for which polynomial interpolation diverges, particularly so near the boundaries of the interplation interval. This behavior is called the Runge phenomenon, and is usually illustrated by means of the following example. Example 9. (Runge phenomenon) Let f (x) for x ∈ [−5, 5]. Interpolate it at equispaced points xj = 10j/N , where j = −N/2, . . . , N/2 and N is even. It is easy to check numerically that the interpolant diverges near the edges of [−5, 5], as N → ∞. See the Trefethen textbook on page 44 for an illustration of the Runge phenomenon. (Figure here) If we had done the same numerical experiment for x ∈ [−1, 1], the interpolant would have converged. This shows that the size of the interval matters. Intuitively, there is divergence when the size of the interval is larger than the ”features”, or characteristic length scale, of the function (here the width of the bump near the origin.) 6 3.2. POLYNOMIAL RULES FOR INTEGRATION The analytical reason for the divergence in the example above is due in no small part to the very large values taken on by πN +1 (x) far away from the origin — in contrast to the relatively small values it takes on near the origin. This is a problem intrinsic to equispaced grids. We will be more quantitative about this issue in the section on Chebyshev interpolants, where a remedy involving non-equispaced grid points will be explained. As a conclusion, polynomial interpolants can be good for small N , and on small intervals, but may fail to converge (quite dramatically) when the interpolation interval is large. 3.2 Polynomial rules for integration Rb In this section, we return to the problem of approximating a u(x)dx by a weighted sum of samples u(xj ), also called a quadrature. The plan is to form interpolants of the data, integrate those interpolants, and deduce corresponding quadrature formulas. We can formulate rules of arbitrarily high order this way, although in practice we almost never go beyond order 4 with polynomial rules. 3.2.1 Polynomial rules Without loss of generality, consider the local interpolants of u(x) formed near the origin, with x0 = 0, x1 = h and x−1 = −h. The rectangle rule does not belong in this section: it is not formed from an interpolant. • The trapezoidal rule, where we approximate u(x) by a line joining (0, u(x0 )) and (h, u(x1 )) in [0, h]. We need 2 derivatives to control the error: u00 (ξ(x)) u(x) = p1 (x) + x(x − h), 2 p1 (x) = u(0)L0 (x) + u(h)L1 (x), h−x x L0 (x) = , L1 (x) = , h h Z h Z h L0 (x)dx = L1 (x)dx = h/2, (areas of triangles) 0 0 Z h 00 | 0 u (ξ(x)) x(x − h)| dx ≤ C max |u00 (ξ)|h3 . ξ 2 7 CHAPTER 3. INTERPOLATION The result is Z h  u(x) dx = h 0 u(0) + u(h) 2  + O(h3 ). As we have seen, the terms combine as Z 1 N −1 X h h u(x) dx = u(x0 ) + h u(xj −1 ) + u(xN ) + O(h2 ). 2 2 0 j=1 • Simpson’s rule, where we approximate u(x) by a parabola through (−h, u(x−1 )), (0, u(x0 )), and (h, u(x1 )) in [−h, h]. We need three derivatives to control the error: u(x) = p2 (x) + u000 (ξ(x)) (x + h)x(x − h), 6 p2 (x) = u(−h)L−1 (x) + u(0)L0 (x) + u(h)L1 (x), x(x − h) (x + h)(x − h) (x + h)x L−1 (x) = , L (x) = , L (x) = , 0 1 2h2 −h2 2h2 Z h Z h Z h L1 (x)dx = h/3, L0 (x)dx = 4h/3, L−1 (x)dx = −h −h Z h −h 000 | 0 u (ξ(x)) (x + h)x(x − h)| dx ≤ C max |u000 (ξ)|h4 . ξ 6 The result is Z h  u(x) dx = h −h u(−h) + 4u(0) + u(h) 3  + O(h4 ). The composite Simpson’s rule is obtained by using this approximation on [0, 2h] ∪ [2h, 4h] ∪ . . . [1 − 2h, 1], adding the terms, and recognizing that the samples at 2nh (except 0 and 1) are represented twice.   Z 1 u(0) + 4u(h) + 2u(2h) + 4u(3h) + . . . + 2u(1 − 2h) + 4u(1 − h) + u(1) u(x)dx = h 3 0 It turns out that the error is in fact O(h5 ) on [−h, h], and O(h4 ) on [0, 1], a result that can be derived by using symmetry considerations (or canceling the terms from Taylor expansions in a tedious way.) For this, we need u to be four times differentiable (the constant in front of h5 involves u0000 .) 8 3.3. POLYNOMIAL RULES FOR DIFFERENTIATION The higher-order variants of polynomial rules, called Newton-Cotes rules, are not very interesting because the Runge phenomenon kicks in again. Also, the weights (like 1,4,2,4,2, etc.) become negative, which leads to unacceptable error magnification if the samples of u are not known exactly. Piecewise spline interpolation is not a good choice for numerical integration either, because of the two leftover degrees of freedom, whose arbitrary choice affects the accuracy of quadrature in an unacceptable manner. (We’ll study splines in a later section.) We’ll return to (useful!) integration rules of arbitrarily high order in the scope of spectral methods. 3.3 Polynomial rules for differentiation A systematic way of deriving finite difference formulas of higher order is to view them as derivatives of a polynomial interpolant passing through a small number of points neighboring xj . For instance (again we use −h, 0, h as reference points without loss of generality): • The forward difference at 0 is obtained from the line joining (0, u(0)) and (h, u(h)): p1 (x) = u(0)L0 (x) + u(h)L1 (x), h−x x , L1 (x) = , h h u(h) − u(0) p01 (0) = . h We already know that u0 (0) − p01 (0) = O(h). L0 (x) = • The centered difference at 0 is obtained from the line joining (−h, u(−h)) and (h, u(h)) (a simple exercise), but it is also obtained from differentiating the parabola passing through the points (−h, u(−h)), (0, u(0)), and (h, u(h)). Indeed, p2 (x) = u(−h)L−1 (x) + u(0)L0 (x) + u(h)L1 (x), L−1 (x) = x(x − h) , 2h2 (x + h)(x − h) (x + h)x , L1 (x) = , 2 −h 2h2 2x − h 2x 2x + h p02 (x) = u(−h) + u(0) 2 + u(h) , 2 2h −h 2h2 L0 (x) = 9 CHAPTER 3. INTERPOLATION u(h) − u(−h) . 2h We already know that u0 (0) − p02 (0) = O(h2 ). p02 (0) = • Other examples can be considered, such as the centered second difference (-1 2 -1), the one-sided first difference (-3 4 -1), etc. Differentiating one-sided interpolation polynomials is a good way to obtain one-sided difference formulas, which comes in handy at the boundaries of the computational domain. The following result establishes that the order of the finite difference formula matches the order of the polynomial being differentiated. Theorem 6. Let f ∈ C N +1 [a, b], {xj , j = 0, . . . , N } some disjoint points, and pN the corresponding interpolation polynomial. Then f 0 (x) − p0N (x) = f (N +1) (ξ) ∗ πN (x), N! ∗ for some ξ ∈ [a, b], and where πN (x) = (x − η1 ) . . . (x − ηN ) for some ηj ∈ [xj −1 , xj ]. The proof of this result is an application of Rolle’s theorem that we leave out. (It is not simply a matter of differentiating the error formula for polynomial interpolation, because we have no guarantee on dξ/dx.) A consequence of this theorem is that the error in computing the deriva∗ tive is a O(hN ) (which comes from a bound on the product πN (x).) It is interesting to notice, at least empirically, that the Runge’s phenomenon is absent when the derivative is evaluated at the center of the interval over which the interpolant is built. 3.4 Piecewise polynomial interpolation The idea of piecewise polynomial interpolation, also called spline interpolation, is to subdivide the interval [a, b] into a large number of subintervals [xj −1 , xj ], and to use low-degree polynomials over each subintervals. This helps avoiding the Runge phenomenon. The price to pay is that the interpolant is no longer a C ∞ function – instead, we lose differentiability at the junctions between subintervals, where the polynomials are made to match. 10 3.4. PIECEWISE POLYNOMIAL INTERPOLATION If we use polynomials of order n, then the interpolant is piecewise C ∞ , and overall C k , with k ≤ n − 1. We could not expect to have k = n, because it would imply that the polynomials are identical over each subinterval. We’ll see two examples: linear splines when n = 1, and cubic splines when n = 3. (We have seen in the homework why it is a bad idea to choose n = 2.) 3.4.1 Linear splines We wish to interpolate a continuous function f (x) of x ∈ [a, b], from the knowledge of f (xj ) at some points xj , j = 0, . . . , N , not necessarily equispaced. Assume that x0 = a and xN = b. The piecewise linear interpolant is build by tracing a straight line between the points (xj−1 , f (xj−1 )) and (xj , f (xj ))); for j = 1, . . . , N the formula is simply sL (x) = xj −1 − x x − xj −1 f (xj −1 ) + f (xj ), xj − xj −1 xj − xj −1 x ∈ [xj −1 , xj ]. In this case we see that sL (x) is a continuous function of x, but that its derivative is not continuous at the junction points, or nodes xj . If the function is at least twice differentiable, then piecewise linear interpolation has second-order accuracy, i.e., an error O(h2 ). Theorem 7. Let f ∈ C 2 [a, b], and let h = maxj=1,...,N (xj − xj−1 ) be the grid diameter. Then h2 kf − sL kL∞ [a,b] ≤ kf 00 kL∞ [a,b] . 8 Proof. Let x ∈ [xj −1 , xj ] for some j = 1, . . . , N . We can apply the basic result of accuracy of polynomial interpolation with n = 1: there exists ξ ∈ [xj−1 , xj ] such that 1 f (x) − sL (x) = f 00 (ξ) (x − xj −1 )(x − xj ), 2 x ∈ [xj−1 , xj ]. Let hj = xj − xj−1 . It is easy to check that the product (x − xj −1 )(x − xj ) x +x takes its maximum value at the midpoint j−12 j , and that the value is h2j /4. We then have |f (x) − sL (x)| ≤ h2j 8 max |f 00 (ξ)|, ξ∈[xj−1 ,xj ] x ∈ [xj−1 , xj ]. The conclusion follows by taking a maximum over j. 11 CHAPTER 3. INTERPOLATION We can express any piecewise linear interpolant as a superposition of “tent” basis functions φk (x): X sL (x) = ck φk (x), k where φk (x) is the piecewise linear function equal to 1 at x = xk , and equal to zero at all the other grid points xj , j 6= k. Or in short, φk (xj ) = δjk . On an equispaced grid xj = jh, an explicit formula is φk (x) = 1 S(1) (x − xk ), h where S(1) (x) = x+ − 2(x − h)+ + (x − 2h)+ , and where x+ denotes the positive part of x (i.e., x if x ≥ 0, and zero if x < 0.) Observe that we can simply take ck = f (xk ) above. The situation will be more complicated when we pass to higher-order polynomials. 3.4.2 Cubic splines Let us now consider the case n = 3 of an interpolant which is a third-order polynomial, i.e., a cubic, on each subinterval. The most we can ask is that the value of the interpolant, its derivative, and its second derivative be continuous at the junction points xj . Definition 5. (Interpolating cubic spline) Let f ∈ C[a, b], and {xj , j = 0, . . . , N } ⊂ [a, b]. An interpolating cubic spline is a function s(x) such that 1. s(xj ) = f (xj ); 2. s(x) is a polynomial of degree 3 over each segment [xj−1 , xj ]; 3. s(x) is globally C 2 , i.e., at each junction point xj , we have the relations + s(x− j ) = s(xj ), 0 + s0 (x− j ) = s (xj ), 00 + s00 (x− j ) = s (xj ), + where the notations s(x− j ) and s(xj ) refer to the adequate limits on the left and on the right. 12 3.4. PIECEWISE POLYNOMIAL INTERPOLATION Let us count the degrees of freedom. A cubic polynomial has 4 coefficients. There are N + 1 points, hence N subintervals, for a total of 4N numbers to be specified. The interpolating conditions s(xj ) = f (xj ) specify two degrees of freedom per polynomial: one value at the left endpoint xj −1 , and one value at the right endpoint xj . That’s 2N conditions. Continuity of s(x) follows automatically. The continuity conditions on s0 , s00 are imposed only at the interior grid points xj for j = 1, . . . , N − 1, so this gives rise to 2(N − 1) additional conditions, for a total of 4N − 2 equations. There is a mismatch between the number of unknowns (4N ) and the number of conditions (4N − 2). Two more degrees of freedom are required to completely specify a cubic spline interpolant, hence the precaution to write “an” interpolant in the definition above, and not “the” interpolant. The most widespread choices for fixing these two degrees of freedom are: • Natural splines: s00 (x0 ) = s00 (xN ) = 0. If s(x) measures the displacement of a beam, a condition of vanishing second derivative corresponds to a free end. • Clamped spline: s0 (x0 ) = p0 , s0 (xN ) = pN , where p0 and pN are specified values that depend on the particular application. If s(x) measures the displacement of a beam, a condition of vanishing derivative corresponds to a horizontal clamped end. • Periodic spline: assuming that s(x0 ) = s(xN ), then we also impose s0 (x0 ) = s0 (xN ), s00 (x0 ) = s00 (xN ). Let us now explain the algorithm most often used for determining a cubic spline interpolant, i.e., the 4N coefficients of the N cubic polynomials, from the knowledge of f (xj ). Let us consider the natural spline. It is advantageous to write a system of equations for the second derivatives at the grid points, that we denote σj = s00 (xj ). Because s(x) is piecewise cubic, we know that s00 (x) is piecewise linear. Let hj = xj − xj−1 . We can write xj − x x − xj −1 s00 (x) = σj −1 + σj , x ∈ [xj−1 , xj ]. hj hj 13 CHAPTER 3. INTERPOLATION Hence (xj − x)3 (x − xj−1 )3 s(x) = σj − 1 + σj + αj (x − xj −1 ) + βj (xj − x). 6hj 6hj We could have written ax + b for the effect of the integration constants in the equation above, but writing it in terms of αj and βj makes the algebra that follows simpler. The two interpolation conditions for [xj −1 , xj ] can be written s(xj −1 ) = f (xj −1 ) σj −1 h2j f (xj −1 ) = + hj βj , 6 ⇒ σj h2j + hj αj . 6 One then isolates αj , βj in this equation; substitutes those values in the 0 + equation for s(x); and evaluates the relation s0 (x− j ) = s (xj ). Given that σ0 = σN = 0, we end up with a system of N − 1 equations in the N − 1 unknowns σ1 , . . . , σN . Skipping the algebra, the end result is   f (xj+1 ) − f (xj ) f (xj ) − f (xj −1 ) − . hj σj −1 + 2(hj+1 +hj )σj + hj+1 σj+1 = 6 hj+1 hj s(xj ) = f (xj ) ⇒ f (xj ) = We are in presence of a tridiagonal system for σj . It can be solved efficiently with Gaussian elimination, yielding a LU decomposition with bidiagonal factors. Unlike in the case of linear splines, there is no way around the fact that a linear system needs to be solved. Notice that the tridiagonal matrix of the system above is diagonally dominant (each diagonal element is strictly greater than the sum of the other elements on the same row, in absolute value), hence it is always invertible. One can check that the interpolation for cubic splines is O(h4 ) well away from the endpoints. This requires an analysis that is too involved for the present set of notes. Finally, like in the linear case, let us consider the question of expanding a cubic spline interpolant as a superposition of basis functions X s(x) = ck φk (x). k There are many ways of choosing the functions φk (x), so let us specify that they should have as small a support as possible, that they should have the 14 3.4. PIECEWISE POLYNOMIAL INTERPOLATION same smoothness C 2 as s(x) itself, and that they should be translates of each other when the grid xj is equispaced with spacing h. The only solution to this problem is the cubic B-spline. Around any interior grid point xk , is supported in the interval [xk−2 , xk+2 ], and is given by the formula 1 φk (x) = 3 S(3) (x − xk−2 ), 4h where 3 S(3) (x) = x+ − 4(x − h)3+ + 6(x − 2h)3+ − 4(x − 3h)3+ + (x − 4h)3+ , and where x+ is the positive part of x. One can check that φk (x) takes the value 1 at xk , 1/4 at xk±1 , zero outside of [xk−2 , xk+2 ], and is C 2 at each junction. It is a bell-shaped curve. it is an interesting exercise to check that it can be obtained as the convolution of two “tent” basis functions of linear interpolation: S(3) (x) = cS(1) (x) ∗ S(1) (x), where c is some constant, and the symbol ∗ denotes convolution: Z f ∗ g(x) = f (y)g(x − y) dy. R Now with cubic B-splines, we cannot put ck = f (xk ) anymore, since 6 δjk . The particular values of ck are the result of solving a linφk (xj ) = ear system, as mentioned above. Again, there is no way around solving a linear system. In particular, if f (xk ) changes at one point, or if we add a datum point (xk , fk ), then the update requires re-computing the whole interpolant. Changing one point has ripple effects throughout the whole interval [a, b]. This behavior is ultimately due to the more significant overlap between neighboring φk in the cubic case than in the linear case. 15