Main Content

La traducción de esta página aún no se ha actualizado a la versión más reciente. Haga clic aquí para ver la última versión en inglés.

Cómo construir splines

Este ejemplo muestra cómo construir splines de varias formas utilizando las funciones para splines de Curve Fitting Toolbox™.

Interpolación

Puede construir una interpolación por splines cúbicos que coincida con la función de coseno en los siguientes sitios x, utilizando el comando csapi.

x = 2*pi*[0 1 .1:.2:.9];
y = cos(x);
cs = csapi(x,y);

A continuación, puede ver el spline de interpolación utilizando fnplt.

fnplt(cs,2);
axis([-1 7 -1.2 1.2])
hold on
plot(x,y,'o')
hold off

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers

Comprobar la interpolación

La función de coseno es 2*pi periódico. ¿Cómo es de efectiva nuestra interpolación por splines cúbicos en ese sentido? Una forma de comprobarlo es calcular la diferencia en la primera derivada en los dos extremos de línea.

diff( fnval( fnder(cs), [0 2*pi] ) )
ans = -0.1375

Para aplicar periodicidad, utilice csape en lugar de csapi.

csp = csape( x, y, 'periodic' );
hold on
fnplt(csp,'g')
hold off

Figure contains an axes object. The axes object contains 3 objects of type line. One or more of the lines displays its values using only markers

Ahora la comprobación da

diff( fnval( fnder(csp), [0 2*pi] ) )
ans = -2.2806e-17

Ahora incluso la segunda derivada coincide en los extremos de línea.

diff( fnval( fnder(csp, 2), [0 2*pi] ) )
ans = -2.2204e-16

La interpolación lineal por tramos para los mismos datos está disponible con spapi. Aquí la añadimos a la gráfica anterior, en color rojo.

pl = spapi(2, x, y);
hold on
fnplt(pl, 'r', 2)
hold off

Figure contains an axes object. The axes object contains 4 objects of type line. One or more of the lines displays its values using only markers

Suavizado

Si los datos son ruidosos, normalmente buscamos aproximar en lugar de interpolar. Por ejemplo, con estos datos

x = linspace(0,2*pi,51);
noisy_y = cos(x) + .2*(rand(size(x))-.5);
plot(x,noisy_y,'x')
axis([-1 7 -1.2 1.2])

Figure contains an axes object. The axes contains a line object which displays its values using only markers.

la interpolación resultante sería serpenteante, como se muestra a continuación en azul.

hold on
fnplt( csapi(x, noisy_y) )
hold off

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers

En cambio, el suavizado con una tolerancia adecuada

tol = (.05)^2*(2*pi)
tol = 0.0157

ofrece una aproximación suavizada, que se muestra a continuación en rojo.

hold on
fnplt( spaps(x, noisy_y,  tol), 'r', 2 )
hold off

Figure contains an axes object. The axes object contains 3 objects of type line. One or more of the lines displays its values using only markers

La aproximación es mucho peor cerca de los extremos del intervalo, y está lejos de ser periódica. Para aplicar periodicidad, aproxime a datos ampliados periódicamente y, a continuación, restrinja la aproximación al intervalo original.

noisy_y([1 end]) = mean( noisy_y([1 end]) );
lx = length(x);
lx2 = round(lx/2);
range = [lx2:lx 2:lx 2:lx2];
sps = spaps([x(lx2:lx)-2*pi x(2:lx) x(2:lx2)+2*pi],noisy_y(range),2*tol);

Así, se obtiene la aproximación más periódica, que se muestra en negro.

hold on
fnplt(sps, [0 2*pi], 'k', 2)
hold off

Figure contains an axes object. The axes object contains 4 objects of type line. One or more of the lines displays its values using only markers

Aproximación por mínimos cuadrados

También podría utilizar la aproximación por mínimos cuadrados a los datos ruidosos con un spline con pocos grados de libertad.

Por ejemplo, puede probar con un spline cúbico con solo cuatro tramos.

spl2 = spap2(4, 4, x, noisy_y);
fnplt(spl2,'b',2);
axis([-1 7 -1.2 1.2])
hold on
plot(x,noisy_y,'x')
hold off

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers

Selección de nudos

Cuando utiliza spapi o spap2, normalmente tiene que especificar un espacio de splines concreto. Se hace especificando una secuencia de nudos y un orden, lo que podría suponer un problema. Sin embargo, cuando hace una interpolación por splines de los datos x,y utilizando un spline de orden k, puede utilizar la función optknt para proporcionar una secuencia de nudos buena, como en el siguiente ejemplo.

k = 5;   % order 5, i.e., we are working with quartic splines
x = 2*pi*sort([0 1 rand(1,10)]);
y = cos(x);
sp = spapi( optknt(x,k), x, y );
fnplt(sp,2,'g');
hold on
plot(x,y,'o')
hold off
axis([-1 7 -1.1 1.1])

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers

Cuando hace una aproximación por mínimos cuadrados, puede utilizar la aproximación actual para determinar una selección de nudos que posiblemente sea mejor con ayuda de newknt. Por ejemplo, la siguiente aproximación a la función exponencial no es tan buena, como demuestra su error, representado en color rojo.

x = linspace(0,10,101);
y = exp(x);
sp0 = spap2( augknt(0:2:10,4), 4, x, y );
plot(x,y-fnval(sp0,x),'r','LineWidth',2)

Figure contains an axes object. The axes object contains an object of type line.

Sin embargo, puede utilizar esa aproximación inicial para crear otra con el mismo número de nudos, pero mejor distribuidos. Su error se representa en negro.

sp1 = spap2( newknt(sp0), 4, x, y );
hold on
plot(x,y-fnval(sp1,x),'k','LineWidth',2)
hold off

Figure contains an axes object. The axes object contains 2 objects of type line.

Datos de malla

Todos los comandos de interpolación por splines y de aproximación de Curve Fitting Toolbox también pueden gestionar datos de malla, en cualquier número de variables.

Por ejemplo, a continuación puede ver una interpolación por splines bicúbicos para la función sombrero mexicano.

x =.0001+(-4:.2:4);
y = -3:.2:3;
[yy,xx] = meshgrid(y,x);
r = pi*sqrt(xx.^2+yy.^2);
z = sin(r)./r;
bcs = csapi({x,y}, z);
fnplt(bcs)
axis([-5 5 -5 5 -.5 1])

Figure contains an axes object. The axes object contains an object of type surface.

A continuación puede ver la aproximación por mínimos cuadrados a valores ruidosos de esa misma función en la misma cuadrícula.

knotsx = augknt(linspace(x(1), x(end), 21), 4);
knotsy = augknt(linspace(y(1), y(end), 15), 4);
bsp2 =  spap2({knotsx,knotsy},[4 4], {x,y},z+.02*(rand(size(z))-.5));
fnplt(bsp2)
axis([-5 5 -5 5 -.5 1])

Figure contains an axes object. The axes object contains an object of type surface.

Curvas

Gestionar los datos de malla resulta fácil porque Curve Fitting Toolbox puede abordar splines con valor vectorial. Así, también resulta más fácil trabajar con curvas paramétricas.

En este caso, por ejemplo, puede ver una aproximación al infinito, obtenida pasando una curva de spline cúbico por los puntos marcados en la siguiente figura.

t = 0:8;
xy = [0 0;1 1; 1.7 0;1 -1;0 0; -1 1; -1.7 0; -1 -1; 0 0].';
infty = csape(t, xy, 'periodic');
fnplt(infty, 2)
axis([-2 2 -1.1 1.1])
hold on
plot(xy(1,:),xy(2,:),'o')
hold off

Figure contains an axes object. The axes object contains 2 objects of type line. One or more of the lines displays its values using only markers

Aquí puede ver la misma curva, pero con movimiento en una tercera dimensión.

roller = csape( t , [ xy ;0 1/2 1 1/2 0 1/2 1 1/2 0], 'periodic');
fnplt( roller , 2, [0 4],'b' )
hold on
fnplt( roller, 2, [4 8], 'r')
plot3(0,0,0,'o')
hold off

Figure contains an axes object. The axes object contains 3 objects of type line. One or more of the lines displays its values using only markers

Las dos mitades de la curva se representan en distintos colores y se marca el origen para ayudarle a visualizar esta curva con dos alas en el espacio.

Superficies

Los splines bivariados de producto tensorial con valores en R^3 dan como resultado superficies. Por ejemplo, a continuación puede ver una buena aproximación a un toro.

x = 0:4;
y = -2:2;
R = 4;
r = 2;
v = zeros(3,5,5);
v(3,:,:) = [0 (R-r)/2 0 (r-R)/2 0].'*[1 1 1 1 1];
v(2,:,:) = [R (r+R)/2 r (r+R)/2 R].'*[0 1 0 -1 0];
v(1,:,:) = [R (r+R)/2 r (r+R)/2 R].'*[1 0 -1 0 1];
dough0 = csape({x,y},v,'periodic');
fnplt(dough0)
axis equal, axis off

Aquí puede ver una corona de valores normales para esa superficie.

nx = 43;
xy = [ones(1,nx); linspace(2,-2,nx)];
points = fnval(dough0,xy)';
ders = fnval(fndir(dough0,eye(2)),xy);
normals = cross(ders(4:6,:),ders(1:3,:));
normals = (normals./repmat(sqrt(sum(normals.*normals)),3,1))';
pn = [points;points+normals];
hold on
for j=1:nx
   plot3(pn([j,j+nx],1),pn([j,j+nx],2),pn([j,j+nx],3))
end
hold off

Finalmente, aquí puede ver su proyección en el plano (x,y).

fnplt(fncmb(dough0, [1 0 0; 0 1 0]))
axis([-5.25 5.25 -4.14 4.14]), axis off

Datos dispersos

También es posible interpolar a valores proporcionados en sitios de datos sin malla en el plano. Considere, por ejemplo, la tarea de aplicar correctamente la unidad cuadrada al disco unidad. Construimos los valores de los datos, marcados como círculos, y los sitios de datos correspondientes, marcados como x. Cada sitio de datos se conecta con su valor asociado mediante una flecha.

n = 64;
t = linspace(0,2*pi,n+1); t(end) = [];
values = [cos(t); sin(t)];
plot(values(1,:),values(2,:),'or')
axis equal, axis off

sites = values./repmat(max(abs(values)),2,1);
hold on
plot(sites(1,:),sites(2,:),'xk')
quiver(sites(1,:),sites(2,:), ...
       values(1,:)-sites(1,:), values(2,:)-sites(2,:))
hold off

A continuación, use tpaps para construir un spline bivariado de thin-plate de interpolación y con valor vectorial.

st = tpaps(sites, values, 1);

En efecto, el spline asigna la unidad cuadrada correctamente (aproximadamente) al disco unidad, como indica su gráfica mediante fnplt. La gráfica muestra la imagen de una cuadrícula espaciada de manera uniforme siguiendo la aplicación de splines en st.

hold on
fnplt(st)
hold off