Algoritmo en MATLAB para trazador cúbico.


Por: Carlos Armando De Castro P.


El siguiente algoritmo recibe un número arbitrario de pares de datos en la forma de una matriz de 2*n, donde las abcisas se encuentran en la primera fila (o renglón) y las ordenadas en la segunda fila de la matriz, devolviendo los coeficientes de los polinomios cúbicos que interpolan a los datos y además entrega la gráfica del trazador cúbico ("spline" de orden 3):

 

 

function [a,b,c,d]=spline3(X)

%Pasos básicos del algoritmo obtenidos del libro Análisis Numérico de
%Richard Burden, 2a. Edición, Grupo Editorial Iberoamérica.

n=length(X(1,:));

for i=1:n;
    a(i)=X(2,i);
end

for i=1:n-1;
    h(i)=X(1,i+1)-X(1,i);
end

for i=2:n-1;
    alfa(i)=3/h(i)*(a(i+1)-a(i))-3/h(i-1)*(a(i)-a(i-1));
end

l(1)=1;
mu(1)=0;
z(1)=0;

for i=2:n-1;
    l(i)=2*(X(1,i+1)-X(1,i-1))-h(i-1)*mu(i-1);
    mu(i)=h(i)/l(i);
    z(i)=(alfa(i)-h(i-1)*z(i-1))/l(i);
end

l(n)=1;
z(n)=0;
c(n)=0;

for i=n-1:-1:1;
    c(i)=z(i)-mu(i)*c(i+1);
    b(i)=(a(i+1)-a(i))/h(i)-h(i)*(c(i+1)+2*c(i))/3;
    d(i)=(c(i+1)-c(i))/(3*h(i));
end

for i=1:n-1;
    x=X(1,i):0.1:X(1,i+1);
    y=a(i)+b(i)*(x-X(1,i))+c(i)*(x-X(1,i)).^2+d(i)*(x-X(1,i)).^3;
    hold on;
    plot(x,y,'b');
end

for i=1:n;
    hold on;
    plot (X(1,i),X(2,i),'*','MarkerEdgeColor','r','LineWidth',1);
    title('Interpolación por "splines" de orden 3.');
end

 

 

Por ejemplo, para los datos {(1,0),(2,3),(3,4),(4,-6),(5,2),(6,4),(7,0),(8,4),(9,3)}, se escribe en el Command Window:

>>X=[1 2 3 4 5 6 7 8 9; 0 3 4 -6 2 4 0 4 3];
>>[m,b]=spline3(X)

Y el programa entrega los resultados:

a = 0 3 4 -6 2 4 0 4 3

b = 2.3806 4.2388 -7.3357 -1.8960 8.9196 -3.7826 0.2107 2.9398

c = 0 1.8582 -13.4326 18.8723 -8.0567 -4.6455 8.6388 -5.9097 0

d = 0.6194 -5.0969 10.7683 -8.9763 1.1371 4.4281 -4.8495 1.9699