代码之家 › 专栏 › 技术社区 › Paul

Scilab中的样条系数

scilab spline interpolation

1

Paul · 技术社区 · 9 年前

我想将分段三次样条拟合到一组大数据中。我认为我不需要B样条曲线,因为我希望样条曲线精确地穿过数据点。我能在Scilab中看到的唯一方法是 splin 和 interp .

然而,我想计算每个样条曲线在节点之间的系数(因为我需要取这些系数并将它们放在不同的软件中)。全部的 花键 给你的是导数。有办法得到三次样条系数吗?或者有没有一种方法可以很容易地从一阶导数中生成系数?

1 回复 | 直到 9 年前

1

user3717023 user3717023 9 年前

是的,可以从您的y值和 splin 每个区间[x(i),x(i+1)]有4个系数要求,4个方程:两端的值,两端的导数。最直接的方法是告诉Scilab为每个子区间求解这个4乘4系统:这不应该比样条曲线本身的求值花费太多时间。下面的程序可以做到这一点。

x = [0,1,2,3,4,5]   // x values
y = [1,0,1,0,1,0]   // y values
d = splin(x,y)      
n = length(x)-1     // number of subintervals
cfs = zeros(4,n)    // matrix to store coefficients in
for i=1:n
    a = x(i)
    b = x(i+1)
    cfs(:,i) = [1,a,a^2,a^3; 1,b,b^2,b^3; 0,1,2*a,3*a^2; 0,1,2*b,3*b^2] \ [y(i);y(i+1);d(i);d(i+1)]
end

前两个方程式 1,a,a^2,a^3; 1,b,b^2,b^3 将多项式的值与y值相关联;另外两个 0,1,2*a,3*a^2; 0,1,2*b,3*b^2 对其衍生物也同样如此。(公式只是x的幂的导数)

上述脚本的输出:

    1.     1.   - 8.6    13.     13.   
  - 3.4  - 3.4    11.  - 10.6  - 10.6  
    3.1    3.1  - 4.1    3.1     3.1   
  - 0.7  - 0.7    0.5  - 0.3   - 0.3

每列有四个系数:例如,样条曲线的第一段是 1-3.4x+3.1x^2-0.7x^3 。由于这是非结样条线 花键 命令,第二块与第一块相同;最后一个与倒数第二个相同。

您可以通过绘制零件来检查这是否正确工作:

for i=1:n  
    t = linspace(x(i),x(i+1))
    plot(t,cfs(:,i)'*[ones(t); t; t.^2; t.^3])
end

也就是说,将形成样条曲线的多项式表示为

 p(x) = y(i) + A*(x-x(i)) + B*(x-x(i))*(x-x(i+1)) + C*(x-x(i))^2*(x-x(i+1))

其中系数很容易逐个找到,而不需要求解线性系统:

A = (y(i+1)-y(i))/(x(i+1)-x(i)) 通过将x(i+1)处的值相等
B = (d(i)-A)/(x(i)-x(i+1)) ,通过在x(i)处求导数
C = (d(i+1)-A-B*(x(i+1)-x(i)))/(x(i+1)-x(i))^2 ,通过在x(i+1)处求导数

当然,这些系数应该与上述适当的多项式一起使用。这是另一个版本

for i=1:n
    A = (y(i+1)-y(i))/(x(i+1)-x(i))
    B = (d(i)-A)/(x(i)-x(i+1))
    C = (d(i+1)-A-B*(x(i+1)-x(i)))/(x(i+1)-x(i))^2
    cfs(:,i) = [y(i);A;B;C]
end
// Again, plot for testing
for i=1:n
    t = linspace(x(i),x(i+1))
    plot(t,cfs(:,i)'*[ones(t); t-x(i); (t-x(i)).*(t-x(i+1)); ((t-x(i)).^2).*(t-x(i+1))])
end