资源描述：
《北航数值分析1-Jacobi法计算矩阵特征值.doc》由会员上传分享，免费在线阅读，更多相关内容在应用文档-天天文库。
1、数值分析2015/11/10准备工作Ø算法设计矩阵特征值的求法有幂法、Jacobi法、QR法等，其中幂法可求得矩阵按模最大的特征值（反幂法可求得按模最小特征值），Jacobi法则可以求得对称阵的所有特征值。分析一：由题目中所给条件λ1≤λ2≤…≤λn，可得出λ1、λn按模并不一定严格小于或大于其他特征值，且即使按模严格小于或大于其他特征值，也极有可能出现
2、λs
3、<λ1
4、<
5、λn
6、或
7、λs
8、<λn
9、<
10、λ1
11、的情况，导致按幂法和反幂法无法求解λ1或λn二者中的一者；分析二：题目要求求解与数μk=λ1+k(λn-λ1)/40最接
12、近的特征值λik(k=1,2,3…39)，这个问题其实可以转换为求A-μk按模最小的特征值的问题，但因为在第一个问题中无法确定能肯定的求得λ1和λn，所以第二个问题暂先搁浅；分析三：cond(A)2=
13、
14、A
15、
16、*
17、
18、A-1
19、
20、=
21、λ
22、max*
23、λ
24、min,这可以用幂法和反幂法求得，det(A)=λ1*λ2*…*λn，这需要求得矩阵A的所有特征值。由以上分析可知，用幂法和反幂法无法完成所有问题的求解，而用Jacobi法求得矩阵所有特征值后可以求解题目中所给的各个问题。所以该题可以用Jacobi法求解。Ø模块设计由Jacobi法
25、的经典4步骤，设计以下模块：矩阵初始化模块voidinit(double*m)迭代模块查找非对角线最大模voidfind_pq(double*m)计算旋转角度voidcalCosSin(double*m)更新矩阵voidrefreshMextrix(double*m)计算最终结果模块voidcalFinal(double*m)Ø数据结构设计由于矩阵是对称阵，上下带宽均为2，所以可以考虑用二维数组压缩存储矩阵上半带或下半带。但由于Jacobi法在迭代过程中会破坏矩阵的形态，所以原来为零的元素可能会变为非零，这就导致原来的二维数
26、组无法存储迭代后的矩阵。基于此的考虑，决定采用一维数组存储整个下三角阵，以此保证迭代的正确进行。完整代码如下（编译环境windows10+visualstudio2010）：11数值分析2015/11/10完整代码//math.cpp:定义控制台应用程序的入口点。//#include"stdafx.h"#include#include#include#defineN501#defineV(N+1)*N/2+1#definee2.718281828459045235360287
27、471352662497757247093699959574966967627724076630353#definea(i)(1.64-0.024*(i))*sin(0.2*(i))-0.64*pow(e,0.1/(i))#defineb0.16#definec-0.064#defineepspow((double)10.0,-12)#definePFbits"%10.5f"#definePFrols5#definePFe%.11e#defineFK39intp;intq;doublecosz;doublesinz;doub
28、leMAX;intkk;//#definePTSpts11数值分析2015/11/10#ifdefPTSvoidPTS(double*m){printf("-----------------------------------------------------------------------");printf("迭代第%d次",kk);for(inti=1;i<=PFrols;i++){for(intj=(i-1)*i/2+1;j<=(i+1)*i/2;j++){printf(PFbits,m[j]);}put
29、char(10);}for(inti=1;i<=PFrols+1;i++){printf("...");}putchar(10);printf("..");printf("..");printf("..");for(inti=1;i<=PFrols+2;i++){printf("...");}putchar(10);}#elsevoidPTS(double*m){}#endifvoidrecounti(inti,int*pp,int*qq){for(intj=0;j<=N-1;j++){if((i-(j+1)*j
30、/2)<=j+1){*pp=j+1;11数值分析2015/11/10*qq=i-(j+1)*j/2;break;}}}voidrefreshMetrix(double*m){intipr,ipc,iqr,iqc;m[(p+1)*p/2]=m[(p+1)*p/2]*pow(cosz,2)+m
