伪谱的定义、性质和计算

1、伪谱的定义、性质和计算近代数学选讲伪谱引言 :二十世纪九十年代以前，研究矩阵的传统工具是特征值(谱) ，它们可以揭示线性和非线性系统的特征，包括稳定性、共振、矩阵迭代的可行性等，因此它们是数学学科的一个重要的标准工具。在计算数学方面，该问题的理论和数值计算也取得了很多成果。然而，在科学和工程应用中，人们经常遇到这样的现象: 根据特征值或谱的性质所作的判断与许多观察的现象或数值结果不相匹配。究其原因，主要是这些问题所包含的矩阵往往是非正规的，甚至是高度非正规的。所以，特征值( 谱 )对分析非正规矩阵是一个不完美的工具。作为谱的自然延伸，伪谱是一个针对非正规系统的新工具。摘要 :本文首先介绍了伪谱

2、的定义及性质，然后介绍了经典的FOVT法来粗略地给出了伪谱范围的矩形界定，之后介绍了伪谱计算的两种方法，即随机扰动法和SVD方法，最后给出了伪谱的一个应用。 关键字 :伪谱定义及性质、矩形界定、伪谱计算记号及说明 :HA文中所有矩阵均为定义在复数域上的阶方阵，表示矩阵的AnI 共轭转置，表示相应阶的单位阵，分别表示复数的实 Re(z)Im(z)z 部，虚部数值，表示以为中心，为半径的闭圆域。 B(z,)z,正文 :伪谱的定义及性质 伪谱的定义:n, nAA假定有矩阵，的谱是指矩阵白特征值的全体，可表示 A,C如下 :,(A),z,C:det(zI,A),0,1(zI,A) 我们知道，当 z,(

3、A) 时，是没有意义的。如果我们定义,1,1|(zI,A)|,， ,|(zI,A)| 。那么当有限而且非常大时，又会如何 , 这A就导致了伪谱最初的一个定义。给定，(A),矩阵的伪谱(),，定义：,1,1,(A),z,C:|(zI,A)|,， (1) ,等价地，伪谱也可以用扰动矩阵的特征值来定义:,(A),z,C:z,(A， E),|E|,， (2) ,AA也就是说，的伪谱是的任何一个扰动矩阵的特征值全体。,AAz,(A) 为矩阵的一个伪特征值，简称的伪特征值。相称,应的，每一个伪特征值有一个特征向量(一般不唯一) ，这样就导致了伪谱的第三种定义 :n,(A),z,C:|(A,zI)v|,，

4、v,C,|v|,1， (3) ,1,1由于如果是2-范数，定义(1)可|(zI,A)|,(,(zI,A)|,|min改为如下形式，即得伪谱的第四种定义:,(A),z,C:,(zI,A),， (4) ,min这种表示在计算机上更为方便。容易证明，伪谱的上述四种定义是等价的。伪谱的性质:, ， ,(A),(,I， ,A) 性质 1( 线性性 ): 对任意，成立。 ,R,|,| 证明 : 分为两步来证 :,(A ， ,I),， ,(A) 先证 : ,1,1,(A ， ,I),z,C:|(zI,(A， ,I)|, 由伪谱的定义(1) 知， ,即有,1,1,(A , ,I),z,C:|(z,)I,A)|

5、,1,1,(A , ,I),(w , ,),C:|(wI,A)|, ,(A)令，则有，z,w,(,A),(A)再证：，|,|,1,1,1由伪谱的定义 知，,(A),z,C:|(zI,A)|,|,|,|z,1,1,(A),z,C:|(I,A)|,即有，,|,1,1 令，则有，,(A),(,v),C:|(vI,A)|,(A)z,v,|,(A),(,I-A)故有，。,|,|0,(A),(A) 性质 2(单调性)：若，则有，并且,(A) , B(0,),(A),证明：分为两部分来证：先证：若则有。，(A),(A)0,E对任意的，由伪谱定义(2)知，存在一个矩阵()，z,(A)|E|,使得，z,(A ,

6、 E)显然，又，则乙,(A)|E|,z,(A, E),从而有，。，(A),(A),0,则有,(A),(A), B(0,)。再证若,，n对任意的乙,(A),由伪谱定义(3)知，存在()，使得|u|,1u,C,一,对任意的满足，即。|(zI,A)u|,z,B(0,)|z|,|zu|,一从而，|(z, z)I,A)u|,|(zI,A)u, zu|,|(zI,A)u|,|zu|,。A性质3:对任意矩阵都有,(A) , ,(A),0,z,C:|z|,其中。,z,(A) , ,z, , z(,(A),z,)证明：，有，由于，即是 z|z,|,A, zI的特征值。z,(A)又由于由伪谱的定义(2)知，|zI

7、|,(A) , ,(A)故,0,(A),(A)矩阵伪谱是复平面上的一个闭子集，当时，有。 0当足够小时,伪谱为其矩阵特征值周围的一个个联通闭区域。我们 ，A称每一个联通闭区域为伪谱的一个联通部分。如果矩阵有个不同mA的特征值，则至多有个不同的伪谱部分。m作为例子，我们取了一个3阶的上三角矩阵，225,A,0i3 , ,00,1,A当依次取值为0.05,0.1,0.3 时，矩阵的伪谱依次有3,2,1个联通,“部分，如图所示，图中仅仅画出的是其相应伪谱的边界曲线。随着，不断增大，区域开始坍塌，不同的联通部分开始合并在一起。不A,|,|难得到，当时，的伪谱必'为复平面上的单联通，,ijmax

8、,(A),ij区域。3251fl. 50田苫'11.5伪诰图伪谱范围的矩形界定虽然减少矩阵的阶数能估计局部的伪谱信息，但不一定可靠。因2,(A), 此为了减少计算量，可以选择尽可能小，并且，从而 ,C, 减少网格的点数。下面我们介绍经典的FOV方法。FOVT 法A通过估计的数值域HnF(A),xAx|x,C,|x|,1，A得到了的伪谱区域的一个闭包，该方法称为FOVJ法。FOV,(A),方法的理论基础是下述定理，它说明可以通过扩大数值域的一个,带域来界定伪谱区域。定理 1: 对于任意的，有,(A),F(A) ， B(0,) 。 ,A证明：根据伪谱的定义，任何x,(A)都是的一个扰动后矩

9、,E 阵的特征值。因而存在矩阵( 满足 |E|,) 和单位向量，使得xHHH|xEx|,(A ， E)x,zx 。进而得到。由于，故z,xAx ， xEx,(A),F(A) ， B(0,)z,F(A)， B(0,) ，即。 ,A 矩阵 F(A) 的数值值域有如下二个性质:HH(A),(A ， A)/2F(H(A),Re(F(A) ，其中 ; 1)HS(A),(A,A)/2F(S(A),i,Im(F(A) ，这里。 2)证明 :(1) 先证 : F(H(A),Re(F(A)y,F(H(A) 对任意的，即存在，其中满足，使得 xx|x|,1H1HHHHxA，x xAx,yxH(A)x,y ，即 H

10、H2HHHHy,xAxy,xAxy,F(A) 若记，显然有且AAA1HRe(y),(y ， y),y,Re(F(A) 从而有， AAAH2再证 Re(F(A),F(H(A):H对任意的，即存在，其中满足，使得， Re(xAx),yxxy,Re(F(A)|x|,1 1HHH1P xAx , xAx,y , 2从而 y,F(H(A)H故有，其中 H(A),(A , A)/2 F(H(A),Re(F(A)。(2)性质 3)的证明与性质2) 的证明同理。上述后两个性质可用于 ( 矩阵界定 ) 伪谱区域 : 令分别为 a,bH(A)的最小，最大特征值，分别为的特征值最小，最大虚部。则 c,dS(A)，进

11、而 F(A),a,b， c,dFOV,(A),a,b， , ， c,d ， , ,这是目前所能得到的一个最好的伪谱界定区域。然而实验结果表明，往往这样界定的区域要远远比 ,(A) 大得多。 ,A例1:矩阵如下定义，用FOVT法求其伪谱界定区域。1111,1111 ？ ,111 ？ 1,A, ？1,111,11,32， 32,3,10 其中让。运行结果 ( 程序 1，见附录 ):险 WU rr*“七 匕，柳g 5由运行结果知，FOVT法得到的伪谱界定区域大致在0.016027,1.984, ,2.5762,2.5762伪谱的计算从上面可知，伪谱有四种等价定义，简单地说计算方法可以不同，但在同一范

12、 数下得到的伪谱集合是一样的。从矩阵伪谱的定义出发，直接计算，可分为随机扰 动法和格点SVD法(又称奇异值分解法)。直接计算思路清晰便于操作，缺点在于计 算量相对比较大。由于集合内有无穷多个点，要全部计算这些点是不实际的也是没 有必要的，有时只要能够把边界刻画出来就可以了。这里给出几种简单的伪谱计算 方法。一、随机扰动法A为了求矩阵的伪谱。首先产生一列随机的复矩阵，(A),N,并且这些矩阵满足条件|E|,然后求出这个矩阵E,(i,1,2,? ,n)ii的特征值，并在复平面上画出，所得到的点集能近似表示(A, E)i,这种方法最简单。，(A),例2:考虑例1中的矩阵，确定其伪谱。值问题本身就比较

13、复杂不容易计算，所以算法不是很通用，不过却可以很好地加深对伪谱概念的理解。、格点SVD?去当矩阵范数取2-范数时，此时我们采用定义（4）来定义伪谱，,伪谱边界就是最小奇异值为的格点连成的曲线。计算并找出该区域基本步骤:E?用MATLAB!机产生一个稠密矩阵，然后令,即可实现。EE?|E|,3扰动为，扰动矩阵个数为50, 100, 150, 200时的伪谱图。,10虽然扰动矩阵个数为200的图比较精确，代价是需要计算很多特征值，而特征内奇异值为的格点，将他们连接起来以确定伪谱的边界。运行结果（程序2,见附录）:忱* F作用】的精*国M L1)将感兴趣的区域划分为份，x(k) , iy(j)为复平

14、面上格点的坐标，其中x(k) , y(j)分别表示对应格点(k,j)的横坐标值和纵坐标值(k,1:M,j,1:L),为虚数 单位；i(zI,A),(x(k), iy(j)I,A2)计算的最小奇异值,(x(k) , iy(j)I,A; min3)通过绘制伪谱边界上的点(即，最小奇异值为的那些点)，就可以得到伪谱边界曲线图。xA例3:考虑例1中的矩阵伪谱，让,10运行结果(程序3,见附录)：(3-a) (3-b)图(3-a)画的是时伪谱的边界图，从图中可以看出随着的，x,4,3,2x,3增大，伪谱所在区域范围扩大(验证了性质2);图(3-b)画的是时伪谱所在的区域，与随机扰动法得到的结果相同。与随

15、机扰动法相比较，SVDt不需要计算特征值，只是计算了奇异值，但是后者的计算要比前者容易的多，也就是说该方法有更好的可行性。虽然随机扰动法和格点SVDt能计算伪谱的范围，但其计算量往往很大。目前还有一些比较好的计算伪谱的方法，如：区域排除法、格点移动法、投影算法、QR分解法等等。伪谱的应用例4:求解下列算子的伪谱：124,Au(x),2u , (cx,dx)u,c,3 , 3i,d,u(,L),0,L,10. 16首先我们用一个切比雪夫点列来离散有限区间并注意边界,L,Ln, 1条件。考虑有个离散点列的情况，让 u,(,L),0j,x,Lcos(),j,0,1,? ,n , 1。点x处的二阶导数

16、差分格式是jjn , 14(u,2u , u)j , jj,11,。(注意：这里 u 表示 u(x)从而u(x),j,1,2,? ,njjj2(x,x)j , j,11将算子离散后得到矩阵：,An, n168, , 24cx,dx,0 ? ? 01122,(x,x)(x,x)2020,816824,cxdx? ?022222,(x,x)(x,x)(x,x)313131,81624,0cx,dx,? ?03322(x,x)(x,x),4242,?,81624000cx,dx,nn,22(x,x)(x,x)n, 1n,1n , 1n,1 ,从理论上说，越大，离散化的程度就越高，矩阵A就可以非nn,

17、 n常逼近原系统。这个矩阵具有众多的特征值，随着矩阵A的规模变n, n大，计算消耗的时间越来越长。n,200这里取，先通过FOVT法大致估计一下伪谱的范围。运行结果：V1 nd Liv艮di IIL-skf oj Im俄区域的左世界是"3仃港噌:其由江世群号i= 3m *M国艳的汇寄W产L53考区域的上透界导h 12S129. rsr|FW方法荐剥的伪谱区城1 60 50051 5由FOVT法得到了伪谱可能存在的大致范围为,这为我们用SVDt确定伪谱白范围提供,747192,32 , ,127531,128128了需要划分的区域。现在，在我们感兴趣的区域内计算伪谱。图,200,30

18、, ,200,3023, 2346, 46(4-a)、(4-b)给出了网格划分分别为与时所得矩阵 v, vA的伪谱的边界曲线。运行每一幅图的程序所花的时间不等。,200 , 200,1,2,6,10,10,? ,10图中，依次从外到内，并且“ *”表示原系统矩阵的特征值RPIMiH i.rw(4-a) (4-b)从这个例子可以看出用FOVT法确定的伪谱界定区域远远大于伪谱的真实区域。参考文献 :1、矩阵拟谱计算的若干研究，加帮平， 2008年 6月 ;2、伪谱的边界曲线及其跟踪算法的步长控制，刘颖，白峰杉，高等学校计算数学学报;3、基于区域排除法和方格移动法的矩阵伪谱计算，周剑，蒋耀林，数值计

19、算与计算机应用 ;4、大规模矩阵伪谱计算的数值方法，孟青云， 2010年 2月 ;附录 :程序 1:A=;n=32;%A矩阵的阶数；%定义A矩阵；for i=1:1:nfor j=1:1:nA(i,j)=0;if i-1>=1A(i,i-1)=-1;endif i=jA(i,j)=1;endif i+1<=nA(i,i+1)=1;endif i+2<=nA(i,i+2)=1;endendendHA=(A+A)./2;% 定义 H(A)矩阵；VHA=eig(HA);%< H(A)的特征值;a=min(VHA);%H(A)特征值的最小值;b=max(VHA);%H(A特征值

20、的最大值;SA=(A-A')./2;% 定义 S(A) 矩阵 ;VSA=eig(SA);%求 S(A)的特征值;imVSA=imag(VSA);%ft S(,)特征值的虚部；c=min(imVSA);%S(A)特征值虚部的最小值 ; d=max(imVSA);%S(A) 特征值虚部的最大值; e=0.001;% 定义 epsilon;y1=c-e;% 区域的下边界;y2=d+e;%区域的上边界；x1=a-e;% 区域的左边界;x2=b+e;%区域的右边界；disp('该区域的左边界是x= 'num2str(x1);disp('该区域的右边界是x= 'nu

21、m2str(x2);disp('该区域的下边界是y= 'num2str(y1);disp('该区域的上边界是y= 'num2str(y2);w=x2-x1;%该区域白宽度；h=y2-y1;%该区域白高度；rectangle('Position',x1,y1,w,h,'LineWidth',2,'LineStyle','-,'EdgeColor','m')xlim(-0.5,2.5)ylim(-3,3)title('FOV 方法得到的伪谱区域')程序 2:A=;X=;Y=;n=32;%定义A的阶数；%定义A 矩阵 ;for i=1:1:nfor j=1:1:nA(i,j)=0;if i-1>=1A(i,i-1)=-1;endif i=jA(i,j)=1;endif i+1<=nA(i,i+1)=1;endif i+2<=nA(i,i+2)=1;endendendm=50;%t义随机矩阵的个数；for k=1:1:mE1=rand(n);% 产生随机矩阵;DE=det(E1);%求随机矩阵的行列式值;e=0.001;%定义 epsil