三次样条插值

打开实时脚本

方法的使用csapi而且csape命令从曲线拟合工具箱™构建三次样条插值。

CSAPI命令

命令

值= csapi(x,y,xx)

返回的值xx的三次样条插值函数(x, y)，使用非结结束条件。该插值函数为分段三次函数，具有间断序列x，它的三次分量连接在一起，形成一个有两个连续导数的函数。“非结”结束条件意味着，在第一个和最后一个内部断裂处，即使是三阶导数也是连续的(直到舍入误差)。

只指定两个数据点就会得到直线插值。

X = [0 1];Y = [2 0];Xx = linspace(0,6,121);情节(xx csapi (x, y, xx),“k -”, x, y,“罗”)标题(“两点插值”）

图中包含一个轴对象。标题为Interpolant to Two Points的axis对象包含2个类型为line的对象。

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

X = [2 3 5];Y = [1 0 4];情节(xx csapi (x, y, xx),“k -”, x, y,“罗”)标题(“三点插值”）

图中包含一个轴对象。标题为Interpolant to Three Points的axes对象包含2个类型为line的对象。

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

X = [1 1.5 2 4.1 5];Y = [1 -1 1 -1 1];情节(xx csapi (x, y, xx),“k -”, x, y,“罗”)标题(五次样条插值算法）

图中包含一个轴对象。标题为“五点三次样条插值”的坐标轴对象包含2个类型为直线的对象。

如何查看CSAPI的输出

这些看起来是很好的插值，但是我们怎么检验呢csapi就像广告上说的那样?

我们已经看到了csapi插值，因为我们绘制了数据点插值器穿过了这些点。但是为了确保我们得到一个三次样条，最好从期望排序的三次样条的数据开始，并检查是否csapi繁殖这个三次样条，也就是说，返回数据所取自的三次样条。

示例:截断幂函数

三次样条函数的一个简单例子是截断的三次幂，即函数

$f （ x ）＝（（ x - x 我）_{+} ）^{3.} ，$

在哪里西其中一个断行符和“+”下标是否表示截断函数，由命令提供subplus：

帮助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])

图中包含一个轴对象。axis对象包含一个line类型的对象。

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

X = 0:6;Y = + (x-2) ^3;值= csapi(x,y,xx);持有在情节(xx,价值观,“k”, x, y,“罗”)举行从标题(插值到(x-2)_+)^3）

图中包含一个轴对象。坐标轴对象的标题I nt e r p o lan t空白t o空白((x - 2) indexOf + baseline)立方基线包含3个类型为line的对象。

插补误差

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

Plot (xx, values - subplus(xx-2).^3)三次样条插值到(x-2)_+)^3的误差）

图中包含一个轴对象。标题为E r r o r空白in空白Cub i C空白S p lin空白i on空白t o空白((x - 2) indexOf + baseline)立方基线包含一个类型为line的对象。

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

Max_y = max(abs(y))

Max_y = 64

一种无法复制的被截断的权力

作为进一步的测试，我们插值一个截断的幂csapi-在0:6处产生的插值不能与之重合。例如，插值样条的第一个内部断裂并不是真正的结，因为csapi使用“非结”条件，因此插值器在该位置有三个连续导数。这意味着我们不应该能够重现以该位置为中心的截断的三次方，因为它的三阶导数在该位置是不连续的。

值= csapi(x,subplus(x-1).^3,xx);Plot (xx, values - subplus(xx-1).^3)非结插值到((x-1)_+)^3的误差）

图中包含一个轴对象。标题为E r r o r blank in blank n o t - a - K n o t blank in t E r p o lan t blank t o blank ((x - 1) indexOf + baseline)立方基线包含一个类型为line的对象。

因为1是第一个内结，所以它对于这个插值器是无效的。

差值高达0.18，但当我们远离1时，差值迅速衰减。这说明了三次样条插值本质上是局部的．

使用ppform代替Values

可以将插值三次样条保留为适合于后续计算或计算其导数或其他操作的形式。这是通过调用来完成的csapi在表格中

Pp = csapi(x,y)

它返回插值函数的ppform。你可以在一些新的点上计算这个形式xx通过命令

值= fnval(pp,xx)

您可以用命令来区分插值器

DPP = fnder(pp)

或者通过命令积分

Ipp = fnint(pp)

它们分别返回导数或积分的ppform。

示例:微分和积分插值

为了显示插值函数的微分，我们画出这个截断幂的导数

$f_{2}^{”} （ x ）＝ 3. （（ x - 2 ）_{+} ）^{2} ，$

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

情节(xx, 3 * subplus (xx-2)。^ 2,“y”，“线宽”，3) pp = csapi(x,subplus(x-2).^3);DPP = fnder(pp);持有在情节(xx fnval(民进党,xx),“k”)举行从标题(插值函数对(x-2)_+)^3的导数）

图中包含一个轴对象。标题为D e r i v a ti v e blank of blank ine r p o lan t blank t o blank ((x - 2) indexOf + baseline)三次基线包含2个类型为line的对象。

同样，更有信息的比较是绘制它们的差异，和前面一样，这并不比舍入误差大。

Plot (xx, fnval(dpp,xx) - 3*subplus(xx-2).^2)插值函数((x-2)_+)^3的导数误差）

图中包含一个轴对象。标题为E r r o r空白in空白D E r i v a ti v E空白of空白in E r p o lan t空白t o空白((x - 2) indexOf + baseline)立方基线包含一个类型为line的对象。

截断幂的二阶导数是

$f_{2}^{” ”} （ x ）＝ 6 （ x - 2 ）_{+}$

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

DDPP = fnder(dpp);Plot (xx, fnval(ddpp,xx) - 6*subplus(xx-2))插值函数((x-2)_+)^3二阶导数的误差）

图中包含一个轴对象。标题为E r r o r空白in空白S E con d空白d E r i v a ti v E空白of空白in E r p o lan t空白t o空白((x - 2) indexOf + baseline)立方基线包含一个类型为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积分的误差）

图中包含一个轴对象。标题为E r r o r空白in空白in t E g r a l空白of空白in E r p o l nt空白t o blank ((x - 2) indexOf + baseline)立方基线包含一个类型为line的对象。

CSAPE命令

就像csapi,csape命令提供给定数据的三次样条插值。但是，它允许各种附加的结束条件。最简单的说法是，

Pp = csape(x,y)

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

例如，再次考虑对函数进行插值

$f_{1} （ x ）＝（（ x - 1 ）_{+} ）^{3.} ，$

哪一个csapi不能很好地繁殖。我们画出由返回的非结插值的误差csapi(黑色)，以及从中获得的插值误差csape(红色)。

Exact = subplus(xx-1).^3;Plot (xx, fnval(csapi(x,subplus(x-1).^3)，xx) -准确，“k”)举行在Plot (xx, fnval(csape(x,subplus(x-1).^3)，xx) -准确，“r”)标题(“非结的错误与拉格朗日结束条件”)({传奇“Not-a-Knot”“拉格朗日”});持有从

图中包含一个轴对象。标题为Error in Not-a-Knot vs. Lagrange End Conditions的axis对象包含2个类型为line的对象。这些物体代表非结，拉格朗日。

在这种情况下，这两个插补之间没有太大的区别。

其他结束条件:'自然'样条插值

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

Pp = csape(x,y，'变分')

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

这一步展示了如何对函数应用“自然”三次样条插值

$f_{2} （ x ）＝（（ x - 2 ）_{+} ）^{3.} ，$

然后画出误差。下面的代码计算“自然”样条插值与替代参数语法等效“变分”论点:使用“第二”指定csape应将极端数据点的二阶导数设置为默认值0。

Pp = csape(x, x-2的次幂)^3，“第二”）;Plot (xx, fnval(pp,xx) - subplus(xx-2).^3)“自然样条插值误差((x-2)_+)^3”）

图中包含一个轴对象。标题为E r r o r空白in空白n a t u r a l '空白S p in空白in r p o l空白t on空白t o空白((x - 2) indexOf + baseline)立方基线包含一个类型为line的对象。

注意右端附近的大误差。这是由于“自然”结束条件隐含地坚持在那里有一个零二阶导数。

其他最终条件:规定二次衍生品

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

Endcond = 6*subplus(x([1 end])-2);

然后我们创建插值，指定端点处的二阶导数与我们刚刚计算的二阶导数值相匹配。我们通过提供来做到这一点endcond (1)对于左端点条件，和endcond (2)对于右边，连同数据值一起。

Pp = csape(x，[endcond(1) subplus(x-2).]^ 3 endcond (2)),“第二”）;Plot (xx, fnval(pp,xx) - subplus(xx-2).^3，“r”)标题(样条插值到(x-1)_+)^3的误差；.．.“当在端点处匹配二阶导数时”]）

图中包含一个轴对象。坐标轴对象与标题E r r o r空白我n空白S p l i n E空白i n t E r p o l t i o n空白t o空白((x - 1) indexOf +基线)立方基线空白空白空白W h E n空白M t c h i n g空白t h E空白2 n d空白d E r我v t v E空白t E n d S空白空白包含一个类型的对象。

其他终点条件:规定坡度

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

Pp = csape(x，[sl,y,sr]，'clamped')

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

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

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

$f_{1} （ x ）＝（（ x - 1 ）_{+} ）^{3.}$

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

因此，通过匹配左边的斜率和右边的曲率，我们期望得到的插值没有错误。

Pp = csape(x，[0次+ (x-1))。^3 30]， [1 2]);Plot (xx, fnval(pp,xx) - subplus(xx-1).^3) title([样条插值到(x-1)_+)^3的误差；.．.“终端条件混合。”]）

图中包含一个轴对象。坐标轴对象与标题E r r o r空白我n空白S p l i n E空白i n t E r p o l t i o n空白t o空白((x - 1) indexOf +基线)立方基线空白空白空白空白空白空白空白空白空白w我t h空白M x E d空白E n d空白C o n d i t i o n S。空白空白空白空白空白空白空白空白空白空白空白空白空白空白空白空白包含一个对象的类型行。

其他结束条件:周期条件

也可以开处方周期结束条件。例如，正弦函数是2*pi周期函数，并且有这些值[0 -1 0 1 0]在现场(π/ 2)* (2:2)．正弦函数和它的周期三次样条插值在这些位置上的差别只有2%。不坏。

X = (pi/2)*(-2:2);Y = [0 -1 0 1 0];Pp = csape(x,y，“周期”）;Xx = linspace(-pi,pi,201);Plot (xx, sin(xx) - fnval(pp,xx)，“x”)标题(周期三次样条插值到sin(x)的误差）

图中包含一个轴对象。标题为Error in Periodic Cubic Spline Interpolation to sin(x)的axis对象包含一个类型为line的对象。

CSAPI或CSAPE未明确涵盖的结束条件

未明确涵盖的任何结束条件csapi或csape可以通过构造插值来处理csape默认边条件，然后向它添加一个适当的标量乘以插值零值和一些边条件。如果有两个“非标准”侧条件要满足，你可能必须先解一个2 × 2线性系统。

例如，假设你想计算三次样条插值年代数据

X = 0:.25:3;Q = @(x) x.*(-1+x.* (-1+x.*x/5));Y = q(x);

并执行这个条件

lambda(s):= a * (Ds)(e) + b * (D²s)(e) = c

的一阶和二阶导数年代在这一点上e．

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

E = x(1);A = 2;B = -3;C = 4;

为了构造满足这个特定条件的插值函数，我们首先用默认的结束条件构造插值函数

Pp1 = csape(x,y);

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

Dp1 = fnder(fnbrk(pp1,1));

此外，我们构造了三次样条插值到零数据值，指定它在的斜率为1e，

Pp0 = csape(x，[1, 0 (size(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);

插补误差图说明了这一点年代拟合四次多项式稍好e而不是插值函数pp1在默认条件下。

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

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

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

(s) = (D³s)(3) = 14.6

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

Pp2 = csape(x，[0, 0 (size(y))，1]，[0,1]);

然后我们求出系数d0而且d2在线性组合中

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

解线性方程组

s = c

Mu (s) = 14.6

请注意pp0而且pp2在所有插值点消失，因此年代将匹配给定的数据为任何选择d0而且d2．

为了好玩，我们使用MATLAB®编码工具来编写一个循环来计算λ(pp_j)而且μ(pp_j),因为j= 0:2。

Dd = 0 (2,3);为j=0:2 j= num2str(j);eval ([“民进党”J' =曾经(pp 'J”),“]);eval ([“ddpp”J=曾经(民进党的J”),“]);eval ([“dd(1,1 +”J“* fnval(民进党)=”J”,e) + b * fnval (ddpp 'J“e);”]);eval ([“dd(2, 1 +”J') = fnval(曾经(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;Yyy = q(xxx);Plot (xxx, yyy - fnval(s,xxx)，“x”)标题(' y = x*(-1+x* (-1+x*x/5))的样条插值错误'）

图中包含一个轴对象。标题为Error in Spline Interpolant to y = x*(-1+x* (-1+x*x/5))的axis对象包含一个类型为line的对象。

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

持有在Plot (xxx, yyy - fnval(csape(x，[-1,y，-7+(4/5)*27]，“夹紧”), xxx),“o”)举行从传奇({“非标准条件”“Endslope条件”})

图中包含一个轴对象。标题为Error in Spline Interpolant to y = x*(-1+x* (-1+x*x/5))的axis对象包含2个类型为line的对象。这些对象代表非标准条件，端坡条件。

误差只在接近终点的地方有所不同(而且差别不大)，证明了这两个事实pp0而且pp2只有在它们各自的终点附近才相当大。

作为最后的检查，我们验证了这一点年代满足3处的三阶导数条件。

fnval(曾经(年代,3),3)

Ans = 14.6000