数值分析计算实习二

数值分析计算实习二一．算法描述：用QR方法求矩阵的特征值和特征向量1.矩阵的拟上三角化实矩阵可以通过构造household矩阵的方法相似变换化为拟上三角矩阵。为节省计算量，使用教材上的算法。对完成拟上三角化的矩阵进行QR分解能减少计算量。2.带双步位移的QR方法QR方法产生的矩阵序列本质上收敛于分块上三角阵。对角块均是一阶或二阶子块，对每个子块求特征值就得到整个矩阵的所有特征值。为加速收敛，采用带双步位移的QR分解法。记最大迭代次数为Max，取为1000，迭代精度为e=1e-12。每次迭代前检查，若下一个子块满足条件（一阶或二阶），输所有特征值，计算完成。若达到最大迭代次数仍未输出所有特征值，则返回错误信息，计算失败。采用书本上的算法可以节省计算量。3.用高斯消元法解特征向量。对每个实特征值𝜆𝑖，解十阶方程组(A−𝜆𝑖𝐼)𝑋=0,因为系数矩阵非满秩的，统一置特征向量的最后一个值为1。二．C语言源代码：//带双步位移的QR方法#includestdio.h#includestdlib.h#includemath.h#definee1e-12//迭代的精度水平#defineMax1000//迭代的最大次数typedefstruct{doubleRe;doubleIm;}complex;voidcreatematrix(double**a);//创建题中要求的矩阵inthessenberg(double**a);//对矩阵进行拟上三角化intmuiltiply(double**a,double**b,double**c,intm);//计算m阶矩阵的乘法a*b=cintprintmatrix(double**a);//打印矩阵voidQRmethod(double**a);//带双步位移的QR分解法voidcheck(double**a);//将满足精度要求可看作零的元素全部置为0voidoutput(complex*L);//输出特征值和特征向量voidcalculate(double**M,double**a,intm);//执行QR方法中迭代计算的步骤voidgauss(doublelambda);//高斯消元voidmaxline(double**ptr,intn,intk);intmain(){double**a;inti;a=(double**)malloc(10*sizeof(double*));for(i=0;i10;i++)a[i]=(double*)malloc(10*sizeof(double));creatematrix(a);hessenberg(a);check(a);printmatrix(a);//打印已完成上三角化的矩阵QRmethod(a);}voidcreatematrix(double**a){inti,j;for(i=0;i10;i++){for(j=0;j10;j++){if(i!=j)a[i][j]=sin(0.5*(i+1)+0.2*(j+1));elsea[i][j]=1.52*cos(i+1+1.2*(j+1));}}}inthessenberg(double**a)//构造的拟上三角阵存储在原来a的地址中{intr,i,j;doublec,d,h,t,u[10],p[10],q[10],w[10];for(r=0;r8;r++){check(a);c=0;d=0;h=0;for(i=r+2;i10;i++)d+=a[i][r]*a[i][r];if(d==0)continue;else{d+=a[r+1][r]*a[r+1][r];d=pow(d,0.5);if(a[r+1][r]!=0)c=-fabs(a[r+1][r])/a[r+1][r]*d;elsec=d;h=c*c-c*a[r+1][r];for(i=0;ir+1;i++)u[i]=0;u[r+1]=a[r+1][r]-c;for(i=r+2;i10;i++)u[i]=a[i][r];for(i=0;i10;i++){p[i]=0;q[i]=0;for(j=0;j10;j++){p[i]+=a[j][i]*u[j]/h;q[i]+=a[i][j]*u[j]/h;}}t=0;for(i=0;i10;i++)t+=p[i]*u[i]/h;for(i=0;i10;i++)w[i]=q[i]-t*u[i];for(i=0;i10;i++){for(j=0;j10;j++)a[i][j]-=(w[i]*u[j]+u[i]*p[j]);}}}}intmuiltiply(double**a,double**b,double**c,intm){double**ctemp;//临时存储计算结果；inti,j,k;ctemp=(double**)malloc(m*sizeof(double*));for(i=0;im;i++)ctemp[i]=(double*)malloc(m*sizeof(double));for(i=0;im;i++){for(j=0;jm;j++){ctemp[i][j]=0;for(k=0;km;k++)ctemp[i][j]+=a[i][k]*b[k][j];}}for(i=0;im;i++){for(j=0;jm;j++)c[i][j]=ctemp[i][j];}for(i=0;im;i++)free(ctemp[i]);free(ctemp);return0;}intprintmatrix(double**a){inti,j;for(i=0;i10;i++){for(j=0;j10;j++)printf(%.12e,a[i][j]);printf(\n);}printf(\n\n\n);return0;}voidQRmethod(double**a){intk,m,i,j,r;doubles,t,det;double**M;complexL[10];//为Q，R，M分配地址M=(double**)malloc(10*sizeof(double*));for(k=0;k10;k++)M[k]=(double*)malloc(10*sizeof(double));m=9;r=0;for(k=0;kMax;k++){if(m==0){L[r].Re=a[m][m];L[r].Im=0;break;}elseif(m0)break;if(fabs(a[m][m-1])e){L[r].Re=a[m][m];L[r].Im=0;m--;r++;}else{if(m==1){det=(a[m][m]+a[m-1][m-1])*(a[m][m]+a[m-1][m-1])-4*(a[m][m]*a[m-1][m-1]-a[m-1][m]*a[m][m-1]);if(det0){L[r].Re=(a[m][m]+a[m-1][m-1])/2+sqrt(det)/2;L[r].Im=0;L[r+1].Re=(a[m][m]+a[m-1][m-1])/2-sqrt(det)/2;L[r+1].Im=0;}else{L[r].Re=(a[m][m]+a[m-1][m-1])/2;L[r].Im=sqrt(-det)/2;L[r+1].Re=(a[m][m]+a[m-1][m-1])/2;L[r+1].Im=-sqrt(-det)/2;}m-=2;r+=2;continue;}elseif(fabs(a[m-1][m-2])e){det=(a[m][m]+a[m-1][m-1])*(a[m][m]+a[m-1][m-1])-4*(a[m][m]*a[m-1][m-1]-a[m-1][m]*a[m][m-1]);if(det0){L[r].Re=(a[m][m]+a[m-1][m-1])/2+sqrt(det)/2;L[r].Im=0;L[r+1].Re=(a[m][m]+a[m-1][m-1])/2-sqrt(det)/2;L[r+1].Im=0;}else{L[r].Re=(a[m][m]+a[m-1][m-1])/2;L[r].Im=sqrt(-det)/2;L[r+1].Re=(a[m][m]+a[m-1][m-1])/2;L[r+1].Im=-sqrt(-det)/2;}m-=2;r+=2;continue;}else{s=a[m-1][m-1]+a[m][m];t=a[m-1][m-1]*a[m][m]-a[m][m-1]*a[m-1][m];muiltiply(a,a,M,m+1);for(i=0;i10;i++){for(j=0;j10;j++)M[i][j]-=s*a[i][j];M[i][i]+=t;}calculate(M,a,m+1);check(a);}}}check(a);printmatrix(a);output(L);}voidcheck(double**a){inti,j;for(i=0;i10;i++){for(j=0;j10;j++){if(a[i][j]e&&a[i][j]-e)a[i][j]=0;}}}voidoutput(complex*L){intr;for(r=0;r10;r++){if(L[r].Im==0){printf(lambda[%d]=%e\n,r+1,L[r].Re);gauss(L[r].Re);}elseprintf(lambda[%d]=%e+i*%e\n,r+1,L[r].Re,L[r].Im);}}voidcalculate(double**M,double**a,intm){intr,i,j;doublec,d,h,t,u[10],v[10],p[10],q[10],w[10];for(r=0;rm-1;r++){check(M);c=0;d=0;h=0;for(i=r+1;im;i++)d+=M[i][r]*M[i][r];if(fabs(d)==0)continue;else{d+=M[r][r]*M[r][r];d=pow(d,0.5);if(M[r][r]!=0)c=-fabs(M[r][r])/M[r][r]*d;elsec=d;h=c*c-c*M[r][r];for(i=0;ir;i++)u[i]=0;u[r]=M[r][r]-c;for(i=r+1;im;i++)u[i]=M[i][r];for(i=0;im;i++){v[i]=0;for(j=0;jm;j++)v[i]+=M[j][i]*u[j]/h;}for(i=0;im;i++){p[i]=0;q[i]=0;for(j=0;jm;j++){M[i][j]-=u[i]*v[j];p[i]+=a[j][i]*u[j]/h;q[i]+=a[i][j]*u[j]/h;}}t=0;for(i=0;im;i++)t+=p[i]*u[i]/h;for(i=0;im;i++)w[i]=q[i]-t*u[i];for(i=0;im;i++){for(j=0;jm;j++){a[i][j]-=(w[i]*u[j]+u[i]*p[j]);}}}}}voidgauss(doublelambda){double**a;double*X;doublem,sigma;inti,j,k,n=10;X=(double*)malloc(10*sizeof(double));a=(double**)malloc(10*sizeof(double*));for(i=0;in;i++)a[i]=(double*)malloc((11)*sizeof(double));creatematrix(a);for(i=0;in;i++){a[i][i]-=lambda;a[i][10]=0;}for(k=0;kn-1;k++){maxline(a,n,k);for(i=k+1;in;i++){m=a