% Matlab code fmgv.m % For "Applied Numerical Linear Algebra", Question 6.16 % Written by James Demmel, Jul 10, 1993 % Modified, Jun 2, 1997 % % Full Multigrid V-Cycle for Poisson's equation on a square grid % with zero Dirichlet boundary conditions. % Algorithm from Brigg's ``Multigrid Tutorial'' % Include zero boundary values in arrays for ease of programming. % Assume dimension n = 2^k + 1 % Inputs: % x = initial guess (n by n matrix with zeros on boundary) % b = right hand side (n by n matrix) % jac1, jac2 = number of weighted Jacobi steps to do before and % after recursive call to mgv % Outputs: % z = improved solution (n by n matrix with zeros on boundary) % % function z=fmgv(x,b,jac1,jac2); % function z=fmgv(x,b,jac1,jac2); [nn,m]=size(b); % Compute residual r=zeros(nn,nn); tmp = b(2:nn-1,2:nn-1) - ( 4*x(2:nn-1,2:nn-1) ... - x(1:nn-2,2:nn-1) - x(3:nn,2:nn-1) ... - x(2:nn-1,1:nn-2) - x(2:nn-1,3:nn) ); r(2:nn-1,2:nn-1) = tmp; % Get right hand side for coarsest grid by repeated restriction of % original right hand side to coarser grids rc=mgvrhs(nn,3,r); % Solve 1 by 1 problem z=zeros(3,3); z(2,2)=rc(2,2)/4; % Assume nn = 2^k+1 k = round(log(nn-1)/log(2)); % Loop from the next-to-coarsest to finest grids for i=2:k, n=2^i+1; m=2^(i-1)+1; % Interpolate coarser solution to this level zstrt=zeros(n,n); zstrt(3:2:n-2,3:2:n-2)=z(2:m-1,2:m-1); zstrt(2:2:n-1,3:2:n-2)=.5*(z(1:m-1,2:m-1)+z(2:m,2:m-1)); zstrt(3:2:n-2,2:2:n-1)=.5*(z(2:m-1,1:m-1)+z(2:m-1,2:m)); zstrt(2:2:n-1,2:2:n-1)=.25*(z(1:m-1,1:m-1)+z(1:m-1,2:m)+ ... z(2:m,1:m-1)+z(2:m,2:m) ); % Get right hand side for current grid by repeated restriction of % original right hand side to coarser grids rhs=mgvrhs(nn,n,r); % Do Multigrid V-cycle z=mgv(zstrt,rhs,jac1,jac2); % mesh(z),title(int2str(n)), pause end % Add final correction to initial guess z=x+z; return