%Programa ejemplo  donde se implementa el método de  colocación  no simétrica 
%para aproximar la solución de la ecuación de Poisson con condiciones de frontera Dirichlet 
%utilizando el kernel capa delagada.
% 
%Entrada : NI  número de nodos interiores
%                 NF número de nodos frontera
%Salida:    Gráfica de la aproximación y solución 
%
%Autor: Daniel A. Cervantes Cabrera
%Supervisión : Pedro González Casanova.

function ucm_thinplate(NI,NF) 
func  = 'phitps4xy'; 
[I F] = gnifa(0,1,0,1,NI,NF); 
[ni t] = size(I); 
[nf t] = size(F); 
  
N = ni + nf; 
  
xn0 = I(1:ni)'; 
yn0 = I(ni+1:2*ni)'; 
 
sol = laplaciano(xn0,yn0); 
 
nft = (nf + 4)/4; 
 
xn1 = F(1:nft-1)'; 
xn2 = F(nft:2*nft-2)'; 
xn3 = F(2*nft-1:3*nft-3)'; 
xn4 = F(3*nft-2:4*nft-4)'; 
 
yn1 = F(4*nft-3:5*nft-5)'; 
yn2 = F(5*nft-4:6*nft-6)'; 
yn3 = F(6*nft-5:7*nft-7)'; 
yn4 = F(7*nft-6:8*nft-8)'; 
 
 
 
sol(ni+1:ni+nft-1)          = g3(xn1,yn1); 
sol(ni+nft:ni+2*nft-2)      = g2(xn2,yn2); 
sol(ni+2*nft-1:ni+3*nft-3)  = g4(xn3,yn3); 
sol(ni+3*nft-2:ni+4*nft-4)  = g1(xn4,yn4); 
sol(ni+4*nft-3:ni+4*nft+2)  = 0.0;  
 
 
x =[xn0 ; xn1; xn2; xn3; xn4]; 
y =[yn0 ; yn1; yn2; yn3; yn4]; 
  
 
 [n1 t] = size(xn1); 
 [n2 t] = size(xn2); 
 [n3 t] = size(xn3); 
 [n4 t] = size(xn4); 
 
 
 for i=1:ni 
    for j=1:N 
      A(i,j) = feval(func,xn0(i),x(j),yn0(i),y(j),1) + feval(func,xn0(i),x(j),yn0(i),y(j),2);  
    end 
 end 
  
 for i=1:ni 
       A(i,N+1) = 2.0; 
       A(i,N+2) = 2.0; 
       A(i,N+3:N+6) = 0.0; 
 end 

for i=1:n1 
  for j=1:N  
    A1(i,j) = feval(func,xn1(i),x(j),yn1(i),y(j),0);  
  end 
end 
 
 for i=1:n1 
     for j=1:6 
       switch j 
           case 1, 
             aij = xn1(i)^2; 
           case 2, 
             aij = yn1(i)^2; 
           case 3, 
             aij = xn1(i)*yn1(i); 
           case 4, 
               aij = xn1(i); 
           case 5, 
               aij = yn1(i); 
           case 6, 
               aij = 1; 
       end 
       P(i,j) = aij; 
     end 
 end 
 
A1 = [A1 P]; 

for i=1:n2 
  for j=1:N 
    A2(i,j) = feval(func,xn2(i),x(j),yn2(i),y(j),0);  
  end 
end 
 
%size(A2)  
 
 for i=1:n2 
     for j=1:6 
       switch j 
           case 1, 
             aij = xn2(i)^2; 
           case 2, 
             aij = yn2(i)^2; 
           case 3, 
             aij = xn2(i)*yn2(i); 
           case 4, 
               aij = xn2(i); 
           case 5, 
               aij = yn2(i); 
           case 6, 
               aij = 1; 
       end 
       P2(i,j) = aij; 
     end 
 end 
 
A2 = [A2 P2]; 

for i=1:n3 
  for j=1:N 
    A3(i,j) = feval(func,xn3(i),x(j),yn3(i),y(j),0);  
  end 
end 
 
%size(A1)  
 
 for i=1:n3 
     for j=1:6 
       switch j 
           case 1, 
             aij = xn3(i)^2; 
           case 2, 
             aij = yn3(i)^2; 
           case 3, 
             aij = xn3(i)*yn3(i); 
           case 4, 
               aij = xn3(i); 
           case 5, 
               aij = yn3(i); 
           case 6, 
               aij = 1; 
       end 
       P3(i,j) = aij; 
     end 
 end 
 
A3 = [A3 P3]; 

for i=1:n4 
  for j=1:N 
    A4(i,j) = feval(func,xn4(i),x(j),yn4(i),y(j),0);  
  end 
end 
 
 for i=1:n4 
     for j=1:6 
       switch j 
           case 1, 
             aij = xn4(i)^2; 
           case 2, 
             aij = yn4(i)^2; 
           case 3, 
             aij = xn4(i)*yn4(i); 
           case 4, 
               aij = xn4(i); 
           case 5, 
               aij = yn4(i); 
           case 6, 
               aij = 1; 
       end 
       P4(i,j) = aij; 
     end 
 end 
 
A4 = [A4 P4]; 

%Momentos 
 
    for i=1:ni 
        A5(1,i) = xn0(i) * xn0(i); 
        A5(2,i) = yn0(i) * yn0(i); 
        A5(3,i) = xn0(i) * yn0(i); 
        A5(4,i) = xn0(i); 
        A5(5,i) = yn0(i); 
        A5(6,i) = 1; 
    end  
 
    j = ni + 1; 
    for i=1:n1 
        A5(1,j) = xn1(i) * xn1(i); 
        A5(2,j) = yn1(i) * yn1(i); 
        A5(3,j) = xn1(i) * yn1(i); 
        A5(4,j) = xn1(i); 
        A5(5,j) = yn1(i); 
        A5(6,j) = 1; 
        j = j + 1; 
    end  
 
    %j = j + n1 + 1; 
    for i=1:n2 
        A5(1,j) = xn2(i) * xn2(i); 
        A5(2,j) = yn2(i) * yn2(i); 
        A5(3,j) = xn2(i) * yn2(i); 
        A5(4,j) = xn2(i); 
        A5(5,j) = yn2(i); 
        A5(6,j) = 1; 
        j = j + 1; 
    end  
     
    %j = j + n2 + 1; 
    for i=1:n3 
        A5(1,j) = xn3(i) * xn3(i); 
        A5(2,j) = yn3(i) * yn3(i); 
        A5(3,j) = xn3(i) * yn3(i); 
        A5(4,j) = xn3(i); 
        A5(5,j) = yn3(i); 
        A5(6,j) = 1; 
        j = j + 1; 
    end  
     
    %j = j + n3 + 1; 
     for i=1:n4 
        A5(1,j) = xn4(i) * xn4(i); 
        A5(2,j) = yn4(i) * yn4(i); 
        A5(3,j) = xn4(i) * yn4(i); 
        A5(4,j) = xn4(i); 
        A5(5,j) = yn4(i); 
        A5(6,j) = 1; 
        j = j + 1; 
     end  
 
     for i =1:6 
         A5(i,N+1:N+6) = 0; 
     end 
     
   
     G = [A ; A1; A2; A3; A4; A5];    
    
    [nt mt] = size(G); 
    fid =  fopen('matriz.txt','w');   
    
   
 
    for i=1:nt 
      for j=1:mt 
        fprintf(fid,'%2.5f ',G(i,j));  
      end 
      fprintf(fid,'\n');  
   end 
 
    fclose(fid); 
 
%Determino la solucion 
lambda = G\(sol); 
 
   
%Gr'afica del la aproximaci'on 
xm = linspace(0,1,101); 
ym = linspace(0,1,101); 
[t,nxym] = size(xm); 
 
hold on 
for i=1:nxym 
    for j=1:nxym 
      zm(i,j) = aprox_thinplate(xm(j),ym(i),x,y,lambda,func); 
    end 
end 
 
for i=1:n1 
    plot(xn1(i),yn1(i),'g.'); 
end 
 
for i=1:n2 
    plot(xn2(i),yn2(i),'g.'); 
end 
 
for i=1:n3 
    plot(xn3(i),yn3(i),'g.'); 
end 
 
for i=1:n4 
    plot(xn4(i),yn4(i),'g.'); 
end 
 
 
for i=1:ni 
   zn0(i) = aprox_thinplate(xn0(i),yn0(i),x,y,lambda,func); 
   plot3(xn0(i),yn0(i),zn0(i),'ro');  
   plot(xn0(i),yn0(i),'r.');  
end 
    
[xm2d,ym2d] = meshgrid(0:.01:1); 
mesh(xm2d,ym2d,zm); 
led = sprintf('Aproximacion.'); 
title(led,'fontsize',18); 
xlabel(sprintf('Condicion %f',cond(G))); 
 
   
%Gr'afica de la soluci'on  
figure 
[xsol,ysol] = meshgrid(0:.01:1); 
zsol = solucion(xsol,ysol); 
mesh(xsol,ysol,zsol); 
led = sprintf('Exacta.'); 
title(led,'fontsize',18);

</pre>'