p29.m

% p29.m - solve Poisson equation on the unit disk
% (compare p16.m and p28.m)
% Laplacian in polar coordinates:
N = 31; [D,r] = cheb(N); N2 = (N-1)/2; D2 = D^2;
D1 = D2(2:N2+1,2:N2+1); D2 = D2(2:N2+1,N:-1:N2+2);
E1 = D(2:N2+1,2:N2+1); E2 = D(2:N2+1,N:-1:N2+2);
M = 40; dt = 2*pi/M; t = dt*(1:M)'; M2 = M/2;
D2t = toeplitz([-pi^2/(3*dt^2)-1/6 ...
.5*(-1).^(2:M)./sin(dt*(1:M-1)/2).^2]);
R = diag(1./r(2:N2+1)); Z = zeros(M2); I = eye(M2);
L = kron(D1+R*E1,eye(M))+kron(D2+R*E2,[Z I;I Z])+kron(R^2,D2t);
% Right-hand side and solution for u:
[rr,tt] = meshgrid(r(2:N2+1),t); rr = rr(:); tt = tt(:);
f = -rr.^2.*sin(tt/2).^4 + sin(6*tt).*cos(tt/2).^2; u = L\f;
% Reshape results onto 2D grid and plot them:
u = reshape(u,M,N2); u = [zeros(M+1,1) u([M 1:M],:)];
[rr,tt] = meshgrid(r(1:N2+1),t([M 1:M]));
[xx,yy] = pol2cart(tt,rr);
clf, subplot('position',[.1 .4 .8 .5])
mesh(xx,yy,u), view(20,40), colormap([0 0 0])
axis([-1 1 -1 1 -.01 .05]), xlabel x, ylabel y, zlabel u