% psa.m - Simple code for 2-norm pseudospectra of given matrix A. % Typically about N/4 times faster than the obvious SVD method. % Comes with no guarantees! - L. N. Trefethen, March 1999. % Set up grid for contour plot: npts = 20; % <- ALTER GRID RESOLUTION s = .8*norm(A,1); xmin = -s; xmax = s; ymin = -s; ymax = s; % <- ALTER AXES x = xmin:(xmax-xmin)/(npts-1):xmax; y = ymin:(ymax-ymin)/(npts-1):ymax; [xx,yy] = meshgrid(x,y); zz = xx + sqrt(-1)*yy; % Compute Schur form and plot eigenvalues: [U,T] = schur(A); if isreal(A), [U,T] = rsf2csf(U,T); end, T = triu(T); eigA = diag(T); hold off, plot(real(eigA),imag(eigA),'.','markersize',15), hold on axis([xmin xmax ymin ymax]), axis square, grid on, drawnow % Reorder Schur decomposition and compress to interesting subspace: select = find(real(eigA)>-250); % <- ALTER SUBSPACE SELECTION n = length(select); for i = 1:n for k = select(i)-1:-1:i G([2 1],[2 1]) = planerot([T(k,k+1) T(k,k)-T(k+1,k+1)]')'; J = k:k+1; T(:,J) = T(:,J)*G; T(J,:) = G'*T(J,:); end end T = triu(T(1:n,1:n)); I = eye(n); % Compute resolvent norms by inverse Lanczos iteration and plot contours: sigmin = Inf*ones(length(y),length(x)); for i = 1:length(y) if isreal(A) & (ymax==-ymin) & (i>length(y)/2) sigmin(i,:) = sigmin(length(y)+1-i,:); else for j = 1:length(x) z = zz(i,j); T1 = z*I-T; T2 = T1'; if real(z)<100 % <- ALTER GRID POINT SELECTION sigold = 0; qold = zeros(n,1); beta = 0; H = []; q = randn(n,1) + sqrt(-1)*randn(n,1); q = q/norm(q); for k = 1:99 v = T1\(T2\q) - beta*qold; alpha = real(q'*v); v = v - alpha*q; beta = norm(v); qold = q; q = v/beta; H(k+1,k) = beta; H(k,k+1) = beta; H(k,k) = alpha; if (alpha>1e100), sig = alpha; else sig = max(eig(H(1:k,1:k))); end if (abs(sigold/sig-1)<.001) | (sig<3 & k>2) break, end sigold = sig; end %text(x(j),y(i),num2str(k)) % <- SHOW ITERATION COUNTS sigmin(i,j) = 1/sqrt(sig); end end end, disp(['finished line ' int2str(i) ' out of ' int2str(length(y))]) end contour(x,y,log10(sigmin+1e-20),-8:-1); % <- ALTER LEVEL LINES