北航数值分析大作业第二次的.docx
《北航数值分析大作业第二次的.docx》由会员分享，可在线阅读，更多相关《北航数值分析大作业第二次的.docx（22页珍藏版）》请在冰点文库上搜索。
北航数值分析大作业第二次的
数值分析
第二次大作业
姓名：
路子威
学号：
S2*******1
题目：
使用带双步位移的QR分解法求矩阵
的全部特征值，并对其中的每一个实特征值求相应的特征向量。
已知：
（i,j=1,2,……,10）
要求：
1.用幂法、反幂法和QR方法求矩阵的特征值时，要求迭代精度水平为e=10-12。
2.打印以下内容：
（1）全部源程序；
（2）矩阵A经过拟上三角化后所得的矩阵A（n-1）；
（3）对矩阵A（n-1）实行QR方法迭代结束后所得的矩阵；
（4）矩阵A的全部特征值λi=（RI,Ii）;（i=1,2,…,10）
（5）A的相应于实特征值的特征向量。
3.采用e型输出实型数，并且至少显示12位有效数字。
算法：
1）输入需要求解的矩阵A
2）对上述生成的矩阵进行拟上三角化
3）对拟上三角化后的矩阵进行QR分解（带双步位移）。
4）用函数characteristic（）求解矩阵A（n-1）即A的所有特征值。
5）用函数characteristicvector（）求解矩阵A（n-1）即A的所有特征向量。
6）输出A的所有特征值λ、A的所有实特征值对应的特征向量、拟上三角矩阵A（n-1）、及其Q、R和R*Q
说明：
为了减少求特征值和特征向量过程中的计算量，在对矩阵进行QR分解前先进行拟上三角化。
拟上三角化之后矩阵中包含有大量的0元素，在进行QR迭代时可以减少大量的0元素乘法运算，节省计算时间。
且上三角矩阵进行QR分解之后，通过RQ得到新的矩阵A仍然是上三角矩阵，因此可以在整个迭代过程中使用。
而带双步位移的QR分解法可以加快收敛速度，使矩阵在更短的迭代次数下完成收敛。
源程序：
#include
#include
#include
#definee0.000000000001/*设置精度水平*/
#defineN10/*设置矩阵阶数*/
#defineL2500/*设置迭代次数*/
/*子程序*/
/*sgn函数*//*保证ci与aii符号相反*/
doublesgn（doublea）
{
doublex;
if（a>e）x=1;
elsex=-1;
returnx;
}
/*求根函数*//*计算双步位移法中的两个位移量，即右下角二阶子式的两个特征值*/
doublepig（doubleb,doublec）
{
doublem;
m=b*b-4*c;
returnm;
}
/*矩阵的拟上三角化*//*只要对A1进行拟上三角化就能保证每个矩阵Ak都为拟上三角矩阵，大大简化计算量*/
voidhessenberg（doubleA[N][N]）
{
inti,j,k;
intm=0;
doubled,c,h,t;
doubleu[N],p[N],q[N],w[N];
for（i=0;i{
for（j=i+2;jif（m==（N-2-i））continue;
for（j=i+1,d=0;jd=sqrt（d）;
c=-1*sgn（A[i+1][i]）*d;
h=c*c-c*A[i+1][i];
for（j=i+2;jfor（j=0;ju[i+1]=A[i+1][i]-c;
for（j=0;j{for（k=i+1,p[j]=0;kfor（j=0;j{for（k=i+1,q[j]=0;kfor（j=0,t=0;jfor（j=0;jfor（j=0;j{
for（k=0;k}
}
}
/*矩阵的QR分解*/
voidQRdiv（doubleA[N][N],doubleQ[N][N],doubleR[N][N]）
{
inti,j,k;
//intm=0;
doubled,c,h;
doubleu[N],w[N],p[N];
for（i=0;i{for（j=0;jfor（i=0;i{for（j=0;jfor（i=0;i{
//for（j=i+1;j//if（m==（N-1-i））continue;
for（j=i,d=0;jd=sqrt（d）;
c=-1*sgn（R[i][i]）*d;
h=c*c-c*R[i][i];
for（j=i+1;jfor（j=0;ju[i]=R[i][i]-c;
for（j=0;jfor（j=0;jfor（j=0;j{for（k=i,p[j]=0;kfor（j=0;j{
for（k=0;k}
}
}
/*求解矩阵的所有特征值*/
voidcharacteristic（doubleA[N][N],doublechaR[N],doublechaI[N]）
{
intk=0,m=N-1;
inti,j;
intl;
doubles,t,x;
doubleM[N][N],B[N][N];
intf=0;
doubled,c,h;
doubleu[N],w[N],p[N];
doubleQ[N][N],R[N][N];
for（l=0;l{
next:
if（m==0）{chaR[0]=A[0][0];chaI[0]=0;break;}
if（fabs（A[m][m-1]）<=e）
{
chaR[m]=A[m][m];chaI[m]=0;
m--;
gotonext;
}
s=A[m-1][m-1]+A[m][m];
t=A[m-1][m-1]*A[m][m]-A[m][m-1]*A[m-1][m];
if（m==1）
{
x=pig（s,t）;
if（x>=e）{x=sqrt（x）;chaR[m]=s/2+x/2;chaR[m-1]=s/2-x/2;chaI[m]=0;chaI[m-1]=0;}
else{x=sqrt（fabs（x））;chaR[m]=s/2;chaR[m-1]=s/2;chaI[m]=x/2;chaI[m-1]=-x/2;}
break;
}
if（fabs（A[m-1][m-2]）<=e）
{
x=pig（s,t）;
if（x>=e）{x=sqrt（x）;chaR[m]=s/2+x/2;chaR[m-1]=s/2-x/2;chaI[m]=0;chaI[m-1]=0;}
else{x=sqrt（fabs（x））;chaR[m]=s/2;chaR[m-1]=s/2;chaI[m]=x/2;chaI[m-1]=-x/2;}
m=m-2;
gotonext;
}
for（i=0;i<=m;i++）
{
for（j=0;j<=m;j++）
{
if（i==j）
{
for（k=0,M[i][j]=0;k<=m;k++）M[i][j]=A[i][k]*A[k][j]+M[i][j];
M[i][j]=M[i][j]-s*A[i][j]+t;
}
else
{
for（k=0,M[i][j]=0;k<=m;k++）M[i][j]=A[i][k]*A[k][j]+M[i][j];
M[i][j]=M[i][j]-s*A[i][j];
}
}
}
/*带双步位移的QR分解中M的分解*/
for（i=0;i<=m;i++）
{for（j=0;j<=m;j++）{if（i==j）Q[i][j]=1;elseQ[i][j]=0;}}
for（i=0;i<=m;i++）
{for（j=0;j<=m;j++）R[i][j]=M[i][j];}
for（i=0;i{
for（j=i+1;j<=m;j++）if（R[j][i]<=e）f=f+1;
if（f==（m-i））continue;
for（j=i,d=0;j<=m;j++）d=d+R[j][i]*R[j][i];
d=sqrt（d）;
c=-1*sgn（R[i][i]）*d;
h=c*c-c*R[i][i];
for（j=i+1;j<=m;j++）u[j]=R[j][i];
for（j=0;ju[i]=R[i][i]-c;
for（j=0;j<=m;j++）
{for（k=0,w[j]=0;k<=m;k++）w[j]=Q[j][k]*u[k]+w[j];}
for（j=0;j<=m;j++）
{for（k=0;k<=m;k++）Q[j][k]=Q[j][k]-w[j]*u[k]/h;}
for（j=0;j<=m;j++）
{for（k=i,p[j]=0;k<=m;k++）p[j]=R[k][j]*u[k]+p[j];p[j]=p[j]/h;}
for（j=0;j<=m;j++）
{
for（k=0;k<=m;k++）R[j][k]=R[j][k]-u[j]*p[k];
}
}
for（j=0;j<=m;j++）
{
for（k=0;k<=m;k++）M[j][k]=Q[j][k];
}
for（i=0;i<=m;i++）
{
for（j=0;j<=m;j++）
{for（k=0,B[i][j]=0;k<=m;k++）B[i][j]=M[k][i]*A[k][j]+B[i][j];}
}
for（i=0;i<=m;i++）
{
for（j=0;j<=m;j++）
{for（k=0,A[i][j]=0;k<=m;k++）A[i][j]=B[i][k]*M[k][j]+A[i][j];}
}
}
}
/*求矩阵的所有特征值的特征向量*/
voideigenvector（doubleV[N][N],doubleT[N]）
{
doubleA[N][N],baoz[N][N],guod[N];
doublec;
inti,j,k,m,t;
intW=0;
for（i=0;ifor（t=0;t<6;t++）
{
for（i=0;ifor（i=0;ifor（i=0;i{
for（j=i;je）{k=j;break;}
for（j=i;jfor（j=i;j{
c=A[j][i];
if（fabs（c）>e）for（m=i;m}
for（j=0;j{
c=A[j][i];
if（j!
=i）{for（m=i;m}
}
V[t][N-1]=-1;
for（i=N-2;i>=0;i--）
{
V[t][i]=A[i][N-1];
}
}
}
/*主函数*/
voidmain（）
{
doublea[N][N],b[N][N],chaR[N],chaI[N];
doubleq[N][N],r[N][N],qr[N][N];
doubleDOG[N];
doublef,g;
inti,j,k;
for（i=0;i{
for（j=0;j{
if（i!
=j）a[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
elsea[i][j]=1.5*cos（（i+1）+1.2*（j+1））;
}
}
hessenberg（a）;
printf（"拟上三角化后A（n-1）:
\n"）;
for（i=0;i{
for（j=0;j{
if（fabs（a[i][j]）printf（"%22.11e",a[i][j]）;
}
printf（"\n"）;
}
printf（"\n"）;
for（i=0;i{
for（j=N-5;j{
if（fabs（a[i][j]）printf（"%22.11e",a[i][j]）;
}
printf（"\n"）;
}
printf（"\n"）;
QRdiv（a,q,r）;
printf（"\n"）;
printf（"\n"）;
printf（"\n"）;
printf（"生成的上三角矩阵R：
\n"）;
for（i=0;i{
for（j=0;jprintf（"\n"）;
}printf（"\n"）;
for（i=0;i{
for（j=N-5;jprintf（"\n"）;
}printf（"\n"）;
printf（"生成的正交矩阵Q：
\n"）;
for（i=0;i{
for（j=0;jprintf（"\n"）;
}printf（"\n"）;
for（i=0;i{
for（j=N-5;jprintf（"\n"）;
}printf（"\n"）;
for（i=0;i{
for（j=0;j}printf（"\n"）;printf（"\n"）;printf（"\n"）;
printf（"生成的R*Q阵：
\n"）;
for（i=0;i{
for（j=0;jprintf（"\n"）;
}printf（"\n"）;printf（"\n"）;printf（"\n"）;
for（i=0;i{
for（j=N-5;jprintf（"\n"）;
}printf（"\n"）;printf（"\n"）;printf（"\n"）;
characteristic（a,chaR,chaI）;
for（i=1;i{
if（i<3）{f=chaR[i];g=chaI[i];chaR[i]=chaR[7+i];chaI[i]=chaI[7+i];chaR[7+i]=f;chaI[7+i]=g;}
if（i==5）{f=chaR[i];g=chaI[i];chaR[i]=chaR[7];chaI[i]=chaI[7];chaR[7]=f;chaI[7]=g;}
}
printf（"矩阵的全部特征值：
\n"）;
for（j=0;j{
if（fabs（chaI[j]）<=e）printf（"λ%2d=%19.11e\n",j+1,chaR[j]）;
elseif（chaI[j]>e）printf（"λ%2d=%18.11e+%18.11ei\n",j+1,chaR[j],chaI[j]）;
elseprintf（"λ%2d=%19.11e%19.11ei\n",j+1,chaR[j],chaI[j]）;
}
printf（"\n"）;printf（"\n"）;printf（"\n"）;
for（i=0;i{
for（j=0;j{
if（i!
=j）a[i][j]=sin（0.5*（i+1）+0.2*（j+1））;
elsea[i][j]=1.5*cos（（i+1）+1.2*（j+1））;
}
}//重新输入矩阵A
for（i=0;i<6;i++）
{
printf（"λ%d对应的特征向量为：
\n",（i+1））;
for（j=0;jprintf（"\n"）;
}
getch（）;
}
计算结果：
拟上三角化后A（n-1）
-0.894521698228-0.0993313649183-1.09983175888-0.7665038709080.170760
114146-1.93488255889-0.08390208705250.913256511314-0.640797700919
0.194673367868
-2.347878362422.37205792161.827998552320.3266556884710.208236058364
2.088987009940.184********9-1.263015266080.67906946685-0.46721508865
01.73595446995-1.16502336748-1.24674444352-0.629822548908-1.98482
0180990.297575006080.633930059659-0.130********70.30403010361
0-3.42119688154e-017-1.29293756392-1.12623922591.19078291192
-1.30877298390.186********70.423673393688-0.1019600826550.194366091451
05.94010094099e-017-3.28119028842e-0171.577711153030.816935
8328160.446153172383-0.0436509254161-0.4665979167190.294123156618
-0.103442111366
0-3.49494913888e-0178.59590318918e-0179.17131065135e-018
-0.772897513499-1.60102824405-0.291268547483-0.2434337858320.673628608451
0.262477290494
0-2.13633528681e-016-3.94299506639e-017-8.80991580387e-017
0-0.729677394636-0.007965456279830.971073910201-0.129896736857
.027*********
0-6.48100629069e-0178.50487510592e-017-5.02328897109e-018
0-5.3832322921e-0170.794553961298-0.4525143454610.504890152758
-0.121121019351
07.034391
升级会员 