% 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)
%
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