MATLAB数值域的绘图演算法.doc

设A为n阶复数矩阵,其数值域(numericalrange)为复数子集合W(A)={x*Ax:,x*x=1}.为人所知,=2时,W(A)为一椭圆盘,其焦距为两个固有值l,m,短轴长为|A|22-|l|2-|m|(A)的绘图演算法,:Theorem1([7])LetAbeann-(A)=ÈW(Axy),x,=[x*Axx*Ayy*Axy*Ay]ÎM2应用Theorem1,只需用BASIC程式即可绘W(A),参阅[7]的BASIC程式:step1随机取两个orthonormal向量x,ystep2绘椭圆W(Axy):Theorem2([6])LetAbeann-(eiqA)isx=l1(q),themaximaleigenvalueofthematrixHq=(eiqA+e-iqA*)/2应用Theorem2绘W(A)边界切线的MATLAB演算法:putethemaximaleigenvalueofHqstep5drawthelinel1(q)e-iqwhichisthesupportinglineofe-iqW(eiqA)=W(A)step6gotostepstep2下面我们以矩阵a=[020;001;003]为例,画50条直线,可以看出一曲线与50条直线相切,这个曲线称为50条直线的包络(envelope),亦即W(A);a=[020;001;003]%a=input('matrixA=')i=sqrt(-1);H=(a+a')/2;K=(a-a')/(2*i);fork=1:50t=(k/50)*pi;Ht=cos(t)*H-sin(t)*K;Kt=cos(t)*K+sin(t)*H;a1=max(real(eig(Ht)));b1=max(real(eig(Kt)));a2=min(real(eig(Ht)));b2=min(real(eig(Kt)));A1=[a1a1a1]';B1=[b1b2b2]';X(:,2*k-1)=cos(t)*A1+sin(t)*B1;Y(:,2*k-1)=-sin(t)*A1+cos(t)*B1;A2=[a1a2a2]';B2=[b1b1b2]';X(:,2*k)=cos(t)*A2+sin(t)*B2;Y(:,2*k)=-sin(t)*A2+cos(t)*B2;end%forplot(X,Y)a=020001003方法三:Theorem3([3])LetAbeann-=xq*AxqliesontheboundaryofW(A),wherexqisaneigenvectorofthematrixHqcorrespondingtotheeigenvaluel1(q)应用Theorem3绘W(A)边界点的MATLAB演算法:putethemaximaleigenvaluel1(q)ofHqanditscorrespondingeigenv