Сплайн-коэффициенты в Scilab

Question

Я хочу поместить кусочек сплайна кусочка в большой набор данных. Я не думаю, что хочу B-сплайны, потому что я хочу, чтобы сплайны проходили точно через точки данных. Единственный способ, которым я могу это сделать в Scilab, - это splin и interp.

Тем не менее, я хочу, чтобы коэффициенты каждой части сплайна между узлами (потому что мне нужно взять эти коэффициенты и поместить их в другую часть программного обеспечения). Все splin дает вам производные. Есть ли способ получить кубические сплайн-коэффициенты? Или есть способ легко получить коэффициенты из первых производных?

Answer 1

Да, можно получить коэффициенты из значений 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 связывающие значения полинома у-значений; другие два 0,1,2*a,3*a^2; 0,1,2*b,3*b^2 делают то же самое для своей производной. (Формулы являются лишь производные степеням х.)

Выход вышеуказанного сценария:

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. Поскольку это не-узел-сплайн, по умолчанию используется команда splin, вторая часть такая же, как и первая; и последний такой же, как и второй.

Вы можете проверить, что это работает правильно, откладывая части:

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)) путем приравнивания значение в точке х (г + 1)
• 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, путем приравнивания производной в точке х (я + 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