spline

pchip

Contents

Data

x = 1:6 y = [16 18 21 17 15 12]

x = 1 2 3 4 5 6 y = 16 18 21 17 15 12

set(0, 'defaultlinelinewidth' ,2) clf plot(x,y, '-o' ) axis([0 7 7.5 25.5]) title( 'plip' )

plip

plip

plip

The PCHIP Family

x = 1:1/64:2; s = 16 + 2*(x-1) + (49/18)*(x-1).^2.*(x-2) - (89/18)*(x-1).*(x-2).^2; p = 16 + 2*(x-1) + (2/5)*(x-1).^2.*(x-2) - (1/2)*(x-1).*(x-2).^2; clf axis([0 3 15 19]) box on line(x,s, 'color' ,[0 3/4 3/4]) line(x,p, 'color' ,[3/4 0 3/4]) line(x(1),s(1), 'marker' , 'o' , 'color' ,[0 0 3/4]) line(x(end),s(end), 'marker' , 'o' , 'color' ,[0 0 3/4])

spline

interpgui

x = 1:6; y = [16 18 21 17 15 12]; interpgui(x,y,3)

sppchip

sppchip

pchip

spline

pchip

pchip

pchip

spline

pchip

spline

interpgui(x,y,3:4)

spline vs. pchip

spline

pchip

spline

pchip

spline

pchip

splinevspchip

Locality

pchip

pchip

pchip

spline

spline

pchip

spline

B-splines

x = 1:8; y = zeros(1,8); y(4) = 1; interpgui(x,y,3:4)

interp1

interp1

method

'linear'

'spline'

'pchip'

'cubic'

'pchip'

pchip

spline

'v5cubic'

interpgui

'v5cubic'

interp1

spline

pchip

x = 1:6; y = [16 18 21 11 15 12]; interpgui_with_v5cubic(x,y,3:5)

Resources

doc curvefit

MATLAB has two different functions for piecewise cubic interpolation,and. Why are there two? How do they compare?Here is the data that I will use in this post.Here is a plot of the data.With line type '-o', the MATLAB plot command plots six 'o's at the six data points and draws straight lines between the points. So I added the titlebecause this is a graph of the. There is a different linear function between each pair of points. Since we want the function to go through the data points, that isthe data, and since two points determine a line, thefunction is unique.A PCHIP, a, is any piecewise cubic polynomial that interpolates the given data, AND has specified derivatives at the interpolation points. Just as two points determine a linear function, two points and two given slopes determine a cubic. The data points are known as "knots". We have the-values at the knots, so in order to get a particular PCHIP, we have to somehow specify the values of the derivative,, at the knots. Consider these two cubic polynomials in $x$ on the interval $1 \le x \le 2$ . These functions are formed by adding cubic terms that vanish at the end points to the linear interpolatant. I'll tell you later where the coefficients of the cubics come from. $$ s(x) = 16 + 2(x-1) + \textstyle{\frac{49}{18}}(x-1)^2(x-2) - \textstyle{\frac{89}{18}}(x-1)(x-2)^2 $$ $$ p(x) = 16 + 2(x-1) + \textstyle{\frac{2}{5}}(x-1)^2(x-2) - \textstyle{\frac{1}{2}}(x-1)(x-2)^2 $$ These functions interpolate the same values at the ends. $$ s(1) = 16, \ \ \ s(2) = 18 $$ $$ p(1) = 16, \ \ \ p(2) = 18 $$ But they have different first derivatives at the ends. In particular, $s'(1)$ is negative and $p'(1)$ is positive. $$ s'(1) = - \textstyle{\frac{53}{18}}, \ s'(2) = \textstyle{\frac{85}{18}} $$ $$ p'(1) = \textstyle{\frac{3}{2}}, \ \ \ p'(2) = \textstyle{\frac{12}{5}} $$ Here's a plot of these two cubic polynomials. The magenta cubic, which is $p(x)$, just climbs steadily from its initial value to its final value. On the other hand, the cyan cubic, which is $s(x)$, starts off heading in the wrong direction, then has to hurry to catch up.If we piece together enough cubics like these to produce a piecewise cubic that interpolates many data points, we have a PCHIP. We could even mix colors and still have a PCHIP. Clearly, we have to be specific when it comes to specifying the slopes. One possibility that might occur to you briefly is to use the slopes of the lines connecting the end points of each segment. But this choice just produces zeros for the coefficients of the cubics and leads back to the piecewise linear interpolant. After all, a linear function is a degenerate cubic. This illustrates the fact that the PCHIP family includes many functions.By far, the most famous member of the PCHIP family is the piecewise cubic spline. All PCHIPs are continuous and have a continuous first derivative. A spline is a PCHIP that is exceptionally smooth, in the sense that its second derivative, and consequently its curvature, also varies continuously. The function derives its name from the flexible wood or plastic strip used to draw smooth curves. Starting about 50 years ago, Carl de Boor developed much of the basic theory of splines. He wrote a widely adopted package of Fortran software, and a widely cited book, for computations involving splines. Later, Carl authored the MATLAB Spline Toolbox. Today, the Spline Toolbox is part of the Curve Fitting Toolbox. When Carl began the development of splines, he was with General Motors Research in Michigan. GM was just starting to use numerically controlled machine tools. It is essential that automobile parts have smooth edges and surfaces. If the hood of a car, say, does not have continuously varying curvature, you can see wrinkles in the reflections in the show room. In the automobile industry, a discontinuous second derivative is known as a "dent". The requirement of a continuous second derivative leads to a set of simultaneous linear equations relating the slopes at the interior knots. The two end points need special treatment, and the default treatment has changed over the years. We now choose the coefficients so that thederivative does not have a jump at the first and last interior knots. Single cubic pieces interpolate the first three, and the last three, data points. This is known as the "not-a-knot" condition. It adds two more equations to set of equations at the interior points. If there areknots, this gives a well-conditioned, almost symmetric, tridiagonal $n$ -by- $n$ linear system to solve for the slopes. The system can be solved by the sparse backslash operator in MATLAB, or by a custom, non-pivoting tridiagonal solver. (Other end conditions for splines are available in the Curve Fitting Toolbox.) As you probably realized, the cyan function $s(x)$ introduced above, is one piece of the spline interpolating our sample data. Here is a graph of the entire function, produced byfrom NCM,I just made up that name,. It stands for. The actual name of the MATLAB function is just. This function is not as smooth as. There may well be jumps in the second derivative. Instead, the function is designed so that it never locallythe data. The slope at each interior point is taken to be a weighted harmonic mean of the slopes of the piecewise linear interpolant. One-sided slope conditions are imposed at the two end points. Theslopes can be computed without solving a linear system.was originally developed by Fred Fritsch and his colleagues at Lawrence Livermore Laboratory around 1980. They described it as "visually pleasing". Dave Kahaner, Steve Nash and I included some of Fred's Fortran subroutines in our 1989 book,. We madepart of MATLAB in the early '90s. Here is a comparison ofandon our data. In this case the spline overshoot on the first subinterval is caused by the not-a-knot end condition. But with more data points, or rapidly varying data points, interior overshoots are possible withHere are eight subplots comparingandon a slightly larger data set. The first two plots show the functions $s(x)$ and $p(x)$. The difference between the functions on the interior intervals is barely noticeable. The next two plots show the first derivatives. You can see that the first derivative of, $s'(x)$, is smooth, while the first derivative of, $p'(x)$, is continuous, but shows "kinks". The third pair of plots are the second derivatives. Thesecond derivative $s''(x)$ is continuous, while thesecond derivative $p''(x)$ has jumps at the knots. The final pair are the third derivatives. 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.is local. The behavior ofon a particular subinterval is determined by only four points, the two data points on either side of that interval.is unaware of the data farther away.is global. The behavior ofon a particular subinterval is determined by all of the data, although the sensitivity to data far away is less than to nearby data. Both behaviors have their advantages and disadvantages. Here is the response to a unit impulse. You can see that the support ofis confined to the two intervals surrounding the impulse, while the support ofextends over the entire domain. (There is an elegant set of basis functions for cubic splines known asthat do have compact support.)Thefunction in MATLAB, has severaloptions. The, andoptions are the same interpolants we have been discussing here. We decided years ago to make theoption the same asbecause we thought the monotonicity property ofwas generally more desirable than the smoothness property of. Theoption is yet another member of the PCHIP family, which has been retained for compatibility with version 5 of MATLAB. It requires the's to be equally spaced. The slope of the v5 cubic at point $x_n$ is $(y_{n+1} - y_{n-1})/2$. The resulting piecewise cubic does not have a continuous second derivative and it does not always preserve shape. Because the abscissa are equally spaced, the v5 cubic can be evaluated quickly by a convolution operation. Here is our example data, modified slightly to exaggerate behavior, andmodified to include theoption of. The v5 cubic is the black curve betweenandA extensive collection of tools for curve and surface fitting, by splines and many other functions, is available in the Curve Fitting Toolbox."NCM",, has more mathematical details. NCM is available online . Here is the interpolation chapter . Here is interpgui . SIAM publishes a print edition . Here are the script splinevspchip.m and the modified version of interpgui interpgui_with_v5cubic.m that I used in this post.