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

第四章数值域的绘图演算法
设A为 n阶复数矩阵, 其数值域(numerical range)为复数子集合
W(A)={x*Ax: , 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-plex matrix. Then

W(A) = È W(Axy),
x,y are .
where Axy = [x*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-plex matrix. Then the right vertical
supporting line of W(eiqA) is x = l1(q), the maximal eigenvalue of the matrix
Hq=(eiqA+e-iqA*)/2
应用 Theorem 2绘W(A)边界切线的MATLAB演算法:
step1 input A
step2 initialize q from 0 to 2p
step3 form Hq
pute 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条直线,可以看出一曲线与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=[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(: