function gs5ptdemo %GS5PTDEMO graphics demonstration of convergence of Gauss-Seidel % % demonstration of Gauss-Seidel for 5-point Laplacian, illustrating % the smoothing property % Douglas N. Arnold, 2009-10-10 niters = 250; % supply f %f = @(x,y) -2*pi^2*sin(pi*x).*sin(pi*y); f = @(x,y) 190*x.^2.*(1-y)-26.4*exp((x-.25).^2+(y-.25).^2); % set mesh size n = 64; 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 %%%%%% Direct solve %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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); fh = f(Xint,Yint); % fh is a mesh function on interior grid points % multiply right hand side by h^2 and adjust for boundary condition gh = h^2*fh; gh(:,1) = gh(:,1) - uh(Iint,1); gh(:,m) = gh(:,m) - uh(Iint,n+1); gh(1,:) = gh(1,:) - uh(1,Jint); gh(m,:) = gh(m,:) - uh(n+1,Jint); % convert the right hand side from a grid function into a vector F = reshape(gh,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') set(gca,'zlim',[-1,1]) l=light('position',[-5,5,10]); lighting gouraud %%%%%% End direct solve %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %%%%%% Gauss-Seidel iteration %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % set up initial iterate (random grid function vanishing on boundary) u = zeros(n+1,n+1); rand('state',0); u(Iint,Jint) = 2*rand(m,m)-1; figure(2) clf uplot=surf(X,Y,u); set(uplot,'facecolor','yellow') set(gca,'zlim',[-1,1]) hold on l=light('position',[-5,5,10]); lighting gouraud % print max norm error fprintf('iteration 0. max norm error: %9.5f\n',max(max(abs(uh-u)))) % iterate, updating plot each time for niter = 1:niters pause for i=2:n for j=2:n u(i,j) = (u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1)-h^2*fh(i-1,j-1))/4; end end clf uplot=surf(X,Y,u); set(uplot,'facecolor','yellow') set(gca,'zlim',[-1,1]) l=light('position',[-5,5,10]); lighting gouraud % print max norm error fprintf('iteration %4d. max norm error: %9.5f\n',niter,max(max(abs(uh-u)))) end