三次样条插值

这个示例说明了如何使用曲线拟合工具箱中的 csapi 和 csape 构造三次样条插值。

光滑数据的三次样条插值

假设您想要插值一些光滑数据，比如

using TyPlot
using TyCurveFitting
using TyMath
rng = MT19937ar(6)
x = (4*pi).*[0;1;rand(rng,15)]
y = sin.(x)

您可以使用下述指令生成三次样条插值

cs = csapi(x,y)

然后使用以下指令绘制样条与数据：

fnplt(cs)
hold("on")
plot(x,y,"o")
legend(["cubic spline","data"])
hold("off")

这会生成以下的图。

光滑数据的三次样条插值

更准确地说，这是带非扭结边界条件的三次样条插值，这意味着它是它是唯一一个在所有内部数据位点的断点，也就是除了最左边以及最右边的点，都具有二阶连续导数的分段三次多项式。这是使用基础数学库 spline 指令 spline(x,y) 生成的同样的插值（格式不同）。

周期数据

正弦函数以 2π 为周期，为了检测您的插值是否具有同样的周期，可以计算，比如它在两个端点的一阶导数的差，

diff(fnval(fnder(cs),[0 4*pi])[:])

ans = 0.27039142633574487

这看上去不是很好。如果您更想要两个端点 0 和 4*pi 的一阶和二阶导相匹配的插值，可以使用 csape，这允许采取不同种类的边界条件，包括周期边界条件。因此，可以使用

pcs = csape(x,y,"periodic")

然后使用

diff(fnval(fnder(pcs),[0 4*pi])[:])

就可以得到端点导数的差为 - 1.1102230246251565e-16。即使二阶导数的差也很小：

diff(fnval(fnder(pcs,2),[0 4*pi])[:])

ans = 1-element Vector{Float64}:

 5.134781488891349e-16

其他端点条件

也可以采用其他的端点条件，比如，

cs = csape(x',[3 y' -4],[1 2])  

提供了以 x[i] 为端点，最左边数据位点的导数为 3，最右端位点的二阶导数为 -4 的三次样条插值。

一般的样条插值

如果您想要使用具有单节点的除了三次样条的样条或者在非断点位点插值，则您需要使用 spapi 函数。其最简单形式可以为 sp = spapi(k,x,y) 其中第一个参数 k 指定插值样条的阶数，这是每个多项式片段的系数的数目，即比多项式片段的度数大一。比如，接下来的图片展示了对您的数据的线性、二次以及四次样条插值。

sp2 = spapi(2,x,y)
fnplt(sp2,2)
hold("on")
sp3 = spapi(3,x,y)
fnplt(sp3,2,"k--")
sp5 = spapi(5,x,y)
fnplt(sp5,2,"r-.")
plot(x,y,"o")
legend(["linear","quadratic","quartic","data"])
hold("off")

同样的光滑数据不同阶数的样条插值

即使是使用 spapi 得到的三次样条插值也与从 csapi 以及 spline 得到的结果不同。为了强调它们的差别，如下计算以及绘制它们的二阶导数：

fnplt(fnder(spapi(4,x,y),2))
hold("on")
fnplt(fnder(csapi(x,y),2),2,"k--")
plot(x,zeros(size(x)),"o")
legend(["from spapi","from csapi","data sites"])
hold("off")

这会生成以下的图。

同样的光滑数据的两个三次样条插值的二阶导

因为三次样条的二阶导为分段线段，端点为样条的断点，您可以清楚地看到 csapi 的断点就是数据位点，而 spapi 不然。

节点选择

实际上说明使用指令 sp = spapi(knots,x,y) 样条插值如何选择断点是存在这个可能性的，至少它使用的断点肯定是来自于提供的序列 knots。举个例子，因为 y 实际上是 sin.(x)，指令

ch = spapi(augknt(x,4,2)[1],[x;x],[y;cos.(x)])

为正弦函数生成了三次 Hermite 插值，这称之为分段三次函数，其中断点为所有的 x[i]，且与正弦函数在这些所有的 x[i] 的函数值与一阶导数值相匹配。这使得这个插值一阶导数连续，但是通常二阶导数断点处存在阶跃。只是这个指令是怎么知道数据值序列 [y;cos.(x)] 哪一部分为函数值，哪一部分为导数值呢？注意到这里的数据位点序列给为 [x;x]，即每个数据位点均出现两次。同样注意到 y[i] 是与 x[i] 的第一次出现相联系的，而 cos(x[i]) 是与 x[i] 的第二次出现相联系的。因此，与数据位点的第一次出现相联系的数据值将被认为是函数值，而与第二次出现相相联系的就会被认为是一阶导数值，如果存在第三次出现的数据位点，则与其相联系的数据值就会被认为是二阶导数值。关于这里用了生成适宜的“节点序列”的指令 augknt 的更多讨论可以参见构造以及处理 B 样条。

平滑

如果数据是含噪的呢？比如，假设给定的数据为

rng = MersenneTwister(1234)
noisy = y + .3 .* (rand(rng,size(x)...).-.5)

则您可能更想要逼近（而非插值）。比如您可能像试试三次平滑样条，这可以由以下指令给出

scs, = csaps(x,noisy)

然后使用以下指令绘制

fnplt(scs,2)
hold("on")
plot(x,noisy,"o")
legend(["smoothing spline","noisy data"])
hold("off")

含噪数据的三次平滑样条

如果您不想要 csaps(x,y) 得到平滑程度，您可以将平滑参数 p 作为第三个可选参数输入。可以在 0 和 1 之间选取任意值。当 p 从 0 变到 1 时，平滑样条也相应地从一个极端，数据的最小二乘直线拟合，变到另一个极端，数据的“自然”三次样条插值。因为 csaps 实际上也会在第二个输出返回平滑参数，您可以做以下实验：

scs, p = csaps(x,noisy)
fnplt(scs,2)
hold("on")
fnplt(csaps(x,noisy,p/2)[1],2,"k--")
fnplt(csaps(x,noisy,(1+p)/2)[1],2,"r:")
plot(x,noisy,"o")
legend(["smoothing spline","more smoothed","less smoothed","noisy data"])
hold("off")

这会生成以下的图。

含噪数据更多或更少地平滑

有时，您可能更想只要在给定数据值以 norm(noisy-fnval(sp,x))^2 <= tol ，其中 tol 为给定的容差，条件下得到最平滑的三次样条 sp。您可以为您定义的容差 tol 使用指令 sp, = spaps(x,noisy,tol) 来生成这样的样条。

最小二乘

如果您更想要最小二乘逼近，您可以使用指令 sp = spap2(knots,k,x,y) 其中必须提供样条节点序列 knots 与阶数 k。

阶数通常选择为 4，这样可以生成一个三次样条。如果您不知道如何选择节点，可以简单指定您想要使用的多项式片段的数目。比如，

sp = spap2(3,4,x,y)

给出一个包含三个多项式片段的三次样条。如果结果误差不是均匀的，你可以使用 newknt 来生成一个更好的节点分布：

sp = spap2(newknt(sp)[1],4,x,y)