% Matlab code testfmgv.m % For "Applied Numerical Linear Algebra", Question 6.16 % Written by James Demmel, Jul 10, 1993, modified Jun 2, 1997. % Simplified Jan, 2009 by Ivar Gustafsson % % Test Full Multigrid V-Cycle code for solving Poisson's equation % on a square grid with zero Dirichlet boundary conditions % Inputs: (run makemgdemo to initialize inputs for demo) % b = right hand side, an n=2^k+1 by n matrix with % explicit zeros on boundary % x = initial solution guess % xtrue = true solution % jac1, jac2 = number of weight Jacobi steps to do before and % after recursive call to mgv % itmax = maximum number of full multigrid cycles to perform % tolres = for stopping criterion, norm of residual smaller than tolres % Outputs: % xnew = improved solution % res = vector of residual after each iteration % costpi = time per iteration per unknown % (should be constant independent of number of unknowns) % plot of approximate solution after each iteration % residual after each iteration % plot of residual(i+1)/residual(i), which should be constant % independent of right hand side and dimension % plot of residual versus iteration numbers % res=[]; [n,m]=size(b); f2=0; figure(1) hold off, clf subplot(1,2,1),mesh(xtrue),title('True Solution') subplot(1,2,2),mesh(b),title('Right Hand Side') disp('hit return to continue'),pause xnew = x; hold off, clf % Loop over all iterations of Full multigrid fprintf(' Full multigrid \n jac1=%3.0f, jac2=%3.0f, tolres=%10.2e \n',jac1,jac2, tolres) fprintf(' Hit return to start on a new cycle! \n \n') fprintf(' it Residual Error norm \n') costpi=0; for i=1:itmax, t1=cputime; % Do a full multigrid cycle xnew=fmgv(xnew,b,jac1,jac2); % Accumulate the cpu-time costpi=costpi+cputime-t1; % Compute and save residual tmp = b(2:n-1,2:n-1) - ( 4*xnew(2:n-1,2:n-1) ... - xnew(1:n-2,2:n-1) - xnew(3:n,2:n-1) ... - xnew(2:n-1,1:n-2) - xnew(2:n-1,3:n) ); t=norm(tmp,1); res=[res;t]; subplot(1,2,1), mesh(xnew), title(['approximate solution ',int2str(i)]) % Choose one of these depending on whether true solution is known % subplot(1,2,2), mesh(tmp), title('residual') subplot(1,2,2), mesh(xnew-xtrue), title('error') % plot(xnew(round(n/2),:)), title('approximate solution') % disp('iteration, residual, error='),i,t,norm(xnew-xtrue) fprintf('%4.0f %10.2e %10.2e \n',i,t,norm(xnew-xtrue)); pause % disp('hit return to continue'),pause % Stop if residual small enough if t