伪谱的定义、性质和计算

近代数学选讲

——伪谱

引言：

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

摘要：

本文首先介绍了伪谱的定义及性质，然后介绍了经典的FOV方法来粗略地给出了伪谱范围的矩形界定，之后介绍了伪谱计算的两种方法，即随机扰动法和SVD方法，最后给出了伪谱的一个应用。

关键字：

伪谱定义及性质、矩形界定、伪谱计算

记号及说明：

文中所有矩阵均为定义在复数域上的n阶方阵，A表示矩阵A的共轭转置，I表示相

H应阶的单位阵，Re(z)，Im(z)分别表示复数z的实部，虚部数值，B(z,)表示以z为中心，为半径的闭圆域。

正文：

伪谱的定义及性质

伪谱的定义：

假定有矩阵ACnn，A的谱是指矩阵A的特征值的全体，可表示如下：

(A){zC:det(zIA)0}

11||(zIA)||。那么(zIA)z(A)我们知道，当时，是没有意义的。如果我们定义1||(zIA)||有限而且非常大时，又会如何？这就导致了伪谱最初的一个定义。给定，矩当

阵A的伪谱（(A)）定义：

(A){zC:||(zIA)1||1}，（1）

等价地，伪谱也可以用扰动矩阵的特征值来定义：

(A){zC:z(AE),||E||}，（2）

也就是说，A的伪谱是A的任何一个扰动矩阵的特征值全体。

称z(A)为矩阵A的一个伪特征值，简称A的伪特征值。相应的，每一个伪特征值有一个特征向量（一般不唯一），这样就导致了伪谱的第三种定义：

(A){zC:||(AzI)v||,vCn,||v||1}，（3）

11||(zIA)||((zIA))||•||min由于如果是2-范数，，定义（1）可改为如下形式，即得

伪谱的第四种定义：

(A){zC:min(zIA)}，（4）

这种表示在计算机上更为方便。

容易证明，伪谱的上述四种定义是等价的。

伪谱的性质：

(A)||||(IA)性质1（线性性）：对任意,R，成立

。

证明：分为两步来证： 先证：(AI)(A)

由伪谱的定义（1）知，

(AI){zC:||(zI(AI))1||1}

即有，

(AI){zC:||((z)IA)1||1}

令zw，则有，

(AI){(w)C:||(wIA)1||1}(A)

再证：

||||(A)(A)

由伪谱的定义（1）知，

||||(A){zC:||(zIA)1||||||11}

即有，

||||(A){zC:||(zIA)1||1}

||||(A){(v)C:||(vIA)1||1}(A)zv令，则有，

故有，

(A)||||(IA)。

性质2（单调性）：若0，则有(A)(A)，并且(A)B(0,)(A)

证明：分为两部分来证：

先证：若0,则有(A)(A)。

对任意的z(A)，由伪谱定义（2）知，存在一个矩阵E（||E||），使得z(AE)，

显然||E||，又z(AE)，则z(A)

从而有，(A)(A)。

再证若0，则有(A)(A)B(0,)。

对任意的z(A)，由伪谱定义（3）知，存在uC（||u||1），使得||(zIA)u||，

n~~~zB(0,)||z||||对任意的满足，即zu||。

~~~从而，||((zz)IA)u||||(zIA)uzu||||(zIA)u||||zu||。

性质3：对任意矩阵A都有

(A)(A),0，其中{zC:|z|}。

z~z((A),~z)~|zI的特征值。，由于z|，即z是A~

证明：

z(A)，有

~|又由于zI|由伪谱的定义（2）知，z(A)

故(A)(A)。

矩阵伪谱是复平面上的一个闭子集，当0时，有0(A)(A)。当足够小时，伪谱为其矩阵特征值周围的一个个联通闭区域。我们称每一个联通闭区域为伪谱的一个联通部分。如果矩阵A有m个不同的特征值，则A至多有m个不同的伪谱部分。

作为例子，我们取了一个3阶的上三角矩阵，

225A0i3001

当依次取值为0.05,0.1,0.3时，矩阵A的伪谱依次有3,2,1个联通部分，如图所示，

图中仅仅画出的是其相应伪谱的边界曲线。随着不断增大，区域开始‘坍塌’，不同的联通部分开始合并在一起。不难得到，当单联通区域。

max|ij|i,j(A)时，A的伪谱必为复平面上的

伪谱范围的矩形界定

虽然减少矩阵的阶数能估计局部的伪谱信息，但不一定可靠。因此为了减少计算量，可以选择C尽可能小，并且(A)，从而减少网格的点数。下面我们介绍经典的FOV

2方法。

FOV方法

通过估计A的数值域

F(A){xHAx|xCn,||x||1}，

得到了A的伪谱区域(A)的一个闭包，该方法称为FOV方法。FOV方法的理论基础是下述定理，它说明可以通过扩大数值域的一个带域来界定伪谱区域。

定理1：对于任意的，有(A)F(A)B(0,)。

证明：根据伪谱的定义，任何x(A)都是A的一个扰动后矩阵的特征值。因而存

HH||E||(AE)xzxzxAxxEx。由于xE在矩阵（满足）和单位向量，使得。进而得到

||xHEx||，故zF(A)B(0,)，即(A)F(A)B(0,)。

矩阵A的数值值域F(A)有如下二个性质：

H1）F(H(A))Re(F(A))，其中H(A)(AA)/2；

HS(A)(AA)/2。 F(S(A))iIm(F(A))2），这里

证明：（1）先证F(H(A))Re(F(A))：

HxyF(H(A))||x||1xxH对任意的，即存在，其中满足，使得H(A)xyH，即

1H[xAxxHAHx]yH2

HHHyxAx且yAF(A) yxAxA若记A，显然有

H从而有，

Re(yA)1H(yAyA)yHRe(F(A))2

再证Re(F(A))F(H(A))：

HRe(xAx)y，即yRe(F(A))||x||1xx对任意的，即存在，其中满足，使得

1H[xAxxHAHx]y2，

从而yF(H(A))

HH(A)(AA)/2。 F(H(A))Re(F(A))故有，，其中

（2）性质3）的证明与性质2）的证明同理。

上述后两个性质可用于（矩阵界定）伪谱区域：令a,b分别为H(A)的最小，最大特征值，c,d分别为S(A)的特征值最小，最大虚部。则F(A)[a,b][c,d]，进而

(A)[a,b][c,d]FOV

这是目前所能得到的一个最好的伪谱界定区域。然而实验结果表明，往往这样界定的区域要远远比(A)大得多。

例1：矩阵A如下定义，用FOV方法求其伪谱界定区域。

11111111111A11111113232

3其中让10。

运行结果（程序1，见附录）：

由运行结果知，FOV方法得到的伪谱界定区域大致在[0.016027,1.984][2.5762,2.5762]

伪谱的计算

从上面可知，伪谱有四种等价定义，简单地说计算方法可以不同，但在同一范数下得到的伪谱集合是一样的。从矩阵伪谱的定义出发，直接计算，可分为随机扰动法和格点SVD法（又称奇异值分解法）。直接计算思路清晰便于操作，缺点在于计算量相对比较大。由于

集合内有无穷多个点，要全部计算这些点是不实际的也是没有必要的，有时只要能够把边界刻画出来就可以了。这里给出几种简单的伪谱计算方法。

一、随机扰动法

为了求矩阵A的伪谱(A)。首先产生一列随机的复矩阵Ei,(i1,2,,n)，并且这些矩阵满足条件||Ei||，然后求出这N个矩阵(AEi)的特征值，并在复平面上画出，所得到的点集能近似表示(A)，这种方法最简单。

例2：考虑例1中的矩阵，确定其伪谱。

EE3ˆ||ˆ||E10E用MATLAB随机产生一个稠密矩阵，然后令即可实现。扰动为，扰

动矩阵个数为50，100，150，200时的伪谱图。

运行结果（程序2，见附录）：

虽然扰动矩阵个数为200的图比较精确，代价是需要计算很多特征值，而特征值问题

本身就比较复杂不容易计算，所以算法不是很通用，不过却可以很好地加深对伪谱概念的理解。

二、格点SVD法

当矩阵范数取2-范数时，此时我们采用定义（4）来定义伪谱，伪谱边界就是最小奇异值为的格点连成的曲线。计算并找出该区域内奇异值为的格点，将他们连接起来以确定伪谱的边界。

基本步骤：

1）将感兴趣的区域划分为ML份，x(k)iy(j)为复平面上格点的坐标，其中x(k)，y(j)分别表示对应格点(k,j)的横坐标值和纵坐标值(k1:M,j1:L)，i为虚数单位；

2）计算(zIA)[(x(k)iy(j))IA]的最小奇异值min[(x(k)iy(j))IA]；

3）通过绘制伪谱边界上的点（即，最小奇异值为的那些点）就可以得到伪谱边界曲线图。

x10A例3：考虑例1中的矩阵伪谱，让

运行结果（程序3，见附录）：

(3-a) (3-b)

图(3-a)画的是x4,3,2时伪谱的边界图，从图中可以看出随着的增大，伪谱所在区域范围扩大（验证了性质2）；图(3-b)画的是x3时伪谱所在的区域，与随机扰动法得到的结果相同。

与随机扰动法相比较，SVD法不需要计算特征值，只是计算了奇异值，但是后者的计算要比前者容易的多，也就是说该方法有更好的可行性。

虽然随机扰动法和格点SVD法能计算伪谱的范围，但其计算量往往很大。目前还有一

些比较好的计算伪谱的方法，如：区域排除法、格点移动法、投影算法、QR分解法等等。

伪谱的应用

例4：求解下列算子的伪谱：

1,u(L)0,L10.16

Au(x)2u(cx2dx4)u,c33i,d首先我们用一个切比雪夫点列来离散有限区间[L,L]并注意边界条件u(L)0。考虑

jxjLcos(),j0,1,,n1xn1n1有个离散点列的情况，让。点j处的二阶导数差分格式是

u(xj)4(uj12ujuj1)(xj1xj1)2,j1,2,,n。（注意：这里

uj表示

u(xj)）从而将算子离散后得到矩

阵：

Ann16824cxdx0011(x2x0)2(x2x0)2816824cxdx022222(x3x1)(x3x1)(x3x1)816240cx3dx3022(xx)(xx)424281624000cxdxnn(xn1xn1)2(xn1xn1)2

从理论上说，n越大，离散化的程度就越高，矩阵Ann就可以非常逼近原系统。这个矩阵具有众多的特征值，随着矩阵Ann的规模变大，计算消耗的时间越来越长。

这里取n200，先通过FOV方法大致估计一下伪谱的范围。

运行结果：

由FOV方法得到了伪谱可能存在的大致范围为[747192,32][127531,128128]，这为我们用SVD法确定伪谱的范围提供了需要划分的区域。

现在，在我们感兴趣的区域[200,30][200,30]内计算伪谱。图(4-a)、(4-b)给出了网格划分vv分别为2323与4646时所得矩阵A200200的伪谱的边界曲线。运行每一幅图的程

12610,10,,10序所花的时间不等。图中，依次从外到内，并且“*”表示原系统矩阵的

特征值。

(4-a) (4-b)

从这个例子可以看出用FOV方法确定的伪谱界定区域远远大于伪谱的真实区域。

参考文献：

1、《矩阵拟谱计算的若干研究》，加帮平，2008年6月；

2、《伪谱的边界曲线及其跟踪算法的步长控制》，刘颖，白峰杉，高等学校计算数学学报；

3、《基于区域排除法和方格移动法的矩阵伪谱计算》，周剑，蒋耀林，数值计算与计算机应用；

4、《大规模矩阵伪谱计算的数值方法》，孟青云，2010年2月；

附录：

程序1：

A=[];

n=32;%A矩阵的阶数；

%定义A矩阵；

for i=1:1:n

for j=1:1:n

A(i,j)=0;

if i-1>=1

A(i,i-1)=-1;

end

if i==j

A(i,j)=1;

end

if i+1<=n

A(i,i+1)=1;

end

if i+2<=n

A(i,i+2)=1;

end

end

end

HA=(A+A')./2;%定义H(A)矩阵；

VHA=eig(HA);%求H(A)的特征值；

a=min(VHA);%H(A)特征值的最小值；

b=max(VHA);%H（A）特征值的最大值；

SA=(A-A')./2;%定义S(A)矩阵；

VSA=eig(SA);%求S(A)的特征值；

imVSA=imag(VSA);%求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= ' num2str(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:n

for j=1:1:n

A(i,j)=0;

if i-1>=1

A(i,i-1)=-1;

end

if i==j

A(i,j)=1;

end

if i+1<=n

A(i,i+1)=1;

end

if i+2<=n

A(i,i+2)=1;

end

end

end

m=50;%定义随机矩阵的个数；

for k=1:1:m

E1=rand(n);%产生随机矩阵；

DE=det(E1);%求随机矩阵的行列式值；

e=0.001;%定义epsilon；

a=e/DE;

E=a.*E1;%得到行列式值为epsilon的随机矩阵；

V=eig(A+E);%求（A+E）的特征值；

reV=real(V);%求（A+E）特征值的实部；

imV=imag(V);%求（A+E）特征值的虚部；

X=[X,reV];

Y=[Y,imV];

end

plot(X,Y,'.m')

title('扰动矩阵个数为50的伪谱图')

程序3：

A=[];

n=32; %矩阵A的阶数

%定义矩阵A;

for k=1:1:n

for j=1:1:n

A(k,j)=0;

if k-1>=1

A(k,k-1)=-1;

end

if k==j

A(k,j)=1;

end

if k+1<=n

A(k,k+1)=1;

end

if k+2<=n

A(k,k+2)=1;

end

end

end

%感兴趣的区域是[0,2]×[-2.5,2.5]；

e=eig(A);

ex=real(e);ey=imag(e);

x=-2.5:0.05:2.5;

y=-2.5:0.05:2.5;

l1=length(x);

l2=length(y);

for j=1:l1

for k=1:l2

svdmin(j,k)=min(svd((x(j)+i*y(k))*eye(n)-A)); % 注意在此之前不要出现i, 否则会用最后一步的值

end

end

z=log10(svdmin);

[C,h]=contour(x,y,z.',-4:-2); %里面的数是epsilon取了对数

clabel(C,h);

hold on

plot(ex,ey,'*') % 画谱点

xlim([-0.5,2.5])

ylim([-3,3])

title('伪谱图')