% Matlab code RayleighContour.m % For "Applied Numerical Linear Algebra", Figure 5.1 % Written by James Demmel, Jun 2, 1997 % % Plots contours of Rayleigh quotient rho(u,A) on the unit sphere % % Input % A = any 3-by-3 symmetric matrix % % Output % Plot of contours of rho(u,A) % [u,d]=eig(A); ds=-sort(-diag(d)); if ds(1)==ds(3), % A is a multiple of the identity disp(['rho(u,A) = ',num2str(ds(1)),' is constant.']) else % choose contour values CON to plot % normalize spectrum to [1 lam 0] ds=ds-ds(3); ds=ds/ds(1); lam = ds(2); if (lam < 1e-13) lam = 1e-13; end if lam == 1, CON = [0 .001 .01 (.1:.1:.9) .99 .999 1]; else omlam=1-lam; CON = [ 0 lam/100 (lam/10:lam/5:lam*9/10) lam ... (lam+omlam/10:omlam/5:lam+omlam*9/10) 1-omlam/100 1 ]; end theta=(0:.025:2*pi); C = cos(theta); S = sin(theta); hold off, clg for subplt=1:2 subplot(1,2,subplt) for c = CON, if (c <= lam) x=sqrt(c)*S; y=(sqrt(c)/sqrt(lam))*C; else t=acos(sqrt(lam/c)); dt = (pi-2*t)/1000; TT = [(t:dt:pi-t),(pi+t:dt:2*pi-t)]; SS = sin(TT); CC = cos(TT); x=sqrt(c)*SS; y=(sqrt(c)/sqrt(lam))*CC; end z=sqrt(1-x.^2-y.^2); R=find(imag(z)==0); x=x(R);y=y(R);z=z(R); if (length(x)>0), rotateP = u*[x;y;z]; xP=rotateP(1,:); yP=rotateP(2,:); zP=rotateP(3,:); rotateM = u*[x;y;-z]; xM=rotateM(1,:); yM=rotateM(2,:); zM=rotateM(3,:); x = [xP,xM]; y = [yP,yM]; z = [zP,zM]; if length(x) > 0 & subplt == 2, vue = view; vue = vue(3,1:3); HS = vue*[x;y;z]; front = find(HS >= 0); x = x(front); y = y(front); z = z(front); end if (c < lam) % Color values less than lam green plot3(x,y,z,'g.') axis('square') hold on elseif (c == lam) % Color values less than lam white plot3(x,y,z,'.w') axis('square') hold on else % Color values greater than lam red plot3(x,y,z,'r.') axis('square') hold on end end end xlabel('x = e_1') ylabel('y = e_2') zlabel('z = e_3') % title('Contour of Rayleigh Quotient. Red = max, Green = min') if subplt == 1, grid title('Transparent') else title('Opaque') end axis([-1 1 -1 1 -1 1]) end end