Scilab 中的样条系数

Spline coefficients in Scilab

我想用分段三次样条拟合大量数据。我不认为我想要 B 样条曲线,因为我希望样条曲线准确地通过数据点。我在 Scilab 中看到的唯一方法是使用 splininterp

但是,我想要结点之间每段样条的系数(因为我需要获取这些系数并将它们放在不同的软件中)。 splin 给你的只是衍生品。有没有办法获得三次样条系数?或者有没有一种方法可以轻松地从一阶导数生成系数?

是的,可以从您拥有的 y-values 和 splin 返回的导数值中获取系数。每个区间 [x(i), x(i+1)] 你有 4 个系数和 4 个方程:两端的值,两端的导数。最直接的方法就是告诉 Scilab 为每个子区间解决这个 4 x 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-values 相关联;另外两个 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。由于这是 not-a-knot 样条,splin 命令的默认模式,第二部分与第一部分相同;最后一个与 second-to-last.

相同

您可以通过绘制片段来检查它是否正常工作:

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