北航数值分析实习第二题.docx

数值分析计算实习报告 第二题 所在班级 A1 班 学生姓名 学生学号 2015年11月 北航数值分析计算实习报告 第 PAGE \* MERGEFORMAT 16 页 算法设计方案 矩阵的输入 A矩阵是一个10乘10的矩阵，同时并没有对称性，所以不能压缩，直接存储即可。 对矩阵进行上三角化 A矩阵的上三角化采用常规方法即可，得到。 QR分解法 对上三角化后的矩阵进行QR分解可以减少计算量。 RQ这个矩阵是与原矩阵A相似的矩阵，与A具有相同的特征值和特征向量。 求解特征值 采用带双步位移的QR分解法对求解特征值即可。 求解实特征值的特征向量 求解特征向量，本质上就是求解方程 的解，其中I是单位向量。直接采用列主元素Gauss消去法求解即可。但是，值得注意的是，对应于一个特征值的特征向量是无穷多的，也就是经过Gauss消去法后，最后一行全是0。因此，求解时候，需要给特征向量某个元素赋值，本题，将最后一个元素统一赋值为1，即可求出对应于某个特征值的特征向量的基础解系，而全部特征向量为，k取不等于0的任意实数。 C++程序 #includestdio.h #includemath.h #includestdlib.h #includetime.h #define n 11 #define err 1e-12 #define L 2500 void caculateA();//计算矩阵A的系数 void hessenberg();//将A拟上三角化 void QR();//对矩阵进行QR分解 void QRshuangbu();//对矩阵进行带双步位移的QR分解 void gauss();//列主元的高斯消元法求解特征向量 int i,j,s,p,k,ik,nR,nC; double A[n][n],q[n][n],r[n][n],rq[n][n],I[n][n]; double P[n],W[n],u[n],Q[n]; double dr,cr,hr,ar,tr; double s1,t,x,tzR[n],tzC[2][n],sum,M[n][n],v[n]; void main() {for(i=1;in;i++) {for(j=1;jn;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()); } void caculateA()//计算矩阵A的系数 { for(i=1;in;i++) {for(j=1;jn;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);} } } void hessenberg()//将A拟上三角化 { for(s=1;sn-2;s++) {for(ar=0.0,i=s+2;in;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;in;i++) u[i]=A[i][s]; for(j=1;jn;j++) {for(P[j]=0.0,i=1;in;i++) P[j]+=A[i][j]*u[i]/hr;} for(tr=0.0,i=1;in;i++) {tr+=P[i]*u[i]/hr;} for(i=1;in;i++) {for(Q[i]=0.0,j=1;jn;j++) Q[i]+=A[i][j]*u[j]/hr;} for(i=1;in;i++) {W[i]=Q[i]-tr*u[i];} for(i=1;in;i++)