Coefficients de spline dans Scilab

Question

Je souhaite adapter une spline cubique par morceaux à un grand nombre de données. Je ne pense pas que je veuille B-splines parce que je veux que les splines traversent exactement les points de données. La seule façon que je peux voir pour le faire dans Scilab est avec splin et interp.

Cependant, je veux des coefficients de chaque morceau de la spline entre les points de noeud (parce que j'ai besoin de prendre ces coefficients et de les mettre dans un logiciel différent). Tous splin vous donne les dérivés. Y a-t-il un moyen d'obtenir les coefficients splines cubiques? Ou existe-t-il un moyen de générer facilement les coefficients à partir des premiers dérivés?

Answer 1

Oui, on peut obtenir des coefficients à partir des valeurs y que vous avez et des valeurs dérivées retournées par splin. Un chaque intervalle [x (i), x (i + 1)] vous avez 4 coefficients à trouver, et 4 équations: des valeurs aux deux extrémités, des dérivées aux deux extrémités. Le plus simple est de dire à Scilab de résoudre ce système 4 par 4 pour chaque sous-intervalle: cela ne devrait pas prendre beaucoup plus de temps que l'évaluation de la spline elle-même. Le programme ci-dessous le fait.

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 

Les deux premières équations 1,a,a^2,a^3; 1,b,b^2,b^3 se rapportent les valeurs de polynôme des valeurs de y; les deux autres 0,1,2*a,3*a^2; 0,1,2*b,3*b^2 font de même pour sa dérivée. (Les formules ne sont que des dérivés de puissances de x.)

La sortie du script ci-dessus:

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 

Chaque colonne a quatre coefficients: par exemple, la première pièce de la spline est 1-3.4x+3.1x^2-0.7x^3. Comme il s'agit d'une spline not-a-knot, le mode par défaut de la commande splin, le second élément est le même que le premier; et le dernier est le même que l'avant-dernier.

Vous pouvez vérifier que cela fonctionne correctement en traçant les pièces:

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

---

Cela dit, il serait plus facile de représenter les polynômes formant la spline sous la forme

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

où les coefficients sont faciles à trouver un par un, sans résoudre un système linéaire:

• A = (y(i+1)-y(i))/(x(i+1)-x(i)) en assimilant la valeur à x (i + 1)
• B = (d(i)-A)/(x(i)-x(i+1)), en mettant en équation dérivée en x (i)
• C = (d(i+1)-A-B*(x(i+1)-x(i)))/(x(i+1)-x(i))^2, en mettant en équation dérivée en x (i + 1)

De Bien sûr, ces coefficients doivent être utilisés avec les polynômes appropriés comme ci-dessus. Voici cette version alternative

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