function mit18086_multigrid
%MIT18086_MULTIGRID
%    Sets up a 1d Poisson test problem and solves it by multigrid.
%    Method uses twogrid recursively using Gauss-Seidel for smoothing
%    and elimination to solve at coarsest level (n<4).
%    Number of pre- and postsmoothing and coarse grid iteration steps
%    can be prescribed.

% 04/2007 by Benjamin Seibold
%            http://www-math.mit.edu/~seibold/
% Feel free to modify for teaching and learning.
%---------------------------------------------------------------------
levels = 5;                                          % size of problem
nu1 = 2;                           % number of presmoothing iterations
nu2 = 2;                          % number of postsmoothing iterations
gamma = 2;   % number of coarse grid iterations (1=V-cycle, 2=W-cycle)
%---------------------------------------------------------------------
n = 2^(levels+2)-1;                            % number of grid points
h = 1/(n+1);
x = (h:h:(1-h))';
f = pi^2*(sin(pi*x)+4^2*sin(pi*4*x)+9^2*sin(pi*9*x));
A = spdiags(ones(n,1)*[-1 2 -1],-1:1,n,n);
b = f*h^2;
uc = A\b;

global t level, t = 0; level = levels;
clf, subplot(2,2,3:4), hold on
u = twogrid(A,b,nu1,nu2,gamma);
hold off, axis tight
subplot(2,2,1), plot(x,u,'b.-',x,uc,'r.--')
title('correct solution and multigrid approximation')
subplot(2,2,2), plot(x,uc-u,'r.-')
title('error')

%=====================================================================

function x = twogrid(A,b,nu1,nu2,gamma,x0)
%TWOGRID
%    Recursive twogrid cycle for 1d Poisson problem.
%    nu1 = number of presmoothing iterations (Gauss-Seidel)
%    nu2 = number of postsmoothing iterations (Gauss-Seidel)
%    gamma = number of coarse grid iterations (1=V-cycle, 2=W-cycle)
%    x0 = starting vector (0 if not prescribed)
global t level
n = length(b);
if n<4
   x = A\b;                          % solve exactly at coarsest level
else
   G = speye(n)-tril(A)\A; cG = tril(A)\b;       % create Gauss-Seidel
   I = spdiags(ones(n-2,1)*[1 2 1],-2:0,n,n-2); % create interpolation
   I = I(:,1:2:end)/2; R = I'/2;            % and restriction matrices
   if nargin<6, x = b*0; else, x = x0; end           % starting vector
   for i = 1:nu1, x = G*x+cG; end                       % presmoothing
   r = b-A*x;                                       % compute residual
   rh = R*r;                        % restrict residual to coarse grid
   t = t+1; level = level-1; plot([t-1 t],[level+1 level],'bo-')
   eh = rh*0;                                        % starting vector
   for i = 1:gamma
      eh = twogrid(R*A*I,rh,nu1,nu2,gamma,eh); % coarse grid iteration
   end
   e = I*eh;                                       % interpolate error
   t = t+1; level = level+1; plot([t-1 t],[level-1 level],'bo-')
   x = x+e;                                          % update solution
   for i = 1:nu2, x = G*x+cG; end                      % postsmoothing
end