%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% multigrid demo on simple -D^2 u = f  u(0)=u(1)=0
%
% finest grid as 2^m+1 points (including boundaries)
%
% requires: buildA.m, coarsen.m, interp.m, jacobi.m 
%
% Same as multigrid.m but with plots.
%
% Fabian Waleffe, 11/15/2000 for Math 801
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

m=5; j=1:2^m+1;

A=buildA(m); f=randn(length(A),1); ue=A\f; %A, f and exact solution

% First relaxation:

u=jacobi(A,f); rh=f-A*u; 

e2hs=[]; r2hs=[];

hold off
plot(j,ue,'r',j,u,'--'), title('EXACT SOLUTION and JACOBI APPROXIMATION'), pause
plot(rh), title('RESIDUAL'), pause


for n=m-1:-1:2,

hold on
  r2h=coarsen(rh); 
  A=buildA(n); e2h=jacobi(A,r2h); rh=r2h-A*e2h; 
  r2hs=[r2h;r2hs]; e2hs=[e2h;e2hs];       % need to save results

  jc=[1:2^(m-n):2^m+1];
  plot(jc,r2h,'--'), title('COARSENED RESIDUAL'), pause
hold off
  plot(jc,rh),  title('RESIDUAL AFTER ITERATIONS'), pause

%  plot(jc,e2h), title('APPROXIMATE ERROR'), pause


end

 r2h=coarsen(rh); A=buildA(1); e2h=A\r2h; % EXACT SOLUTION ON COARSEST GRID

  jc=[1:2^(m-1):2^m+1];
  hold on
  plot(jc,r2h,'--'), title('COARSEST RESIDUAL'), pause
  hold off
  plot(jc,e2h), title('COARSEST  ERROR'), pause

i2=0;
for n=2:m-1,

  eh=interp(e2h);  

  i1=i2+1;i2=i2+(2^n+1); eh=e2hs(i1:i2)+eh;   % correction step

  A=buildA(n); rh=r2hs(i1:i2); e2h=jacobi(A,rh,eh);

  jc=[1:2^(m-n):2^m+1];
  plot(jc,e2hs(i1:i2),'--'), title('APPROXIMATE ERROR'), pause
  hold on 
  plot(jc,eh,'m'),   title('CORRECTED ERROR'), pause
  plot(jc,e2h),  title('ERROR AFTER ITERATIONS'), pause
  hold off

end

hold off
plot(j,ue,'r',j,u,'--'), title('EXACT SOLUTION and JACOBI APPROXIMATION'), pause

eh=interp(e2h); u=u+eh; 
hold on
plot(j,u,'k--'), title('CORRECTED APPROXIMATION'), pause
A=buildA(m); u=jacobi(A,f,u); %u=jacobi(A,f,u,3,10); 

hold on
plot(j,u), title('FINAL RELAXATION')
hold off
