Scilab 中的样条系数
Spline coefficients in Scilab
我想用分段三次样条拟合大量数据。我不认为我想要 B 样条曲线,因为我希望样条曲线准确地通过数据点。我在 Scilab 中看到的唯一方法是使用 splin
和 interp
。
但是,我想要结点之间每段样条的系数(因为我需要获取这些系数并将它们放在不同的软件中)。 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
我想用分段三次样条拟合大量数据。我不认为我想要 B 样条曲线,因为我希望样条曲线准确地通过数据点。我在 Scilab 中看到的唯一方法是使用 splin
和 interp
。
但是,我想要结点之间每段样条的系数(因为我需要获取这些系数并将它们放在不同的软件中)。 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