scipy - 使用 scipy 从结和系数创建样条

Question

我正在尝试从论文中重现一个函数，该函数仅根据样条结和系数指定。在stackoverflow 上找到这个后，给定一个 scipy 插值对象，从它的结和系数，我可以重新创建 scipy 插值。但是，对于论文中指定的功能，该方法失败了。要重现 scipy 插值，我可以这样做：

using PyCall, PyPlot, Random
Random.seed!(5)
sp = pyimport("scipy.interpolate")
x = LinRange(0,1,50)
y = (0.9 .+ 0.1rand(length(x))).*sin.(2*pi*(x.-0.5))
t = collect(x[2:2:end-1]) # knots
s1 = sp.LSQUnivariateSpline(x, y, t)
x2 = LinRange(0, 1, 200) # new x-grid
y2 = s1(x2) # evaluate spline on that new grid
figure()
plot(x, y, label="original")
plot(x2, y2, label="interp", color="k")
knots = s1.get_knots()
c = s1.get_coeffs()
newknots(knots, k) = vcat(fill(knots[1],k),knots,fill(knots[end],k)) # func for boundary knots of order k
forscipyknots = newknots(knots, 3)
s2 = sp.BSpline(forscipyknots, c, 3)
y3 = s2(x2)
plot(x2,y3,"--r", label="reconstructed \nfrom knots and coeff")
legend()

它按预期提供以下内容：

在尝试重现具有指定和应该产生的 knots = [.4,.4,.4,.4,.7]系数的函数（下图）时：c = [2,-5,5,2,-3,-1,2]

使用下面的代码和上面的结和系数：

knots = [.4,.4,.4,.4,.7]
c = [2,-5,5,2,-3,-1,2]
forscipyknots = newknots(knots, 3)
s2 = sp.BSpline(forscipyknots, c, 3)
figure()
plot(x2, s2(x2))

我得到以下（下）。我确定我弄乱了边界结 - 我该如何解决这个问题？

标签： scipyjuliainterpolation

Answer 1

简短的回答：

内结序列t=[0.4,0.4,0.4,0.4,0.7]和参数c=[2,-5,5,2,-3,-1,2]不允许构造样条，示例包含错误（稍后会详细介绍）。最好的办法是移除其中一个0.4结并构造一个二次（二次）样条曲线，如下所示

tt = [0.0,0.0,0.0,0.4,0.4,0.4,0.7,1.0,1.0,1.0]
c  = [2,-5,5,2,-3,-1,2]
s2 = BSpline(tt,c,2)

这将产生以下图表

长答案：

示例 3中的节点序列仅包含内部节点，因此，您需要添加边界节点。由于您要在区间上评估样条，[0,1]因此完整的结序列需要覆盖点 0 和 1。最简单的方法是将 0 添加到序列的开头，将 1 添加到序列的末尾，并根据需要根据需要复制它们的样条。三次（三次）样条曲线需要四个边界节点（即四个零和四个一），二次样条曲线需要三个边界节点（三个零和三个一）。

然而，有一个问题。三次样条曲线需要 9 个参数，而示例 3 只给您 7 个参数。因此，您不能由此构造三次样条曲线。使用给定的七个参数，您可以构建二次样条，但是，问题在于对于二次样条，每个点在内结序列中最多只能出现 3 次。并且 0.4 出现四次（这表明三次样条）。因此，您所能做的就是移除其中一个 0.4 节并构造一个二阶样条，如上面的简短答案中所示。

现在我将解释你做错了什么。在第一个示例中，您使用 获得了现有样条曲线的结序列knots = s1.get_knots()，这给了您knots=[0,0.02,0.04,...,0.98,1]。该序列包含边界节点 0 和 1（尽管只有一次）。因此，为了构造三次样条，您将这三个样条中的每一个都复制了 3 次以获得forscipyknots = [0,0,0,0,0.02,0.04,...,0.98,1,1,1,1]. 到目前为止，一切都很好。

然而，在示例 3中，节点序列不包含边界点。正如您之前所做的那样，您最终将 0.4 和 0.7 节复制了 3 次，结果为forscipyknots = [0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.7,0.7,0.7,0.7]. 您不能在此序列上构造样条曲线，由此产生的任何东西都不是样条曲线。相反，您需要的是forscipyknots = [0.0,0.0,0.0,0.0,0.4,0.4,0.4,0.4,0.7,1.0,1.0,1.0,1.0]（这不会起作用，因为您没有足够的系数；但是您可以自己尝试一下，例如c = [1,2,-5,5,2,-3,-1,2,1]）。为此，您需要将 0 添加到数组的开头，将 1 添加到数组的末尾，然后才使用您的newknots函数。

举个例子，三次样条曲线可能看起来像这样

tt = [0.0,0.0,0.0,0.0,0.4,0.4,0.4,0.4,0.7,1.0,1.0,1.0,1.0]
c  = [1,2,-5,5,2,-3,-1,2,1]
s2 = BSpline(tt,c,3)