% Matlab code mgv.m % For "Applied Numerical Linear Algebra", Question 6.16 % Written by James Demmel, Jul 10, 1993 % Modified, Jun 2, 1997 % % 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=mgv(x,b,jac1,jac2) [n,m]=size(b); if n == 3, % Solve 1 by 1 problem explicitly z=zeros(3,3); z(2,2)=b(2,2)/4; else % Perform jac1 steps of weighted Jacobi with weight w w=2/3; for i=1:jac1; x(2:n-1,2:n-1)=(1-w)*x(2:n-1,2:n-1) ... + (w/4)*( x(1:n-2,2:n-1) + x(3:n,2:n-1) ... + x(2:n-1,1:n-2) + x(2:n-1,3:n) ... + b(2:n-1,2:n-1)); end, % Compute residual on current grid r=zeros(n,n); tmp = b(2:n-1,2:n-1) - ( 4*x(2:n-1,2:n-1) ... - x(1:n-2,2:n-1) - x(3:n,2:n-1) ... - x(2:n-1,1:n-2) - x(2:n-1,3:n) ); r(2:n-1,2:n-1) = tmp; % Restrict residual to coarse grid m=(n+1)/2; rhat=zeros(m,m); rhat(2:m-1,2:m-1)=.25*r(3:2:n-2,3:2:n-2) ... + .125*(r(2:2:n-3,3:2:n-2)+r(4:2:n-1,3:2:n-2) ... + r(3:2:n-2,2:2:n-3)+r(3:2:n-2,4:2:n-1) ) ... + .0625*(r(2:2:n-3,2:2:n-3)+r(4:2:n-1,2:2:n-3) ... + r(2:2:n-3,4:2:n-1)+r(4:2:n-1,4:2:n-1) ); % Solve recursively with zero initial guess % Multiply residual by 4 for coarser grid xhat = mgv(zeros(size(rhat)),4*rhat,jac1,jac2); % Interpolate coarse solution to fine grid, and add to x xcor=zeros(size(x)); xcor(3:2:n-2,3:2:n-2)=xhat(2:m-1,2:m-1); xcor(2:2:n-1,3:2:n-2)=.5*(xhat(1:m-1,2:m-1)+xhat(2:m,2:m-1)); xcor(3:2:n-2,2:2:n-1)=.5*(xhat(2:m-1,1:m-1)+xhat(2:m-1,2:m)); xcor(2:2:n-1,2:2:n-1)=.25*(xhat(1:m-1,1:m-1)+xhat(1:m-1,2:m)+ ... xhat(2:m,1:m-1)+xhat(2:m,2:m) ); z=x+xcor; % Perform jac2 steps of weighted Jacobi with weight w for i=1:jac2; z(2:n-1,2:n-1)=(1-w)*z(2:n-1,2:n-1) ... + (w/4)*( z(1:n-2,2:n-1) + z(3:n,2:n-1) ... + z(2:n-1,1:n-2) + z(2:n-1,3:n) ... + b(2:n-1,2:n-1)); end, end return %end