% twogrid lecture, May 27th, 2015 clear, close N = 100; xc = linspace(0,1,N+2); xf = linspace(0,1,2*N+3); hc = xc(2)-xc(1); hf = xf(2)-xf(1); % coarse and fine grid sine functions yc = sin(2*pi*xc); yf = sin(2*pi*xf); %plot(xc,yc,'ro',xf,yf,'b-*') %break % linear interpolation matrix I = zeros(2*N+1,N); v = [1;2;1]; for j=1:N I((2*j-1):(2*j+1),j) = v; end I = 0.5*I; % linear interpolation y = yc(2:end-1); y_int = [0;I*y';0]; hold on %plot(xf,yf,'b-*',xf,y_int,'g-s') %break % full weighting R = 0.5*I'; yf = yf(2:end-1); y_res = [0;R*yf';0]; hold off %plot(xc,yc,'b-*',xc,y_res,'g-s') %break % Poisson problem on the interior of fine grid with 3-point stencil n = 2*N+1; e = ones(n,1); Af = -spdiags([e,-2*e,e],-1:1,n,n)/hf^2; %spy(Af) %break bf = yf'; ur = Af\bf; plot(xf,[0;ur;0],'b-') %break % coarse grid matrix Ac = R*Af*I; % preconditioners for cg M1 = speye(n,n); M2 = M1; %M1 = diag(diag(Af)); M2 = speye(n,n); %M1 = ichol(Af); M2 = M1'; % cg iterations k = 5; u0 = rand(size(bf)); u = pcg(Af,bf,[],k,M1,M2,u0); err1 = sqrt((u-ur)'*Af*(u-ur)); hold on plot(xf,[0;u;0],'r*') title(err1) %break % coarse grid correction rc = R*(bf - Af*u); %ec = Ac\rc; ec = pcg(Ac,rc,[],10,[],[],zeros(size(rc))); ef = I*ec; % cg iterations u = pcg(Af,bf,[],k,M1,M2,u+ef); err2 = sqrt((u-ur)'*Af*(u-ur)); plot(xf,[0;u;0],'go') title(err2)