第8章 矩陣特征值問題的數(shù)值解法

第8章矩陣特征值問題的數(shù)值解法8.1引言8.2冪法及反冪法8.3雅可比方法8.4QR算法8.1引言

工程技術(shù)中有多種振動(dòng)問題，如橋梁或建筑物的振動(dòng)，機(jī)械零件、飛機(jī)機(jī)翼的振動(dòng)，及一些穩(wěn)定性分析和相關(guān)分析在數(shù)學(xué)上都可轉(zhuǎn)化為求矩陣特征值與特征向量的問題.

下面先復(fù)習(xí)一些矩陣的特征值和特征向量的基礎(chǔ)知識(shí).

定義1⑴已知n階矩陣A=(aij)，則稱為A的特征多項(xiàng)式.一般有n個(gè)根(實(shí)的或復(fù)的，復(fù)根按重?cái)?shù)計(jì)算)稱為A的特征值.用λ(A)表示A的所有特征值的集合.

A的特征方程⑵設(shè)λ為A的特征值，相應(yīng)的齊次方程組

注：當(dāng)A為實(shí)矩陣時(shí)，

(λ)=0為實(shí)系數(shù)n次代數(shù)方程，其復(fù)根是共軛成對(duì)出現(xiàn).的非零解x稱為矩陣A的對(duì)應(yīng)于λ的特征向量.

定理1

設(shè)λ為A∈Rn×n的特征值,且Ax=λx(x0)，則有⑵λ-p為A-pI的特征值，即(A-pI)x=(λ-p)x

；⑴cλ為的cA特征值(c≠0為常數(shù))；

下面敘述有關(guān)特征值的一些結(jié)論：⑶λk為Ak的特征值，即Akx=λkx

；⑷設(shè)A為非奇異矩陣，那么λ≠0,且λ-1為A-1的特征值，即A-1x=λ-1x.

定理2

設(shè)λi(i=1,2,,n)為n階矩陣A=(aij)的特征值，則有⑴稱為A的跡；⑵

定理3設(shè)A∈Rn×n，則有

定義2

設(shè)n階矩陣A=(aij)，令

下面討論矩陣特征值界的估計(jì).⑴；⑵集合稱為復(fù)平面上以aii為圓心，以ri為半徑的n階矩陣A的n個(gè)Gerschgorin圓盤.

定理4(Gerschgorin圓盤定理)

特別地，如果A的一個(gè)圓盤Di是與其它圓盤分離(即孤立圓盤)，則Di中精確地包含A的一個(gè)特征值.⑴設(shè)n階矩陣A＝(aij)，則A的每一個(gè)特征值必屬于下面某個(gè)圓盤之中⑵如果A有m個(gè)圓盤組成一個(gè)連通的并集S，且S與余下n-m個(gè)圓盤是分離的，則S內(nèi)恰包含A的m個(gè)特征值.或者說A的特征值都在n個(gè)圓盤的并集中.

證明只就⑴給出證明.設(shè)λ為A的特征值，即Ax=λx，其中x=(x1,x2,,xn)T0.或

記，考慮Ax=λx的第k個(gè)方程，即于是即這說明，A的每一個(gè)特征值必位于A的一個(gè)圓盤中，并且相應(yīng)的特征值λ一定位于第k個(gè)圓盤中(其中k是對(duì)應(yīng)特征向量x絕對(duì)值最大的分量的下標(biāo)).

例2

估計(jì)矩陣A的特征值范圍，其中

解矩陣A的3個(gè)圓盤為

由定理8，可知A的3個(gè)特征值位于3個(gè)圓盤的并集中，由于D1是孤立圓盤，所以D1內(nèi)恰好包含A的一個(gè)特征值λ1(為實(shí)特征值)，即A的其它兩個(gè)特征值λ2,λ3包含在D2,D3的并集中.

定義3設(shè)A∈Rn×n為對(duì)稱矩陣，對(duì)于任一非零向量x，稱為對(duì)應(yīng)于向量x的瑞利(Rayleigh)商.

定理5設(shè)A∈Rn×n為對(duì)稱矩陣(其特征值次序記為λ1≥λ2≥≥λn)，則1.(對(duì)任何非零x∈Rn)；2.；3..

證明只證1，關(guān)于2,3自己作練習(xí).

由于A為實(shí)對(duì)稱矩陣，可將λ1,λ2,,λn

對(duì)應(yīng)的特征向量x1,x2,,xn

正交規(guī)范化，則有(xi,xj)=δij，設(shè)x0為Rn中任一向量，則有于是從而1成立.結(jié)論1說明瑞利商必位于λn和λ1之間.

關(guān)于計(jì)算矩陣A的特征值問題，當(dāng)n＝2,3時(shí)，我們還可按行列式展開的辦法求(λ)=0的根.但當(dāng)n較大時(shí)，如果按展開行列式的辦法，首先求出(λ)的系數(shù)，再求(λ)的根，工作量就非常大，用這種辦法求矩陣的特征值是不切實(shí)際的，由此需要研究求A的特征值及特征向量的數(shù)值解法.

本章將介紹一些計(jì)算機(jī)上常用的兩類方法，一類是冪法及反冪法（迭代法），另一類是正交相似變換的方法（變換法）.

冪法與反冪法都是求實(shí)矩陣的特征值和特征向量的向量迭代法，所不同的是冪法是計(jì)算矩陣的主特征值(矩陣按模最大的特征值稱為主特征值，其模就是該矩陣的譜半徑)和相應(yīng)特征向量的一種向量迭代法，而反冪法則是計(jì)算非奇異(可逆)矩陣按模最小的特征值和相應(yīng)特征向量的一種向量迭代法.下面分別介紹冪法與反冪法.8.2冪法及反冪法現(xiàn)討論求λ1及x1的方法.

設(shè)實(shí)矩陣A=(aij)有一個(gè)完全的特征向量組，即A有n個(gè)線性無關(guān)的特征向量，設(shè)矩陣A的特征值為λ1,λ2,,λn,相應(yīng)的特征向量為x1,x2,,xn.已知A的主特征值λ1是實(shí)根，且滿足條件

8.2.1冪法（又稱乘冪法）

冪法的基本思想是:任取非零的初始向量v0

,由矩陣A構(gòu)造一向量序列{vk}稱為迭代向量，由假設(shè)，v0可唯一表示為于是其中由假設(shè)故從而為λ1的特征向量.所以當(dāng)k充分大時(shí)，有即為矩陣A的對(duì)應(yīng)特征值1的一個(gè)近似特征向量.用(vk)i

表示vk的第i個(gè)分量，則當(dāng)k充分大時(shí)，有即為A的主特征值1的近似值.由于

迭代公式實(shí)質(zhì)上是由矩陣A的乘冪Ak與非零向量v0相乘來構(gòu)造向量序列{vk}={Akv0}，從而計(jì)算主特征值λ1及其對(duì)應(yīng)的特征向量，這就是冪法的思想.的收斂速度由比值來確定，r越小收斂越快，但當(dāng)r≈1時(shí)收斂可能很慢.

定理6設(shè)A∈Rn×n有n個(gè)線性無關(guān)的特征向量，主特征值λ1滿足條件|λ1|>|λ2|≥≥|λn|，則對(duì)任何非零向量v0(a10)，冪法的算式成立.又設(shè)A有n個(gè)線性無關(guān)的特征向量，λ1對(duì)應(yīng)的r個(gè)線性無關(guān)的特征向量為x1,x2,,xr，則由(2.2)式有

如果A的主特征值為實(shí)的重根,即λ1=λ2==λr,且|λr|>|λr+1|≥≥|λn|，為A的特征向量，這說明當(dāng)A的主特征值是實(shí)的重根時(shí)，定理6的結(jié)論還是正確的.

應(yīng)用冪法計(jì)算A的主特征值λ1及其對(duì)應(yīng)的特征向量時(shí)，如果|λ1|>1(或|λ1|<1)，迭代向量

vk的各個(gè)不等于零的分量將隨k→∞

而趨向于無窮(或趨向于零)，這樣在計(jì)算機(jī)實(shí)現(xiàn)時(shí)就可能“溢出”.為克服這個(gè)缺點(diǎn)，就需要將迭代向量加以規(guī)范化.

設(shè)有一向量v0，將其規(guī)范化得向量為其中max(v)表示v的絕對(duì)值最大的分量.即如果有則max(v)=vq，且q為所有絕對(duì)值最大的分量中的最小下標(biāo).在定理6的條件下冪法可這樣進(jìn)行：任取一初始向量v00(a10)，構(gòu)造規(guī)范化向量序列為由(2.3)式收斂速度由比值r=|λ2/λ1|確定.總結(jié)上述結(jié)論，有同理，可得到

定理7設(shè)A∈Rn×n有n個(gè)線性無關(guān)的特征向量，主特征值λ1滿足|λ1|>|λ2|≥≥|λn|，則對(duì)任意非零初始向量v0=u0(a10)，有冪法計(jì)算公式為則有⑴⑵

例1

用冪法計(jì)算矩陣的主特征值與其對(duì)應(yīng)的特征向量.

解取v0=u0=(0,0,1)T

,則直到k=8時(shí)的計(jì)算結(jié)果見下表k12,4,1,40.5,1,0.2524.5,9,7.7590.5,1,0.861135.7222,11.4444,8.36111.44440.5,1,0.736045.4621,10.9223,8.230610.92230.5,1,0.753655.5075,11.0142,8.257611.01420.5,1,0.749465.4987,10.9974,8.249410.99740.5,1,0.750175.5002,11.0005,8.250111.00050.5,1,0.750085.5000,11.0000,8.250011.00000.5,1,0.7500從而8.2.2冪法的加速方法1、原點(diǎn)平移法

由前面討論知道，應(yīng)用冪法計(jì)算A的主特征值的收斂速度主要由比值r=|λ2/λ1|來決定，但當(dāng)r接近于1時(shí)，收斂可能很慢.這時(shí)，一個(gè)補(bǔ)救辦法是采用加速收斂的方法.其中p為參數(shù)，設(shè)A的特征值為i，則對(duì)矩陣B的特征值為i-p

，而且A,B的特征向量相同.

引進(jìn)矩陣B=A-pI.

如果要計(jì)算A的主特征值1,只要選擇合適的數(shù)p，使1-p為矩陣B=A-pI

的主特征值，且那么，對(duì)矩陣B=A-pI應(yīng)用冪法求其主特征值1-p,收斂速度將會(huì)加快.這種通過求B=A-pI的主特征值和特征向量，而得到A的主特征值和特征向量的方法叫原點(diǎn)平移法.對(duì)于A的特征值的某種分布，它是十分有效的.

例4

設(shè)A∈R4×4有特征值比值r=|λ2/λ1|≈0.9.做變換B=A-12I(p=12),則B的特征值為應(yīng)用冪法計(jì)算B的主特征值μ1的收斂速度的比值為

雖然常常能夠選擇有利的p值,使冪法得到加速,但設(shè)計(jì)一個(gè)自動(dòng)選擇適當(dāng)參數(shù)p的過程是困難的.

定理

設(shè)A∈Rn×n為對(duì)稱矩陣，特征值滿足對(duì)應(yīng)的特征向量xi滿足(xi,xj)=δij

(單位正交向量)

，應(yīng)用冪法公式(2.9)計(jì)算A的主特征值1，則規(guī)范化向量uk的瑞利商給出1的較好的近似值為

由此可見，R(uk)比μk更快的收斂于1.2、瑞利商加速

證明由(2.8)式及得

冪法的瑞利商加速迭代公式可以寫為其中A為n階實(shí)對(duì)稱矩陣.

對(duì)給定的誤差限，當(dāng)|μ

k–μk-1|<時(shí)，取近似值反冪法

反冪法是用于求非奇異矩陣A的按模最小的特征值和對(duì)應(yīng)特征向量的方法.而結(jié)合原點(diǎn)平移法的反冪法則可以求矩陣A的任何一個(gè)具有先了解的特征值和對(duì)應(yīng)的特征向量。設(shè)矩陣A非奇異,其特征值i(i=1,2,,n),滿足其相應(yīng)的特征向量x1,x2,,xn線性無關(guān)，則A-1

的特征值為1/i

，對(duì)應(yīng)的特征向量仍為xi

(i=1,2,,n).此時(shí),A-1的特征值滿足因此,對(duì)A-1應(yīng)用冪法,可求出其主特征值1/n

μ

k

和特征向量

xn

uk.從而求得A的按模最小特征值

n

1/μk

和對(duì)應(yīng)的特征向量

xn

uk

，這種求A-1的方法就稱為反冪法.為了避免求A-1,可通過解線性方程組Avk=uk-1得到vk，采用LU分解法，即先對(duì)A進(jìn)行LU分解A=LU,此時(shí)反冪法的迭代公式為

反冪法的迭代公式為對(duì)給定的誤差，當(dāng)|μk–μk-1|<

時(shí)，得顯然，反冪法的收斂速度取決于比值，比值越小，收斂越快.

定理8設(shè)A∈Rn×n為非奇異矩陣，且有n個(gè)線性無關(guān)的特征向量，其對(duì)應(yīng)的特征值滿足

|1|≥|2|≥≥|n-2|>|n|>0，則對(duì)任意非零初始向量u0(an0)

，由反冪法計(jì)算公式構(gòu)造的向量序列{vk}，{uk}滿足⑴⑵

在反冪法中也可以用原點(diǎn)平移法加速迭代過程，或求其它特征值與其對(duì)應(yīng)的特征向量.

如果矩陣(A-pI)-1存在，顯然其特征值為對(duì)應(yīng)的特征向量仍然是x1,x2,,xn.

如果p是A的特征值j的一個(gè)近似值，且設(shè)j與其它特征值是分離的，即就是說1/(j-p)是矩陣(A-pI)-1的主特征值，可用反冪法應(yīng)用于矩陣(A-pI)-1計(jì)算j及對(duì)應(yīng)的特征向量.現(xiàn)對(duì)矩陣(A-pI)-1應(yīng)用反冪法，得到如下迭代公式

定理9設(shè)A∈Rn×n有n個(gè)線性無關(guān)的特征向量，矩陣A的特征值及對(duì)應(yīng)的特征向量分別記為i

及xi(i=1,2,,n)，而p為j的近似值，(A-pI)-1存在，且⑴⑵則對(duì)任意非零初始向量u0(aj0)

，由反冪法計(jì)算公式(2.12)構(gòu)造的向量序列{vk}，{uk}滿足且收斂速度為

由該定理知，對(duì)A-pI(其中p≈j)應(yīng)用反冪法，可用來計(jì)算特征向量xj，只要選擇p是j的一個(gè)較好的近似且特征值分離情況較好，一般r很小，常常只要迭代一二次就可完成特征向量的計(jì)算.反冪法迭代公式中的vk是通過解方程組求得的,為了節(jié)省工作量,可以先將A-pI進(jìn)行三角分解于是求vk相對(duì)于解兩個(gè)三角形方程組實(shí)驗(yàn)表明,按下述方法選擇u0是較好的:選u0使用回代求解(2.13)即得v1,然后再按公式(2.12)進(jìn)行迭代.例5

求矩陣A最接近于1.9的特征值和相應(yīng)的特征向量取作迭代，結(jié)果如表：8.3Jacobi方法Jacobi

方法的基本思路矩陣的旋轉(zhuǎn)變換Jacobi

迭代法收斂定理

說明解：

例6：

8.4QR算法8.4.1.

化矩陣為Hessenberg形8.4.2QR算法及其收斂性（1）鏡面反射矩陣（2）化一般實(shí)矩陣為上Hessenberg矩陣（3）化對(duì)稱矩陣為三對(duì)角矩陣

當(dāng)A為實(shí)矩陣，如果限制用正交相似變換，由于A有復(fù)的特征值，A不能用正交相似變換約化為上三角陣.用正交相似變換能約化到什么程度呢？(Schur定理)

設(shè)A∈Rn×n，則存在酉矩陣U使其中rii(i=1,2,,n)為A的特征值.

下面給出理論上有關(guān)通過酉相似變換及正交變換可以約化一般矩陣A到什么程度的問題.其中Rii(i=1,2,,m)為一階或二階方陣，且每個(gè)一階Rii是A的實(shí)特征值，每個(gè)二階對(duì)角塊Rii的兩個(gè)特征值是A的兩個(gè)共軛復(fù)特征值.

(實(shí)Schur分解)

設(shè)A∈Rn×n，則存在正交矩陣Q使63（1）鏡面反射矩陣8.4.1.

化矩陣為Hessenberg形鏡面反射矩陣的意義是“成批”消去向量的非零元素。例：解:#(1)用初等反射矩陣作正交相似變換約化一般實(shí)矩陣A為上Hessenberg矩陣.化矩陣為上Hessenberg矩陣討論兩個(gè)問題(2)用初等反射矩陣作正交相似變換約化對(duì)稱矩陣A為對(duì)稱三對(duì)角矩陣.

于是，求原矩陣特征值問題，就轉(zhuǎn)化為求上Hessenberg矩陣或?qū)ΨQ三對(duì)角矩陣的特征值問題.（2）化一般實(shí)矩陣為上Hessenberg矩陣

設(shè)A∈Rn×n，下面來說明，可選擇初等反射矩陣U1,U2,,Un-2使A經(jīng)正交相似變換約化為一個(gè)上Hessenberg陣.

(1)設(shè)其中c1=(a21,,an1)T∈Rn-1

,不妨設(shè)c1≠0，否則這一步不需要約化.于是,可選擇初等反射陣令使其中則其中

(2)第k步約化：重復(fù)上述過程，設(shè)對(duì)A已完成第1步,,第k-1步正交相似變換，即有或且其中為k階上Hessenberg陣，設(shè)ck≠0,于是可選擇初等反射陣Rk使其中，Rk計(jì)算公式為令則其中為k+1階上Hessenberg陣，第k步約化只需計(jì)算及(當(dāng)A為對(duì)稱矩陣時(shí)，只需要計(jì)算).

(3)重復(fù)上述過程，則有

定理11(Household約化矩陣為上Hessenberg陣)設(shè)A∈Rn×n則存在初等反射矩陣U1,U2,,Un-2

使為上Hessenberg矩陣.總結(jié)上述結(jié)論，有

例7用Household方法將矩陣約化為上Hessenberg陣.

解選取初等反射陣R1使,其中c1=(2,4)T.

(1)計(jì)算R1:則有

(2)約化計(jì)算:則得到上Hessenberg陣為（3）化對(duì)稱矩陣為對(duì)稱三對(duì)角矩陣推論(Household約化對(duì)稱矩陣為對(duì)稱三對(duì)角矩陣)設(shè)A∈Rn×n為對(duì)稱矩陣，則存在初等反射矩陣U1,U2,,Un-2使為對(duì)稱三對(duì)角矩陣.8.4.2QR算法及其收斂性

QR算法是一種變換方法，是計(jì)算一般矩陣(中小型矩陣)全部特征值問題的最有效方法之一.(1)上Hessenberg矩陣的全部特征值問題；(2)計(jì)算對(duì)稱三對(duì)角矩陣的全部特征值問題，

目前QR算法主要用來計(jì)算：且QR算法具有收斂快，算法穩(wěn)定等特點(diǎn).

對(duì)于一般矩陣A∈Rn×n(或?qū)ΨQ矩陣)，則首先用Household方法將A化為上Hessenberg陣B(或?qū)ΨQ三對(duì)角陣)，然后再用QR算法計(jì)算B的全部特征值.

設(shè)A∈Rn×n，且A對(duì)進(jìn)行QR分解，即A=QR，其中R為上三角陣,Q為正交陣,于是可得到一新矩陣B=RQ=QTAQ.顯然，B是由A經(jīng)過正交相似變換得到，因此B與A的特征值相同.再對(duì)B進(jìn)行QR分解，又可得一新矩陣,重復(fù)這一過程可得到矩陣序列：

設(shè)A=A1

將A1進(jìn)行QR分解A1=Q1R1

作矩陣A2=R1Q1=Q1TR1Q1

QR算法，就是利用矩陣的QR分解，按上述遞推法則構(gòu)造矩陣序列{Ak}的過程.只要A為非奇異矩陣，則由QR算法就完全確定{Ak}.

定理13(基本QR算法)設(shè)A=A1∈Rn×n,構(gòu)造QR算法:于是

證明(1),(2)顯然，現(xiàn)證(3).用歸納法，顯然，當(dāng)k=1時(shí)有,

設(shè)Ak-1有分解式

定理14(QR算法的收斂性)設(shè)A=(aij)∈Rn×n,(1)如果A的特征值滿足:(2)A有標(biāo)準(zhǔn)型A=XDX-1,其中D=diag(1,2,n),且設(shè)X-1有三角分解X-1=LU(L為單位下三角陣，U為上三角陣)，則由QR算法產(chǎn)生的{Ak}本質(zhì)上收斂于上三角矩陣，即若記Ak=(aij(k))，則(1)(2)當(dāng)i>j時(shí)，當(dāng)i<j時(shí)aij(k)的極限不一定存在.

推論如果對(duì)稱矩陣A滿足定理14的條件,則由QR算法產(chǎn)生的{Ak}收斂于對(duì)角陣D=diag(1,2,n).

設(shè)A∈Rn×n，且A有完備的特征向量集合，如果A的等模特征值中只有實(shí)重特征值或多重共軛復(fù)特征值，則由QR算法產(chǎn)生的{Ak}本質(zhì)收斂于分塊上三角矩陣(對(duì)角塊為一階和二階子塊)，且對(duì)角塊中每一個(gè)2×2子塊給出A的一對(duì)共軛復(fù)特征值，每一個(gè)一階對(duì)角子塊給出A的實(shí)特征值，即