% elliptic linear systems, lecture May 20th, 2015 clear, close % five-point star for 33x33 uniform square mesh n = 33; G = numgrid('S',n); %spy(G) %break A = delsq(G); %spy(A) %break b = ones(size(A,1),1); % direct solver for reference solution u = A\b; N = 50; x0 = rand(size(b)); % Jacobi x = x0; D = diag(diag(A)); L = tril(A,-1); for k = 1:N x = x + D\(b-A*x); e_j(k) = sqrt((x-u)'*A*(x-u)); end % plotting the solution %Ux = G; %Ux(G>0) = full(x(G(G>0))); %mesh(Ux) %plotting the reference solution %U = G; %U(G>0) = full(u(G(G>0))); %figure(2) %mesh(U) %break % damped Jacobi x = x0; omega = 2/3; for k = 1:N x = x + omega*(D\(b-A*x)); e_dj(k) = sqrt((x-u)'*A*(x-u)); end % Gauss-Seidel x = x0; for k=1:N x = x + (D+L)\(b-A*x); e_gs(k) = sqrt((x-u)'*A*(x-u)); end % symmetric Gauss-Seidel x = x0; for k=1:N x = x + (D+L)\(b-A*x); x = x + (D+L')\(b-A*x); e_sgs(k) = sqrt((x-u)'*A*(x-u)); end % CG c_cg = cond(A); for k = 1:N x = pcg(A,b,[],k,[],[],x0); e_cg(k) = sqrt((x-u)'*A*(x-u)); end % plotting energy errors plot(1:N,e_j,'r-*',1:N,e_dj,'b-d',1:N,e_gs,'m-*',... 1:N,e_gs,'c-d', 1:N,e_cg,'g-s') legend('Jacobi','damped Jacobi','GS','symm. GS','CG') %break % PCG: Jacobi M1 = D; M2 = speye(size(A)); M = M1*M2; c_cg = cond(A), c_j = cond(M\A) for k = 1:N x = pcg(A,b,[],k,M1,M2,x0); e_cgj(k) = sqrt((x-u)'*A*(x-u)); end % PCG: SGS M1 = D+L; M2 = D\(D+L'); M = M1*M2; c_g = cond(A), c_gs = cond(M\A), break for k=1:N x = pcg(A,b,[],k,M1,M2,x0); e_cgg(k) = sqrt((x-u)'*A*(x-u)); end % PCG: ichol M1 = ichol(A); M2 = M1'; for k = 1:N x = pcg(A,b,[],k,M1,M2,x0); e_cgc(k) = sqrt((x-u)'*A*(x-u)); end semilogy(1:N,e_cg,1:N,e_cgj,1:N,e_cgg,1:N,e_cgc) legend('CG','PCG:Jacobi','PCG:SGS','PCG:ichol')