北航数值分析实习第二题.docx
- 文档编号:80439
- 上传时间:2022-10-02
- 格式:DOCX
- 页数:17
- 大小:267.83KB
北航数值分析实习第二题.docx
《北航数值分析实习第二题.docx》由会员分享,可在线阅读,更多相关《北航数值分析实习第二题.docx(17页珍藏版)》请在冰豆网上搜索。
数值分析计算实习报告
第二题
所在班级
A1班
学生姓名
学生学号
2015年11月
北航数值分析计算实习报告 第16页
1算法设计方案
1.1矩阵的输入
A矩阵是一个10乘10的矩阵,同时并没有对称性,所以不能压缩,直接存储即可。
1.2对矩阵进行上三角化
A矩阵的上三角化采用常规方法即可,得到。
1.3QR分解法
对上三角化后的矩阵进行QR分解可以减少计算量。
RQ这个矩阵是与原矩阵A相似的矩阵,与A具有相同的特征值和特征向量。
1.4求解特征值
采用带双步位移的QR分解法对求解特征值即可。
1.5求解实特征值的特征向量
求解特征向量,本质上就是求解方程
的解,其中I是单位向量。
直接采用列主元素Gauss消去法求解即可。
但是,值得注意的是,对应于一个特征值的特征向量是无穷多的,也就是经过Gauss消去法后,最后一行全是0。
因此,求解时候,需要给特征向量某个元素赋值,本题,将最后一个元素统一赋值为1,即可求出对应于某个特征值的特征向量的基础解系,而全部特征向量为,k取不等于0的任意实数。
2C++程序
#include
#include
#include
#include
#definen11
#defineerr1e-12
#defineL2500
voidcaculateA();//计算矩阵A的系数
voidhessenberg();//将A拟上三角化
voidQR();//对矩阵进行QR分解
voidQRshuangbu();//对矩阵进行带双步位移的QR分解
voidgauss();//列主元的高斯消元法求解特征向量
inti,j,s,p,k,ik,nR,nC;
doubleA[n][n],q[n][n],r[n][n],rq[n][n],I[n][n];
doubleP[n],W[n],u[n],Q[n];
doubledr,cr,hr,ar,tr;
doubles1,t,x,tzR[n],tzC[2][n],sum,M[n][n],v[n];
voidmain()
{for(i=1;i {for(j=1;j {if(i==j) {I[i][j]=1;q[i][j]=1;} else{I[i][j]=0;q[i][j]=0;}}} caculateA(); hessenberg(); QR(); QRshuangbu(); gauss(); printf("运行时间为%d\n",clock()); } voidcaculateA()//计算矩阵A的系数 { for(i=1;i {for(j=1;j {if(j! =i) A[i][j]=sin(0.5*i+0.2*j); else A[i][j]=1.52*cos(i+1.2*j);} } } voidhessenberg()//将A拟上三角化 { for(s=1;s {for(ar=0.0,i=s+2;i ar+=A[i][s]*A[i][s]; if(ar<=err) continue; else{ ar+=A[s+1][s]*A[s+1][s]; dr=sqrt(ar); if(A[s+1][s]>0) cr=-dr; else cr=dr; hr=cr*cr-cr*A[s+1][s]; for(i=1;i<=s;i++) u[i]=0.0; u[s+1]=A[s+1][s]-cr; for(i=s+2;i u[i]=A[i][s]; for(j=1;j {for(P[j]=0.0,i=1;i P[j]+=A[i][j]*u[i]/hr;} for(tr=0.0,i=1;i {tr+=P[i]*u[i]/hr;} for(i=1;i {for(Q[i]=0.0,j=1;j Q[i]+=A[i][j]*u[j]/hr;} for(i=1;i {W[i]=Q[i]-tr*u[i];} for(i=1;i {for(j=1;j A[i][j]-=(W[i]*u[j]+u[i]*P[j]);} } } for(i=1;i {for(j=1;j if(fabs(A[i][j]) A[i][j]=0.0;} printf("数值分析计算实习第二次作业\n"); printf("A1班\n"); printf("\n"); printf("矩阵A经过拟上三角化后得到的矩阵为: \n"); for(i=1;i { for(j=1;j printf("%1.12e,",A[i][j]); printf("\n");} } voidQR()//对矩阵进行QR分解 { doubleu[n],w[n],F[n]; for(s=1;s for(p=1;p r[s][p]=A[s][p]; for(s=1;s { for(dr=0.0,i=s;i dr+=r[i][s]*r[i][s]; dr=sqrt(dr); if(dr<=err) continue; else {if(A[s][s]>0)cr=-dr; elsecr=dr; hr=cr*cr-cr*r[s][s]; for(i=1;i u[i]=0; u[s]=r[s][s]-cr; for(i=s+1;i u[i]=r[i][s]; for(i=1;i {for(F[i]=0.0,j=1;j F[i]+=r[j][i]*u[j]/hr;} for(i=1;i {for(j=1;j r[i][j]=r[i][j]-u[i]*F[j];} for(i=1;i {for(w[i]=0.0,j=1;j w[i]+=q[i][j]*u[j];} for(i=1;i {for(j=1;j q[i][j]-=w[i]*u[j]/hr;} } } for(i=1;i {for(j=1;j {for(rq[i][j]=0.0,s=1;s rq[i][j]+=r[i][s]*q[s][j];}} printf("\n"); printf("A经过QR分解后得到的正交矩阵Q为: \n"); for(i=1;i {for(j=1;j if(fabs(q[i][j]) q[i][j]=0.0;} for(i=1;i {for(j=1;j printf("%1.12e,",q[i][j]); printf("\n");} printf("\n"); printf("A经过QR分解后得到的上三角矩阵R为: \n"); for(i=1;i {for(j=1;j if(fabs(r[i][j]) r[i][j]=0.0;} for(i=1;i {for(j=1;j printf("%1.12e,",r[i][j]); printf("\n");} printf("\n"); printf("正交矩阵Q乘以上三角矩阵R得到的RQ为: \n"); for(i=1;i {for(j=1;j if(fabs(rq[i][j]) rq[i][j]=0.0;} for(i=1;i {for(j=1;j printf("%1.12e,",rq[i][j]); printf("\n");} } voidQRshuangbu()//对矩阵进行带双步位移的QR分解 { intK=1,m=10; nR=0,nC=0; loop3: if(fabs(A[m][m-1])<=err) {nR++; tzR[nR]=A[m][m]; m--; gotoloop4; } elsegotoloop5; loop4: if(m==1) { nR++; tzR[nR]=A[1][1]; gotoloop9; } elseif(m==2) { s1=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]; x=s1*s1-4*t; if(x>=0) { nR++; tzR[nR]=(s1+sqrt(x))/2; nR++; tzR[nR]=(s1-sqrt(x))/2; } else { nC++; tzC[0][nC]=s1/2; tzC[1][nC]=sqrt(-x)/2; nC++; tzC[0][nC]=s1/2; tzC[1][nC]=-sqrt(-x)/2; } gotoloop9; } elsegotoloop3; loop5: { if(fabs(A[m-1][m-2])<=err) { s1=A[m
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北航 数值 分析 实习 第二
![提示](https://static.bdocx.com/images/bang_tan.gif)