Splines Cubicos

Programa de Splines Cúbicos

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

José Jeremías Caballero

Servicios de Programación en MatLab
https://programacioncaballero.blogspot.com/
Lima - Perú

Para mayor informacion, contactarse a jjcc94@hotmail.com