数值分析第二次大作业.doc

《数值分析第二次大作业.doc》由会员上传分享，免费在线阅读，更多相关内容在行业资料-天天文库。

1、数值分析第二次SY1107214王宽一、题目分析使用带双步位移的QR分解法求矩阵的全部特使用带双步位移的QR分解法求矩阵的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：(i,j=1,2,……,10)二、算法设计1.按照题目给出的矩阵定义对矩阵A赋初值：对应的函数为initA()；2.为了减少求特征值和特征向量过程中的计算量，在对矩阵进行QR分解前先进行拟上三角化：对应的函数为nssj()；3.对拟上三角化后的矩阵A使用带双步位移的QR分解法逐次迭代（最大迭代次数L=500），逐个求出其特征值，对应的函数为tezhengzhi()；这个过程中包含着两个子程序：QR

2、()和qmk()，分别用来对矩阵Mk进行QR分解并得到Ak+1和计算mk的值；4.使用带原点平移的反幂法求出其对应的特征向量，对应的函数为：fmifa()，这个过程中求解线性方程时用到列主元的高斯消元法，对应的函数为：gauss(double);根据数值分析课本的相关知识，步骤3中带双步位移的QR分解法的流程图如下：第18页数值分析第二次SY1107214王宽本次作业所使用的编译环境为：Visualc++6.0第18页数值分析第二次SY1107214王宽一、程序源代码#include#includedoubleA[10][10],rr[10]

3、[10],qq[10][10],rq[10][10],uk[10][10];doublexy[10],x[10]={0};//反幂法中doubletzr[10]={0},tzi[10]={0};//定义矩阵的特征值数组,r实部、i虚部//对矩阵A进行初始化，赋值voidinitA(){inti,j;for(i=1;i<=10;i++){for(j=1;j<=10;j++)A[i-1][j-1]=sin(0.5*i+0.2*j);}for(i=1;i<=10;i++)A[i-1][i-1]=cos(i+1.2*i)*1.5;}//定义函数对矩阵A进行拟上三角化voidnssj(

4、){intr,i,j;intsgn(doublea);doubled,c,h,t,tr,u[10],p[10],q[10],w[10];doublesum=0,tmp=0;//开始循环for(r=0;r<8;r++){t=0;for(i=r+2;i<10;i++)t=t+(A[i][r]==0);if(t==8-r)continue;else{sum=0;for(i=r+1;i<10;i++)sum=sum+A[i][r]*A[i][r];第18页数值分析第二次SY1107214王宽d=sqrt(sum);c=-1*sgn(A[r+1][r])*d;h=c*c-c*A[r+1

5、][r];//step(3)for(i=0;i<10;i++)u[i]=0;for(i=r+2;i<10;i++)u[i]=A[i][r];u[r+1]=A[r+1][r]-c;//step(4)for(i=0;i<10;i++){tmp=0;sum=0;for(j=0;j<10;j++){tmp=tmp+A[j][i]*u[j];sum=sum+A[i][j]*u[j];}p[i]=tmp/h;q[i]=sum/h;}sum=0;for(i=0;i<10;i++)sum=sum+p[i]*u[i];tr=sum/h;for(i=0;i<10;i++)w[i]=q[i]-tr

6、*u[i];for(i=0;i<10;i++){for(j=0;j<10;j++){A[i][j]=A[i][j]-w[i]*u[j]-u[i]*p[j];}}}}//以上的A[][]存在一定的误差，消除误差for(i=0;i<10;i++){for(j=0;j<10;j++)if(-1.0e-12=0)return1;elsereturn-1;}//计算A的特征值的函数（课本P6311步）intte

7、zhengzhi(){intk=1,m=9,n=9;intflag=0;doublems,mt,mk[10][10];//step5doubledetdk,fb,temp,s1_r,s1_i,s2_r,s2_i;voidqmk(doublemk[10][10],intm,doublems,doublemt);//声明求mk的函数voidQR(doublemk[10][10],intm);//声明对mk进行qr分解的函数while(1){if(-1e-12