《MATLAB数值域的绘图演算法.doc》由会员分享，可在线阅读，更多相关《MATLAB数值域的绘图演算法.doc（4页珍藏版）》请在金锄头文库上搜索。
1、第四章 数值域的绘图演算法 设A为 n阶复数矩阵, 其数值域(numerical range)为复数子集合 W(A)=x*Ax: xCn, x*x=1.为人所知, 此集合为包含A的固有值的凸集合. 当n=2时, W(A)为一椭圆盘, 其焦距为两个固有值l, m, 短轴长为 |A|22-|l|2- |m|2 之根号. 本章将讨论W(A)的绘图演算法, 重点是方法二与方法三的MATLAB演算法.方法一 :Theorem 1 (7) Let A be an n-square complex matrix. Then W(A) = W(Axy), x,y are o.n. where Axy = x*
2、Ax x*Ay y*Ax y*Ay M2应用 Theorem 1, 只需用BASIC程式即可绘W(A), 参阅7的BASIC程式: step1 随机取两个orthonormal向量 x, y step2 绘椭圆 W(Axy) step3 goto step1重复直到图形成型.方法二 :Theorem 2 (6) Let A be an n-square complex matrix. Then the right vertical supporting line of W(eiqA) is x = l1(q), the maximal eigenvalue of the matrix Hq=(
3、eiqA+e-iqA*)/2应用 Theorem 2绘W(A)边界切线的MATLAB演算法: step1 input A step2 initialize q from 0 to 2p step3 form Hq step4 compute the maximal eigenvalue of Hq step5 draw the line l1(q)e-iq which is the supporting line of e-iq W(eiqA)=W(A) step6 goto step step2 下面我们以矩阵 a= 0 2 0 ; 0 0 1 ; 0 0 3 为例 , 画50条直线,可以看
4、出一曲线与50条直线相切,这个曲线称为50条直线的包络(envelope), 亦即W(A)的边界. clear; a=0 2 0;0 0 1;0 0 3 %a=input(matrix A=) i=sqrt(-1); H=(a+a)/2;K=(a-a)/(2*i); for k=1:50 t=(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=a1 a1 a1;B1=
5、b1 b2 b2; X(:,2*k-1)=cos(t)*A1+sin(t)*B1; Y(:,2*k-1)=-sin(t)*A1+cos(t)*B1; A2=a1 a2 a2;B2=b1 b1 b2; X(:,2*k)=cos(t)*A2+sin(t)*B2; Y(:,2*k)=-sin(t)*A2+cos(t)*B2; end %for plot(X,Y) a = 0 2 0 0 0 1 0 0 3 方法三 :Theorem 3 (3) Let A be an n-square complex matrix. Then the point pq = xq*Axq lies on the bou
6、ndary of W(A), where xq is an eigenvector of the matrix Hq corresponding to the eigenvalue l1(q)应用 Theorem 3 绘W(A)边界点的MATLAB演算法: step1 input A step2 initialize q from 0 to 2p step3 form Hq step4 compute the maximal eigenvalue l1(q)of Hq and its corresponding eigenvector xq step5 draw the point pq =
7、xq*Axq step6 goto step step2 同样以矩阵a=0 2 0; 0 0 1; 0 0 3为例,下图画50个W(A)的边界点点, 相邻两点连一线, 则此多边行逼近W(A)的边界. clg;clc;clear; axis; axis(square) i=sqrt(-1); %A=input(enter a matrix A= ); %m=input(enter steps m= ); A=0 2 0;0 0 1;0 0 3 m=50; p=; for j=1:m+1 theta=(j-1)*2*pi/m; H=(exp(i*theta)*A+exp(-i*theta)*A)/
8、2; X,L=eig(H); Y,Inde=max(diag(real(L); E=X(:,Inde)/norm(X(:,Inde); S=E*A*E; p(j)=S; end %forplot(real(p),imag(p) A = 0 2 0 0 0 1 0 0 3 方法三的演算法可用来估计W(A)的面积, 一般而言, 平面n凸多边行依反时针n个点p1,p2,pn, 则凸多边行的面积公式为 _ _ _ (Imaginary(p1 p2 + p2 p3 + pn p1)/2 clg;clc;clear; i=sqrt(-1); %A=input(enter a matrix A= ); %m
9、=input(enter steps m= ); A=0 2 0;0 0 1;0 0 3 m=50; p=; for j=1:m+1 theta=(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(:,Inde); S=E*A*E; p(j)=S; end %for area=0 ; for r=1:(m-1) area=area+imag(p(r)*conj(p(r+1)*.5; end %for area=area+imag(p(m)*conj(p(1)*.5; disp(The estimated area of W(A) is ); area A = 0 2 0 0 0 1 0 0 3The estimated area of W(A) is area = 6.7654 第四章 第 1 頁
