北航数值分析大作业第二题精解详解.docx
- 文档编号:25752788
- 上传时间:2023-06-12
- 格式:DOCX
- 页数:16
- 大小:205.89KB
北航数值分析大作业第二题精解详解.docx
《北航数值分析大作业第二题精解详解.docx》由会员分享,可在线阅读,更多相关《北航数值分析大作业第二题精解详解.docx(16页珍藏版)》请在冰豆网上搜索。
北航数值分析大作业第二题精解详解
目标:
使用带双步位移的QR分解法求矩阵
的全部特征值,并对其中的每一个实特征值求相应的特征向量。
已知:
(i,j=1,2,……,10)
开始
输入矩阵A
用函数“nishangsanjiadiv”将矩阵A拟上三角化即为A(n-1)
用函数“characteristic”求解矩阵A(n-1)即A的所有特征值
用函数“qrdiv”将A(n-1)QR分解
用函数“characteristicvector”求解矩阵A(n-1)即A的所有实特征值对应的特征向量
输出A的所有特征值λ、A的所有实特征值对应的特征向量、拟上三角矩阵A(n-1)、及其Q、R和R*Q
算法:
以上是程序运作的逻辑,其中具体的函数的算法,大部分都是数值分析课本上的逻辑,在这里特别写出矩阵A的实特征值对应的一个特征向量的求法:
由于
,因此在经过选主元的高斯消元以后,
即
的最后一行必然为零,左上方变为n-1阶单位矩阵
,右上方变为n-1阶向量
,然后令
,则
。
这样即求出所有A所有实特征值对应的一个特征向量。
#include
#include
#include
#defineN10
#defineE1.0e-12
#defineMAX10000
//以下是符号函数
doublesgn(doublea)
{
doublez;
if(a>E)z=1;
elsez=-1;
returnz;
}
//以下是矩阵的拟三角分解
voidnishangsanjiaodiv(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 } } }//矩阵的QR分解 //以下是二次多项式求根 doubleroot(doubleb,doublec) { doublem; m=b*b-4*c; returnm; }//二次多项式求根 //以下是求解矩阵的所有特征值 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=root(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=root(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]; } } } //以下是M的QR分解 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]; } //以上是M的QR分解 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]; doubleshiyan[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)); } } nishangsanjiaodiv(a); printf("矩阵A的拟上三角分解: \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"); qrdiv(a,q,r);printf("\n");printf("\n");printf("\n"); printf("拟上三角矩阵A的QR分解: \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("矩阵A所有特征值: \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 eigenvector(a,chaR);printf("相应实特征值对应的特征向量: \n"); for(i=0;i<6;i++) { printf("λ%d的一个特征向量为: \n",(i+1)); for(j=0;j printf("\n"); } getch(); }
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 作业 第二 题精解 详解
![提示](https://static.bdocx.com/images/bang_tan.gif)