This preview shows page 1. Sign up to view the full content.
Unformatted text preview: Chapter 3 Interpolation Interpolation is the process of defining a function that takes on specified values at specified points. This chapter concentrates on two closely related interpolants: the piecewise cubic spline and the shape-preserving piecewise cubic named "pchip." 3.1 The Interpolating Polynomial We all know that two points determine a straight line. More precisely, any two points in the plane, (x1 , y1 ) and (x2 , y2 ), with x1 = x2 , determine a unique firstdegree polynomial in x whose graph passes through the two points. There are many different formulas for the polynomial, but they all lead to the same straight line graph. This generalizes to more than two points. Given n points in the plane, (xk , yk ), k = 1, . . . , n, with distinct xk 's, there is a unique polynomial in x of degree less than n whose graph passes through the points. It is easiest to remember that n, the number of data points, is also the number of coefficients, although some of the leading coefficients might be zero, so the degree might actually be less than n - 1. Again, there are many different formulas for the polynomial, but they all define the same function. This polynomial is called the interpolating polynomial because it exactly reproduces the given data: P (xk ) = yk , k = 1, . . . , n. Later, we examine other polynomials, of lower degree, that only approximate the data. They are not interpolating polynomials. The most compact representation of the interpolating polynomial is the Lagrange form x - xj yk . P (x) = xk - xj
k j=k February 15, 2008 1 2 Chapter 3. Interpolation There are n terms in the sum and n - 1 terms in each product, so this expression defines a polynomial of degree at most n - 1. If P (x) is evaluated at x = xk , all the products except the kth are zero. Furthermore, the kth product is equal to one, so the sum is equal to yk and the interpolation conditions are satisfied. For example, consider the following data set. x = 0:3; y = [-5 -6 The command disp([x; displays 0 -5 1 -6 2 -1 3 16 y]) -1 16]; The Lagrangian form of the polynomial interpolating these data is P (x) = (x - 1)(x - 2)(x - 3) x(x - 2)(x - 3) (-5) + (-6) (-6) (2) x(x - 1)(x - 3) x(x - 1)(x - 2) + (-1) + (16). (-2) (6) We can see that each term is of degree three, so the entire sum has degree at most three. Because the leading term does not vanish, the degree is actually three. Moreover, if we plug in x = 0, 1, 2, or 3, three of the terms vanish and the fourth produces the corresponding value from the data set. Polynomials are not usually represented in their Lagrangian form. More frequently, they are written as something like x3 - 2x - 5. The simple powers of x are called monomials, and this form of a polynomial is said to be using the power form. The coefficients of an interpolating polynomial using its power form, P (x) = c1 xn-1 + c2 xn-2 + + cn-1 x + cn , can, in principle, be computed by solving a system of simultaneous linear equations n-1 c1 y1 x1 xn-2 x1 1 1 xn-1 xn-2 x2 1 c2 y2 2 2 . = . . 1 . . . . xn-1 xn-2 xn 1 n n cn yn The matrix V of this linear system is known as a Vandermonde matrix. Its elements are vk,j = xn-j . k 3.1. The Interpolating Polynomial 3 The columns of a Vandermonde matrix are sometimes written in the opposite order, but polynomial coefficient vectors in Matlab always have the highest power first. The Matlab function vander generates Vandermonde matrices. For our example data set, V = vander(x) generates V = 0 1 8 27 Then c = V\y' computes the coefficients. c = 1.0000 0.0000 -2.0000 -5.0000 In fact, the example data were generated from the polynomial x3 - 2x - 5. Exercise 3.6 asks you to show that Vandermonde matrices are nonsingular if the points xk are distinct. But Exercise 3.18 asks you to show that a Vandermonde matrix can be very badly conditioned. Consequently, using the power form and the Vandermonde matrix is a satisfactory technique for problems involving a few well-spaced and well-scaled data points. But as a general-purpose approach, it is dangerous. In this chapter, we describe several Matlab functions that implement various interpolation algorithms. All of them have the calling sequence v = interp(x,y,u) The first two input arguments, x and y, are vectors of the same length that define the interpolating points. The third input argument, u, is a vector of points where the function is to be evaluated. The output v is the same length as u and has elements v(k)=interp(x,y,u(k)) Our first such interpolation function, polyinterp, is based on the Lagrange form. The code uses Matlab array operations to evaluate the polynomial at all the components of u simultaneously. 0 1 4 9 0 1 2 3 1 1 1 1 4 function v = polyinterp(x,y,u) n = length(x); v = zeros(size(u)); for k = 1:n w = ones(size(u)); for j = [1:k-1 k+1:n] w = (u-x(j))./(x(k)-x(j)).*w; end v = v + w*y(k); end Chapter 3. Interpolation To illustrate polyinterp, create a vector of densely spaced evaluation points. u = -.25:.01:3.25; Then v = polyinterp(x,y,u); plot(x,y,'o',u,v,'-') creates Figure 3.1.
25 20 15 10 5 0 -5 -10 -0.5 0 0.5 1 1.5 2 2.5 3 3.5 Figure 3.1. polyinterp. The polyinterp function also works correctly with symbolic variables. For example, create symx = sym('x') Then evaluate and display the symbolic form of the interpolating polynomial with 3.2. Piecewise Linear Interpolation P = polyinterp(x,y,symx) pretty(P) which produces -5 (-1/3 x + 1)(-1/2 x + 1)(-x + 1) - 6 (-1/2 x + 3/2)(-x + 2)x -1/2 (-x + 3)(x - 1)x + 16/3 (x - 2)(1/2 x - 1/2)x 5 This expression is a rearrangement of the Lagrange form of the interpolating polynomial. Simplifying the Lagrange form with P = simplify(P) changes P to the power form P = x^3-2*x-5 Here is another example, with a data set that is used by the other methods in this chapter. x = 1:6; y = [16 18 21 17 15 12]; disp([x; y]) u = .75:.05:6.25; v = polyinterp(x,y,u); plot(x,y,'o',u,v,'-'); produces 1 16 2 18 3 21 4 17 5 15 6 12 and Figure 3.2. Already in this example, with only six nicely spaced points, we can begin to see the primary difficulty with full-degree polynomial interpolation. In between the data points, especially in the first and last subintervals, the function shows excessive variation. It overshoots the changes in the data values. As a result, fulldegree polynomial interpolation is hardly ever used for data and curve fitting. Its primary application is in the derivation of other numerical methods. 3.2 Piecewise Linear Interpolation You can create a simple picture of the data set from the last section by plotting the data twice, once with circles at the data points and once with straight lines connecting the points. The following statements produce Figure 3.3. x = 1:6; y = [16 18 21 17 15 12]; plot(x,y,'o',x,y,'-'); 6
Full degree polynomial interpolation 22 Chapter 3. Interpolation 20 18 16 14 12 10 1 2 3 4 5 6 Figure 3.2. Full-degree polynomial interpolation. To generate the lines, the Matlab graphics routines use piecewise linear interpolation. The algorithm sets the stage for more sophisticated algorithms. Three quantities are involved. The interval index k must be determined so that The local variable, s, is given by The first divided difference is xk x < xk+1 . s = x - xk . k = yk+1 - yk . xk+1 - xk yk+1 - yk xk+1 - xk With these quantities in hand, the interpolant is L(x) = yk + (x - xk ) = yk + sk . This is clearly a linear function that passes through (xk , yk ) and (xk+1 , yk+1 ). The points xk are sometimes called breakpoints or breaks. The piecewise linear interpolant L(x) is a continuous function of x, but its first derivative, L (x), is not continuous. The derivative has a constant value, k , on each subinterval and jumps at the breakpoints. Piecewise linear interpolation is implemented in piecelin.m. The input u can be a vector of points where the interpolant is to be evaluated, so the index k is actually a vector of indices. Read this code carefully to see how k is computed. 3.2. Piecewise Linear Interpolation
Piecewise linear interpolation 22 7 20 18 16 14 12 10 1 2 3 4 5 6 Figure 3.3. Piecewise linear interpolation. function v = piecelin(x,y,u) %PIECELIN Piecewise linear interpolation. % v = piecelin(x,y,u) finds the piecewise linear L(x) % with L(x(j)) = y(j) and returns v(k) = L(u(k)). % First divided difference delta = diff(y)./diff(x); % Find subinterval indices k so that x(k) <= u < x(k+1) n = length(x); k = ones(size(u)); for j = 2:n-1 k(x(j) <= u) = j; end % Evaluate interpolant s = u - x(k); v = y(k) + s.*delta(k); 8 Chapter 3. Interpolation 3.3 Piecewise Cubic Hermite Interpolation Many of the most effective interpolation techniques are based on piecewise cubic polynomials. Let hk denote the length of the kth subinterval: hk = xk+1 - xk . Then the first divided difference, k , is given by k = yk+1 - yk . hk Let dk denote the slope of the interpolant at xk : dk = P (xk ). For the piecewise linear interpolant, dk = k-1 or k , but this is not necessarily true for higher order interpolants. Consider the following function on the interval xk x xk+1 , expressed in terms of local variables s = x - xk and h = hk : P (x) = 3hs2 - 2s3 h3 - 3hs2 + 2s3 yk+1 + yk h3 h3 2 2 s (s - h) s(s - h) + dk+1 + dk . 2 h h2 This is a cubic polynomial in s, and hence in x, that satisfies four interpolation conditions, two on function values and two on the possibly unknown derivative values: P (xk ) = yk , P (xk+1 ) = yk+1 , P (xk ) = dk , P (xk+1 ) = dk+1 . Functions that satisfy interpolation conditions on derivatives are known as Hermite or osculatory interpolants, because of the higher order contact at the interpolation sites. ("Osculari" means "to kiss" in Latin.) If we happen to know both function values and first derivative values at a set of data points, then piecewise cubic Hermite interpolation can reproduce those data. But if we are not given the derivative values, we need to define the slopes dk somehow. Of the many possible ways to do this, we will describe two, which Matlab calls pchip and spline. 3.4 Shape-Preserving Piecewise Cubic The acronym pchip abbreviates "piecewise cubic Hermite interpolating polynomial." Although it is fun to say, the name does not specify which of the many possible interpolants is actually being used. In fact, spline interpolants are also piecewise cubic Hermite interpolating polynomials, but with different slopes. Our 3.4. Shape-Preserving Piecewise Cubic
Shape-preserving Hermite interpolation 22 9 20 18 16 14 12 10 1 2 3 4 5 6 Figure 3.4. Shape-preserving piecewise cubic Hermite interpolation. particular pchip is a shape-preserving, "visually pleasing" interpolant that was introduced into Matlab fairly recently. It is based on an old Fortran program by Fritsch and Carlson  that is described by Kahaner, Moler, and Nash . Figure 3.4 shows how pchip interpolates our sample data. The key idea is to determine the slopes dk so that the function values do not overshoot the data values, at least locally. If k and k-1 have opposite signs or if either of them is zero, then xk is a discrete local minimum or maximum, so we set dk = 0. This is illustrated in the first half of Figure 3.5. The lower solid line is the piecewise linear interpolant. Its slopes on either side of the breakpoint in the center have opposite signs. Consequently, the dashed line has slope zero. The curved line is the shape-preserving interpolant, formed from two different cubics. The two cubics interpolate the center value and their derivatives are both zero there. But there is a jump in the second derivative at the breakpoint. If k and k-1 have the same sign and the two intervals have the same length, then dk is taken to be the harmonic mean of the two discrete slopes: 1 1 1 1 = + . dk 2 k-1 k In other words, at the breakpoint, the reciprocal slope of the Hermite interpolant is the average of the reciprocal slopes of the piecewise linear interpolant on either 10 Chapter 3. Interpolation side. This is shown in the other half of Figure 3.5. At the breakpoint, the reciprocal slope of the piecewise linear interpolant changes from 1 to 5. The reciprocal slope of the dashed line is 3, the average of 1 and 5. The shape-preserving interpolant is formed from the 2 cubics that interpolate the center value and that have slope equal to 1/3 there. Again, there is a jump in the second derivative at the breakpoint. Figure 3.5. Slopes for pchip. If k and k-1 have the same sign, but the two intervals have different lengths, then dk is a weighted harmonic mean, with weights determined by the lengths of the two intervals: w1 + w2 w1 w2 = + , dk k-1 k where w1 = 2hk + hk-1 , w2 = hk + 2hk-1 . This defines the pchip slopes at interior breakpoints, but the slopes d1 and dn at either end of the data interval are determined by a slightly different, one-sided analysis. The details are in pchiptx.m. 3.5 Cubic Spline Our other piecewise cubic interpolating function is a cubic spline. The term "spline" refers to an instrument used in drafting. It is a thin, flexible wooden or plastic tool that is passed through given data points and defines a smooth curve in between. The physical spline minimizes potential energy subject to the interpolation constraints. The corresponding mathematical spline must have a continuous second derivative and satisfy the same interpolation constraints. The breakpoints of a spline are also referred to as its knots. The world of splines extends far beyond the basic one-dimensional, cubic, interpolatory spline we are describing here. There are multidimensional, high-order, variable knot, and approximating splines. A valuable expository and reference text for both the mathematics and the software is A Practical Guide to Splines by Carl 3.5. Cubic Spline 11 de Boor . De Boor is also the author of the spline function and the Spline Toolbox for Matlab. Figure 3.6 shows how spline interpolates our sample data.
Spline interpolation 22 20 18 16 14 12 10 1 2 3 4 5 6 Figure 3.6. Cubic spline interpolation. The first derivative P (x) of our piecewise cubic function is defined by different formulas on either side of a knot xk . Both formulas yield the same value dk at the knots, so P (x) is continuous. On the kth subinterval, the second derivative is a linear function of s = x-xk : P (x) = If x = xk , s = 0 and (6h - 12s)k + (6s - 2h)dk+1 + (6s - 4h)dk . h2 P (xk +) = 6k - 2dk+1 - 4dk . hk The plus sign in xk + indicates that this is a one-sided derivative. If x = xk+1 , s = hk and -6k + 4dk+1 + 2dk P (xk+1 -) = . hk On the (k - 1)st interval, P (x) is given by a similar formula involving k-1 , dk , and dk-1 . At the knot xk , P (xk -) = -6k-1 + 4dk + 2dk-1 . hk-1 Requiring P (x) to be continuous at x = xk means that P (xk +) = P (xk -). 12 This leads to the condition Chapter 3. Interpolation hk dk-1 + 2(hk-1 + hk )dk + hk-1 dk+1 = 3(hk k-1 + hk-1 k ). If the knots are equally spaced, so that hk does not depend on k, this becomes dk-1 + 4dk + dk+1 = 3k-1 + 3k . Like our other interpolants, the slopes dk of a spline are closely related to the differences k . In the spline case, they are a kind of running average of the k 's. The preceding approach can be applied at each interior knot xk , k = 2, . . . , n- 1, to give n - 2 equations involving the n unknowns dk . As with pchip, a different approach must be used near the ends of the interval. One effective strategy is known as "not-a-knot." The idea is to use a single cubic on the first two subintervals, x1 x x3 , and on the last two subintervals, xn-2 x xn . In effect, x2 and xn-1 are not knots. If the knots are equally spaced with all hk = 1, this leads to d1 + 2d2 = and 2dn-1 + dn = 5 1 1 + 2 2 2 1 5 n-2 + n-1 . 2 2 The details if the spacing is not equal to one are in splinetx.m. With the two end conditions included, we have n linear equations in n unknowns: Ad = r. The vector of unknown slopes is d1 d2 d = . . . . dn . 2(hn-2 + hn-1 ) hn-2 hn-1 + hn-2 hn-2 .. The coefficient matrix A is tridiagonal: h h2 + h1 2 h1 h2 2(h1 + h2 ) h3 2(h2 + h3 ) A= .. . The right-hand side is h2 .. . hn-1 . r = 3 h r1 h2 1 + h1 2 h3 2 + h2 3 . . .
n-1 n-2 + hn-2 n-1 rn . 3.6. pchiptx, splinetx 13 The two values r1 and rn are associated with the end conditions. If the knots are equally spaced with all hk = 1, the coefficient matrix is quite simple: 1 2 1 4 1 1 4 1 . A= .. .. .. . . . 1 4 1 2 1 The right-hand side is . r = 3 + n-1 n-2 1 5 6 n-2 + 6 n-1 5 6 1 1 + 1 2 6 + 2 2 + 3 . . . In our textbook function, splinetx, the linear equations defining the slopes are solved with the tridisolve function introduced in Chapter 2, Linear Equations. In the spline functions distributed with Matlab and the Spline Toolbox, the slopes are computed by the Matlab backslash operator d = A\r Because most of the elements of A are zero, it is appropriate to store A in a sparse data structure. The backslash operator can then take advantage of the tridiagonal structure and solve the linear equations in time and storage proportional to n, the number of data points. Figure 3.7 compares the spline interpolant, s(x), with the pchip interpolant, p(x). The difference between the functions themselves is barely noticeable. The first derivative of spline, s (x), is smooth, while the first derivative of pchip, p (x), is continuous, but shows "kinks." The spline second derivative s (x) is continuous, while the pchip second derivative p (x) has jumps at the knots. Because both functions are piecewise cubics, their third derivatives, s (x) and p (x), are piecewise constant. The fact that s (x) takes on the same values in the first two intervals and the last two intervals reflects the "not-a-knot" spline end conditions. 3.6 pchiptx, splinetx The M-files pchiptx and splinetx are both based on piecewise cubic Hermite interpolation. On the kth interval, this is P (x) = 3hs2 - 2s3 h3 - 3hs2 + 2s3 yk+1 + yk h3 h3 2 2 s (s - h) s(s - h) + dk+1 + dk , 2 h h2 14
spline 21 18 15 12 2 5 0 -5 2 20 0 -20 30 0 -30 2 4 6 8 4 6 8 20 0 -20 30 0 -30 2 4 6 8 5 0 -5 2 21 18 15 12 2 Chapter 3. Interpolation
pchip 4 6 8 4 6 8 2 4 6 8 2 4 6 8 4 6 8 Figure 3.7. The spline and pchip interpolants, and their first three derivatives. where s = x - xk and h = hk . The two functions differ in the way they compute the slopes, dk . Once the slopes have been computed, the interpolant can be efficiently evaluated using the power form with the local variable s: P (x) = yk + sdk + s2 ck + s3 bk , where the coefficients of the quadratic and cubic terms are ck = 3k - 2dk - dk+1 , h dk - 2k + dk+1 bk = . h2 Here is the first portion of code for pchiptx. It calls an internal subfunction to compute the slopes, then computes the other coefficients, finds a vector of interval indices, and evaluates the interpolant. After the preamble, this part of the code for splinetx is the same. function v = pchiptx(x,y,u) %PCHIPTX Textbook piecewise cubic Hermite interpolation. % v = pchiptx(x,y,u) finds the shape-preserving piecewise % P(x), with P(x(j)) = y(j), and returns v(k) = P(u(k)). % % See PCHIP, SPLINETX. 3.6. pchiptx, splinetx 15 % First derivatives h = diff(x); delta = diff(y)./h; d = pchipslopes(h,delta); % Piecewise polynomial coefficients n = length(x); c = (3*delta - 2*d(1:n-1) - d(2:n))./h; b = (d(1:n-1) - 2*delta + d(2:n))./h.^2; % Find subinterval indices k so that x(k) <= u < x(k+1) k = ones(size(u)); for j = 2:n-1 k(x(j) <= u) = j; end % Evaluate interpolant s = u - x(k); v = y(k) + s.*(d(k) + s.*(c(k) + s.*b(k))); The code for computing the pchip slopes uses the weighted harmonic mean at interior breakpoints and a one-sided formula at the endpoints. function d = pchipslopes(h,delta) % PCHIPSLOPES Slopes for shape-preserving Hermite cubic % pchipslopes(h,delta) computes d(k) = P'(x(k)). % % % % % % Slopes at interior points delta = diff(y)./diff(x). d(k) = 0 if delta(k-1) and delta(k) have opposites signs or either is zero. d(k) = weighted harmonic mean of delta(k-1) and delta(k) if they have the same sign. n = length(h)+1; d = zeros(size(h)); k = find(sign(delta(1:n-2)).*sign(delta(2:n-1))>0)+1; w1 = 2*h(k)+h(k-1); w2 = h(k)+2*h(k-1); d(k) = (w1+w2)./(w1./delta(k-1) + w2./delta(k)); % Slopes at endpoints 16 Chapter 3. Interpolation d(1) = pchipend(h(1),h(2),delta(1),delta(2)); d(n) = pchipend(h(n-1),h(n-2),delta(n-1),delta(n-2)); function d = pchipend(h1,h2,del1,del2) % Noncentered, shape-preserving, three-point formula. d = ((2*h1+h2)*del1 - h1*del2)/(h1+h2); if sign(d) ~= sign(del1) d = 0; elseif (sign(del1)~=sign(del2))&(abs(d)>abs(3*del1)) d = 3*del1; end The splinetx M-file computes the slopes by setting up and solving a tridiagonal system of simultaneous linear equations. function d = splineslopes(h,delta); % SPLINESLOPES Slopes for cubic spline interpolation. % splineslopes(h,delta) computes d(k) = S'(x(k)). % Uses not-a-knot end conditions. % Diagonals of tridiagonal system n = length(h)+1; a = zeros(size(h)); b = a; c = a; r = a; a(1:n-2) = h(2:n-1); a(n-1) = h(n-2)+h(n-1); b(1) = h(2); b(2:n-1) = 2*(h(2:n-1)+h(1:n-2)); b(n) = h(n-2); c(1) = h(1)+h(2); c(2:n-1) = h(1:n-2); % Right-hand side r(1) = ((h(1)+2*c(1))*h(2)*delta(1)+ ... h(1)^2*delta(2))/c(1); r(2:n-1) = 3*(h(2:n-1).*delta(1:n-2)+ ... h(1:n-2).*delta(2:n-1)); r(n) = (h(n-1)^2*delta(n-2)+ ... (2*a(n-1)+h(n-1))*h(n-2)*delta(n-1))/a(n-1); % Solve tridiagonal linear system d = tridisolve(a,b,c,r); 3.7. interpgui
Piecewise linear interpolation Full degree polynomial interpolation 17 Shape-preserving Hermite interpolation Spline interpolation Figure 3.8. Four interpolants. 3.7 interpgui Figure 3.8 illustrates the tradeoff between smoothness and a somewhat subjective property that we might call local monotonicity or shape preservation. The piecewise linear interpolant is at one extreme. It has hardly any smoothness. It is continuous, but there are jumps in its first derivative. On the other hand, it preserves the local monotonicity of the data. It never overshoots the data and it is increasing, decreasing, or constant on the same intervals as the data. The full-degree polynomial interpolant is at the other extreme. It is infinitely differentiable. But it often fails to preserve shape, particularly near the ends of the interval. The pchip and spline interpolants are in between these two extremes. The spline is smoother than pchip. The spline has two continuous derivatives, while pchip has only one. A discontinuous second derivative implies discontinuous curvature. The human eye can detect large jumps in curvature in graphs and in mechanical parts made by numerically controlled machine tools. On the other hand, pchip is guaranteed to preserve shape, but the spline might not. The M-file interpgui allows you to experiment with the four interpolants discussed in this chapter: piecewise linear interpolant, full-degree interpolating polynomial, 18
Interpolation Chapter 3. Interpolation 24 22 20 18 16 14 12 10 8 0 1 2 3 4 5 6 Figure 3.9. interpgui. piecewise cubic spline, shape-preserving piecewise cubic. The program can be initialized in several different ways: With no arguments, interpgui starts with 8 zeros. With a scalar argument, interpgui(n) starts with n zeros. With one vector argument, interpgui(y) starts with equally spaced x's. With two arguments, interpgui(x,y) starts with a plot of y versus x. After initialization, the interpolation points can be varied with the mouse. If x has been specified, it remains fixed. Figure 3.9 is the initial plot generated by our example data. Exercises
3.1. Reproduce Figure 3.8, with four subplots showing the four interpolants discussed in this chapter. 3.2. Tom and Ben are twin boys born on October 27, 2001. Here is a table of their weights, in pounds and ounces, over their first few months. Exercises % W = [10 11 12 12 01 01 03 Date 27 2001 19 2001 03 2001 20 2001 09 2002 23 2002 06 2002 Tom 5 10 7 4 8 12 10 14 12 13 14 8 16 10 Ben 4 8 5 11 6 4 8 7 10 3 12 0 13 10]; 19 You can use datenum to convert the date in the first three columns to a serial date number measuring time in days. t = datenum(W(:,[3 1 2])); Make a plot of their weights versus time, with circles at the data points and the pchip interpolating curve in between. Use datetick to relabel the time axis. Include a title and a legend. The result should look something like Figure 3.10.
Twins' weights 18 16 14 12 Pounds 10 8 6 4 Nov01 Dec01 Jan02 Feb02 Tom Ben Mar02 Figure 3.10. Twins' weights. 3.3. (a) Interpolate these data by each of the four interpolants discussed in this chapter: piecelin, polyinterp, splinetx, and pchiptx. Plot the results for -1 x 1. 20 x -1.00 -0.96 -0.65 0.10 0.40 1.00 y -1.0000 -0.1512 0.3860 0.4802 0.8838 1.0000 Chapter 3. Interpolation (b) What are the values of each of the four interpolants at x = -0.3? Which of these values do you prefer? Why? (c) The data were actually generated from a low-degree polynomial with integer coefficients. What is that polynomial? 3.4. Make a plot of your hand. Start with figure('position',get(0,'screensize')) axes('position',[0 0 1 1]) [x,y] = ginput; Place your hand on the computer screen. Use the mouse to select a few dozen points outlining your hand. Terminate the ginput with a carriage return. You might find it easier to trace your hand on a piece of paper and then put the paper on the computer screen. You should be able to see the ginput cursor through the paper. (Save these data. We will refer to them in other exercises later in this book.) Figure 3.11. A hand. Exercises 21 Now think of x and y as two functions of an independent variable that goes from one to the number of points you collected. You can interpolate both functions on a finer grid and plot the result with n = length(x); s = (1:n)'; t = (1:.05:n)'; u = splinetx(s,x,t); v = splinetx(s,y,t); clf reset plot(x,y,'.',u,v,'-'); Do the same thing with pchiptx. Which do you prefer? Figure 3.11 is the plot of my hand. Can you tell if it was done with splinetx or pchiptx? 3.5. The previous exercise uses the data index number as the independent variable for two-dimensional parametric interpolation. This exercise uses, instead, the angle from polar coordinates. In order to do this, the data must be centered so that they lie on a curve that is starlike with respect to the origin, that is, every ray emanating from the origin meets the data only once. This means that you must be able to find values x0 and y0 so that the Matlab statements x = x - x0 y = y - y0 theta = atan2(y,x) r = sqrt(x.^2 + y.^2) plot(theta,r) produce a set of points that can be interpolated with a single-valued function, r = r(). For the data obtained by sampling the outline of your hand, the point (x0 , y0 ) is located near the base of your palm. See the small circle in Figure 3.11. Furthermore, in order to use splinetx and pchiptx, it is also necessary to order the data so that theta is monotonically increasing. Choose a subsampling increment, delta, and let t = (theta(1):delta:theta(end))'; p = pchiptx(theta,r,t); s = splinetx(theta,r,t); Examine two plots: plot(theta,r,'o',t,[p s],'-') and plot(x,y,'o',p.*cos(t),p.*sin(t),'-',... s.*cos(t),s.*sin(t),'-') Compare this approach with the one used in the previous exercise. Which do you prefer? Why? 22 Chapter 3. Interpolation 3.6. This exercise requires the Symbolic Toolbox. (a) What does vandal(n) compute and how does it compute it? (b) Under what conditions on x is the matrix vander(x) nonsingular? 3.7. Prove that the interpolating polynomial is unique. That is, if P (x) and Q(x) are two polynomials with degree less than n that agree at n distinct points, then they agree at all points. 3.8. Give convincing arguments that each of the following descriptions defines the same polynomial, the Chebyshev polynomial of degree five, T5 (x). Your arguments can involve analytic proofs, symbolic computation, numeric computation, or all three. Two of the representations involve the golden ratio, 1+ 5 = . 2 (a) Power form basis: (b) Relation to trigonometric functions: (c) Horner representation: (d) Lagrange form: T5 (x) = ((((16x + 0)x - 20)x + 0)x + 5)x + 0. x1 , x6 = 1, x2 , x5 = /2, x3 , x4 = ( - 1)/2, T5 (x) = 16x5 - 20x3 + 5x. T5 (x) = cos (5 cos-1 x). (e) Factored representation: yk = (-1)k , k = 1, . . . , 6, x - xj yk . T5 (x) = xk - xj
k j=k z1 , z5 = z2 , z4 = z3 = 0, (2 + )/4, (3 - )/4, (x - zk ). T5 (x) = 16 (f) Three-term recurrence: 5 1 T0 (x) = 1, T1 (x) = x, Tn (x) = 2xTn-1 (x) - Tn-2 (x) for n = 2, . . . , 5. Exercises 23 3.9. The M-file rungeinterp.m provides an experiment with a famous polynomial interpolation problem due to Carl Runge. Let f (x) = 1 , 1 + 25x2 and let Pn (x) denote the polynomial of degree n - 1 that interpolates f (x) at n equally spaced points on the interval -1 x 1. Runge asked whether, as n increases, Pn (x) converges to f (x). The answer is yes for some x, but no for others. (a) For what x does Pn (x) f (x) as n ? (b) Change the distribution of the interpolation points so that they are not equally spaced. How does this affect convergence? Can you find a distribution so that Pn (x) f (x) for all x in the interval? 3.10. We skipped from piecewise linear to piecewise cubic interpolation. How far can you get with the development of piecewise quadratic interpolation? 3.11. Modify splinetx and pchiptx so that, if called with two output arguments, they produce both the value of the interpolant and its first derivative. That is, [v,vprime] = pchiptx(x,y,u) and [v,vprime] = splinetx(x,y,u) compute P (u) and P (u). 3.12. Modify splinetx and pchiptx so that, if called with only two input arguments, they produce PP, the piecewise polynomial structure produced by the standard Matlab functions spline and pchip and used by ppval. 3.13. (a) Create two functions perpchip and perspline by modifying pchiptx and splinetx to replace the one-sided and not-a-knot end conditions with periodic boundary conditions. This requires that the given data have yn = y1 and that the resulting interpolant be periodic. In other words, for all x, P (x + ) = P (x), where = xn - x1 . The algorithms for both pchip and spline involve calculations with yk , hk , and k to produce slopes dk . With the periodicity assumption, all of these quantities become periodic functions, with period n - 1, of the subscript k. In other words, for all k, 24 yk hk k dk = yk+n-1 , = hk+n-1 , = k+n-1 , = dk+n-1 . Chapter 3. Interpolation This makes it possible to use the same calculations at the endpoints that are used at the interior points in the nonperiodic situation. The special case code for the end conditions can be eliminated and the resulting M-files are actually much simpler. For example, the slopes dk for pchip with equally spaced points are given by dk = 0 if sign(k-1 ) = sign(k ), 1 1 1 1 = + if sign(k-1 ) = sign(k ). dk 2 k-1 k With periodicity, these formulas can also be used at the endpoints where k = 1 and k = n because 0 = n-1 and n = 1 . For spline, the slopes satisfy a system of simultaneous linear equations for k = 2, . . . , n - 1: hk dk-1 + 2(hk-1 + hk )dk + hk-1 dk+1 = 3(hk k-1 + hk-1 k ). With periodicity, this becomes h1 dn-1 + 2(hn-1 + h1 )d1 + hn-1 d2 = 3(h1 n-1 + hn-1 1 ) at k = 1 and hn dn-1 + 2(hn-1 + h1 )dn + hn-1 d2 = 3(h1 n-1 + hn-1 1 ) at k = n. The resulting matrix has two nonzero elements outside the tridiagonal structure. The additional nonzero elements in the first and last rows are A1,n-1 = h1 and An,2 = hn-1 . (b) Demonstrate that your new functions work correctly on x = 0:pi/4:2*pi; y = cos(x); u = 0:pi/50:2*pi; v = your_function(x,y,u); plot(x,y,'o',u,v,'-') (c) Once you have perchip and perspline, you can use the NCM M-file interp2dgui to investigate closed-curve interpolation in two dimensions. You should find that the periodic boundary conditions do a better job of reproducing symmetries of closed curves in the plane. Exercises 3.14. (a) Modify splinetx so that it forms the full tridiagonal matrix A = diag(a,-1) + diag(b,0) + diag(c,1) 25 and uses backslash to compute the slopes. (b) Monitor condest(A) as the spline knots are varied with interpgui. What happens if two of the knots approach each other? Find a data set that makes condest(A) large. 3.15. Modify pchiptx so that it uses a weighted average of the slopes instead of the weighted harmonic mean. 3.16. (a) Consider x = -1:1/3:1 interpgui(1-x.^2) Which, if any, of the four interpolants linear, spline, pchip, and polynomial are the same? Why? (b) Same questions for interpgui(1-x.^4) 3.17. Why does interpgui(4) show only three graphs, not four, no matter where you move the points? 3.18. (a) If you want to interpolate census data on the interval 1900 t 2000 with a polynomial, P (t) = c1 t10 + c2 t9 + + c10 t + c11 , you might be tempted to use the Vandermonde matrix generated by t = 1900:10:2000 V = vander(t) Why is this a really bad idea? (b) Investigate centering and scaling the independent variable. Plot some data, pull down the Tools menu on the figure window, select Basic Fitting, and find the check box about centering and scaling. What does this check box do? (c) Replace the variable t with t- . ~ This leads to a modified polynomial P (s). How are its coefficients related to those of P (t)? What happens to the Vandermonde matrix? What values of and lead to a reasonably well conditioned Vandermonde matrix? One possibility is s= mu = mean(t) sigma = std(t) but are there better values? 26 Chapter 3. Interpolation Bibliography
 C. de Boor, A Practical Guide to Splines, Springer-Verlag, New York, 1978.  F. N. Fritsch and R. E. Carlson, Monotone Piecewise Cubic Interpolation, SIAM Journal on Numerical Analysis, 17 (1980), pp. 238-246.  D. Kahaner, C. Moler, and S. Nash, Numerical Methods and Software, PrenticeHall, Englewood Cliffs, NJ, 1989. 27 ...
View Full Document
- Fall '11