function mit18086_smoothing %MIT18086_SMOOTHING % Sets up a 1d Poisson problem with an oscillatory right hand % side and applies various iterative solvers to it. % Observe the decay rate of various error frequencies: % Jacobi smooths medium frequencies best, but fails for % grid-scale oscillations. Weighted-Jacobi remedies this % problem at the expense of a slower decay rate for small % frequencies. Gauss-Seidel performs well as a smoother. % SOR smooths best, but it is also more expensive in practice. % Observe how the error gets smoothed very quickly, but the % actual convergence to the correct solution is very slow. % 04/2008 by Benjamin Seibold % Feel free to modify for teaching and learning. %------------------------------------------------------------------ b = 11; % number of oscillations n = 50; % number of grid points %------------------------------------------------------------------ % set up 1d Poisson problem with oscillatory right hand side h = 1/(n+1); x = (h:h:(1-h))'; sx = 2*pi*b*x.^2; f = 1e2*sx.*sin(sx); A = spdiags(ones(n,1)*[-1 2 -1],-1:1,n,n); b = f*h^2; uc = A\b; I = speye(n); D = diag(diag(A)); L = tril(A)-D; %------------------------------------------------------------------ % different preconditioners name = {'Jacobi','weighted Jacobi','Gauss-Seidel','SOR'}; P = {D,3/2*D,D+L,D+1.8*L}; %------------------------------------------------------------------ np = length(P); for k = 1:np M{k} = I-P{k}\A; c{k} = P{k}\b; u{k} = zeros(n,1); end for i = 1:100 clf for k = 1:np, u{k} = M{k}*u{k}+c{k}; end subplot(1,2,1) plot(x,uc,'k.-') hold all, for k = 1:np, plot(x,u{k},'.-'), end, hold off title(sprintf('solution after %d steps',i)) subplot(1,2,2) hold all, for k = 1:np, plot(x,uc-u{k},'.-'), end, hold off legend(name) title(sprintf('error after %d steps',i)) drawnow end