如何构造样条

---

这个示例展示了如何使用曲线拟合工具箱中的样条函数以各种方式构造样条。

插值

可以使用命令在位点 x 构造与余弦函数匹配的三次样条插值。

using TyPlot
using TyBase
using TyCurveFitting
using TyMath
x = 2 * pi * [0; 1; 0.1:0.2:0.9...]'
y = cos.(x)
cs = csapi(x, y)

使用 fnplt 查看插值样条。

fnplt(cs, 2)
axis([-1 7 -1.2 1.2])
hold("on")
plot(x, y, "o")
hold("off")

检查插值

余弦函数以 2*pi 为周期，那么我们的三次样条插值在这个层面上效果怎么样呢？一种检查的方法是计算其在两个端点处一阶导数的差。

diff(fnval(fnder(cs), [0 2 * pi]), dims=2)

ans = 1×1 Matrix{Float64}:

 -0.13747690654157635

如果需要这种周期性 请使用 caspe 而不是 csapi。

csp = csape(x, y, "periodic")
hold("on")
fnplt(csp, "g")
hold("off")

然后再检查，

diff(fnval(fnder(csp), [0 2 * pi]), dims=2)

ans = 1×1 Matrix{Float64}:

 1.9798345746301846e-17

现在就算是端点处的二阶导数都是匹配的。

diff(fnval(fnder(csp, 2), [0 2 * pi]), dims=2)

1×1 Matrix{Float64}:

 4.440892098500626e-16

对相同数据的分段线性插值可通过 spapi 获得。在这里，我们将其添加到上一个图中，以红色显示。

pl = spapi(2, x, y)
hold("on")
fnplt(pl, "r", 2)
hold("off")

平滑

如果数据含噪，您通常需要近似而不是插值。 例如，有了这些数据

x = LinRange(0, 2 * pi, 51)
rng = MersenneTwister(1234)
noisy_y = cos.(x) .+ 0.2 * (rand(rng, size(x, 1), 1) .- 0.5)
plot(x, noisy_y, "x")
axis([-1 7 -1.2 1.2])

插值将给出如下蓝色所示的摆动插值。

hold("on")
fnplt(csapi(x, noisy_y))
hold("off")

相反，以适当的公差进行平滑

tol = (0.05)^2 * (2 * pi)

tol = 0.015707963267948967

给出一个平滑的近似，如下图红色所示。

hold("on")
fnplt(spaps(x, noisy_y, tol)[1], "r", 2)
hold("off")

这个近似在接近区间末端的很差，而且远非周期。 为了保持周期性，先对周期延拓的数据进行近似，然后再将近似限制在原始区间。

noisy_y[[1 end]] .= mean(noisy_y[[1 end]])
lx = length(x)
lx2 = Int.(round(lx / 2))
range = Int.([lx2:lx... 2:lx... 2:lx2...])
sps, = spaps([x[lx2:lx] .- 2 * pi; x[2:lx]; x[2:lx2] .+ 2 * pi], noisy_y[range], 2 * tol)

这给出了更接近周期性的近似值，以黑色显示。

hold("on")
fnplt(sps, [0, 2 * pi], "k", 2)
hold("off")

最小二乘逼近

或者，您可以通过具有少自由度的样条对噪声数据使用最小二乘逼近。

例如，您可以尝试仅包含四个片段的三次样条。

spl2 = spap2(4, 4, x, noisy_y)
fnplt(spl2, "b", 2)
axis([-1 7 -1.2 1.2])
hold("on")
plot(x, noisy_y, "x")
hold("off")

节点选择

使用 spapi 或 spap2 时，通常必须指定特定的样条空间。 这是通过指定节点序列和阶数来完成的，这可能会有点问题。 但是，当使用 k 阶样条对 x,y 数据进行样条插值时，您可以使用函数 optknt 来提供良好的节点序列，如下例所示。

k = 5   # order 5, i.e., we are working with quartic splines
x = 2 * pi * sort([0 1 rand(1, 10)], dims=2)
y = cos.(x)
sp = spapi(optknt(x, k)', x, y)
fnplt(sp, 2, "g")
hold("on")
plot(x, y, "o")
hold("off")
axis([-1 7 -1.1 1.1])

在进行最小二乘逼近时，您可以使用当前逼近来借助 newknt 确定可能更好的节点选择。 例如，以下对指数函数的近似并不是那么好，这从它用红色绘制的误差中可以看出来。

x = LinRange(0, 10, 101)
y = exp.(x)
sp0 = spap2(augknt([0:2:10...], 4)[1], 4, x, y)
plot([x...], vec(y .- vec(fnval(sp0, x))), "r", linewidth=2)

但是，您可以使用该初始近似值来创建另一个具有相同节点数目但分布更好的近似值。 它的误差用黑色绘制。

sp1 = spap2(newknt(sp0)[1], 4, x, y)
hold("on")
plot([x...], y .- vec(fnval(sp1, x)), "k", linewidth=2)
hold("off")

网格数据

曲线拟合工具箱中的所有样条插值和近似命令也可以处理任意数量变量的网格数据。

例如，这里是墨西哥帽子函数的双三次样条插值。

x = 0.0001 .+ (-4:0.2:4)
y = -3:0.2:3
yy, xx = meshgrid2(y, x)
r = pi * sqrt.(xx .^ 2 + yy .^ 2)
z = sin.(r) ./ r
bcs = csapi((x, y), z)
fnplt(bcs)
axis([-5 5 -5 5 -0.5 1])

这是同一网格上同一函数的含噪数据的最小二乘近似值。

knotsx = augknt([LinRange(x[1], x[end], 21)...], 4)[1]
knotsy = augknt([LinRange(y[1], y[end], 15)...], 4)[1]
bsp2 = spap2((knotsx, knotsy), [4, 4], (x, y), z .+ 0.02 * (rand(size(z, 1), size(z, 2)) .- 0.5))
fnplt(bsp2)
axis([-5 5 -5 5 -0.5 1])

曲线

网格数据可以被轻松处理，因为曲线拟合工具箱可以处理向量值样条曲线。这也使得使用参数曲线变得容易。

例如，这里是一个近似无穷大，通过将三次样条曲线穿过下图中标记的点而获得。

t = 0:8
xy = transpose([0 0; 1 1; 1.7 0; 1 -1; 0 0; -1 1; -1.7 0; -1 -1; 0 0])
infty = csape(t, xy, "periodic")
fnplt(infty, 2)
axis([-2 2 -1.1 1.1])
hold("on")
plot(xy[1, :], xy[2, :], "o")
hold("off")

这是相同的曲线，但在三维中存在一个运动。

roller = csape(t, [xy; 0 1/2 1 1/2 0 1/2 1 1/2 0], "periodic")
fnplt(roller, 2, [0, 4], "b")
hold("on")
fnplt(roller, 2, [4, 8], "r")
plot3(0, 0, 0, "o")
hold("off")

曲线的两半用不同的颜色绘制，并标记了原点，以帮助可视化这条双翼空间曲线。

曲面

具有 R^3 值的二元张量积样条曲线给出曲面。 例如，这是一个很好的甜甜圈近似值。

x = 0:4
y = -2:2
R = 4
r = 2
v = zeros(3, 5, 5)
v[3, :, :] = transpose([0 (R - r) / 2 0 (r - R) / 2 0]) * [1 1 1 1 1]
v[2, :, :] = transpose([R (r + R) / 2 r (r + R) / 2 R]) * [0 1 0 -1 0]
v[1, :, :] = transpose([R (r + R) / 2 r (r + R) / 2 R]) * [1 0 -1 0 1]
dough0 = csape((x, y), v, "periodic")
fnplt(dough0)
axis("equal")
axis("off")

这是该表面的法线冠。

nx = 43
xy = [ones(1, nx); [LinRange(2, -2, nx)...]']
points = fnval(dough0, xy)'
ders = fnval(fndir(dough0, eye(2)), xy);
normals = mapslices(x -> cross(x[1:3], x[4:6]), [ders[4:6, :]; ders[1:3, :]], dims=1)
normals = (normals ./ repeat(sqrt.(sum(normals .* normals, dims=1)), 3, 1))'
pn = [points; points + normals]
hold("on")
for j = 1:nx
    plot3(pn[[j, j + nx], 1], pn[[j, j + nx], 2], pn[[j, j + nx], 3])
end
hold("off")

最后，这是它在 (x,y) 平面上的投影。

fnplt(fncmb(dough0, [1 0 0; 0 1 0]))
axis([-5.25 5.25 -4.14 4.14])
axis("off")

散点数据

也可以插值到平面中未网格化的数据位点给出的值。例如，考虑将单位正方形平滑映射到单位圆盘。我们构建了数据值，标记为圆圈，以及相应的数据位点，标记为 x。 每个数据位点都通过一个箭头与其对应值相连。

n = 64
t = LinRange(0, 2 * pi, n + 1)
t = t[1:end-1]'
values = [cos.(t); sin.(t)]
plot(values[1, :], values[2, :], "or")
axis("equal")
axis("off")

sites = values ./ repeat(mapslices(maximum, abs.(values), dims=1), 2, 1)
hold("on")
plot(sites[1, :], sites[2, :], "xk")
quiver(sites[1, :], sites[2, :], values[1, :] - sites[1, :], values[2, :] - sites[2, :], length=0.9)
hold("off")

然后使用 tpaps 构造一个二元插值向量值薄板样条。

st = tpaps(sites, values, 1)[1]

样条曲线确实将单位正方形平滑地（近似地）映射到单位圆盘，正如其通过 fnplt 绘制的图所示。 该图显示了均匀间隔的方形网格的图像.

hold("on")
fnplt(st)
hold("off")