北航数值分析1Jacobi法计算矩阵特征值.docx
《北航数值分析1Jacobi法计算矩阵特征值.docx》由会员分享，可在线阅读，更多相关《北航数值分析1Jacobi法计算矩阵特征值.docx（15页珍藏版）》请在冰点文库上搜索。
北航数值分析1Jacobi法计算矩阵特征值
准备工作
Ø算法设计
矩阵特征值的求法有幂法、Jacobi法、QR法等，其中幂法可求得矩阵按模最大的特征值（反幂法可求得按模最小特征值），Jacobi法则可以求得对称阵的所有特征值。
分析一：
由题目中所给条件λ1≤λ2≤…≤λn，可得出λ1、λn按模并不一定严格小于或大于其他特征值，且即使按模严格小于或大于其他特征值，也极有可能出现|λs|<λ1|<|λn|或|λs|<λn|<|λ1|的情况，导致按幂法和反幂法无法求解λ1或λn二者中的一者；
分析二：
题目要求求解与数μk=λ1+k（λn-λ1）/40最接近的特征值λik（k=1,2,3…39），这个问题其实可以转换为求A-μk按模最小的特征值的问题，但因为在第一个问题中无法确定能肯定的求得λ1和λn，所以第二个问题暂先搁浅；
分析三：
cond（A）2=||A||*||A-1||=|λ|max*|λ|min,这可以用幂法和反幂法求得，det（A）=λ1*λ2*…*λn，这需要求得矩阵A的所有特征值。
由以上分析可知，用幂法和反幂法无法完成所有问题的求解，而用Jacobi法求得矩阵所有特征值后可以求解题目中所给的各个问题。
所以该题可以用Jacobi法求解。
Ø模块设计
由Jacobi法的经典4步骤，设计以下模块：
矩阵初始化模块
voidinit（double*m）
迭代模块
查找非对角线最大模
voidfind_pq（double*m）
计算旋转角度
voidcalCosSin（double*m）
更新矩阵
voidrefreshMextrix（double*m）
计算最终结果模块
voidcalFinal（double*m）
Ø数据结构设计
由于矩阵是对称阵，上下带宽均为2，所以可以考虑用二维数组压缩存储矩阵上半带或下半带。
但由于Jacobi法在迭代过程中会破坏矩阵的形态，所以原来为零的元素可能会变为非零，这就导致原来的二维数组无法存储迭代后的矩阵。
基于此的考虑，决定采用一维数组存储整个下三角阵，以此保证迭代的正确进行。
完整代码如下（编译环境windows10+visualstudio2010）：
完整代码
//math.cpp:
定义控制台应用程序的入口点。
//
#include"stdafx.h"
#include
#include
#include
#defineN501
#defineV（N+1）*N/2+1
#definee2.6630353
#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
#defineFK39
intp;
intq;
doublecosz;
doublesinz;
doubleMAX;
intkk;
//#definePTSpts
#ifdefPTS
voidPTS（double*m）
{
printf（"-----------------------------------------------------------------------\n"）;
printf（"迭代第%d次\n",kk）;
for（inti=1;i<=PFrols;i++）
{
for（intj=（i-1）*i/2+1;j<=（i+1）*i/2;j++）
{
printf（PFbits,m[j]）;
}
putchar（10）;
}
for（inti=1;i<=PFrols+1;i++）
{
printf（"..."）;
}
putchar（10）;
printf（"..\n"）;
printf（"..\n"）;
printf（"..\n"）;
for（inti=1;i<=PFrols+2;i++）
{
printf（"..."）;
}
putchar（10）;
}
#else
voidPTS（double*m）
{
}
#endif
voidrecounti（inti,int*pp,int*qq）
{
for（intj=0;j<=N-1;j++）
{
if（（i-（j+1）*j/2）<=j+1）
{
*pp=j+1;
*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[（q+1）*q/2]*pow（sinz,2）+2*m[（p-1）*p/2+q]*cosz*sinz;
m[（q+1）*q/2]=m[（p+1）*p/2]*pow（sinz,2）+m[（q+1）*q/2]*pow（cosz,2）-2*m[（p-1）*p/2+q]*cosz*sinz;
for（inti=1;i<=N;i++）
{
if（i!
=p&&i!
=q）
{
if（i>p）
{
ipr=i;
ipc=p;
}
else
{
ipr=p;
ipc=i;
}
if（i>q）
{
iqr=i;
iqc=q;
}
else
{
iqr=q;
iqc=i;
}
m[（ipr-1）*ipr/2+ipc]=m[（ipr-1）*ipr/2+ipc]*cosz+m[（iqr-1）*iqr/2+iqc]*sinz;
m[（iqr-1）*iqr/2+iqc]=-m[（ipr-1）*ipr/2+ipc]*sinz+m[（iqr-1）*iqr/2+iqc]*cosz;
}
}
m[（p-1）*p/2+q]=0;
PTS（m）;
}
//
voidcalCosSin（double*m）
{
doubleapp=m[（p+1）*p/2];
doubleaqq=m[（q+1）*q/2];
doubleapq=m[（p-1）*p/2+q];
cosz=cos（atan（2*apq/（app-aqq））/2）;
sinz=sin（atan（2*apq/（app-aqq））/2）;
}
//
voidfind_pq（double*m）
{
doublemax=0.0;
intpp=0;
intqq=0;
for（inti=1;i<=V;i++）
{
if（fabs（m[i]）>max）
{
recounti（i,&pp,&qq）;
if（pp!
=qq）
{
max=fabs（m[i]）;
p=pp;
q=qq;
}
}
}
MAX=max;
}
voidinit（double*m）
{
for（inti=1;i<=N;i++）
m[（i+1）*i/2]=a（i）;
for（inti=2;i<=N;i++）
m[（i-1）*i/2+i-1]=b;
for（inti=3;i<=N;i++）
m[（i-1）*i/2+i-2]=c;
PTS（m）;
}
voidcalFinal（double*m）
{
printf（"---------------------------------------------------------------------------------------------------\n"）;
printf（"结果输出:
\n\n"）;
doubleconda;
doubledeta=1.0;
doubleminlumda=pow（（double）10.0,12）;
doublemaxlumda=pow（（double）10.0,-12）;
doubleabsminlumda=pow（（double）10.0,12）;
for（inti=1;i<=N;i++）
{
if（m[（i+1）*i/2]>maxlumda）
maxlumda=m[（i+1）*i/2];
if（m[（i+1）*i/2]minlumda=m[（i+1）*i/2];
if（fabs（m[（i+1）*i/2]）absminlumda=fabs（m[（i+1）*i/2]）;
deta*=m[（i+1）*i/2];
}
if（fabs（minlumda）conda=fabs（maxlumda）/absminlumda;
else
conda=fabs（minlumda）/absminlumda;
printf（"Lumda
（1）=%.11eLumda（%d）=%.11eLumda（s）=%.11e\n",minlumda,N,maxlumda,absminlumda）;
printf（"Cond（A）=%.11e\n",conda）;
printf（"Det（A）=%.11e\n\n",deta）;
for（inti=1;i<=FK;i++）
{
doublemuk=minlumda+i*（maxlumda-minlumda）/40;
doublelumdak=0.0;
doubletempabsmin=pow（（double）10.0,12）;
for（intj=1;j<=N;j++）
{
if（fabs（muk-m[（j+1）*j/2]）{
lumdak=m[（j+1）*j/2];
tempabsmin=fabs（muk-m[（j+1）*j/2]）;
}
}
printf（"Lumda（i%d）=%.11e",i,lumdak）;
if（i%3==0）
putchar（10）;
}
putchar（10）;
printf（"------------------------------------------------------------------------------------------------------\n"）;
putchar（10）;
putchar（10）;
}
int_tmain（intargc,_TCHAR*argv[]）
{
doublem[（N+1）*N/2+1]={0.0};
kk=0;
MAX=1.0;
time_tt0,t1;
t0=time（&t0）;
init（m）;
#ifndefPTS
printf（"正在计算...\n\n"）;
#endif
while（true）
{
kk++;
find_pq（m）;
if（MAX#ifdefPTS
printf（"p=%dq=%d|max|=%e\n",p,q,MAX）;
printf（"-----------------------------------------------------------------------\n\n"）;
#endif
calCosSin（m）;
refreshMetrix（m）;
}
#ifdefPTS
printf（"p=%dq=%d|max|=%e\n",p,q,MAX）;
printf（"-----------------------------------------------------------------------\n\n"）;
#endif
printf（"矩阵最终形态...\n"）;
for（inti=1;i<=PFrols;i++）
{
for（intj=（i-1）*i/2+1;j<=（i+1）*i/2;j++）
{
printf（PFbits,m[j]）;
}
putchar（10）;
}
for（inti=1;i<=PFrols+1;i++）
{
printf（"..."）;
}
putchar（10）;
printf（"..\n"）;
printf（"..\n"）;
printf（"..\n"）;
for（inti=1;i<=PFrols+2;i++）
{
printf（"..."）;
}
putchar（10）;
t1=time（&t1）;
#ifdefPTS
printf（"计算并输出用时%.2f秒\n\n",difftime（t1,t0））;
#else
printf（"迭代次数%d,计算用时%.2f秒\n\n",kk,difftime（t1,t0））;
#endif
calFinal（m）;
return0;
}
运行结果如下：
中间运行状态如下：
结果分析
Ø有效性分析
1.由输出结果可见矩阵经过21840次迭代后，非对角元全部为零或接近于零；
2.代码中有定义预编译宏//#definePTS控制程序运行过程是否输出中间迭代结果，如果输出中间迭代结果，可以发现对角元素在迭代的后期变化非常小，达到收敛的效果；
3.算法在多次运行中基本可以在45秒左右完成计算（酷睿i5双核处理器，10G存，64位windows10操作系统）。
由此三条可得出结论：
算法是有效的。
Ø时空效率分析
1.算法的空间复杂度为o（N2），与矩阵的维数相关，当矩阵维数较高时，算法可能在空间的开销上比较大，目前没有想到更好的数据结构来节省空间开销。
2.算法的时间复杂度主要取决于矩阵的收敛速度，而每次迭代的时间复杂度为o（N2），主要耗费在顺序查找非对角线最大元和数组下标分解为矩阵坐标这两个子算法上。
目前可以考虑在查找非对角线最大元的算法上做一些改进：
在每次查找的过程中记录下上次查找到的第二大元素（第一大元素会被清零），然后在矩阵刷新之后的查找过程中，用上次查找到的第二大元素和刷新过的元素做比较得出本次迭代的最大元和第二大元。
Ø数值角度分析
1.计算过程中，涉及到矩阵变元的量均采用double型变量存储，满足精度要求；
2.由计算结果可见矩阵的条件数较大，可能会由于舍入误差的影响对结果造成较大扭曲。
升级会员 