% lap5pt.m -- solve Dirichlet's problem for the Poisson equation %LAP5PT 5 point Laplacian difference method for Poisson equation % % Laplacian_h u = f on interior grid points % u = g on boundary grid points % % using 5-point Laplacian. Domain is unit square. % Douglas N. Arnold, adapted from poisson.m of Randy Leveque % supply f and, if known, exact solution u as functions f = @(x,y) -2*pi^2*sin(pi*x).*sin(pi*y); u = @(x,y) sin(pi*x).*sin(pi*y); %f = @(x,y) .5*((x-3/4).^2+(y-1/3).^2<1/256)-((x-1/3).^2+(y-3/4).^2 <1/256); % set mesh size n = 8; h = 1/n; m = n-1; x = linspace(0,1,n+1); % x grid points including boundaries y = linspace(0,1,n+1); % y grid points including boundaries [X,Y] = meshgrid(x,y); % 2d arrays of x,y values X = X'; % transpose so that X(i,j),Y(i,j) are Y = Y'; % coordinates of (i,j) grid point, i.e., ((i-1)h,(j-1)h) Iint = 2:n; % indices of interior points in x Jint = 2:n; % indices of interior points in y Xint = X(Iint,Jint); % interior points Yint = Y(Iint,Jint); fh = f(Xint,Yint); % fh is a mesh function on interior grid points % set uh on the boundary to given values (interior values will be ignored) uh = zeros(n+1,n+1); % can use true solution if it is known % form 5-point Laplacian matrix (times h^2) I = speye(m); e = ones(m,1); A = spdiags([e,-4*e,e],[-1,0,1],m,m); J = spdiags([e,e],[-1,1],m,m); Lh = kron(I,A) + kron(J,I); % multiply right hand side by h^2 and adjust for boundary condition fh = h^2*fh; fh(:,1) = fh(:,1) - uh(Iint,1); fh(:,m) = fh(:,m) - uh(Iint,n+1); fh(1,:) = fh(1,:) - uh(1,Jint); fh(m,:) = fh(m,:) - uh(n+1,Jint); % convert the right hand side from a grid function into a vector F = reshape(fh,m*m,1); % Solve the linear system for the grid function uh on the interior, written % as vector U = Lh\F; % convert from vector into grid function use this for uh at the interior % grid points (uh is already given at the boundary grid points) uh(Iint,Jint) = reshape(U,m,m); % plot results: figure(1) clf uhplot=surf(X,Y,uh); set(uhplot,'facecolor','yellow') l=light('position',[-5,5,10]); lighting gouraud % skip the rest if the true solution is not known return utrue = u(X,Y); % utrue is the solution at all the grid points % print error err = max(max(abs(uh-utrue))); fprintf('Absolute error in max norm = %10.3e\n',err) fprintf('Relative error in max norm = %10.3e\n',err/max(max(utrue))) % plot true solution on fine mesh for comparison figure(2) clf xf = linspace(0,1,101); % x grid points including boundaries yf = linspace(0,1,101); % y grid points including boundaries [Xf,Yf] = meshgrid(xf,yf); % 2d arrays of x,y values Xf = Xf'; % transpose so that X(i,j),Y(i,j) are Yf = Yf'; % coordinates of (i,j) point uplot=surf(Xf,Yf,u(Xf,Yf)); set(uplot,'edgecolor','none','facecolor','yellow') l=light('position',[-5,5,10]); lighting gouraud