偽譜的定義、性質(zhì)和計(jì)算

1、偽譜的定義、性質(zhì)和計(jì)算近代數(shù)學(xué)選講偽譜引言 :二十世紀(jì)九十年代以前，研究矩陣的傳統(tǒng)工具是特征值(譜) ，它們可以揭示線性和非線性系統(tǒng)的特征，包括穩(wěn)定性、共振、矩陣迭代的可行性等，因此它們是數(shù)學(xué)學(xué)科的一個(gè)重要的標(biāo)準(zhǔn)工具。在計(jì)算數(shù)學(xué)方面，該問題的理論和數(shù)值計(jì)算也取得了很多成果。然而，在科學(xué)和工程應(yīng)用中，人們經(jīng)常遇到這樣的現(xiàn)象: 根據(jù)特征值或譜的性質(zhì)所作的判斷與許多觀察的現(xiàn)象或數(shù)值結(jié)果不相匹配。究其原因，主要是這些問題所包含的矩陣往往是非正規(guī)的，甚至是高度非正規(guī)的。所以，特征值( 譜 )對(duì)分析非正規(guī)矩陣是一個(gè)不完美的工具。作為譜的自然延伸，偽譜是一個(gè)針對(duì)非正規(guī)系統(tǒng)的新工具。摘要 :本文首先介紹了偽譜

2、的定義及性質(zhì)，然后介紹了經(jīng)典的FOVT法來粗略地給出了偽譜范圍的矩形界定，之后介紹了偽譜計(jì)算的兩種方法，即隨機(jī)擾動(dòng)法和SVD方法，最后給出了偽譜的一個(gè)應(yīng)用。 關(guān)鍵字 :偽譜定義及性質(zhì)、矩形界定、偽譜計(jì)算記號(hào)及說明 :HA文中所有矩陣均為定義在復(fù)數(shù)域上的階方陣，表示矩陣的AnI 共軛轉(zhuǎn)置，表示相應(yīng)階的單位陣，分別表示復(fù)數(shù)的實(shí) Re(z)Im(z)z 部，虛部數(shù)值，表示以為中心，為半徑的閉圓域。 B(z,)z,正文 :偽譜的定義及性質(zhì) 偽譜的定義:n, nAA假定有矩陣，的譜是指矩陣白特征值的全體，可表示 A,C如下 :,(A),z,C:det(zI,A),0,1(zI,A) 我們知道，當(dāng) z,(

3、A) 時(shí)，是沒有意義的。如果我們定義,1,1|(zI,A)|,， ,|(zI,A)| 。那么當(dāng)有限而且非常大時(shí)，又會(huì)如何 , 這A就導(dǎo)致了偽譜最初的一個(gè)定義。給定，(A),矩陣的偽譜(),，定義：,1,1,(A),z,C:|(zI,A)|,， (1) ,等價(jià)地，偽譜也可以用擾動(dòng)矩陣的特征值來定義:,(A),z,C:z,(A， E),|E|,， (2) ,AA也就是說，的偽譜是的任何一個(gè)擾動(dòng)矩陣的特征值全體。,AAz,(A) 為矩陣的一個(gè)偽特征值，簡(jiǎn)稱的偽特征值。相稱,應(yīng)的，每一個(gè)偽特征值有一個(gè)特征向量(一般不唯一) ，這樣就導(dǎo)致了偽譜的第三種定義 :n,(A),z,C:|(A,zI)v|,，

4、v,C,|v|,1， (3) ,1,1由于如果是2-范數(shù)，定義(1)可|(zI,A)|,(,(zI,A)|,|min改為如下形式，即得偽譜的第四種定義:,(A),z,C:,(zI,A),， (4) ,min這種表示在計(jì)算機(jī)上更為方便。容易證明，偽譜的上述四種定義是等價(jià)的。偽譜的性質(zhì):, ， ,(A),(,I， ,A) 性質(zhì) 1( 線性性 ): 對(duì)任意，成立。 ,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) 性質(zhì) 2(單調(diào)性)：若，則有，并且,(A) , B(0,),(A),證明：分為兩部分來證：先證：若則有。，(A),(A)0,E對(duì)任意的，由偽譜定義(2)知，存在一個(gè)矩陣()，z,(A)|E|,使得，z,(A ,

6、 E)顯然，又，則乙,(A)|E|,z,(A, E),從而有，。，(A),(A),0,則有,(A),(A), B(0,)。再證若,，n對(duì)任意的乙,(A),由偽譜定義(3)知，存在()，使得|u|,1u,C,一,對(duì)任意的滿足，即。|(zI,A)u|,z,B(0,)|z|,|zu|,一從而，|(z, z)I,A)u|,|(zI,A)u, zu|,|(zI,A)u|,|zu|,。A性質(zhì)3:對(duì)任意矩陣都有,(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)矩陣偽譜是復(fù)平面上的一個(gè)閉子集，當(dāng)時(shí)，有。 0當(dāng)足夠小時(shí),偽譜為其矩陣特征值周圍的一個(gè)個(gè)聯(lián)通閉區(qū)域。我們 ，A稱每一個(gè)聯(lián)通閉區(qū)域?yàn)閭巫V的一個(gè)聯(lián)通部分。如果矩陣有個(gè)不同mA的特征值，則至多有個(gè)不同的偽譜部分。m作為例子，我們?nèi)×艘粋€(gè)3階的上三角矩陣，225,A,0i3 , ,00,1,A當(dāng)依次取值為0.05,0.1,0.3 時(shí)，矩陣的偽譜依次有3,2,1個(gè)聯(lián)通,“部分，如圖所示，圖中僅僅畫出的是其相應(yīng)偽譜的邊界曲線。隨著，不斷增大，區(qū)域開始坍塌，不同的聯(lián)通部分開始合并在一起。不A,|,|難得到，當(dāng)時(shí)，的偽譜必'為復(fù)平面上的單聯(lián)通，,ijmax

8、,(A),ij區(qū)域。3251fl. 50田苫'11.5偽誥圖偽譜范圍的矩形界定雖然減少矩陣的階數(shù)能估計(jì)局部的偽譜信息，但不一定可靠。因2,(A), 此為了減少計(jì)算量，可以選擇盡可能小，并且，從而 ,C, 減少網(wǎng)格的點(diǎn)數(shù)。下面我們介紹經(jīng)典的FOV方法。FOVT 法A通過估計(jì)的數(shù)值域HnF(A),xAx|x,C,|x|,1，A得到了的偽譜區(qū)域的一個(gè)閉包，該方法稱為FOVJ法。FOV,(A),方法的理論基礎(chǔ)是下述定理，它說明可以通過擴(kuò)大數(shù)值域的一個(gè),帶域來界定偽譜區(qū)域。定理 1: 對(duì)于任意的，有,(A),F(A) ， B(0,) 。 ,A證明：根據(jù)偽譜的定義，任何x,(A)都是的一個(gè)擾動(dòng)后矩

9、,E 陣的特征值。因而存在矩陣( 滿足 |E|,) 和單位向量，使得xHHH|xEx|,(A ， E)x,zx 。進(jìn)而得到。由于，故z,xAx ， xEx,(A),F(A) ， B(0,)z,F(A)， B(0,) ，即。 ,A 矩陣 F(A) 的數(shù)值值域有如下二個(gè)性質(zhì):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) 對(duì)任意的，即存在，其中滿足，使得 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對(duì)任意的，即存在，其中滿足，使得， 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)性質(zhì) 3)的證明與性質(zhì)2) 的證明同理。上述后兩個(gè)性質(zhì)可用于 ( 矩陣界定 ) 偽譜區(qū)域 : 令分別為 a,bH(A)的最小，最大特征值，分別為的特征值最小，最大虛部。則 c,dS(A)，進(jìn)

11、而 F(A),a,b， c,dFOV,(A),a,b， , ， c,d ， , ,這是目前所能得到的一個(gè)最好的偽譜界定區(qū)域。然而實(shí)驗(yàn)結(jié)果表明，往往這樣界定的區(qū)域要遠(yuǎn)遠(yuǎn)比 ,(A) 大得多。 ,A例1:矩陣如下定義，用FOVT法求其偽譜界定區(qū)域。1111,1111 ？ ,111 ？ 1,A, ？1,111,11,32， 32,3,10 其中讓。運(yùn)行結(jié)果 ( 程序 1，見附錄 ):險(xiǎn) WU rr*“七 匕，柳g 5由運(yùn)行結(jié)果知，F(xiàn)OVT法得到的偽譜界定區(qū)域大致在0.016027,1.984, ,2.5762,2.5762偽譜的計(jì)算從上面可知，偽譜有四種等價(jià)定義，簡(jiǎn)單地說計(jì)算方法可以不同，但在同一范

12、 數(shù)下得到的偽譜集合是一樣的。從矩陣偽譜的定義出發(fā)，直接計(jì)算，可分為隨機(jī)擾 動(dòng)法和格點(diǎn)SVD法(又稱奇異值分解法)。直接計(jì)算思路清晰便于操作，缺點(diǎn)在于計(jì) 算量相對(duì)比較大。由于集合內(nèi)有無窮多個(gè)點(diǎn)，要全部計(jì)算這些點(diǎn)是不實(shí)際的也是沒 有必要的，有時(shí)只要能夠把邊界刻畫出來就可以了。這里給出幾種簡(jiǎn)單的偽譜計(jì)算 方法。一、隨機(jī)擾動(dòng)法A為了求矩陣的偽譜。首先產(chǎn)生一列隨機(jī)的復(fù)矩陣，(A),N,并且這些矩陣滿足條件|E|,然后求出這個(gè)矩陣E,(i,1,2,? ,n)ii的特征值，并在復(fù)平面上畫出，所得到的點(diǎn)集能近似表示(A, E)i,這種方法最簡(jiǎn)單。，(A),例2:考慮例1中的矩陣，確定其偽譜。值問題本身就比較

13、復(fù)雜不容易計(jì)算，所以算法不是很通用，不過卻可以很好地加深對(duì)偽譜概念的理解。、格點(diǎn)SVD?去當(dāng)矩陣范數(shù)取2-范數(shù)時(shí)，此時(shí)我們采用定義（4）來定義偽譜，,偽譜邊界就是最小奇異值為的格點(diǎn)連成的曲線。計(jì)算并找出該區(qū)域基本步驟:E?用MATLAB!機(jī)產(chǎn)生一個(gè)稠密矩陣，然后令,即可實(shí)現(xiàn)。EE?|E|,3擾動(dòng)為，擾動(dòng)矩陣個(gè)數(shù)為50, 100, 150, 200時(shí)的偽譜圖。,10雖然擾動(dòng)矩陣個(gè)數(shù)為200的圖比較精確，代價(jià)是需要計(jì)算很多特征值，而特征內(nèi)奇異值為的格點(diǎn)，將他們連接起來以確定偽譜的邊界。運(yùn)行結(jié)果（程序2,見附錄）:忱* F作用】的精*國(guó)M L1)將感興趣的區(qū)域劃分為份，x(k) , iy(j)為復(fù)平

14、面上格點(diǎn)的坐標(biāo)，其中x(k) , y(j)分別表示對(duì)應(yīng)格點(diǎn)(k,j)的橫坐標(biāo)值和縱坐標(biāo)值(k,1:M,j,1:L),為虛數(shù) 單位；i(zI,A),(x(k), iy(j)I,A2)計(jì)算的最小奇異值,(x(k) , iy(j)I,A; min3)通過繪制偽譜邊界上的點(diǎn)(即，最小奇異值為的那些點(diǎn))，就可以得到偽譜邊界曲線圖。xA例3:考慮例1中的矩陣偽譜，讓,10運(yùn)行結(jié)果(程序3,見附錄)：(3-a) (3-b)圖(3-a)畫的是時(shí)偽譜的邊界圖，從圖中可以看出隨著的，x,4,3,2x,3增大，偽譜所在區(qū)域范圍擴(kuò)大(驗(yàn)證了性質(zhì)2);圖(3-b)畫的是時(shí)偽譜所在的區(qū)域，與隨機(jī)擾動(dòng)法得到的結(jié)果相同。與隨

15、機(jī)擾動(dòng)法相比較，SVDt不需要計(jì)算特征值，只是計(jì)算了奇異值，但是后者的計(jì)算要比前者容易的多，也就是說該方法有更好的可行性。雖然隨機(jī)擾動(dòng)法和格點(diǎn)SVDt能計(jì)算偽譜的范圍，但其計(jì)算量往往很大。目前還有一些比較好的計(jì)算偽譜的方法，如：區(qū)域排除法、格點(diǎn)移動(dòng)法、投影算法、QR分解法等等。偽譜的應(yīng)用例4:求解下列算子的偽譜：124,Au(x),2u , (cx,dx)u,c,3 , 3i,d,u(,L),0,L,10. 16首先我們用一個(gè)切比雪夫點(diǎn)列來離散有限區(qū)間并注意邊界,L,Ln, 1條件。考慮有個(gè)離散點(diǎn)列的情況，讓 u,(,L),0j,x,Lcos(),j,0,1,? ,n , 1。點(diǎn)x處的二階導(dǎo)數(shù)

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常逼近原系統(tǒng)。這個(gè)矩陣具有眾多的特征值，隨著矩陣A的規(guī)模變n, n大，計(jì)算消耗的時(shí)間越來越長(zhǎng)。n,200這里取，先通過FOVT法大致估計(jì)一下偽譜的范圍。運(yùn)行結(jié)果：V1 nd Liv艮di IIL-skf oj Im俄區(qū)域的左世界是"3仃港噌:其由江世群號(hào)i= 3m *M國(guó)艷的匯寄W產(chǎn)L53考區(qū)域的上透界導(dǎo)h 12S129. rsr|FW方法薦剝的偽譜區(qū)城1 60 50051 5由FOVT法得到了偽譜可能存在的大致范圍為,這為我們用SVDt確定偽譜白范圍提供,747192,32 , ,127531,128128了需要?jiǎng)澐值膮^(qū)域。現(xiàn)在，在我們感興趣的區(qū)域內(nèi)計(jì)算偽譜。圖,200,30

18、, ,200,3023, 2346, 46(4-a)、(4-b)給出了網(wǎng)格劃分分別為與時(shí)所得矩陣 v, vA的偽譜的邊界曲線。運(yùn)行每一幅圖的程序所花的時(shí)間不等。,200 , 200,1,2,6,10,10,? ,10圖中，依次從外到內(nèi)，并且“ *”表示原系統(tǒng)矩陣的特征值RPIMiH i.rw(4-a) (4-b)從這個(gè)例子可以看出用FOVT法確定的偽譜界定區(qū)域遠(yuǎn)遠(yuǎn)大于偽譜的真實(shí)區(qū)域。參考文獻(xiàn) :1、矩陣擬譜計(jì)算的若干研究，加幫平， 2008年 6月 ;2、偽譜的邊界曲線及其跟蹤算法的步長(zhǎng)控制，劉穎，白峰杉，高等學(xué)校計(jì)算數(shù)學(xué)學(xué)報(bào);3、基于區(qū)域排除法和方格移動(dòng)法的矩陣偽譜計(jì)算，周劍，蔣耀林，數(shù)值計(jì)

19、算與計(jì)算機(jī)應(yīng)用 ;4、大規(guī)模矩陣偽譜計(jì)算的數(shù)值方法，孟青云， 2010年 2月 ;附錄 :程序 1:A=;n=32;%A矩陣的階數(shù)；%定義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;% 區(qū)域的下邊界;y2=d+e;%區(qū)域的上邊界；x1=a-e;% 區(qū)域的左邊界;x2=b+e;%區(qū)域的右邊界；disp('該區(qū)域的左邊界是x= 'num2str(x1);disp('該區(qū)域的右邊界是x= 'nu

21、m2str(x2);disp('該區(qū)域的下邊界是y= 'num2str(y1);disp('該區(qū)域的上邊界是y= 'num2str(y2);w=x2-x1;%該區(qū)域白寬度；h=y2-y1;%該區(qū)域白高度；rectangle('Position',x1,y1,w,h,'LineWidth',2,'LineStyle','-,'EdgeColor','m')xlim(-0.5,2.5)ylim(-3,3)title('FOV 方法得到的偽譜區(qū)域')程序 2:A=;X=;Y=;n=32;%定義A的階數(shù)；%定義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義隨機(jī)矩陣的個(gè)數(shù)；for k=1:1:mE1=rand(n);% 產(chǎn)生隨機(jī)矩陣;DE=det(E1);%求隨機(jī)矩陣的行列式值;e=0.001;%定義 epsil