% % This script contains a code that illustrates convergence % behavior for Jacobi, Gauss-Seidel, and SOR iterations. % % Example 0 illustrates behavior of all three on a well % behaved problem - asympotic convergence rate % is predictive % % Example 1 illustrates behavior of SOR on a highly % non-normal problem - asympotic convergence rate % is not predictive at all % % Example 2 illustrates behavior of Gauss Seidel on a flow % problem - with vs. against the direction of the % flow. Asympotic convergence rate is predictive % with the flow but not against it % %---------------------------------------------------------- % % D.C. Sorensen % 11 Sep 02 % % Slight modification of code from M. Embree (10 Sep 02) % %%%%%%%%%%%%%%%%%%%%%% Example 0 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % PDE type matrix for 1-d diffusion reaction equation % % alpha = 4; n = 100; A = alpha*eye(n) - diag(ones(n-1,1),-1) - diag(ones(n-1,1),1); D = diag(diag(A)); E = -tril(A,-1); F = -triu(A,1); % % Take b = 0, so that Ax=b has solution x=0. % Start with x_0 = 1e-2*ones(N,1); % % Apply Jacobi iteration % G = D\(E + F); x0 = 1e-2*ones(n,1); xk = x0; maxit = 400; resvec = [norm(A*x0)]; for j=1:maxit xk = G*xk; resvec = [resvec;norm(A*xk)]; if resvec(end) < 1e-15, break, end end figure(1), clf semilogy([0:length(resvec)-1],resvec, 'g-', 'linewidth',2) rho = max(abs(eig(G))); hold SpecRadJac = max(abs(eig(G))) % % Apply Gauss-Seidel iteration % G = (D - E)\F; x0 = 1e-2*ones(n,1); xk = x0; maxit = 400; resvec = [norm(A*x0)]; for j=1:maxit xk = G*xk; resvec = [resvec;norm(A*xk)]; if resvec(end) < 1e-15, break, end end semilogy([0:length(resvec)-1],resvec, 'k-', 'linewidth',2) SpecRadGS = max(abs(eig(G))) % % Apply SOR iteration with w = 1.5 % w = 1.5; G = (D - w*E)\((1-w)*D + w*F); x0 = 1e-2*ones(n,1); xk = x0; maxit = 400; resvec = [norm(A*x0)]; for j=1:maxit xk = G*xk; resvec = [resvec;norm(A*xk)]; if resvec(end) < 1e-15, break, end end semilogy([0:length(resvec)-1],resvec, 'b-', 'linewidth',2) SpecRadSOR = max(abs(eig(G))) % % Apply SOR iteration with w = optimal % w = 2/(1 + sqrt(1 - rho*rho)); G = (D - w*E)\((1-w)*D + w*F); x0 = 1e-2*ones(n,1); xk = x0; maxit = 400; resvec = [norm(A*x0)]; for j=1:maxit xk = G*xk; resvec = [resvec;norm(A*xk)]; if resvec(end) < 1e-15, break, end end semilogy([0:length(resvec)-1],resvec, 'r-', 'linewidth',2) xlabel('iteration, k', 'fontsize', 18) ylabel('residual norm, ||r_k||', 'fontsize', 18) legend('Jacobi','G-S','SOR','SOR-opt') hold SpecRadOpt = max(abs(eig(G))) disp('any key to continue') pause %%%%%%%%%%%%%%%%%%%%%% Example 1 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % (from Quarteroni, Sacco, Saleri; matrix credited to an % NPL tech report by Hammarling & Wilkinson, 1976) % N.B. spectral radius of iteration matrix is 0.5. % % Apply SOR iteration with w = 1.5 % n = 100; A = 1.5*eye(n) + diag(ones(n-1,1),-1); w = 1.5; D = diag(diag(A)); E = -tril(A,-1); F = -triu(A,1); G = (D - w*E)\((1-w)*D + w*F); r = max(abs(eig(G))); t = 1; % % Take b = 0, so that Ax=b has solution x=0. % Start with x_0 = 1e-2*ones(N,1); x0 = 1e-2*ones(n,1); xk = x0; maxit = 400; resvec = [norm(A*x0)]; atote = [t]; for j=1:maxit xk = G*xk; resvec = [resvec;norm(A*xk)]; t = t*r; atote = [atote; t]; if resvec(end) < 1e-15, break, end end figure(2), clf semilogy([0:length(resvec)-1],resvec, 'b-', 'linewidth',2) xlabel('iteration, k', 'fontsize', 18) ylabel('residual norm, ||r_k||', 'fontsize', 18) adjust = resvec(end)/atote(end); atote = atote*adjust; hold; semilogy([310:length(resvec)-1],atote(311:end), 'g:', 'linewidth',2) legend('actual','predicted') hold; specradG = r Con_ratio = resvec(382)/resvec(381) %break disp('any key to continue') pause %%%%%%%%%%%%%%%%%%%%%% Example 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % % (from Elman and Chernesky, SINUM 30 (1993), 1268-1290) % 1-d convection diffusion problem % -nu * u''(x) + gamma*u'(x) = f(x) % for x \in [0,1] with u(0) = alpha, u(1) = beta, nu>0, gamma>0. % % We will solve this system using Gauss-Seidel. We will choose gamma >0 % and in this case the "flow" of G-S is in the direction of the % fluid flow (with the wind). % % We will then choose gamma < 0 and agian solve the system using Gauss-Seidel. % In this case the "flow" of G-S is opposite to the direction of the % fluid flow (against the wind). % % Another way to look at this is to keep gamma > 0 but order the nodes % (grid points) in standard order (with the wind) and then in reverse % order (against the wind). % The spectral radius is the same if we label the nodes "with the wind" % (i.e., from left to right) or "against the wind", but there is a % transient plateau if we order against the wind. (For the usual left-to-right % ordering of the unkowns, "with the wind" corresponds to normal Gauss-Seidel, % while "against the wind" is Gauss-Seidel with the roles of upper and lower % triangular parts swapped. n = 64; nu = 1; gamma = 3*n/2; w = 1; % % With Wind: Take gamma positive ... flow is to right % a = -nu-gamma/(2*n); b = 2*nu; c = -nu+gamma/(2*n); A = a*diag(ones(n-1,1),-1) + ... b*eye(n) + ... c*diag(ones(n-1,1),1); Gdown = tril(A)\(-triu(A,1)); % downwind (with wind) G-S iteration matrix % % Against Wind: Take gamma negative ... flow is to left % a = -nu+gamma/(2*n); b = 2*nu; c = -nu-gamma/(2*n); A = a*diag(ones(n-1,1),-1) + ... b*eye(n) + ... c*diag(ones(n-1,1),1); Gup = tril(A)\(-triu(A,1)); % upwind (against wind) G-S iteration matrix % Solve Ax=0 b = zeros(n,1); x0 = ones(n,1); maxit = 150; resdown = [norm(A*x0)]; r = max(abs(eig(Gdown))); t = 1; atote = [t]; xk = x0; for j=1:maxit xk = Gdown*xk; resdown = [resdown;norm(A*xk)]; t = t*r; atote = [atote; t]; if resdown(end) < 1e-15, break, end end adjust = resdown(end)/atote(end); atote = atote*adjust; resup = [norm(A*x0)]; xk = x0; for j=1:maxit xk = Gup*xk; hold; resup = [resup;norm(A*xk)]; if resup(end) < 1e-15, break, end end hold; figure(3), clf semilogy([0:length(resdown)-1],resdown, 'b-', 'linewidth',2); hold on semilogy([40:length(resdown)-1],atote(41:end), 'g:', 'linewidth',2) semilogy([0:length(resup)-1],resup, 'r-', 'linewidth',2 ); hold; axis([0 140 1e-15 10]) legend('downwind','predicted','upwind') xlabel('iteration, k', 'fontsize', 18) hold; ylabel('residual norm, ||r_k||', 'fontsize', 18) specrad_Gup = max(abs(eig(Gup))) specrad_Gdown = max(abs(eig(Gdown)))