雅克比法求矩陣特征值特征向量

1、C語言課程設(shè)計報告課程名稱： 計算機綜合課程設(shè)計 學(xué) 院： 土木工程學(xué)院 設(shè)計題目： 矩陣特征值分解 級 別： B 學(xué)生姓名： 學(xué) 號： 同組學(xué)生： 無 學(xué) 號： 無 指導(dǎo)教師： 2012年 9 月 5 日C語言課程設(shè)計任務(wù)書（以下要求需寫入設(shè)計報告書）學(xué)生選題說明: 以所發(fā)課程設(shè)計要求為準，請同學(xué)們仔細閱讀； 本任務(wù)書提供的設(shè)計案例僅供選題參考；也可自選，但難易程度需難度相當； 鼓勵結(jié)合本專業(yè)（土木工程、力學(xué)）知識進行選題，編制程序解決專業(yè)實際問題。 限2人選的題目可由1-2人完成（A級）；限1人選的題目只能由1人單獨完成（B級）； 設(shè)計總體要求: 采用模塊化程序設(shè)計； 鼓勵可視化編程； 源

2、程序中應(yīng)有足夠的注釋； 學(xué)生可自行增加新功能模塊（視情況可另外加分）； 必須上機調(diào)試通過； 注重算法運用，優(yōu)化存儲效率與運算效率； 需提交源程序(含有注釋)及相關(guān)文件(數(shù)據(jù)或數(shù)據(jù)庫文件)； (cpp文件、txt或dat文件等) 提交設(shè)計報告書，具體要求見以下說明。設(shè)計報告格式：目錄1 課程設(shè)計任務(wù)書（功能簡介、課程設(shè)計要求）；2 系統(tǒng)設(shè)計（包括總體結(jié)構(gòu)、模塊、功能等，輔以程序設(shè)計組成框圖、流程圖解釋）；3 模塊設(shè)計（主要模塊功能、源代碼、注釋（如函數(shù)功能、入口及出口參數(shù)說明，函數(shù)調(diào)用關(guān)系描述等）；4 調(diào)試及測試：（調(diào)試方法，測試結(jié)果的分析與討論，截屏、正確性分析）；5 設(shè)計總結(jié)：（編程中遇到的

3、問題及解決方法）；6 心得體會及致謝；參考文獻 1.課程設(shè)計任務(wù)書功能簡介：a） 輸入一個對稱正方矩陣A，從文本文件讀入；b） 對矩陣A進行特征值分解，將分解結(jié)果：即U矩陣、S矩陣輸出至文本文件；c） 將最小特征值及對應(yīng)的特征向量輸出至文本文件；d） 驗證其分解結(jié)果是否正確。提示：A=USUT,具體算法可參考相關(guān)文獻。功能說明：矩陣特征值分解被廣泛運用于土木工程問題的數(shù)值計算中，如可用于計算結(jié)構(gòu)自振頻率與自振周期、結(jié)構(gòu)特征屈曲問題等。注：以三階對稱矩陣為例總體結(jié)構(gòu) 2.系統(tǒng)設(shè)計遞推求矩陣U遞推求矩陣S找矩陣A中非對角元素中的最大值，并記下其位置從文本文件中讀入數(shù)組A向屏幕和txt文件輸入矩陣S

4、求最小特征值及其對應(yīng)特征向量，并輸出到屏幕和txt文件中向屏幕和txt文件輸入矩陣Utezheng函數(shù)main函數(shù)3. 模塊設(shè)計#include#include#includeint main()FILE *fp;int tezheng(double *a,int n,double *s,double *u,double eps,int itmax); /函數(shù)調(diào)用聲明int i,j,p,itmax=1000; /itmax為最大循環(huán)次數(shù)double eps=1e-7,s33,u33; /eps為元素精度，s為對角矩陣S，u為矩陣U double a9;/a為待分解矩陣Ai=tezheng(a,

5、3,s,u,eps,1000);if(i0) /i對應(yīng)函數(shù)中的返回值itif(fp=fopen(juzhen.txt,w)=NULL) /打開待輸入txt文件printf(無法打開文件.n);return;printf(U矩陣為:n); /下幾句分別向屏幕和txt文件輸入矩陣Ufprintf(fp,U矩陣為:n);for(i=0;i3;i+)for(j=0;j3;j+)printf(%10.6f,uij);fprintf(fp,%10.6f,uij);printf(n);fprintf(fp,n);printf(S對角矩陣為:n); /下幾句分別向屏幕和txt文件輸入矩陣Sfprintf(fp

6、,S對角矩陣為:n);for(i=0;i3;i+)for(j=0;j3;j+)printf(%10.6f,sij);fprintf(fp,%10.6f,sij); printf(n);fprintf(fp,n); p=0;for(i=0;i3;i+)/下面幾句為求最小特征值及其對應(yīng)特征向量，并輸出到屏幕和txt文件中if(siis00)p=i;printf(最小特征值為:%10.6fn,spp);fprintf(fp,最小特征值為:%10.6fn,spp);printf(對應(yīng)特征向量為:n);fprintf(fp,對應(yīng)特征向量為:n);for(i=0;i3;i+)printf(%10.6fn,

7、uip);fprintf(fp,%10.6fn,uip);int tezheng(double *a,int n,double s33,double u33,double eps,int itmax)FILE *A;char line1000; /存放從文件juzhen A.txt讀入的數(shù)據(jù)char *k= ;int i,j,p,q,it;double sint,cost,sin2t,cos2t,d,tmp,t1,t2,t3,fm;it=0;if(A=fopen(juzhen A.txt,r)=NULL)printf(無法打開文件n);return; fgets(line,1000,A); /

8、從文件指針A指向的文件，即juzhen A.txt中讀入字符數(shù)據(jù)到數(shù)組line中 strtok(line,k); /原型char *strtok(char *s, const char *delim);說明：strtok()用來將字符串分割成一個個片段。參數(shù)s指向欲分割的字符串，參數(shù)delim則為分割字符串，當strtok()在參數(shù)s的字符串中發(fā)現(xiàn)到參數(shù)delim的分割字符時則會將該字符改為0 字符。在第一次調(diào)用時，strtok()必需給予參數(shù)s字符串，往后的調(diào)用則將參數(shù)s設(shè)置成NULL。每次調(diào)用成功則返回被分割出片段的指針。 for(i=0;in*n;i+) ai = atof(strtok

9、(NULL,k); /atof函數(shù)功能: 把字符串轉(zhuǎn)換成浮點數(shù) fclose(A); /釋放A指針for(i=0;in;i+) /以下幾句初始化矩陣Uuii=1.0;for(j=0;jn;j+)if(i!=j)uij=0.0;while(ititmax)it+;d=0.0;for(i=1;in;i+) /該for循環(huán)找矩陣A中非對角元素中的最大值，并記下其位置for(j=0;jd)d=tmp;p=i;q=j;if(deps) /該if條件句用數(shù)組元素精度控制循環(huán)結(jié)束與否for(i=0;in;i+)for(j=0;jn;j+)sij=ai*n+j;return (it);sint=2*ap*n+

10、q; /以下幾句為遞推公式,求S矩陣主對角元素 sin2t=2*ap*n+q/(sqrt(4*ap*n+q*ap*n+q+(aq*n+q-ap*n+p)*(aq*n+q-ap*n+p);if(aq*n+q-ap*n+p0.0)sin2t=-sin2t;cos2t=sqrt(1.0-sin2t*sin2t);sint=sin2t/(sqrt(2.0*(1.0+cos2t);cost=sqrt(1.0-sint*sint);tmp=ap*n+p;t1=tmp*cost*cost;t2=aq*n+q*cost*cost;t3=ap*n+q*sin2t;ap*n+p=t1+aq*n+q-t2-t3;a

11、q*n+q=tmp-t1+t2+t3;ap*n+q=0.0; /置該非對角元素為0.0，下次循環(huán)最大便不是它了aq*n+p=0.0; /同上for(j=0;jn;j+) /以下兩for語句求S矩陣非主對角元素if(j!=p)&(j!=q)tmp=ap*n+j;ap*n+j=tmp*cost-aq*n+j*sint;aq*n+j=tmp*sint+aq*n+j*cost;for(i=0;in;i+)if(i!=p)&(i!=q)ai*n+p=ap*n+i;ai*n+q=aq*n+i;for(i=0;in;i+) /該for語句遞推求矩陣U fm=uip;uip=fm*cost-uiq*sint;

12、uiq=fm*sint+uiq*cost;return (0);4. 調(diào)試及測試 1.0 2.0 3.0矩陣A= 2.0 2.0 5.0 3.0 5.0 1.01.屏幕輸出如下：2. 文本文件輸出見文件“juzhen.txt”。3. 結(jié)果正確性分析利用數(shù)學(xué)軟件Mathematica 8.0計算特征值（即矩陣S主對角線元素）及對應(yīng)特征向量組U如下截圖。（Eigenvalues 為計算矩陣特征值函數(shù)，Eigenvectors 為計算矩陣特征向量函數(shù)）5. 設(shè)計總結(jié)拿到這個題目，我首先想到用解方程組的方法來求解矩陣U和S,但是后來發(fā)現(xiàn)：假如這是一個三階矩陣，那么通過A=USUT可以得到含9個方程的方

13、程組，而矩陣U和S中共12個未知數(shù)，故通過這種方法建模被否定了。后查閱資料，理解了這是一個矩陣特征分解的問題，故轉(zhuǎn)化為兩部分求解：特征值和特征向量?？梢酝ㄟ^雅克比變換來求解。從txt文本文件中讀入數(shù)據(jù)又是一大障礙，本以為存在文本文件中的浮點數(shù)據(jù)只需通過浮點格式字符串就能將其讀出，結(jié)果運行后屏幕顯示一團亂。后查閱資料得知，文本文件中的數(shù)據(jù)都是字符數(shù)據(jù)。故先是將其中的字符數(shù)據(jù)讀出，然后用字符與浮點之間的轉(zhuǎn)換函數(shù)還原浮點數(shù)據(jù)。最后碰到的一個問題是主函數(shù)與被調(diào)用函數(shù)之間的銜接問題，就是一個指針所指的實參和形參，要么都是一維數(shù)組，要么都是二維數(shù)組，否則就會出錯。6. 心得體會及致謝起先聽說我們這是編程課，感覺很激動，心想上學(xué)期所學(xué)內(nèi)容終于有用武之地。但拿到題目后，我像立刻被潑了一盆冷水似的，感覺這個題目無從下手。咬著牙掙扎了兩天仍然是毫無頭緒，于是乎我想到了放棄，找到老師表達想換題目的想法。起初對老師的拒絕感到很不爽，可老師建議我去查閱資料，先把矩陣論中沒學(xué)過的只是給搞懂然后再去編程。我先查閱了一些關(guān)于算法的書籍，結(jié)果其中的程序看不懂，故又借了關(guān)于矩陣的書一起研究。經(jīng)過一個多星期的努力，