matlab数值域的绘图演算法

《matlab数值域的绘图演算法》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、第四章数值域的绘图演算法设A为n阶复数矩阵,其数值域(numericalrange)为复数子集合W(A)={x*Ax:xÎCn,x*x=1}.为人所知,此集合为包含A的固有值的凸集合.当n=2时,W(A)为一椭圆盘,其焦距为两个固有值l,m,短轴长为

2、A

3、22-

4、l

5、2-

6、m

7、2之根号.本章将讨论W(A)的绘图演算法,重点是方法二与方法三的MATLAB演算法.方法一:Theorem1([7])LetAbeann-squarecomplexmatrix.ThenW(A)=ÈW(Axy),x,yareo.n.whereAxy=[x*Axx*Ayy*

8、Axy*Ay]ÎM2应用Theorem1,只需用BASIC程式即可绘W(A),参阅[7]的BASIC程式:step1随机取两个orthonormal向量x,ystep2绘椭圆W(Axy)step3gotostep1重复直到图形成型.方法二:Theorem2([6])LetAbeann-squarecomplexmatrix.ThentherightverticalsupportinglineofW(eiqA)isx=l1(q),themaximaleigenvalueofthematrixHq=(eiqA+e-iqA*)/2应用Theorem2

9、绘W(A)边界切线的MATLAB演算法:step1inputAstep2initializeqfrom0to2pstep3formHqstep4computethemaximaleigenvalueofHqstep5drawthelinel1(q)e-iqwhichisthesupportinglineofe-iqW(eiqA)=W(A)step6gotostepstep2第四章第4頁下面我们以矩阵a=[020;001;003]为例,画50条直线,可以看出一曲线与50条直线相切,这个曲线称为50条直线的包络(envelope),亦即W(A)的边

10、界.clear;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)*

11、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第四章第4頁方法三:Theorem3([3])LetAbeann-squarecomplexmatrix.Thenthepointpq=xq*AxqliesontheboundaryofW(A),wherexqisaneigenvect

12、orofthematrixHqcorrespondingtotheeigenvaluel1(q)应用Theorem3绘W(A)边界点的MATLAB演算法:step1inputAstep2initializeqfrom0to2pstep3formHqstep4computethemaximaleigenvaluel1(q)ofHqanditscorrespondingeigenvectorxqstep5drawthepointpq=xq*Axqstep6gotostepstep2同样以矩阵a=[020;001;003]为例,下图画50个W(A)的

13、边界点点,相邻两点连一线,则此多边行逼近W(A)的边界.clg;clc;clear;axis;axis('square')i=sqrt(-1);%A=input('enteramatrixA=');%m=input('enterstepsm=');A=[020;001;003]m=50;p=[];forj=1:m+1theta=(j-1)*2*pi/m;H=(exp(i*theta)*A+exp(-i*theta)*A')/2;[X,L]=eig(H);[Y,Inde]=max(diag(real(L)));E=X(:,Inde)/norm(X

14、(:,Inde));S=E'*A*E;p(j)=S;end%forplot(real(p),imag(p))A=020001003第四章第4頁方法三的