数值分析第二次大作业

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

1、-、题目分析使用带双步位移的QR分解法求矩阵A=[6/..]10x10的全部特使用带双步位移的QR分解法求矩阵的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：_pin(0.5/+0.27)(zVj)Uij—tl.5cos(z+1.2j)(/=/)(i,j=l,2,,10)算法设计1.按照题目给出的矩阵定义对矩阵A赋初值：对应的函数为initA()；2.为了减少求特征值和特征向量过程中的计算量，在对矩阵进行QR分解前先进行拟上三角化：对应的函数为nssj()；3.对拟上三角化后的矩阵A使用带双步位移的QR分解法逐次迭代(最大迭代

2、次数L=500),逐个求出其特征值，对应的函数为tezhengzhi()；这个过程中包含着两个子程序：QR()和qmk(),分别用来对矩阵Mk进行QR分解并得到Ak+i和计算mk的值；4.使用带原点平移的反幕法求岀其对应的特征向量，对应的函数为:fmifaO,这个过程中求解线性方程时用到列主元的高斯消元法，对应的函数为：gauss(double);根据数值分析课木的相关知识，步骤3中带双步位移的QR分解法的流程图如下:所使用的编译环境为：VisualC++6.0本次作业#include#includedou

3、bleA[10][10],rr[10][10],qq[10][10]jq[10][10],uk[10][10];doublexy[10],x[10]={0};〃反幕法中doubletzr[10]={0},tzi[10]={0};〃定义矩阵的特征值数组J实部、i虚部〃对矩阵A进行初始化，赋值voidinitA()inti,j;for(i=1;i<=10;i++)for(j=1;j<=10;j++)A[i-l][j-l]=sin(0.5*i+0.2*j);}for(i=1;i<=10;i++)A[i-l][i-l]=cos(i+1.2*i)*1

4、.5;}〃定义函数对矩阵A进行拟上三角化voidnssj(){intrJJ;intsgn(doublea);doubled,c,hXtr,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]==O);讦(t二二8-r)continue;elsesum=0;for(i=r+l;i<10;i++)sum=sum+A[i][r]*A[i][r];d=sqrt(sum);c=-l*sgn(A[r+l][r]

5、)*d;h=c*c-c*A[r+l][r];//step⑶for(i=0;i<10;i卄)u[i]二0;for(i=r+2;i<10;i++)u[i]=A[i][r];u[r+l]=A[r+l][r]-c;//step⑷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

6、;for(i=0;i<10;i++)w[i]=q[i]-tr*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]*pU];}}}}〃以上的A[][]存在一定的误差，消除误差for(i=0;i<10;i++){for(j=0;j<10;j++)if(-1.0e-12=0)return1;elsereturn-1;}〃计算A的特征

7、值的函数(课本P6311步)inttezhengzhi(){intk=l,m=9,n=9;intflag=0;doublems,mt,mk[10][10];//step5doubledetdk,fb,tempzsl_r/slJ,s2_r,s2J;voidqmk(doublemk[10][10],intm,doublems,doublemt);〃声明求mk的函数voidQR(doublemk[10][10],intm);〃声明对mk进行qr分解的函数while(1){if(-le-12

8、tzr[n]=A[m][m];m-;flag=l;}//step4if(flag){flag=0;if(m==0){tzr[n]=A[m][m];return0;}elseif(