{"id":196,"date":"2012-07-16T13:39:15","date_gmt":"2012-07-16T18:39:15","guid":{"rendered":"https:\/\/blogs.mathworks.com\/cleve\/?p=196"},"modified":"2017-02-24T20:10:29","modified_gmt":"2017-02-25T01:10:29","slug":"splines-and-pchips","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/cleve\/2012\/07\/16\/splines-and-pchips\/","title":{"rendered":"Splines and Pchips"},"content":{"rendered":"&nbsp;\r\n<div class=\"content\"><!--introduction-->MATLAB has two different functions for piecewise cubic interpolation, <tt>spline<\/tt> and <tt>pchip<\/tt>. Why are there two? How do they compare?\r\n\r\n<!--\/introduction-->\r\n<h3>Contents<\/h3>\r\n<div>\r\n<ul>\r\n \t<li><a href=\"#59fac6b0-93a3-4c9b-a3a2-bcdfb3253b0b\">Data<\/a><\/li>\r\n \t<li><a href=\"#6eb05117-4ff1-435f-81d4-c61ee1f715f8\">plip<\/a><\/li>\r\n \t<li><a href=\"#ee54c20b-0ecc-4ac9-8986-8a0774a1763f\">The PCHIP Family<\/a><\/li>\r\n \t<li><a href=\"#6bd949e0-f5f7-45da-bdbb-70ffb9959cfc\">spline<\/a><\/li>\r\n \t<li><a href=\"#e832779f-f543-4521-8a9d-9da4347692e9\">sppchip<\/a><\/li>\r\n \t<li><a href=\"#98ccb1df-b614-41d4-b1b5-e090a87e0d46\">spline vs. pchip<\/a><\/li>\r\n \t<li><a href=\"#6f1d49c3-6007-44d9-8d8d-e03477c67673\">Locality<\/a><\/li>\r\n \t<li><a href=\"#f3bace9d-40da-4119-9814-95d285141d37\">interp1<\/a><\/li>\r\n \t<li><a href=\"#59fe8852-238e-4351-9285-8f9a17018a89\">Resources<\/a><\/li>\r\n<\/ul>\r\n<\/div>\r\n<h4>Data<a name=\"59fac6b0-93a3-4c9b-a3a2-bcdfb3253b0b\"><\/a><\/h4>\r\nHere is the data that I will use in this post.\r\n<pre class=\"codeinput\">x = 1:6\r\ny = [16 18 21 17 15 12]\r\n<\/pre>\r\n<pre class=\"codeoutput\">x =\r\n\r\n     1     2     3     4     5     6\r\n\r\n\r\ny =\r\n\r\n    16    18    21    17    15    12\r\n\r\n<\/pre>\r\nHere is a plot of the data.\r\n<pre class=\"codeinput\">set(0,<span class=\"string\">'defaultlinelinewidth'<\/span>,2)\r\nclf\r\nplot(x,y,<span class=\"string\">'-o'<\/span>)\r\naxis([0 7 7.5 25.5])\r\ntitle(<span class=\"string\">'plip'<\/span>)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/pchips_01.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>plip<a name=\"6eb05117-4ff1-435f-81d4-c61ee1f715f8\"><\/a><\/h4>\r\nWith 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 title <tt>plip<\/tt> because this is a graph of the <i>piecewise linear interpolating polynomial<\/i>. There is a different linear function between each pair of points. Since we want the function to go through the data points, that is <i>interpolate<\/i> the data, and since two points determine a line, the <tt>plip<\/tt> function is unique.\r\n<h4>The PCHIP Family<a name=\"ee54c20b-0ecc-4ac9-8986-8a0774a1763f\"><\/a><\/h4>\r\nA PCHIP, a <i>Piecewise Cubic Hermite Interpolating Polynomial<\/i>, 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 <i>y<\/i>-values at the knots, so in order to get a particular PCHIP, we have to somehow specify the values of the derivative, <i>y'<\/i>, at the knots.\r\n\r\nConsider 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.\r\n\r\n$$ s(x) = 16 + 2(x-1) + \\textstyle{\\frac{49}{18}}(x-1)^2(x-2)\r\n- \\textstyle{\\frac{89}{18}}(x-1)(x-2)^2 $$\r\n\r\n$$ p(x) = 16 + 2(x-1) + \\textstyle{\\frac{2}{5}}(x-1)^2(x-2)\r\n- \\textstyle{\\frac{1}{2}}(x-1)(x-2)^2 $$\r\n\r\nThese functions interpolate the same values at the ends.\r\n\r\n$$ s(1) = 16, \\ \\ \\ s(2) = 18 $$\r\n\r\n$$ p(1) = 16, \\ \\ \\ p(2) = 18 $$\r\n\r\nBut they have different first derivatives at the ends. In particular, $s'(1)$ is negative and $p'(1)$ is positive.\r\n\r\n$$ s'(1) = - \\textstyle{\\frac{53}{18}}, \\ s'(2) = \\textstyle{\\frac{85}{18}} $$\r\n\r\n$$ p'(1) = \\textstyle{\\frac{3}{2}}, \\ \\ \\ p'(2) = \\textstyle{\\frac{12}{5}} $$\r\n\r\nHere'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.\r\n<pre class=\"codeinput\">x = 1:1\/64:2;\r\ns = 16 + 2*(x-1) + (49\/18)*(x-1).^2.*(x-2) - (89\/18)*(x-1).*(x-2).^2;\r\np = 16 + 2*(x-1) + (2\/5)*(x-1).^2.*(x-2) - (1\/2)*(x-1).*(x-2).^2;\r\n\r\nclf\r\naxis([0 3 15 19])\r\nbox <span class=\"string\">on<\/span>\r\nline(x,s,<span class=\"string\">'color'<\/span>,[0 3\/4 3\/4])\r\nline(x,p,<span class=\"string\">'color'<\/span>,[3\/4 0 3\/4])\r\nline(x(1),s(1),<span class=\"string\">'marker'<\/span>,<span class=\"string\">'o'<\/span>,<span class=\"string\">'color'<\/span>,[0 0 3\/4])\r\nline(x(end),s(end),<span class=\"string\">'marker'<\/span>,<span class=\"string\">'o'<\/span>,<span class=\"string\">'color'<\/span>,[0 0 3\/4])\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/pchips_02.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n\r\nIf 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.\r\n\r\nOne 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.\r\n<h4>spline<a name=\"6bd949e0-f5f7-45da-bdbb-70ffb9959cfc\"><\/a><\/h4>\r\nBy 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.\r\n\r\nStarting 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.\r\n\r\nWhen 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\".\r\n\r\nThe 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 the <i>third<\/i> derivative 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 are <i>n<\/i> knots, 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.)\r\n\r\nAs 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 by <tt>interpgui<\/tt> from NCM, <i>Numerical Computing with MATLAB<\/i>.\r\n<pre class=\"codeinput\">x = 1:6;\r\ny = [16 18 21 17 15 12];\r\ninterpgui(x,y,3)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/pchips_03.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>sppchip<a name=\"e832779f-f543-4521-8a9d-9da4347692e9\"><\/a><\/h4>\r\nI just made up that name, <tt>sppchip<\/tt>. It stands for <i>shape preserving piecewise cubic Hermite interpolating polynomial<\/i>. The actual name of the MATLAB function is just <tt>pchip<\/tt>. This function is not as smooth as <tt>spline<\/tt>. There may well be jumps in the second derivative. Instead, the function is designed so that it never locally <i>overshoots<\/i> the 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. The <tt>pchip<\/tt> slopes can be computed without solving a linear system.\r\n\r\n<tt>pchip<\/tt> 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, <i>Numerical Methods and Software<\/i>. We made <tt>pchip<\/tt> part of MATLAB in the early '90s.\r\n\r\nHere is a comparison of <tt>spline<\/tt> and <tt>pchip<\/tt> on 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 with <tt>spline<\/tt>.\r\n<pre class=\"codeinput\">interpgui(x,y,3:4)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/pchips_04.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>spline vs. pchip<a name=\"98ccb1df-b614-41d4-b1b5-e090a87e0d46\"><\/a><\/h4>\r\nHere are eight subplots comparing <tt>spline<\/tt> and <tt>pchip<\/tt> on 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 <tt>spline<\/tt>, $s'(x)$, is smooth, while the first derivative of <tt>pchip<\/tt>, $p'(x)$, is continuous, but shows \"kinks\". The third pair of plots are the second derivatives. The <tt>spline<\/tt> second derivative $s''(x)$ is continuous, while the <tt>pchip<\/tt> second 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.\r\n<pre class=\"codeinput\">splinevspchip\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/pchips_05.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>Locality<a name=\"6f1d49c3-6007-44d9-8d8d-e03477c67673\"><\/a><\/h4>\r\n<tt>pchip<\/tt> is local. The behavior of <tt>pchip<\/tt> on a particular subinterval is determined by only four points, the two data points on either side of that interval. <tt>pchip<\/tt> is unaware of the data farther away. <tt>spline<\/tt> is global. The behavior of <tt>spline<\/tt> on 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.\r\n\r\nHere is the response to a unit impulse. You can see that the support of <tt>pchip<\/tt> is confined to the two intervals surrounding the impulse, while the support of <tt>spline<\/tt> extends over the entire domain. (There is an elegant set of basis functions for cubic splines known as <tt>B-splines<\/tt> that do have compact support.)\r\n<pre class=\"codeinput\">x = 1:8;\r\ny = zeros(1,8);\r\ny(4) = 1;\r\ninterpgui(x,y,3:4)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/pchips_06.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>interp1<a name=\"f3bace9d-40da-4119-9814-95d285141d37\"><\/a><\/h4>\r\nThe <tt>interp1<\/tt> function in MATLAB, has several <tt>method<\/tt> options. The <tt>'linear'<\/tt>, <tt>'spline'<\/tt>, and <tt>'pchip'<\/tt> options are the same interpolants we have been discussing here. We decided years ago to make the <tt>'cubic'<\/tt> option the same as <tt>'pchip'<\/tt> because we thought the monotonicity property of <tt>pchip<\/tt> was generally more desirable than the smoothness property of <tt>spline<\/tt>.\r\n\r\nThe <tt>'v5cubic'<\/tt> option is yet another member of the PCHIP family, which has been retained for compatibility with version 5 of MATLAB. It requires the <i>x<\/i>'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.\r\n\r\nHere is our example data, modified slightly to exaggerate behavior, and <tt>interpgui<\/tt> modified to include the <tt>'v5cubic'<\/tt> option of <tt>interp1<\/tt>. The v5 cubic is the black curve between <tt>spline<\/tt> and <tt>pchip<\/tt>.\r\n<pre class=\"codeinput\">x = 1:6;\r\ny = [16 18 21 11 15 12];\r\ninterpgui_with_v5cubic(x,y,3:5)\r\n<\/pre>\r\n<img decoding=\"async\" src=\"https:\/\/blogs.mathworks.com\/images\/cleve\/pchips_07.png\" alt=\"\" hspace=\"5\" vspace=\"5\" \/>\r\n<h4>Resources<a name=\"59fe8852-238e-4351-9285-8f9a17018a89\"><\/a><\/h4>\r\nA extensive collection of tools for curve and surface fitting, by splines and many other functions, is available in the Curve Fitting Toolbox.\r\n<pre class=\"codeinput\">doc <span class=\"string\">curvefit<\/span>\r\n<\/pre>\r\n\"NCM\", <i>Numerical Computing with MATLAB<\/i>, has more mathematical details. NCM is available <a href=\"https:\/\/www.mathworks.com\/moler\/index_ncm.html\">online<\/a>. Here is the <a href=\"https:\/\/www.mathworks.com\/content\/dam\/mathworks\/mathworks-dot-com\/moler\/interp.pdf\">interpolation chapter<\/a>. Here is <a href=\"https:\/\/www.mathworks.com\/content\/dam\/mathworks\/mathworks-dot-com\/moler\/ncm\/interpgui.m\">interpgui<\/a>. SIAM publishes a <a href=\"https:\/\/bookstore.siam.org\/OT87\">print edition<\/a>.\r\n\r\nHere are the script <a href=\"https:\/\/blogs.mathworks.com\/images\/cleve\/splinevspchip.m\">splinevspchip.m<\/a> and the modified version of interpgui <a href=\"https:\/\/blogs.mathworks.com\/images\/cleve\/interpgui_with_v5cubic.m\">interpgui_with_v5cubic.m<\/a> that I used in this post.\r\n\r\n<script>\/\/ <![CDATA[\r\nfunction grabCode_1add2164eaeb45548380e6801ab008b6() {\r\n        \/\/ Remember the title so we can use it in the new page\r\n        title = document.title;\r\n\r\n        \/\/ Break up these strings so that their presence\r\n        \/\/ in the Javascript doesn't mess up the search for\r\n        \/\/ the MATLAB code.\r\n        t1='1add2164eaeb45548380e6801ab008b6 ' + '##### ' + 'SOURCE BEGIN' + ' #####';\r\n        t2='##### ' + 'SOURCE END' + ' #####' + ' 1add2164eaeb45548380e6801ab008b6';\r\n    \r\n        b=document.getElementsByTagName('body')[0];\r\n        i1=b.innerHTML.indexOf(t1)+t1.length;\r\n        i2=b.innerHTML.indexOf(t2);\r\n \r\n        code_string = b.innerHTML.substring(i1, i2);\r\n        code_string = code_string.replace(\/REPLACE_WITH_DASH_DASH\/g,'--');\r\n\r\n        \/\/ Use \/x3C\/g instead of the less-than character to avoid errors \r\n        \/\/ in the XML parser.\r\n        \/\/ Use '\\x26#60;' instead of '<' so that the XML parser\r\n        \/\/ doesn't go ahead and substitute the less-than character. \r\n        code_string = code_string.replace(\/\\x3C\/g, '\\x26#60;');\r\n\r\n        copyright = 'Copyright 2012 The MathWorks, Inc.';\r\n\r\n        w = window.open();\r\n        d = w.document;\r\n        d.write('\r\n\r\n<pre>\\n');\r\n        d.write(code_string);\r\n\r\n        \/\/ Add copyright line at the bottom if specified.\r\n        if (copyright.length > 0) {\r\n            d.writeln('');\r\n            d.writeln('%%');\r\n            if (copyright.length > 0) {\r\n                d.writeln('% _' + copyright + '_');\r\n            }\r\n        }\r\n\r\n        d.write('<\/pre>\r\n\r\n\r\n\\n');\r\n\r\n        d.title = title + ' (MATLAB code)';\r\n        d.close();\r\n    }\r\n\/\/ ]]><\/script>\r\n<p style=\"text-align: right; font-size: xx-small; font-weight: lighter; font-style: italic; color: gray;\">\r\n<a><span style=\"font-size: x-small; font-style: italic;\">Get\r\nthe MATLAB code<noscript>(requires JavaScript)<\/noscript><\/span><\/a>\r\n\r\nPublished with MATLAB\u00ae 7.14<\/p>\r\n<p class=\"footer\">\r\nPublished with MATLAB\u00ae 7.14<\/p>\r\n\r\n<\/div>\r\n<!--\r\n1add2164eaeb45548380e6801ab008b6 ##### SOURCE BEGIN #####\r\n%% Splines and Pchips\r\n% MATLAB has two different functions for piecewise cubic interpolation,\r\n% |spline| and |pchip|.  Why are there two? How do they compare?\r\n\r\n%% Data\r\n% Here is the data that I will use in this post.\r\n\r\nx = 1:6\r\ny = [16 18 21 17 15 12]\r\n\r\n%%\r\n% Here is a plot of the data.\r\n\r\nset(0,'defaultlinelinewidth',2)\r\nclf\r\nplot(x,y,'-o')\r\naxis([0 7 7.5 25.5])\r\ntitle('plip')\r\n\r\n%% plip\r\n% With line type '-o', the MATLAB plot command plots six 'o's at the six\r\n% data points and draws straight lines between the points.\r\n% So I added the title |plip| because this is a graph of the\r\n% _piecewise linear interpolating polynomial_.\r\n% There is a different linear function between each pair of points.\r\n% Since we want the function to go through the data points,\r\n% that is _interpolate_ the data, and since two points determine a line,\r\n% the |plip| function is unique.\r\n\r\n%% The PCHIP Family\r\n% A PCHIP, a _Piecewise Cubic Hermite Interpolating Polynomial_,\r\n% is any piecewise cubic polynomial that interpolates the given data,\r\n% AND has specified derivatives at the interpolation points.\r\n% Just as two points determine a linear function, two points and\r\n% two given slopes determine a cubic.\r\n% The data points are known as \"knots\".\r\n% We have the _y_-values at the knots, so in order to get a particular PCHIP,\r\n% we have to somehow specify the values of the derivative, _y'_, at the knots.\r\n\r\n%%\r\n% Consider these two cubic polynomials in $x$ on the interval $1 \\le x \\le 2$ .\r\n% These functions are formed by adding cubic terms that vanish at\r\n% the end points to the linear interpolatant.\r\n% I'll tell you later where the coefficients of the cubics come from.\r\n%\r\n% $$ s(x) = 16 + 2(x-1) + \\textstyle{\\frac{49}{18}}(x-1)^2(x-2)\r\n%    - \\textstyle{\\frac{89}{18}}(x-1)(x-2)^2 $$\r\n%\r\n% $$ p(x) = 16 + 2(x-1) + \\textstyle{\\frac{2}{5}}(x-1)^2(x-2)\r\n%    - \\textstyle{\\frac{1}{2}}(x-1)(x-2)^2 $$\r\n%\r\n% These functions interpolate the same values at the ends.\r\n%\r\n% $$ s(1) = 16, \\ \\ \\ s(2) = 18 $$\r\n%\r\n% $$ p(1) = 16, \\ \\ \\ p(2) = 18 $$\r\n%\r\n% But they have different first derivatives at the ends.\r\n% In particular, $s'(1)$ is negative and $p'(1)$ is positive.\r\n%\r\n% $$ s'(1) = - \\textstyle{\\frac{53}{18}}, \\ s'(2) = \\textstyle{\\frac{85}{18}} $$\r\n%\r\n% $$ p'(1) = \\textstyle{\\frac{3}{2}}, \\ \\ \\ p'(2) = \\textstyle{\\frac{12}{5}} $$\r\n%\r\n\r\n%%\r\n% Here's a plot of these two cubic polynomials.\r\n% The magenta cubic, which is $p(x)$, just climbs steadily from its\r\n% initial value to its final value.\r\n% On the other hand, the cyan cubic, which is $s(x)$, starts off\r\n% heading in the wrong direction, then has to hurry to catch up.\r\n\r\nx = 1:1\/64:2;\r\ns = 16 + 2*(x-1) + (49\/18)*(x-1).^2.*(x-2) - (89\/18)*(x-1).*(x-2).^2;\r\np = 16 + 2*(x-1) + (2\/5)*(x-1).^2.*(x-2) - (1\/2)*(x-1).*(x-2).^2;\r\n\r\nclf\r\naxis([0 3 15 19])\r\nbox on\r\nline(x,s,'color',[0 3\/4 3\/4])\r\nline(x,p,'color',[3\/4 0 3\/4])\r\nline(x(1),s(1),'marker','o','color',[0 0 3\/4])\r\nline(x(end),s(end),'marker','o','color',[0 0 3\/4])\r\n\r\n%%\r\n% If we piece together enough cubics like these to produce a\r\n% piecewise cubic that interpolates many data points, we have a PCHIP.\r\n% We could even mix colors and still have a PCHIP.\r\n% Clearly, we have to be specific when it comes to specifying the slopes.\r\n\r\n%%\r\n% One possibility that might occur to you briefly is to use the\r\n% slopes of the lines connecting the end points of each segment.\r\n% But this choice just produces zeros for the coefficients of the cubics\r\n% and leads back to the piecewise linear interpolant.\r\n% After all, a linear function is a degenerate cubic.\r\n% This illustrates the fact that the PCHIP family includes many functions.\r\n\r\n%% spline\r\n% By far, the most famous member of the PCHIP family is the piecewise\r\n% cubic spline.\r\n% All PCHIPs are continuous and have a continuous first derivative.\r\n% A spline is a PCHIP that is exceptionally smooth, in the sense that its\r\n% second derivative, and consequently its curvature, also varies continuously.\r\n% The function derives its name from the flexible wood or plastic strip\r\n% used to draw smooth curves.\r\n\r\n%%\r\n% Starting about 50 years ago, Carl de Boor developed much of the\r\n% basic theory of splines.  He wrote a widely adopted package of Fortran\r\n% software, and a widely cited book, for computations involving splines.\r\n% Later, Carl authored the MATLAB Spline Toolbox.\r\n% Today, the Spline Toolbox is part of the Curve Fitting Toolbox.\r\n\r\n%%\r\n% When Carl began the development of splines, he was with\r\n% General Motors Research in Michigan.\r\n% GM was just starting to use numerically controlled machine tools.\r\n% It is essential that automobile parts have smooth edges and surfaces.\r\n% If the hood of a car, say, does not have continuously varying curvature,\r\n% you can see wrinkles in the reflections in the show room.\r\n% In the automobile industry, a discontinuous second derivative is\r\n% known as a \"dent\".\r\n\r\n%%\r\n% The requirement of a continuous second derivative leads to a set of\r\n% simultaneous linear equations relating the slopes at the interior knots.\r\n% The two end points need special treatment, and the default treatment\r\n% has changed over the years.\r\n% We now choose the coefficients so that the _third_ derivative does not have\r\n% a jump at the first and last interior knots.\r\n% Single cubic pieces interpolate the first three, and the last three,\r\n% data points.  This is known as the \"not-a-knot\" condition.\r\n% It adds two more equations to set of equations at the interior points.\r\n% If there are _n_ knots, this gives a well-conditioned, almost symmetric,\r\n% tridiagonal $n$ -by- $n$ linear system to solve for the slopes.\r\n% The system can be solved by the sparse backslash operator in MATLAB,\r\n% or by a custom, non-pivoting tridiagonal solver.\r\n% (Other end conditions for splines are available in the Curve Fitting Toolbox.)\r\n\r\n%%\r\n% As you probably realized, the cyan function $s(x)$ introduced\r\n% above, is one piece of the spline interpolating our sample data.\r\n% Here is a graph of the entire function, produced by |interpgui|\r\n% from NCM, _Numerical Computing with MATLAB_.\r\n\r\nx = 1:6;\r\ny = [16 18 21 17 15 12];\r\ninterpgui(x,y,3)\r\n\r\n%% sppchip\r\n% I just made up that name, |sppchip|.\r\n% It stands for _shape preserving piecewise cubic Hermite interpolating\r\n% polynomial_.  The actual name of the MATLAB function is just |pchip|.\r\n% This function is not as smooth as |spline|.\r\n% There may well be jumps in the second derivative.\r\n% Instead, the function is designed so that it never locally _overshoots_\r\n% the data.\r\n% The slope at each interior point is taken to be a weighted harmonic\r\n% mean of the slopes of the piecewise linear interpolant.\r\n% One-sided slope conditions are imposed at the two end points.\r\n% The |pchip| slopes can be computed without solving a linear system.\r\n\r\n%%\r\n% |pchip| was originally developed by Fred Fritsch and his colleagues\r\n% at Lawrence Livermore Laboratory around 1980.\r\n% They described it as \"visually pleasing\".\r\n% Dave Kahaner, Steve Nash and I included some of Fred's Fortran\r\n% subroutines in our 1989 book, _Numerical Methods and Software_.\r\n% We made |pchip| part of MATLAB in the early '90s.\r\n\r\n%%\r\n% Here is a comparison of |spline| and |pchip| on our data.\r\n% In this case the spline overshoot on the first subinterval\r\n% is caused by the not-a-knot end condition.\r\n% But with more data points, or rapidly varying data points,\r\n% interior overshoots are possible with |spline|.\r\n\r\ninterpgui(x,y,3:4)\r\n\r\n%% spline vs. pchip\r\n% Here are eight subplots comparing |spline| and |pchip| on a slightly\r\n% larger data set.\r\n% The first two plots show the functions $s(x)$ and $p(x)$.\r\n% The difference between the functions on the interior intervals\r\n% is barely noticeable.\r\n% The next two plots show the first derivatives.\r\n% You can see that the first derivative of |spline|, $s'(x)$, is smooth,\r\n% while the first derivative of |pchip|, $p'(x)$, is continuous,\r\n% but shows \"kinks\".\r\n% The third pair of plots are the second derivatives.\r\n% The |spline| second derivative $s''(x)$ is continuous,\r\n% while the |pchip| second derivative $p''(x)$ has jumps at the knots.\r\n% The final pair are the third derivatives.\r\n% Because both functions are piecewise cubics, their third derivatives,\r\n% $s'''(x)$ and $p'''(x)$, are piecewise constant.\r\n% The fact that $s'''(x)$ takes on the same values in the first two\r\n% intervals and the last two intervals reflects the\r\n% \"not-a-knot\" spline end conditions.\r\n\r\nsplinevspchip\r\n\r\n%% Locality\r\n% |pchip| is local.  The behavior of |pchip| on a particular subinterval\r\n% is determined by only four points, the two data points on either side\r\n% of that interval.  |pchip| is unaware of the data farther away.\r\n% |spline| is global.  The behavior of |spline| on a particular subinterval\r\n% is determined by all of the data, although the sensitivity to data far away\r\n% is less than to nearby data.\r\n% Both behaviors have their advantages and disadvantages.\r\n\r\n%%\r\n% Here is the response to a unit impulse.\r\n% You can see that the support of |pchip| is confined to the two\r\n% intervals surrounding the impulse, while the support of |spline|\r\n% extends over the entire domain.\r\n% (There is an  elegant set of basis functions for cubic splines known as\r\n% |B-splines| that do have compact support.)\r\n\r\nx = 1:8;\r\ny = zeros(1,8);\r\ny(4) = 1;\r\ninterpgui(x,y,3:4)\r\n\r\n%% interp1\r\n% The |interp1| function in MATLAB, has several |method| options.\r\n% The |'linear'|, |'spline'|, and |'pchip'| options are the same\r\n% interpolants we have been discussing here.\r\n% We decided years ago to make the |'cubic'| option the same as |'pchip'|\r\n% because we thought the monotonicity property of |pchip| was generally\r\n% more desirable than the smoothness property of |spline|.\r\n\r\n%%\r\n% The |'v5cubic'| option is yet another member of the PCHIP family,\r\n% which has been retained for compatibility with version 5 of MATLAB.\r\n% It requires the _x_'s to be equally spaced.\r\n% The slope of the v5 cubic at point $x_n$ is $(y_{n+1} - y_{n-1})\/2$.\r\n% The resulting piecewise cubic does not have a continuous second\r\n% derivative and it does not always preserve shape.\r\n% Because the abscissa are equally spaced, the v5 cubic can be\r\n% evaluated quickly by a convolution operation.\r\n\r\n%%\r\n% Here is our example data, modified slightly to exaggerate behavior,\r\n% and |interpgui| modified to include the |'v5cubic'| option of |interp1|.\r\n% The v5 cubic is the black curve between |spline| and |pchip|.\r\n\r\nx = 1:6;\r\ny = [16 18 21 11 15 12];\r\ninterpgui_with_v5cubic(x,y,3:5)\r\n\r\n%% Resources\r\n% A extensive collection of tools for curve and surface fitting, by\r\n% splines and many other functions, is available in the Curve Fitting Toolbox.\r\ndoc curvefit\r\n\r\n%%\r\n% \"NCM\", _Numerical Computing with MATLAB_, has more mathematical details.\r\n% NCM is available <https:\/\/www.mathworks.com\/moler\/index_ncm.html online>.\r\n% Here is the <https:\/\/www.mathworks.com\/moler.htmlinterp.pdf interpolation chapter>.\r\n% Here is <https:\/\/www.mathworks.com\/moler.htmlncm\/interpgui.m interpgui>.\r\n% SIAM publishes a <http:\/\/ec-securehost.com\/SIAM\/ot87.html print edition>.\r\n\r\n%%\r\n% Here are the script <https:\/\/blogs.mathworks.com\/images\/cleve\/splinevspchip.m % splinevspchip.m> and the modified version of interpgui\r\n% <https:\/\/blogs.mathworks.com\/images\/cleve\/interpgui_with_v5cubic.m % interpgui_with_v5cubic.m> that I used in this post.\r\n\r\n##### SOURCE END ##### 1add2164eaeb45548380e6801ab008b6\r\n-->","protected":false},"excerpt":{"rendered":"<!--introduction-->MATLAB has two different functions for piecewise cubic interpolation, <tt>spline<\/tt> and <tt>pchip<\/tt>. Why are there two? How do they compare?\r\n\r\n<!--\/introduction-->... <a class=\"read-more\" href=\"https:\/\/blogs.mathworks.com\/cleve\/2012\/07\/16\/splines-and-pchips\/\">read more >><\/a><\/p>","protected":false},"author":78,"featured_media":0,"comment_status":"open","ping_status":"closed","sticky":false,"template":"","format":"standard","meta":[],"categories":[11],"tags":[],"_links":{"self":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/196"}],"collection":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts"}],"about":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/types\/post"}],"author":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/users\/78"}],"replies":[{"embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/comments?post=196"}],"version-history":[{"count":19,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/196\/revisions"}],"predecessor-version":[{"id":2356,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/posts\/196\/revisions\/2356"}],"wp:attachment":[{"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/media?parent=196"}],"wp:term":[{"taxonomy":"category","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/categories?post=196"},{"taxonomy":"post_tag","embeddable":true,"href":"https:\/\/blogs.mathworks.com\/cleve\/wp-json\/wp\/v2\/tags?post=196"}],"curies":[{"name":"wp","href":"https:\/\/api.w.org\/{rel}","templated":true}]}}