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