Exemplo 2. Equação de Poisson

% IFPB 30/11/2017
% EQUAÇÃO DE POISSON COMK CONDIÇÕES DE CONTORNO DE DIRICHLET
% Uxx - Uyy = 4
% Condições de contorno
% U(x,1)=x^2+1 e U(x,-1)=x^2+1, U(1,y)=y^2+1, U(-1,y)=y^2+1
% Aproximação no quadrado [-1;1]x[-1;1]
% Solução analítica u(x,y)=x^2+y^2

close all, clear all, clc;

% malha computacional bidiminensional
n=100;  h=2/(n-1); x=-1:h:1; y=x;

% Condições de contorno
u=zeros(n);
u(1,:)=y.^2+1; u(n,:)=y.^2+1;
u(:,1)=x.^2+1; u(:,n)=x.^2+1;

% Equação de diferenças
for iter=1:1000,    for i=2:n-1,    for j=2:n-1
            u(i,j)=0.25*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-h^2*4);
        end, end, end

fig1=figure(1); set(fig1,'Position',[1 33 1024 662],'Color',[1 1 1]);
surfc(u), shading interp
label(1)=xlabel('x'); label(2)=ylabel('y');
label(3)=zlabel('u(x,y)');
label(4)=title('Equação de Poisson: Uxx - Uyy = f');
label(5)=legend('Diferenças Finitas',1);
set(label, 'FontSize',16); drawnow,pause(1)

fig2=figure(2); set(fig2,'Position',[1 33 1024 662],'Color',[1 1 1]);
[x1,y1] = meshgrid(x,y); [px,py] = gradient(u,h,h);
[c,h] = contour(x1,y1,u);
cl=clabel(c,h); set(cl,'FontSize', 14);
colorbar; set(h,'Linewidth', 2.5);
hold on, s=1:10:100;
quiver(x(s),y(s),-px(s,s),-py(s,s)), hold off,
label(1)=xlabel('x'); label(2)=ylabel('y');
label(3)=title('Equação de Poisson: Uxx - Uyy = f');
label(4)=legend('Diferenças Finitas',1);
set(label, 'FontSize',16); drawnow