使用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))

我得到的是下面的结果 我确定我把边界结搞乱了--如何解决这个问题？

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)

这就产生了以下图表

长答案。

结号顺序是： 例三 只包含内结，因此，你需要添加边界结。由于你想在区间上评估花键，所以你需要添加边界结。[0,1] 最简单的方法是在序列的开始和结束处分别加上0和1，并根据花键的度数进行必要的复制。三次元（三度）花键需要四个边界结（即四个零和四个一），二次元花键需要三个边界结（三个零和三个一）。

然而，有一个问题。三次方花键需要9个参数，而例3只给了你7个参数，因此，你不能用这个来构造三次方花键。用给定的7个参数，你可以构造一个二次花键，但是，问题是，对于二次花键，每个点在内结序列中最多只能出现三次。而0.4出现了四次（这就说明是三次花键）。因此，你能做的就是去掉其中一个0.4的结，然后构建一个二度花键，就像上面的简答一样。

现在我将解释你做错了什么。在第一个例子中，你从一个现有的花键中获得了结点序列，使用了 knots = s1.get_knots()这给了你 knots=[0,0.02,0.04,...,0.98,1]. 这个序列包含边界结0和1（虽然只有一次）。因此，要构造一个立方花键，你要把这三个序列分别复制三次，得到的是 forscipyknots = [0,0,0,0,0.02,0.04,...,0.98,1,1,1,1]. 到目前为止还不错。

在 例三然而，结点序列并不包含边界点。由于你和之前一样，你最终将0.4和0.7的结重复了三次，结果是 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)