数值分析计算实习二.docx
《数值分析计算实习二.docx》由会员分享，可在线阅读，更多相关《数值分析计算实习二.docx（19页珍藏版）》请在冰豆网上搜索。
数值分析计算实习二
题目：
试求矩阵
的全部特征值，并对其中的每个实特征值求相应的特征向量，已知
说明：
1.要求迭代的精度水平为
。
2.打印以下内容：
（1）采用带双步位移的QR分解法，说明算法设计方案和思路；
（2）全部源程序；
（3）矩阵
经过拟上三角化后所得的矩阵
；
（4）对矩阵
实行QR方法迭代结束后所得的矩阵；
（5）矩阵A的全部特征值
其中
，
。
若
是实数，则令
=0；
（6）
的相应于实特征值的特征向量；
（7）讨论：
发现的现象与遇到的问题。
3.采用e型数输出实型数，并且至少显示12位有效数字。
算法设计方案：
a）声明二维10×10数组空间，对其进行赋值得到矩阵
；
b）将矩阵
拟上三角化得到矩阵
，具体算法见课本P61-62；
c）对矩阵
进行双步位移QR分解迭代，具体算法大致同课本P63-65，但有改进。
改进部分为：
先判断迭代矩阵右下角是否为二阶分块矩阵，若是，再求解其一对特征值。
d）根据求得的实特征值计算对应的特征向量。
因为特征值已知，只需要求解以下线性方程组的一非零解即可：
为了方便起见，直接利用QR分解函数分解矩阵B，则原方程组等价于
然后利用直接法解出方程组即可。
其中矩阵R右下角元素必为0，则使对应行未知数为1，然后回代即可。
当矩阵R对应特征值的重数为2时，在矩阵R中除了右下角元素为0外，必有另一个对角元素为0，则令最后行对应未知数为0，该行对应未知数为1.当重数大于2时依次类推，即可求出对应于特定重数特征值的特征向量。
结果：
1、
拟上三角矩阵
如下。
2、迭代最终矩阵如下，迭代次数为14次。
3、矩阵所有特征值如下。
从上至下10个特征值，其中左列为实部，右列为虚部，可知有两对复特征值。
4、对应于实特征值的特征向量如下
问题发现与讨论
1、双步位移QR迭代过程中迭代矩阵
始终为拟上三角矩阵。
2、双步位移QR方法与基本QR方法比较
基本QR方法迭代过程中迭代矩阵
始终为拟上三角矩阵。
下图为基本QR方法得到的分块上三角矩阵：
下图为基本QR方法得到的特征值：
在基本QR迭代中发现，当矩阵收敛到分块上三角矩阵后，再次迭代它仍旧为分块上三角矩阵，然而矩阵中的元素改变了，不过块的阶数和位置不会改变。
基本QR迭代收敛速度为545次。
可以看出虽然两者收敛到的分块上三角矩阵不同，但得到的最终特征值是一样的。
也从而可初步验证所求得得特征值的正确性。
收敛速度上基本QR方法明显慢于双步位移QR方法。
源程序：
#include"stdafx.h"
#include"stdio.h"
#include"math.h"
#include"stdlib.h"
#defineN10
#defineEp1e-12
int_tmain（intargc,_TCHAR*argv[]）
{
intar0（doublea[N][N],ints,intn）;//判断矩阵多少位后是0
intsgn（doublex）;//0为1的符号函数
voidatra（doublea[N][N],intn）;//矩阵转置
voidaadd（doublea[N],doubleb[N],doublec[N],intn,intco）;//向量加法
voidaaadd（doublea[N][N],doubleb[N][N],doublec[N][N],intn,intco）;//矩阵加法
voidnamul（doublea[N],doublex,doubleb[N],intn）;//向量数乘
voidnaamul（doublea[N][N],doublex,doubleb[N][N],intn）;//矩阵数乘
voidamul（doublea[N],doubleb[N],doublec[N][N],intn）;//向量乘法
doubleaneiji（doublea[N],doubleb[N],intn）;//向量内积
voidaamul（doublea[N],doubleb[N][N],doublec[N],intn）;//矩阵向量乘法
voidaaamul（doublea[N][N],doubleb[N][N],doublec[N][N],intn）;//矩阵矩阵乘法
voidsol2（doublea11,doublea12,doublea21,doublea22,double（*p1）[2],double（*p2）[2]）;//计算二阶矩阵的特征值
voidqrs（doublea[N][N],doubleq[N][N]）;
voidGauss（doublea[N][N],doubley[N][N],int*p）;
doubleA[N][N];
doublenmd[N][2]={0};
for（inti=0;i{
for（intj=0;j{
if（i==j）A[i][j]=1.52*cos（i+1+1.2*（j+1））;
elseA[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
}
}
for（intr=0;r{
if（ar0（A,r+2,r）==0）continue;
else
{
doubled=0;
for（inti=r+1;idoubledr=sqrt（d）;
doublecr=-sgn（A[r+1][r]）*dr;
doublehr=pow（cr,2）-cr*A[r+1][r];
doubleur[N]={0};
ur[r+1]=A[r+1][r]-cr;
for（inti=r+2;idoublepr[N];
atra（A,N）;
aamul（ur,A,pr,N）;
namul（pr,1/hr,pr,N）;
doubleqr[N];
atra（A,N）;
aamul（ur,A,qr,N）;
namul（qr,1/hr,qr,N）;
doubletr;
tr=aneiji（pr,ur,N）/hr;
doublewr[N];
namul（ur,-tr,wr,N）;
aadd（qr,wr,wr,N,1）;
doubleB[N][N];
amul（wr,ur,B,N）;
aaadd（A,B,A,N,-1）;
amul（ur,pr,B,N）;
aaadd（A,B,A,N,-1）;//求一步操作后的矩阵Ar
}
}
for（inti=0;i{
for（intj=0;j{
printf（"%.11lf",A[i][j]）;
if（j==N-1）printf（"\n"）;
}
}
printf（"\n"）;//打印拟上三角矩阵Ar
doubleAqr[N][N]={0};
aaadd（Aqr,A,Aqr,N,1）;
doubleqqr[N][N]={0};
for（inti=0;i<10000;i++）
{
doubleAqr1[N][N]={0};
aaadd（Aqr,Aqr1,Aqr1,N,1）;
qrs（Aqr,qqr）;
doubleB[N][N]={0};
aaamul（Aqr,qqr,B,N）;
for（ints=0;s{
for（intk=0;kAqr[s][k]=B[s][k];
}
ints=0;
for（intt=0;t{
if（fabs（Aqr[t+1][t]）}
if（s==N-2）break;
}
for（inti=0;i{
for（intj=0;j{
printf（"%.11lf",Aqr[i][j]）;
if（j==N-1）printf（"\n"）;
}
}
//打印A（n-1）进行qr迭代后的最终结果
intm=N-1;
intk=0;
while（m>-1）
{
if（（fabs（A[m][m-1]）{
nmd[m][0]=A[m][m];
nmd[m][1]=0;
m=m-1;
continue;
}
elseif（（fabs（A[m-1][m-2]）{
sol2（A[m-1][m-1],A[m-1][m],A[m][m-1],A[m][m],（nmd+m）,（nmd+（m-1）））;//求解二阶矩阵的特征值
m=m-2;
continue;
}
else
{
doubles=A[m-1][m-1]+A[m][m];
doublet=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m];
doubleM[N][N]={0};
doubleD[N][N]={0};
aaamul（A,A,M,m+1）;
naamul（A,s,D,m+1）;
aaadd（M,D,M,m+1,-1）;
for（inti=0;iM[i][i]=M[i][i]+t;
for（intr=0;r{
if（ar0（A,r+1,r）==0）continue;
else
{
doubled=0;
for（inti=r;id=d+pow（M[i][r],2）;
doubledr=sqrt（d）;
doublecr=-sgn（M[r][r]）*dr;
doublehr=pow（cr,2）-cr*M[r][r];
doubleur[N]={0};
ur[r]=M[r][r]-cr;
for（inti=r+1;idoublevr[N];
atra（M,m+1）;
aamul（ur,M,vr,m+1）;
namul（vr,1/hr,vr,m+1）;
amul（ur,vr,D,m+1）;
atra（M,m+1）;
aaadd（M,D,M,m+1,-1）;
doublepr[N];
atra（A,m+1）;
aamul（ur,A,pr,m+1）;
namul（pr,1/hr,pr,m+1）;
doubleqr[N];
atra（A,m+1）;
aamul（ur,A,qr,m+1）;
namul（qr,1/hr,qr,m+1）;
doubletr;
tr=aneiji（pr,ur,m+1）/hr;
doublewr[N];
namul（ur,-tr,wr,m+1）;
aadd（qr,wr,wr,m+1,1）;
amul（wr,ur,D,m+1）;
aaadd（A,D,A,m+1,-1）;
amul（ur,pr,D,m+1）;
aaadd（A,D,A,m+1,-1）;
}
}
}
k++;
if（k>10000）//算法失败判定
{
printf（"算法失败\n"）;
break;
}
}
for（inti=0;i{
for（intj=0;j{
printf（"%.11lf",A[i][j]）;
if（j==N-1）printf（"\n"）;
}
}
for（inti=0;i{
printf（"%1.11e",nmd[i][0]）;
printf（"%1.11e\n",nmd[i][1]）;
}
for（inti=0;i{
if（fabs（nmd[i][1]）>Ep）continue;
for（intj=i+1;j{
if（fabs（nmd[j][1]）>Ep）continue;
elseif（fabs（nmd[j][0]-nmd[i][0]）}
}
for（intk=0;k{
if（fabs（nmd[k][1]）>Ep）continue;
for（inti=0;i{
for（intj=0;j{
if（i==j）A[i][j]=1.52*cos（i+1+1.2*（j+1））-nmd[k][0];
elseA[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
}
}
doubleq[N][N];
doubley[N][N];
intp[1]={0};
qrs（A,q）;//对A-I*nmd进行QR分解得到新矩阵R
Gauss（A,y,p）;//齐次不满秩的高斯消去法,求非0根
printf（"对应于λ%d实特征值的特征向量为\n",k+1）;
for（intj=0;j{
for（inti=0;iprintf（"%1.12e\n",y[i][j]）;
}
}
return0;
}
intar0（doublea[N][N],ints,intn）
{
for（inti=s;i{
if（fabs（a[i][n]）>Ep）
{
return1;
break;
}
}
return0;
}
intsgn（doublex）
{
if（x>0）return1;
elsereturn-1;
}
voidatra（doublea[N][N],intn）
{
doublex;
for（inti=0;i{
for（intj=i+1;j{
x=a[i][j];
a[i][j]=a[j][i];
a[j][i]=x;
}
}
}
voidaadd（doublea[N],doubleb[N],doublec[N],intn,intco）
{
if（co==1）
for（inti=0;ielse
{
if（co==-1）
for（inti=0;ielse
printf（"向量加减选择错误\n"）;
}
}
voidaaadd（doublea[N][N],doubleb[N][N],doublec[N][N],intn,intco）
{
if（co==1）
{
for（inti=0;i{
for（intj=0;j{
c[i][j]=a[i][j]+b[i][j];
}
}
}
else
{
if（co==-1）
{
for（inti=0;i{
for（intj=0;j{
c[i][j]=a[i][j]-b[i][j];
}
}
}
else
printf（"矩阵加减选择错误\n"）;
}
}
voidnamul（doublea[N],doublex,doubleb[N],intn）
{
for（inti=0;ib[i]=x*a[i];
}
voidnaamul（doublea[N][N],doublex,doubleb[N][N],intn）
{
for（inti=0;ifor（intj=0;jb[i][j]=x*a[i][j];
}
voidamul（doublea[N],doubleb[N],doublec[N][N],intn）
{
for（inti=0;ifor（intj=0;jc[i][j]=a[i]*b[j];
}
doubleaneiji（doublea[N],doubleb[N],intn）
{
doublet=0;
for（inti=0;i{
t=t+a[i]*b[i];
}
returnt;
}
voidaamul（doublea[N],doubleb[N][N],doublec[N],intn）
{
for（inti=0;i{
c[i]=0;
for（intj=0;j{
c[i]=c[i]+b[i][j]*a[j];
}
}
}
voidaaamul（doublea[N][N],doubleb[N][N],doublec[N][N],intn）
{
for（inti=0;i{
for（intj=0;j{
c[i][j]=0;
for（ints=0;sc[i][j]=c[i][j]+a[i][s]*b[s][j];
}
}
}
voidsol2（doublea11,doublea12,doublea21,doublea22,double（*p1）[2],double（*p2）[2]）
{
doubledeta;
deta=a11*a11+a22*a22-2*a11*a22+4*a12*a21;
if（deta<0）
{
（*p1）[0]=（a11+a22）/2;
（*p2）[0]=（a11+a22）/2;
（*p1）[1]=sqrt（0-deta）/2;
（*p2）[1]=0-sqrt（0-deta）/2;
}
else
{
（*p1）[0]=（a11+a22）/2+sqrt（deta）/2;
（*p1）[1]=0;
（*p2）[0]=（a11+a22）/2-sqrt（deta）/2;
（*p2）[1]=0;
}
}
voidqrs（doublea[N][N],doubleq[N][N]）
{
for（inti=0;i{
for（intj=0;j{
if（i==j）q[i][j]=1;
elseq[i][j]=0;
}
}
for（intr=0;r{
if（ar0（a,r+1,r）==0）continue;
else
{
doubled=0,dr;
doublecr;
doublehr;
for（inti=r
升级会员 