function v = mg5pt(f,u) % MG5PT multigrid iteration for 5 point Laplacian % % Usage: v = mg5pt(f,u) % % f is a vector giving the right hand side % u is the initial guess, a grid function defined on the n x n interior grid points % output v is the improved solution after 1 V-cycle iteration % % size of vectors must be n = 2^p-1 with p an integer % % To solve the 5 point Laplacian, we choose an initial guess for u (e.g., zero), and % iterate several times: u=zeros(size(f)); for niter=1:niters; u=mg5pt(f,u); end % Douglas N. Arnold, 2009-10-15 n = length(u); if n == 1 v = f/16; % n=1 when h=1/2, so 4/h^2=16 else h = 1/(n+1); n2 = (n-1)/2; % presmooth with Gauss-Seidel % add zero values at boundary grid points U = [zeros(1,n+2);zeros(n,1),u,zeros(n,1);zeros(1,n+2)]; for i=2:n+1 for j=2:n+1 U(i,j) = (U(i+1,j)+U(i-1,j)+U(i,j+1)+U(i,j-1)+h^2*f(i-1,j-1))/4; end end u = U(2:n+1,2:n+1); % compute residual r = f -(4*u-U(1:n,2:n+1)-U(3:n+2,2:n+1)-U(2:n+1,1:n)-U(2:n+1,3:n+2))/h^2; % transfer to coarser mesh r2 = ( 4*r(2:2:n-1,2:2:n-1) ... +2*(r(1:2:n-2,2:2:n-1)+r(3:2:n,2:2:n-1)+r(2:2:n-1,1:2:n-2)+r(2:2:n-1,3:2:n)) ... +r(1:2:n-2,1:2:n-2)+r(1:2:n-2,3:2:n)+r(3:2:n,1:2:n-2)+r(3:2:n,3:2:n) )/16; % recursively apply multigrid iteration on coarse mesh e2 = mg5pt(r2,zeros(size(r2))); % transfer solution to fine mesh % add boundary values E2 = [zeros(1,n2+2);zeros(n2,1),e2,zeros(n2,1);zeros(1,n2+2)]; e = zeros(n,n); e(2:2:n-1,2:2:n-1) = e2; e(1:2:n,2:2:n-1) = (E2(1:n2+1,2:n2+1)+E2(2:n2+2,2:n2+1))/2; e(2:2:n-1,1:2:n) = (E2(2:n2+1,1:n2+1)+E2(2:n2+1,2:n2+2))/2; e(1:2:n,1:2:n) = (E2(1:n2+1,1:n2+1)+E2(1:n2+1,2:n2+2)+E2(2:n2+2,1:n2+1)+E2(2:n2+2,2:n2+2))/4; % correct v = u + e; % postsmooth with Gauss-Seidel in reverse order % add zero valves at boundary grid points V = [zeros(1,n+2);zeros(n,1),v,zeros(n,1);zeros(1,n+2)]; V(2:n+1,2:n+1) = v; for i=n+1:-1:2 for j=n+1:-1:2 V(i,j) = (V(i+1,j)+V(i-1,j)+V(i,j+1)+V(i,j-1)+h^2*f(i-1,j-1))/4; end end v = V(2:n+1,2:n+1); end