北航-数值分析第二次大作业(带双步位移的QR方法)

一、算法设计方案：按题目要求，本程序运用带双步位移的QR方法求解给定矩阵的特征值，并对每一实特征值，求解其相应的特征向量。总体思路：1）初始化矩阵首先需要将需要求解的矩阵输入程序。为了防止矩阵在后面的计算中被破坏保存A[][]。2）对给定的矩阵进行拟上三角化为了尽量减少计算量，提高程序的运行效率，在对矩阵进行QR分解之前，先进行拟上三角化。由于矩阵的QR分解不改变矩阵的结构，所以具有拟上三角形状的矩阵的QR分解可以减少大量的计算量。这里用函数voidQuasiTriangularization()来实现，函数形参为double型N维方阵doublea[][N]。3）对拟上三角化后的矩阵进行QR分解对拟上三角化的矩阵进行QR分解会大大减小计算量。用子程序voidQR_decomposition()来实现，将Q、R设为形参，然后将计算出来的结果传入Q和R，然后求出RQ乘积。4）对拟上三角化后的矩阵进行带双步位移的QR分解为了加速收敛，对QR分解引入双步位移，适当选取位移量，可以避免进行复数运算。为了进一步减少计算量，在每次进行QR分解之前，先判断是否可以直接得到矩阵的一个特征值或者通过简单的运算得到矩阵的一对特征值。若可以，则得到特征值，同时对矩阵进行降阶处理；若不可以，则进行QR分解。这里用函数intTwoStepDisplacement_QR()来实现。这是用来存储计算得到的特征值的二维数组。考虑到特征值可能为复数，因此将所有特征值均当成复数处理。此函数中，QR分解部分用子函数voidQR_decompositionMk()实现。这里形参有三个，分别用来传递引入双步位移后的Mk阵，A矩阵，以及当前目标矩阵的维数m。5）计算特征向量得到特征值后，计算实特征值相应的特征向量。这里判断所得特征值的虚数部分是否为零。求实特征值的特征向量采用求解相应的方程组((A-λI)x=0)的方法。因此先初始化矩阵Array，计算(A-λI)，再求解方程组。方程组的求解采用列主元的高斯消去法，由函数intGauss_column(doublea[][N],doubleb[],doubleX[])实现。由于此给定矩阵的特殊性，其没有重根，所有对应于每一特征值只有一个特征向量，因此可以用高斯消去法求解此奇异的线性方程组。首先用高斯消去将矩阵(A-λI)化为上三角阵，其最后一行必定全为零。因此在反代的过程中，只要令x[]的最后一个元素为“1”，即可得到方程组的一个基础解系，从而也就是一个特征向量。6）输出有关结果根据题目要求，需要输出拟上三角化后的矩阵、QR分解结束后的矩阵、矩阵全部特征值及对应实特征值的特征向量。由于输出结果要求都要保留12位有效数字，所以将结果输出到文件result.txt中更有利于数据的打印。程序中构造了两个输出函数专门来解决不同输出结果的问题，voidprint_lambda(complexlambda[]);voidprint_matrix(doublemat[][N])。二、程序源代码：#includestdafx.h#includestdlib.h#includemath.h#defineL100//定义最大迭代次数#defineN10//定义矩阵大小#defineerr1e-12//定义误差限//定义一个结构体来表示复数typedefstructcomplex{doublerpart;doubleipart;};FILE*pReadFile;//主函数int_tmain(intargc,_TCHAR*argv[]){inti,j,t;doubley[N],X[N],a[N][N],A[N][N],B[N][N],Q[N][N],R[N][N],RQ[N][N];structcomplexlambda[N];//声明要调用的函数voidQuasiTriangularization(doublea[][N]);voidQR_decomposition(doubleA[][N],doubleQ[][N],doubleR[][N]);voidQR_decompositionMk(doublemk[][N],intm,doublea[][N]);voidprint_lambda(complexlambda[]);voidprint_matrix(doublemat[][N]);voidmulti_matrix(doublea[][N],doubleb[][N],doublec[][N]);intTwoStepDisplacement_QR(doublea[][N],complexlambda[]);intGauss_column(doublea[][N],doubleb[],doubleX[]);//矩阵初始化for(i=0;iN;i++){for(j=0;jN;j++){if(i==j){a[i][j]=1.5*cos((i+1)+1.2*(j+1));A[i][j]=a[i][j];}else{a[i][j]=sin(0.5*(i+1)+0.2*(j+1));A[i][j]=a[i][j];}}}for(i=0;iN;i++){y[i]=0;}//对矩阵a[][]拟上三角化QuasiTriangularization(a);//打开文件result.txtpReadFile=fopen(result.txt,w);//打印结果到文件result.txt中fprintf(pReadFile,拟上三角化之后的矩阵A[%d][%d]:\n,N,N);//printf(拟上三角化之后的矩阵A[%d][%d]:\n,N,N);print_matrix(a);//对拟上三角化后的矩阵a[][]，进行QR分解QR_decomposition(a,Q,R);fprintf(pReadFile,Q[%d][%d]:\n,N,N);//printf(Q[%d][%d]:\n,N,N);print_matrix(Q);fprintf(pReadFile,R[%d][%d]:\n,N,N);//printf(R[%d][%d]:\n,N,N);print_matrix(R);multi_matrix(R,Q,RQ);fprintf(pReadFile,RQ[%d][%d]:\n,N,N);//printf(RQ[%d][%d]:\n,N,N);print_matrix(RQ);//双步位移QR分解求全部特征值TwoStepDisplacement_QR(a,lambda);//利用校正的列主元素高斯消元法求每个实特征值对应的特征向量for(t=0;tN;t++){if(lambda[t].ipart==0){for(i=0;iN;i++){for(j=0;jN;j++){if(i==j)B[i][j]=A[i][j]-lambda[1].rpart;elseB[i][j]=A[i][j];}}Gauss_column(B,y,X);fprintf(pReadFile,\n矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\n,lambda[t].rpart,N);//printf(\n矩阵与特征值λ=%.12e对应的特征向量为X[%d]：\n,lambda[t].rpart,N);for(i=0;iN;i++){fprintf(pReadFile,%.12e,i,X[i]);//printf(X[%d]:%.12e,i,X[i]);}}}fclose(pReadFile);scanf(%d,&i);return0;}//主函数/*************************************矩阵的拟上三角化输入n阶方阵a[][]，将a[][]拟上三角化无返回值**************************************/voidQuasiTriangularization(doublea[][N]){intr,i,j;doubletr,hr,cr,dr,sum;doubleur[N],pr[N],qr[N],wr[N];for(r=0;rN-2;r++){//判断a[i][r](i=r+2,r+3,...,n-2)是否全为零sum=0;//变量sum使用前清零for(i=r+2;iN;i++){sum=sum||a[i][r];}//如果不是全部都是零，则计算if(sum){//计算drsum=0;for(i=r+1;iN;i++){sum+=a[i][r]*a[i][r];}dr=sqrt(sum);//计算crif(a[r+1][r]0)cr=-dr;elsecr=dr;//计算hrhr=cr*cr-cr*a[r+1][r];//取向量ur[]for(i=0;iN;i++){if(ir+1)ur[i]=0;elseif(i==r+1)ur[i]=a[i][r]-cr;elseur[i]=a[i][r];}//计算向量qr[]for(i=0;iN;i++){sum=0;for(j=r+1;jN;j++)sum+=a[i][j]*ur[j];qr[i]=sum/hr;}//计算向量pr[]for(i=0;iN;i++){sum=0;for(j=r+1;jN;j++)sum+=a[j][i]*ur[j];pr[i]=sum/hr;}//计算trsum=0;for(i=r+1;iN;i++)sum+=pr[i]*ur[i];tr=sum/hr;//计算wr[]for(i=0;iN;i++){if(ir+1)wr[i]=qr[i];elsewr[i]=qr[i]-tr*ur[i];}//计算新产生的矩阵a[][]for(i=0;iN;i++)for(j=0;jN;j++)a[i][j]=a[i][j]-wr[i]*ur[j]-ur[i]*pr[j];}}}/*************************************矩阵的QR分解将A[][]分解为Q[][]*R[][]无返回值**************************************/voidQR_decomposition(doubleA[][N],doubleQ[][N],doubleR[][N]){intr,i,j;doubletr,hr,cr,dr,sum;doubleur[N],pr[N],wr[N];//取矩阵R1[][]为A[][]for(i=0;iN;i++)for(j=0;jN;j++)R[i][j]=A[i][j];//取矩阵Q1[][]为单位矩阵for(i=0;iN;i++)for(j=0;jN;j++){if(i==j)Q[i][j]=1;elseQ[i][j]=0;}//for(r=0;rN-1;r++){//判断R[i][r](i=r+1,r+2,...,N-1)是否全为零sum=0;//变量sum使用前清零for(i=r+1;iN;i++){sum=sum||R[i][r];}if(sum){//计算drsum=0;for(i=r;iN;i++){sum+=R[i][r]*R[i][r];}dr=sqrt(sum);//计算crif(R[r][r]0)cr=-dr;elsecr=dr;//计算hrhr=cr*cr-cr*R[r][r];//取向量ur[]for(i=0;iN;i++){if(ir)ur[i]=0;elseif(i==r)ur[i]=R[i][r]-cr;elseur[i]=R[i][r];}//计算向量wr[]for(i=0;iN;i++){sum=0;for(j=r;jN;j++)sum+=Q[i][j]*ur[j];wr[i]=sum;}//计算新的矩阵Qr+1[][],存储在Q[][]里面for(i=0;iN;i++)for(j=0;jN;j++)Q[i][j]=Q[i][j]-wr[i]*ur[j]/hr