%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>'