北航数值分析大作业第二次的.docx
- 文档编号:2724634
- 上传时间:2023-05-04
- 格式:DOCX
- 页数:22
- 大小:26.32KB
北航数值分析大作业第二次的.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;j if(m==(N-2-i))continue; for(j=i+1,d=0;j d=sqrt(d); c=-1*sgn(A[i+1][i])*d; h=c*c-c*A[i+1][i]; for(j=i+2;j for(j=0;j u[i+1]=A[i+1][i]-c; for(j=0;j {for(k=i+1,p[j]=0;k for(j=0;j {for(k=i+1,q[j]=0;k for(j=0,t=0;j for(j=0;j for(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;j for(i=0;i {for(j=0;j for(i=0;i { //for(j=i+1;j //if(m==(N-1-i))continue; for(j=i,d=0;j d=sqrt(d); c=-1*sgn(R[i][i])*d; h=c*c-c*R[i][i]; for(j=i+1;j for(j=0;j u[i]=R[i][i]-c; for(j=0;j for(j=0;j for(j=0;j {for(k=i,p[j]=0;k for(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;j u[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;i for(t=0;t<6;t++) { for(i=0;i for(i=0;i for(i=0;i { for(j=i;j for(j=i;j for(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;j printf("\n"); }printf("\n"); for(i=0;i { for(j=N-5;j printf("\n"); }printf("\n"); printf("生成的正交矩阵Q: \n"); for(i=0;i { for(j=0;j printf("\n"); }printf("\n"); for(i=0;i { for(j=N-5;j printf("\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;j printf("\n"); }printf("\n");printf("\n");printf("\n"); for(i=0;i { for(j=N-5;j printf("\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;j printf("\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
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业 第二次
![提示](https://static.bingdoc.com/images/bang_tan.gif)