Interpolación a Partir de Placa Delgada: Códigos en Matlab
 Programa Principal y Rutinas Auxiliares (comprimido)
 

 

Programa Principal

%Programa ejemplo del método de interpolación multivariada con funciones base radiales 
%usando el kernel r^4*log(r) y la s(u,v) = 5/4 + cos(5.4y)/6+6(3x-1)^2.
% 
%Parámetros: 
%Entrada     : Entero N - Número de nodos aleatorios para hacer la interpolación. 
%                   por ejemplo  mirbf_thinplate(100) 
%Salida      :  Vector de aproximación de la función s(u,v) en los N puntos aleatorios.  
%
%Autor       :  Daniel A. Cervantes Cabrera, danielcc@rv.unam.mx. 
%Supervisión : Pedro González Casanova, pedrogc@dgsca2.unam.mx.

 

function z_aprox = mirbf_thinplate(N)
    
%kernel r^4*log(r)
func  = 'phitps4';
 
 
%Generación de números aleatorios en el intervalo 
[x y] = gna(0,1,0,1,N);
myc   = [];
 
%Vector  solucion.
sol = solucion(x,y); 
sol(N+1:N+6) = 0;
 
 
 
%Creo la matriz 
 for i=1:N
    for j=1:N
      r= sqrt((x(i)-x(j))*(x(i)-x(j))+(y(i)-y(j))*(y(i)-y(j)));       
      A(i,j) = feval(func,r,myc,0); 
    end
 end
 
 for i=1:N
     for j=1:6
       switch j
           case 1,
             aij = x(i)^2;
           case 2,
             aij = y(i)^2;
           case 3,
             aij = x(i)*y(i);
           case 4,
               aij = x(i);
           case 5,
               aij = y(i);
           case 6,
               aij = 1;
       end
       P(i,j) = aij;
     end
 end
 
 
 
O = zeros(6,6);
 
M =[A P; P' O];
 
 
%Resuelvo el sistema  algebráico
lambda = M\(sol');
 
%Gr'afica del interpolante 
xm = linspace(0,1,101);
ym = linspace(0,1,101);
[t,nxym] = size(xm);
 
for i=1:nxym
    for j=1:nxym
        zm(i,j) = interp(xm(j),ym(i),x,y,lambda,func);
    end
end
 
hold on
[xm2d,ym2d] = meshgrid(0:.01:1);
mesh(xm2d,ym2d,zm);
led = sprintf('Interpolacion.');
title(led,'fontsize',18);
xlabel(sprintf('Condicion %f',cond(M)));
 
%Grafica de los puntos aleatorios
for i=1:N
    za(i) = solucion(x(i),y(i));
    plot3(x(i),y(i),za(i),'ro'); hold on
    plot(x(i),y(i),'r.');
end
 
 
%Gr'afica de la soluci'on.
figure
[xsol,ysol] = meshgrid(0:.01:1);
zsol = solucion(xsol,ysol);
mesh(xsol,ysol,zsol);
 
led = sprintf('Solucion.');
title(led,'fontsize',18);