方法求矩陣全部特征值

1、數(shù) 值 分 析課程設(shè)計QR方法求矩陣全部特征值問題復(fù)述用算法求矩陣特征值：（i） （ii）要求：（1） 根據(jù)算法原理編制求（i）與（ii）中矩陣全部特征值的程序并輸出計算結(jié)果（要求誤差）（2） 直接用現(xiàn)有的數(shù)學(xué)軟件求（i），（ii）的全部特征值，并與（1）的結(jié)果比較。問題分析 QR方法是目前公認(rèn)為計算矩陣全部特征值的最有效的方法。它適用于任一種實矩。QR方法的原理是利用矩陣的正交分解產(chǎn)生一個與矩陣A相似的矩陣迭代序列，這個序列將收斂于一個上三角陣或擬上三角陣，從而求得原矩陣A的全部特征值。QR是一個迭代算法，每一步迭代都要進(jìn)行QR分解，再作逆序的矩陣乘法。要是對矩陣A直接用QR方法，計算量太大

2、，效率不高。因此，一個實際的QR方法計算過程包括兩個步驟，首先要對矩陣A用初等相似變換約化為上Hessenberg矩陣，再用QR方法求上Hessenberg矩陣的全部特征值。示意如下： 對B矩陣的約化只需將每列次對角線上的元素約化為0。因此常常用平面旋轉(zhuǎn)陣（Givens變換陣）來進(jìn)行約化。1、 QR方法原理及收斂性設(shè)已實現(xiàn)了QR分解，記其中是正交陣，是上三角陣。因為，用對作正交相似變換有可改寫為 顯然只是的QR分解因子陣的逆序相乘，而且與原矩陣有相同的特征值。對矩陣再重復(fù)以上過程并繼續(xù)下去，可以得到一個與原矩陣A有相同特征值的矩陣序列。其過程如下： 記 對作正交分解 作矩陣 ， 對作正交分解

3、作矩陣 ，重復(fù)以上過程可得一般的形式為對作正交分解 構(gòu)成矩陣序列 （k=1,2）A 從矩陣A開始得到一個矩陣序列這個矩陣序列中每一個矩陣都與原矩陣相似，即都有與A相同的特征值。這個矩陣序列在實質(zhì)上收斂于依次以為對角元的上三角陣。具體可以表示為其中 的極限不一定存在。2、 用正交相似變換約化矩陣為上Hessenberg陣 用Householder變換可以將一個向量指定的某個分量以下的各分量變?yōu)?。我們只要求消掉A的次對角線以下的元素，即將A約化為上Hessenberg陣。為了使變換前后矩陣的特征值不變，需要用Householder矩陣對A作相似變換，即用正交陣同時左乘和右乘A時，原來已變?yōu)?的元

4、素不再改變。若設(shè)是Householder矩陣，用它對A的第一列元素的變換示意如下： 依次對A的各列進(jìn)行類似的變換，一共要進(jìn)行次變換，最終可以得到一個與原矩陣A有相同特征值的上Hessenberg陣。以上約化A為上Hessenberg陣的過程可以用一系列Householder矩陣來實現(xiàn)。其中，對于每一個有經(jīng)過步約化就可得到一個上Hessenberg陣（A的第列不需要約化）3、 Hessenberg陣的QR算法 設(shè)矩陣，其特征值都是實數(shù)。若已用Householder變換約化為上Hessenberg陣對已得到的上Hessenberg陣可用QR變換，經(jīng)過迭代過程約化為上三角形矩陣以求出A的特征值。只要

5、A的特征值是實數(shù)，將B約化為上三角形矩陣總是可以做到的。 由于B矩陣結(jié)構(gòu)上的特點，對B矩陣的約化只需將每列次對角線上的元素約化為0.可用平面旋轉(zhuǎn)陣（Givens變換陣）來進(jìn)行約化。n階方陣為平面旋轉(zhuǎn)陣。還是一個非對稱的正交陣，有，也是一個平面旋轉(zhuǎn)陣。有以下幾種作用： 1、左乘向量只改變x 的第i個和第j個分量。現(xiàn)構(gòu)造對x作變換使的結(jié)果將x的第j個分量約化為0。令，有調(diào)整角可使。若記，按下式選取于是有。 2、對非零的n維向量x連續(xù)左乘，可將x的第i+1到第n個分量都約化為零；即其中 3、用對矩陣A作變換得到的結(jié)論是左乘A只改變A的第i,j 行；右乘A只改變A的第i,j列； 用對A作正交相似變換改

6、變了A的i行和j行以及i列和j列。用一系列連續(xù)左乘矩陣A，可以將矩陣A化為上三角陣。 數(shù)據(jù)結(jié)構(gòu)描述主要數(shù)據(jù)成員說明double Ann 存放矩陣Adouble Qnn,存放QR分解式的正交矩陣Qdouble Rnn存放QR分解式的上三角陣Rdouble pnnGivens矩陣pdouble InnN階單位陣double Vnn存放Q矩陣的轉(zhuǎn)置double Tnn初等反射陣Tdouble eps精度double max最大值double det存放行列式的值int count存放迭代次數(shù)主要函數(shù)成員說明double Det(double Lnn)用高斯列主元方法求行列式int Non_singu

7、larMatrix(double Lnn)判斷是否是非奇異矩陣void Disp(double Hnn)輸出矩陣int IsZero(double a,int j)判斷數(shù)組是否全為0int sgn(double y)符號函數(shù)void Hessenberg(double Ann)將矩陣化為上Hessenberg矩陣int IsHessenberg(double Enn)判斷是否是上Hessenberg矩陣void QRAlgorithm(double Ann)QR算法求特征值void SeekEigenvalue(double Ann)判斷是否滿足QR算法條件，滿足則進(jìn)行QR方法求特征值算法的描

8、述（流程圖）  源程序C程序為：#include<stdio.h>#include<math.h>#define n 3#define eps 1E-5double Det(double Lnn)/用高斯列主元方法求行列式double det=1,t;for(int k=0;k<n-1;k+)double max=Lkk;int Ik=k;for(int j=k+1;j<n;j+)if(max<Ljk)max=Ljk;Ik=j;if(max=0) return 0;if(Ik!=k)for(int j=k;j<n;j+)t=LIkj;L

9、Ikj=Lkj;Lkj=t;det*=-1;for(int i=k+1;i<n;i+)t=Lik/Lkk;Lik=t;for(int j=k+1;j<n;j+)Lij=Lij-t*Lkj;det*=Lkk;if(Ln-1n-1=0) return 0;elseprintf(" 矩陣A的行列式為:%.5fn",det*Ln-1n-1);return det*Ln-1n-1;int Non_singularMatrix(double Lnn)/判斷是否是非奇異矩陣printf(" 判斷矩陣A的行列式是否為0?n");if(Det(L)!=0) r

10、eturn 1;return 0;void Disp(double Hnn)/輸出矩陣for(int i=0;i<n;i+)for(int j=0;j<n;j+)printf("%.6f ",Hij);printf("n");int IsZero(double a,int j)/判斷數(shù)組是否全為0for(int i=0;i<j;i+)if(ai!=0) return 0;return 1;int sgn(double y)/符號函數(shù)if(y>0) return 1;else if(y=0) return 0;else return

11、 -1;void Hessenberg(double Ann)/將矩陣化為上Hessenberg矩陣printf("原矩陣為:n");Disp(A);/輸出原矩陣double Tnn,Bnn,Cnn,cn-1,vn-1,un-1,Rn-1n-1,In-1n-1,t,w,s;int i,j,k,m;for(k=0;k<n-2;k+)for(i=0;i<n-k-1;i+)for(j=0;j<n-k-1;j+)if(i=j) Iij=1;else Iij=0;/定義單位陣Idouble max=fabs(Ak+1k);/求最大值for(i=0;i<n-k-

12、1;i+)if(max<fabs(Ai+k+1k)max=fabs(Ai+k+1k);for(i=0;i<n-k-1;i+)ci=Ai+k+1k/max;/標(biāo)準(zhǔn)化數(shù)組if(IsZero(c,n-k-1)/判斷數(shù)組是否全為0continue;/數(shù)組為0，則這一步不需要約化for(i=0,t=0.0;i<n-k-1;i+)t+=ci*ci;vk=sgn(Ak+1k)*sqrt(t);/u0=c0+vk;for(j=1;j<n-k-1;j+)uj=cj;/求矩陣uw=vk*(c0+vk);for(i=0;i<n-k-1;i+)for(j=0;j<n-k-1;j+)

13、Rij=Iij-ui*uj/w;for(i=0;i<n;i+)for(j=0;j<n;j+)if(i=j) Tij=1;else Tij=0;/定義矩陣T為單位陣for(i=0;i<n-k-1;i+)for(j=0;j<n-k-1;j+)Ti+k+1j+k+1=Rij;/初等反射陣Tfor(i=0;i<n;i+) /矩陣T左乘矩陣Afor(j=0;j<n;j+)for(m=0,s=0.0;m<n;m+)s+=Tim*Amj;Bij=s;for(i=0;i<n;i+) /矩陣T右乘矩陣Afor(j=0;j<n;j+)for(m=0,s=0.0

14、;m<n;m+)s+=Bim*Tmj;Cij=s;for(i=0;i<n;i+)for(j=0;j<n;j+)Aij=Cij;printf("原矩陣化成上Hessenberg矩陣后為:n");Disp(A);int IsHessenberg(double Enn)/判斷是否是上Hessenberg矩陣for(int i=2;i<n;i+)for(int j=0;j+1<i;j+)if(Eij!=0)return 0;return 1;int count=1;void QRAlgorithm(double Ann)/QR算法求特征值int i,j

15、,k,m;double pnn,Qnn,Rnn,Fnn,Vnn,c,s,t,y,max;for(i=0;i<n;i+)/賦Q為單位陣for(int j=0;j<n;j+)if(i!=j) Qij=0;else Qij=1;for(m=0;m<n-1;m+)/對矩陣A進(jìn)行QR分解for(i=0;i<n;i+)/賦p為單位陣for(j=0;j<n;j+)if(i!=j) pij=0;else pij=1;c=Amm/sqrt(Amm*Amm+Am+1m*Am+1m);s=Am+1m/sqrt(Amm*Amm+Am+1m*Am+1m);pmm=c;pmm+1=s;pm+

16、1m=-s;pm+1m+1=c;/p為Givens矩陣for(i=0;i<n;i+)/p矩陣左乘W矩陣,p矩陣左乘Q矩陣for(j=0;j<n;j+)t=0;y=0;for(k=0;k<n;k+)t+=pik*Akj;y+=pik*Qkj;Rij=t;Fij=y;for(i=0;i<n;i+)/A,R賦新值for(j=0;j<n;j+)Aij=Rij;Qij=Fij;for(i=0;i<n;i+)/Q轉(zhuǎn)置for(j=0;j<n;j+)Vij=Qji;for(i=0;i<n;i+)for(j=0;j<n;j+)Qij=Vij;for(i=0;

17、i<n;i+)/矩陣A為矩陣R和矩陣Q的乘積for(j=0;j<n;j+)t=0;for(k=0;k<n;k+)t+=Rik*Qkj;Aij=t;count+;max=fabs(A10);/求A矩陣下三角的絕對值最大值for(i=1;i<n;i+)for(j=0;j<i;j+)if(fabs(Aij)>max)max=fabs(Aij);if(max<eps)/滿足精度條件，輸出Q矩陣，R矩陣以及特征值printf("在精度%.1e下,QR算法迭代次數(shù)為:%d次n輸出矩陣A%d:n",eps,count-1,count);Disp(

18、A);/輸出A矩陣printf("輸出矩陣Q:n");Disp(Q);printf("輸出矩陣R:n");Disp(R);printf("輸出A矩陣的全部特征值:");for(i=0;i<n;i+)if(i%3=0) printf("n");printf("a%d=%.6ft",i+1,Aii);printf("n");return;QRAlgorithm(A);void SeekEigenvalue(double Ann)/判斷是否滿足QR算法條件，滿足則進(jìn)行QR方法求特征值double Znn;for(int i=0;i<n;i+)for(int j=0;j<n;j+)Zij=Aij;printf("判斷矩陣A是否是非奇異矩陣?n");if(Non_singularMatrix(Z)=0)printf("矩陣A不是非奇異矩陣!不符合分解條件!n");return;elseprintf("矩陣A是非奇異矩陣!符合分解條件!n");printf("判斷矩陣A是否是上Hessenberg矩陣?