function envelope(A,t,m) % The function 'envelope(A,t,m)' plots the envelope, % using 't' rotations and an 'm'-by-'m' partition, % and the eigenvalues of a square matrix 'A'. % % Authors: Panayiotis Psarrakos and Michael Tsatsomeros. figure, tic, hold on % The command tic-toc computes the execution time. M=moviein(t); for j=1:t, cshell(A,m,(2*pi/t)*(j-1)), M(:,j)=getframe; end sp(A), hold off, toc, xlabel('Real Axis'), ylabel('Imaginary Axis') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function cshell(A,N,theta) % The function 'cshell(A,N,theta)' plots the shell of a % square matrix 'A', with respect to the angle 'theta', % using an 'N'-by-'N' partition. if nargin<2, k=500; theta=0; end [n,mx]=size(A); B=exp(i*theta)*A; R=1.2*norm(A,2); xmin=-R; xmax=R; ymin=-R; ymax=R; tolx=(xmax-xmin)/N; toly=(ymax-ymin)/N; hp=(B+B')/2; shp=(B-B')/2; % The hermitian part of 'B'. [S,T]=eig(hp); % Solve the eigenproblem of 'hp'. % The eigenvalues of the hermitian part 'hp' are sorted % in ascending order and the 2 largest eigenvalues 'd1' % and 'd2' of 'hp' are found. [sa,I]=sort(real(diag(T))); d1=sa(n); d2=sa(n-1); % A normalized eigenvector of 'hp' corresponding to 'd1'. y=S(:,I(n)); y=y/norm(y,2); % Find 'v', that is, the quantity v(B). aux=shp*y; v=norm(aux,2)^2; % Find 'u', that is, the quantity u(B), but avoid % unnecessary computations when u = 0. if isreal(B), u=0; else, u=imag(y'*aux); end % Compute the values of the two-variable polynomial % at the partition points. for I=1:N, xx=(xmin+(I-1)*tolx); for J=1:N yy=(ymin+(J-1)*toly); z=exp(i*theta)*(xx+i*yy); xxx=real(z); yyy=imag(z); X(I)=xx; Y(J)=yy; Z(I,J)=(d2-xxx)*((d1-xxx)^2+(u-yyy)^2)+(d1-xxx)*(v-u^2); end, end % Plot the shell (cubic curve). contour(X,Y,exp(-i*theta)*Z',[0 0],'b') %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function sp(A) % The function 'sp(A)' plots the eigenvalues of a square % matrix 'A' as red +'s. [n,m]=size(A); hold on, h=eig(A); for j=1:n, plot(real(h(j)),imag(h(j)),'r+'), end