%%
http://www.lawebdelprogramador.com/foros/Matlab/1327174-Duda_spline.html
% Hola,
estoy haciendo un .m para obtener y dibujar las splines. Primero
%introduces
los puntos (x,y) que tu quieres.
%
% El
siguiente código te permite obtener los coeficientes para cada spline:
% De esta manera, se obtiene
%
% [ A1 B1 C1 D1; A2 B2 C2 D2;............; AN BN CN
DN]
%
% El
problema, es que aún teniendo esta matriz, no sé cómo hacer para
%representar
las funciones spline (correspondiente a cada intervalo).
%Primero
hay que introducir el vecto X e Y
clear all
X=[40 45 50
55 60 65
70];
Y=[390 340 290 250 210 180
160];
N=length(X)-1;
H=diff(X);
E=diff(Y)./H;
diagsupinf=H(2:N-1);
diagprinc=2*(H(1:N-1)+H(2:N));
U=6*diff(E);
%
restricciones del spline natural
Minicial=0;
Mfinal=0;
%
construccion de la matriz de los coeficientes para el calculo de los
%
momentos
A=diag(diagprinc)+diag(diagsupinf,1)+diag(diagsupinf,-1);
%
construccion del vector independiente para el calculo de los momentos
B=U';
B(1)=B(1)-H(1)*Minicial; B(N-1)=B(N-1)-H(N)*Mfinal;
%
resolvemos el sistema y hallamos los momentos
M=A\B;
%
añadimos los dos momentos naturales extremos y lo escribimos como vector
% fila
M=[Minicial,M',Mfinal];
%
calculo de los coeficientes del polinomio cubico i-esimo en potencias
% de (x-x_i)
for i=1:N
S(i,1)=(M(i+1)-M(i))/(6*H(i));
S(i,2)=M(i)/2;
S(i,3)=E(i)-H(i)*(M(i+1)+2*M(i))/6;
S(i,4)=Y(i);
end
%
creamos los polinomios para cada
intervalo
syms x
for i=1:N
s=0;
for j=1:size(S,2)
s=s+expand(S(i,j)*(x-X(j))^(4-j));
end
Es(i,:)=s;
end
%
Es(1)--->f1---> Intervalo X [40,45)
% Es(2)--->f2---> Intervalo X [45,50)
% Es(3)--->f3---> Intervalo X [50,55)
% Es(4)--->f4---> Intervalo X [55,60)
% Es(5)--->f5---> Intervalo X [60,65)
% Es(6)--->f6---> Intervalo X [65,70)
%dibujamos los polinomios
figure(gcf)
hold on
c='rgbcmk';
for i=1:length(X)-1
F=vectorize(inline(Es(i)));
plot(linspace(X(i),X(i+1)),F(linspace(X(i),X(i+1))),c(i),'linewidth',5)
%fplot(F,[X(i) X(i+1)],'r',);
end
plot(X,Y,'linewidth',5,'color',rand(1,3))
axis([X(1) X(end) -200 500])
grid on
hold off
|