Cómo construir splines

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

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

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

Suavizado

Si los datos son con ruido, 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])

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

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

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

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

Aproximación por mínimos cuadrados

También podría utilizar la aproximación por mínimos cuadrados a los datos con ruido 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

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

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)

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

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 un número ilimitado 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])

A continuación puede ver la aproximación por mínimos cuadrados a valores con ruido 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])

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

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

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