MATLAB有两种不同的分段三次插值函数，样条和pchip．为什么有两个?它们是如何比较的?

内容

数据

以下是我将在这篇文章中使用的数据。

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

X = 12 3 4 5 6 y = 16 18 21 17 15 12

这是一张数据图。

集(0,“defaultlinelinewidth”，2) CLF plot(x,y，“o”)轴([0 7 7.5 25.5])标题(“plip”）

plip

使用线类型'-o'， MATLAB plot命令在六个数据点处绘制六个'o'，并在点之间绘制直线。所以我添加了标题plip因为这是分段线性插值多项式．每对点之间有一个不同的线性函数。因为我们想让函数遍历这些数据点插入数据，因为两点决定了一条直线，所以plip功能是唯一的。

PCHIP系列

一个PCHIP，一个分段三次赫米特插值多项式是插值给定数据的任意分段三次多项式，并且在插值点有指定的导数。正如两点决定一个线性函数一样，两点和两个给定的斜率决定一个三次函数。这些数据点被称为“结点”。我们有y-值，所以为了得到一个特定的PCHIP，我们必须以某种方式指定导数的值，y '打结。考虑$x$上的这两个三次多项式在区间$1 \le x \le 2$上。这些函数是通过向线性插值函数中添加在端点处消失的三次项而形成的。我稍后会告诉你们立方的系数是怎么来的。$ $ s (x) = 16 + 2 (x - 1) + \ textstyle{\压裂{49}{18}}(x - 1) ^ 2 (x - 2) - \ textstyle{\压裂{89}{18}}(x - 1) (x - 2) ^ 2 $ $ $ $ p (x) = 16 + 2 (x - 1) + \ textstyle{\压裂{2}{5}}(x - 1) ^ 2 (x - 2) - \ textstyle{\压裂{1}{2}}(x - 1) (x - 2) ^ 2 $ $这些功能末端插入相同的值。$$ s(1) = 16， \ \ \ s(2) = 18 $$ $$ p(1) = 16， \ \ \ p(2) = 18 $$但它们在末端有不同的一阶导数。特别地，$s'(1)$是负的，$p'(1)$是正的。$ $ s’(1)= - \ textstyle{\压裂{53}{18}},\ s (2) = \ textstyle{\压裂{85}{18}}$ $ $ $ p’(1)= \ textstyle{\压裂{3}{2 }}, \ \ \ p ' (2) = \ textstyle{\压裂{12}{5}}$ $这是一块这两个三次多项式。洋红色的立方，也就是p(x)，从初始值稳步上升到最终值。另一方面，青色的立方体，也就是$s(x)$，开始的方向是错误的，然后必须尽快赶上。

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轴([0 3 15 19])框在线(x,年代,“颜色”，[0 3/4 3/4]) line(x,p，“颜色”，[3/4 0 3/4]) line(x(1)，s(1)，“标记”，“o”，“颜色”，[0 0 3/4]) line(x(end)，s(end)，“标记”，“o”，“颜色”，[0 0 3/4]

如果我们把足够多的像这样的立方体拼接在一起，产生一个分段立方体，插入许多数据点，我们就有了一个PCHIP。我们甚至可以混合颜色，仍然有PCHIP。很明显，在确定斜率时，我们必须明确。你可能会想到的一种简单的方法是使用连接每个线段端点的直线的斜率。但是这个选择只会让三次函数的系数为零然后回到分段线性插值。毕竟，线性函数是退化的三次函数。这说明PCHIP系列包含许多功能。

样条

到目前为止，PCHIP家族中最著名的成员是分段三次样条。所有pchip都是连续的，并且具有连续的一阶导数。样条是一种非常光滑的PCHIP，因为它的二阶导数和曲率也是连续变化的。该功能得名于用于绘制光滑曲线的柔性木材或塑料条。从大约50年前开始，卡尔·德·布尔发展了许多样条的基本理论。他写了一个被广泛采用的Fortran软件包和一本被广泛引用的书，用于涉及样条曲线的计算。后来，Carl编写了MATLAB样条工具箱。现在，样条工具箱是曲线拟合工具箱的一部分。当卡尔开始花键的发展，他是通用汽车研究在密歇根州。通用汽车刚刚开始使用数控机床。 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 the第三导数在第一个和最后一个内部结点处没有跳跃。单立方块插值前三个数据点，后三个数据点。这就是所谓的“不打结”情况。它给内部点的方程组增加了两个方程。如果有的话n结，这给出了一个条件良好的，几乎对称的，三对角线n × n线性方程组来解斜率。该系统可以通过MATLAB中的稀疏反斜杠算子或自定义的非旋转三对角线求解器来求解。(其他样条的末端条件可在曲线拟合工具箱中获得。)您可能已经意识到，上面介绍的青色函数$s(x)$是样条插值样例数据的一部分。这是整个函数的图像，由interpgui不合格品,MATLAB数值计算．

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

sppchip

我刚编了个名字sppchip．它代表的是保形分段三次埃尔米特插值多项式．MATLAB函数的实际名称是pchip．这个函数不像样条．二阶导数可能会有跳跃。相反，该函数的设计使其永远不会在本地发生过激的数据。每个内部点的斜率取为分段线性插值函数斜率的加权调和平均值。在两个端点处施加单侧斜率条件。的pchip不用解线性方程组就可以计算斜率。pchip最初是由劳伦斯利弗莫尔实验室的弗雷德·弗里奇和他的同事在1980年左右发明的。他们将其描述为“视觉上的愉悦”。戴夫·卡哈纳、史蒂夫·纳什和我在1989年出版的书中收录了弗雷德的一些Fortran子程序，数值方法与软件．我们做了pchip是90年代早期MATLAB的一部分。这里有一个比较样条和pchip根据我们的数据。在这种情况下，第一个子区间上的样条超调是由非结结束条件引起的。但是随着数据点的增多，或者数据点的快速变化，内部超调就有可能发生样条．

interpgui (x, y, 3:4)

样条vs. pchip

下面是八个比较的子图样条和pchip在稍微大一点的数据集上。前两个图显示了函数$s(x)$和$p(x)$。内部间隔上的函数之间的差异几乎不明显。下面两个图显示了一阶导数。你可以看到一阶导数样条$s'(x)$是光滑的，而pchip， $p'(x)$，是连续的，但显示“扭结”。第三对图是二阶导数。的样条二阶导数$s''(x)$是连续的，而pchip二阶导数p''(x)在结点处跳跃。最后一对是三阶导数。因为两个函数都是分段三次函数，它们的三阶导数$s'' (x)$和$p'' (x)$都是分段常数。$s'' (x)$在前两个区间和后两个区间取相同的值反映了“非结”样条结束条件。

splinevspchip

位置

pchip是本地的。的行为pchip在一个特定的子区间上，仅由四个点决定，即该区间两侧的两个数据点。pchip不知道远处的数据。样条是全球性的。的行为样条在一个特定的子区间上，由所有数据决定，尽管对远数据的敏感性低于对近数据的敏感性。这两种行为都有其优点和缺点。这是对单位脉冲的响应。你可以看到支持金宝apppchip是局限于两个间隔周围的冲量，而支撑金宝app样条扩展到整个域。(有一组优雅的三次样条的基函数叫做b样它们确实有紧密的支持。)金宝app

X = 1:8;Y = 0 (1,8);Y (4) = 1;interpgui (x, y, 3:4)

interp1

的interp1函数，有几个方法选项。的“线性”，样条的,“pchip”选项是我们在这里讨论的相同的插值。我们几年前就决定“立方”选项与“pchip”因为我们认为的单调性pchip通常比的平滑性更可取样条．的“v5cubic”option是PCHIP家族的另一个成员，它保留了与MATLAB 5版的兼容性。它需要x是等距的。v5三次函数在x_n点处的斜率为(y_{n+1} - y_{n-1})/2。得到的分段三次函数没有连续的二阶导数，也不总是保持形状。因为横坐标是等距的，所以v5立方可以通过卷积运算快速求出。下面是我们的示例数据，稍作修改以夸大行为interpgui修改后包括“v5cubic”选择interp1．v5的三次是黑色曲线样条和pchip．

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

资源

曲线拟合工具箱中有大量的曲线和曲面拟合工具，包括样条和许多其他函数。

医生curvefit

“不合格品”,MATLAB数值计算，有更多的数学细节。NCM可用在线．这是插值章．这是interpgui．SIAM出版了印刷版．下面是脚本splinevspchip.m以及修改后的interpguiinterpgui_with_v5cubic.m我在这篇文章中使用的。

获取MATLAB代码使用MATLAB®7.14发布

使用MATLAB®7.14发布