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

数值分析计算实习报告第二题所在班级A1 班学生姓名学生学号2015年11月算法设计方案矩阵的输入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 2500void caculateA();//计算矩阵A的系数void hessenberg();//将A拟上三角化void QR();//对矩阵进行QR分解void QRshuangbu();//对矩阵进行带双步位移的QR分解void gauss();//列主元的高斯消元法求解特征向量inti,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];doubledr,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;elsecr=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++) {for(j=1;jn;j++) A[i][j]-=(W[i]*u[j]+u[i]*P[j]);} } }for(i=1;in;i++) {for(j=1;jn;j++)if(fabs(A[i][j])err) A[i][j]=0.0;}printf(数值分析计算实习第二次作业\n);printf( A1班\n);printf(\n);printf(矩阵A经过拟上三角化后得到的矩阵为:\n);for(i=1;in;i++) {for(j=1;jn;j++) printf(%1.12e,,A[i][j]);printf(\n);}}void QR()//对矩阵进行