三次样条插值

此示例显示如何使用csapi和csape命令从曲线拟合工具箱™构建三次样条插值。

CSAPI命令

命令

值=csapi（x，y，xx）

返回位于的值xx的三次样条插值到给定数据(x、 y)，使用非结结束条件。这个插值是一个带断点序列的分段三次函数x，其立方段连接在一起，形成一个具有两个连续导数的函数。“不结”结束条件意味着，在第一个和最后一个内部断裂处，即使是第三个导数也是连续的（直到舍入误差）。

仅指定两个数据点将生成一条直线插值。

x=[01]；y=[20]；xx=linspace（0,6121）；plot（xx，csapi（x，y，xx），“k -”，x，y，“罗”)头衔("插值到两点")

图中包含一个坐标轴。标题为插值到两点的轴包含两个类型为line的对象。

指定三个数据点可以得到一个抛物线。

x=[2 3 5]；y=[1 0 4]；绘图（xx，csapi（x，y，xx），“k -”，x，y，“罗”)头衔(“插入到三点”)

图中包含一个轴。标题插入到三个点的轴包含两个线型对象。

更一般地说，四个或更多的数据点给出一个三次样条曲线。

x=[11.524.15]；y=[1-11-11]；绘图（xx，csapi（x，y，xx），“k -”，x，y，“罗”)头衔(“三次样条插值到五点”)

图中包含一个轴。标题为三次样条插值到五个点的轴包含两个线型对象。

如何检查CSAPI的输出

这些看起来是很好的插值，但是我们怎么检查呢csapi执行广告吗?

我们已经看到了csapi插值，因为我们绘制了数据点，插值正好穿过这些点。但为了确保我们得到一个三次样条曲线，最好从预期排序的三次样条曲线开始，检查csapi复制这个三次样条，也就是说，返回了数据来自的三次样条。

示例:截断的幂函数

要检查的三次样条函数的一个简单示例是截断的三次方，即函数

$F (x) = ((x - x 我)_{+})^{3.},$

在哪里西是否有一个中断和“+”下标表示截断函数，由命令提供子脉冲:

帮助子脉冲

SUBPLUS正部分.x，如果x>=0 y=SUBPLUS（x）：=（x）{+}=，0，如果x<=0返回x的正部分。用于计算截断幂。

对于特定的选择，截短的三次幂图如下图所示西=2.．不出所料，它在2的左边是0，在2的右边是(x-2)^3。

情节(xx subplus (xx-2)。^ 3,“y”,“线宽”3)轴([0,-10,70])

图中包含一个坐标轴。轴包含一个线型对象。

现在，我们在数据点0:6处插值这条特殊的三次样条曲线，并在样条曲线顶部用黑色绘制插值曲线。

x=0:6；y=subplus（x-2）。^3；value=csapi（x，y，xx）；保持在绘图（xx，数值，“k”，x，y，“罗”)举行从头衔('插入剂到（（x-2）+^3')

图中包含一个坐标轴。标题为Interpolant to ((x-2)_+)^3的轴包含3个类型为line的对象。

插值误差

在比较两个函数时，绘制它们的差异通常会提供更多的信息。

绘图（xx，数值-子绘图（xx-2）。^3）标题(三次样条插值到((x-2)_+)^3的误差)

图中包含一个坐标轴。标题为误差三次样条插值到((x-2)_+)^3的轴包含一个类型为线的对象。

要将它们的差异大小放入上下文中，还可以计算最大数据值。这表明这个错误并不比不可避免的舍入错误更糟。

最大值=最大值（绝对值（y））

最大值y=64

无法复制的被截断的电源

作为进一步的检验，我们插入一个截断的幂csapi-produced interpolant at the sites 0:6 cannot coincide with it。例如，插值样条的第一个内部中断并不是一个真正的结csapi使用“不结”条件，因此插值点在该点处有三个连续导数。这意味着我们不能重现以该点为中心的截断三次方，因为其三次导数在该点处不连续。

数值=csapi（x，子图形（x-1）。^3，xx）；绘图（xx，数值-子图形（xx-1）。^3）标题((x-1)_+)^3的非结插值错误)

图中包含一个轴。插入到（（x-1）+^3的非节点中出现标题错误的轴包含一个line类型的对象。

因为1是第一个内结，它对这个内插不活跃。

这个差值有0.18那么大，但是当我们离开1时，这个差值会迅速衰减。这说明三次样条插值基本上是局部的．

使用ppform而不是值

可以将插值三次样条保留为适合后续求值、计算其导数或其他操作的形式csapi在形式上

pp=csapi（x，y）

它返回插入植物的ppform。您可以在一些新点上计算此表单xx由命令

值= fnval (pp、xx)

你可以通过命令来区分插值

民进党=曾经(pp)

或者通过命令对它进行积分

ipp=fnint（pp）

分别返回导数或积分的平形式。

示例：区分和集成插值

为了显示插值函数的微分，我们绘制了截断幂的导数

$F_{2.}^{'} (x) = 3. ((x - 2.)_{+})^{2.},$

(还是黄色的)然后在上面，插值函数对截断的三次幂函数的导数(还是黑色的)

情节(xx, 3 * subplus (xx-2)。^ 2,“y”,“线宽”pp = csapi(x,subplus(x-2).^3);民进党=曾经(pp);持有在地块（xx，fnval（dpp，xx），“k”)举行从头衔('插入剂对（（x-2）\+^3）的导数')

图中包含一个坐标轴。标题为“插值到((x-2)_+)^3的导数”的轴包含两个类型为line的对象。

同样，信息量更大的比较是绘制它们之间的差异，与之前一样，这并不比舍入误差大。

Plot (xx, fnval(dpp,xx) - 3*subplus(xx-2).^2) title(插值对((x-2)_+)^3求导的误差)

图中包含一个坐标轴。标题为“插值到((x-2)_+)^3的导数误差”的轴包含一个类型为line的对象。

截短后的幂次的二阶导数是

$F_{2.}^{''} (x) = 6. (x - 2.)_{+}$

这个函数和原始函数的插值函数的二阶导数之间的差值图显示现在有跳跃，但它们仍然在舍入误差之内。

ddpp =曾经(民进党);地块(xx, fnval(ddpp,xx) - 6*subplus(xx-2))((x-2)_+)^3插值二阶导数的误差)

图中包含一个坐标轴。标题为“插值到((x-2)_+)^3的二阶导数误差”的轴包含一个类型为line的对象。

截断幂的积分为

$F_{2.} (x) = ((x - 2.)_{+})^{4.} / 4. ．$

该函数与原始函数的插值积分之间的差值图再次表明，误差在舍入误差范围内。

ipp = fnint (pp);Plot (xx, fnval(ipp,xx) - subplus(xx-2).^4/4) title(((x-2)_+)^3的插值积分误差)

图中包含一个坐标轴。题目为((x-2)_+)^3积分误差的轴包含一个类型为line的对象。

CSAPE命令

喜欢csapi这个csape命令提供三次样条插值到给定数据。但是，它允许各种附加的结束条件。最简单的版本,

pp=csape（x，y）

使用拉格朗日结束条件，这是使用的非节点条件的常见替代方案csapi．csape不直接返回插值的值，而只返回它的ppform。

例如，再次考虑函数的插值

$F_{1.} (x) = ((x - 1.)_{+})^{3.},$

哪一个csapi不能很好地复制。我们绘制非结插值的误差csapi（黑色），以及从csape（红色）。

确切= subplus (xx-1)。^ 3;* * * * * * * * * * * * * * * * *“k”)举行在* * * * * * * * * * * * * * * *“r”)头衔(“非结与拉格朗日端点条件下的错误”)传奇({“Not-a-Knot”“拉格朗日”})；持有从

图中包含一个坐标轴。标题为“非结与拉格朗日端点条件误差”的轴包含两个线型对象。这些物体代表非结，拉格朗日。

在这种情况下，两种插入剂之间没有太大区别。

其他结束条件:“自然的”样条插值

这个csape命令还提供了为插值三次样条指定其他几种结束条件的方法。例如，命令

pp=csape（x，y，'变分'）

使用所谓的“自然”结束条件。这意味着二阶导数在两个极值处为零。

此步骤说明如何将“自然”三次样条插值应用于函数

$F_{2.} (x) = ((x - 2.)_{+})^{3.},$

下面的代码使用与“变分”字符串参数等效的替代参数语法计算“自然”样条插值：使用字符串“second”指定csape应将极端数据站点的二阶导数设置为默认值0。

页= csape (x, subplus (x - 2) ^ 3,“第二”)；地块（xx，fnval（pp，xx）-子地块（xx-2）。^3）标题('到（（x-2）\+^3的“自然”样条插值中存在错误')

图中包含一个轴。“自然”样条插值到（（x-2）+^3中标题错误的轴包含一个line类型的对象。

注意右边的大错误。这是由于“自然的”终端条件隐含地坚持二阶导数为零。

其他端部条件:规定的二级衍生物

我们也可以明确地使用正确的二阶导数来得到一个小的误差。首先，我们计算截断的幂在端点处的正确的二阶导数值。

endcond=6*子脉冲（x（[1结束]）-2）；

然后我们创建一个插值函数，指定端点处的二阶导数与我们刚刚计算的二阶导数值相匹配。我们通过提供来做到这一点结束秒（1）对于左端点条件，以及结束秒（2）对于右侧，以及数据值。

pp=csape（x，[endcond（1）子脉冲（x-2）。^3 endcond（2）]，“第二”）;Plot (xx, fnval(pp,xx)) - subplus(xx-2).^3，“r”)标题('样条曲线插值到（（x-1）+^3时出错';．..“当在末端匹配二阶导数时”])

图中包含一个坐标轴。当在端点匹配二阶导数时，标题为“样条插值误差(x-1)_+)^3”的轴包含一个类型为line的对象。

其他结束条件:规定坡度

csape也允许端点的规范山坡上.这是夹紧(或者,完整的)三次样条插值。该声明

pp=csape（x，[sl，y，sr]，'clipped'）

创建数据的三次样条插值(x,Y)也有斜率sl在最左边的数据点和斜坡老在最右边的数据站点。

其他结束条件:混合结束条件

甚至可以混合这些条件。例如，我们经常使用的截断幂函数

$F_{1.} (x) = ((x - 1.)_{+})^{3.}$

斜率为0x=0，二阶导数为30x=6(最后一个数据站点)。

因此，通过匹配左端点的斜率和右端点的曲率，我们期望在得到的插值中没有误差。

pp=csape（x，[0子像素（x-1）。^3 30]，[1 2]）；地块（xx，fnval（pp，xx）-子地块（xx-1）。^3）标题(['样条曲线插值到（（x-1）+^3时出错';．..“具有混合的最终条件。”])

图中包含一个坐标轴。样条插值误差轴为((x-1)_+)^3与混合结束条件。包含类型为line的对象。

其他结束条件:周期条件

也可以开处方周期性结束条件。例如，正弦函数是2*pi周期性函数，其值为[0 -1 0 1 0]在现场(π/ 2)* (2:2)．正弦函数和它的周期三次样条插值在这些点之间的差异仅为2%。不坏。

x =(π/ 2)* (2:2);Y = [0 -1 0 1 0];页= csape (x, y,“周期性”)；xx=linspace（-pi，pi，201）；plot（xx，sin（xx）-fnval（pp，xx），“x”)头衔(周期三次样条插值到sin(x)的误差)

图中包含一个坐标轴。标题为“到sin(x)的周期三次样条插值误差”的轴包含一个类型为line的对象。

CSAPI或cape未显式覆盖的结束条件

未显式覆盖的任何结束条件csapi或csape可以通过构造插值来处理csape默认的边条件，然后向其中添加插值到零值的适当标量倍数和一些边条件。如果需要满足两个“非标准”边条件，则可能必须首先解一个2乘2的线性系统。

例如，假设您想要计算三次样条插值s对数据的修改

x=0:25:3；q=@（x）x.*（-1+x.*（-1+x.*x/5））；y=q（x）；

并强制执行该条件

lambda（s）：=a*（Ds）（e）+b*（D^2 s）（e）=c

的一阶和二阶导数s当时E．

数据是由一个四次多项式产生的，恰好满足这个边条件与特定的参数

e = x (1);= 2;b = 3;c = 4;

为了构造满足此特定条件的插值，我们首先使用默认的结束条件构造插值

pp1 = csape (x, y);

它的第一个多项式部分的一阶导数。

dp1 =曾经(fnbrk (pp1 1));

此外，我们将三次样条插值构造为零数据值，并指定其斜率为1E,

pp0=csape（x，[1，零（大小（y）），0]，[1,0]）；

以及构造它的第一个多项式的一阶导数。

dp0=fnder（fnbrk（pp0,1））；

然后我们计算λ两者皆有pp1和pp0,

lam1=a*fnval（dp1，e）+b*fnval（fnder（dp1，e）；lam0=a*fnval（dp0，e）+b*fnval（fnder（dp0，e）；

并构造正确的线性组合pp1和pp0得到三次样条

S:= pp1 + ((c - lambda(pp1))/lambda(pp0)) * pp0

这确实满足了所需的条件，以及右端点的默认结束条件。这个线性组合是通过fncmb．

s = fncmb (pp0 (c-lam1) / lam0 pp1);

插值误差图显示：s拟合四次多项式的效果稍好一些E而不是插入剂pp1在默认条件下可以。

xx=（.3）：.05:7；yy=q（xx）；图（xx，fnval（pp1，xx）-yy，“x”)举行在图（xx，fnval（s，xx）-yy，“哦”)举行从传奇({“默认条件”“非标准条件”},“位置”,“本身”)

图中包含一个轴。轴包含2个类型为line的对象。这些对象表示默认条件、非标准条件。

如果我们想强制执行这个条件

mu（s）：=（D^3 s）（3）=14.6

在插值的三阶导数上(四次函数满足这个条件)，然后我们构造一个额外的三次样条插值到零值，并且在左端点的一阶导数为零，因此确定是独立的pp0．

pp2 = csape (x,[0, 0(大小(y)), 1], [0,1]);

然后我们找到系数d0和d2在线性组合中

s:=pp1+d0*pp0+d2*pp2

这就解出了线性方程组

lambda（s）=c

μ(s) = 14.6

请注意,这两个pp0和pp2在所有插值点消失，因此s将匹配给定的数据为任何选择d0和d2．

为了好玩，我们使用MATLAB®编码工具编写一个循环来进行计算λ(pp_j)和μ(pp_j)对于J=0:2.

dd=零（2,3）；对于j=0:2j=num2str（j）；评估([“民进党”J'=fnder（pp'J');'])；评估([“ddpp”J'=fnder（民进党'J');'])；评估(['dd（1,1+'J'）=a*fnval（民进党'J“，e）+b*fnval（ddpp”J“，e）；”])；评估(['dd（2,1+'J'）=fnval（fnder（ddpp'J”),3),“]）;结束

给定的值λ和μ对于pp0,pp1,pp2，然后解出定义正确线性组合的系数。

d = dd (: [1,3]) \ ([c; 14.6] dd (:, 2));s = fncmb (fncmb (pp0, d (1) pp2 d (2)), pp1);xxx = 0: .05:3;多= q (xxx);Plot (xxx, yyy - fnval(s,xxx)，“x”)头衔('样条插值到y=x*（-1+x*（-1+x*x/5））的误差')

图中包含轴。样条线插值中标题错误为y=x*（-1+x*（-1+x*x/5））的轴包含线型对象。

为了保证，我们将此误差与此函数的完整三次样条插值得到的误差进行比较:

持有在图（xxx，yyy-fnval）（csape（x，[-1，y，-7+（4/5）*27]，“夹紧”)，xxx），“哦”)举行从传奇({“非标准条件”“端坡条件”})

图中包含一个轴。样条线插值中标题错误为y=x*（-1+x*（-1+x*x/5））的轴包含2个类型为line的对象。这些对象表示非标准条件、端点坡度条件。

错误的差别（不是很大）只是在接近终点时，这证明了这两种错误pp0和pp2只有在它们各自的端点附近才会相当大。

作为最后的检查，我们验证一下s在3处满足所需的三阶导数条件。

fnval（fnder（s，3），3）

ans = 14.6000