北航数值分析第一次大作业(幂法反幂法).doc

一、问题分析及算法描述1.问题的提出：（1）用幂法、反幂法求矩阵A=aij2020的按摸最大和最小特征值，并求出相应的特征向量。其中aij=sin0.5i+0.2j ij1.5cosi+1.2j i=j要求：迭代精度达到10-12。（2）用带双步位移的QR法求上述的全部特征值，并求出每一个实特征值相应的特征向量。2.算法的描述：(1)幂法幂法主要用于计算矩阵的按摸为最大的特征值和相应的特征向量。其迭代格式为：任取非零向量u0=h1(0)，hn(0)Thr(k-1)=max1jnhrk-1 yk-1=uk-1hrk-1 uk=Ayk-1=h1k，hnkT k=sgnhrk-1hrk (k=1，2，)终止迭代的控制选用。幂法的使用条件为nn实矩阵A具有n个线性无关的特征向量x1，x2，xn，其相应的特征值1，2，n满足不等式123n或1=2=m1m+1m+2n幂法收敛速度与比值21或m+11有关，比值越小，收敛速度越快。(2)反幂法反幂法用于计算nn实矩阵A按摸最小的特征值，其迭代格式为：任取非零向量u0Rn k-1=uk-1Tuk-1 yk-1=uk-1k-1 Auk=yk-1 k=yk-1uk k=1，2， 每迭代一次都要求解一次线性方程组Auk=yk-1。当k足够大时，n1k，yk-1可近似的作为矩阵A的属于n的特征向量。比值nn-1越小，收敛的越快。反幂法要求矩阵A非奇异。(3)带双步位移的QR分解法QR方法适用于计算一般实矩阵的全部特征值，尤其适用于计算中小型实矩阵的全部特征值。本算例中采用带双步位移的QR方法，可加速收敛，其迭代格式为：A1=A(n-1) Ak-s1(k)I=QkRk 对Ak-s1kI作QR分解 Ak+1=RkQk+s1kI Ak+1-s2kI=Qk+1Rk+1 (对Ak+1-s2kI作QR分解)Ak+2=Rk+1Qk+1+s2kI (k=1，3，5，) 二、计算结果及分析1.计算结果：(1)幂法：初始条件：最大迭代次数L=1000；向量u0=1，0，0计算结果：第1次迭代结果：最大特征值：0.00000e+000第2次迭代结果：最大特征值：2.48910e+000 相对误差：1.00000e+000第3次迭代结果：最大特征值：1.67719e+000 相对误差：4.84085e-001第4次迭代结果：最大特征值：-2.10960e+000 相对误差：1.79503e+000第5次迭代结果：最大特征值：-6.13203e-001 相对误差：2.44030e+000第794次迭代结果：最大特征值：-1.97638e+000 相对误差：4.30074e-011第795次迭代结果：最大特征值：-1.97638e+000 相对误差：3.04354e-013*最终迭代结果*特征值：-1.97638e+000 相对误差：3.04354e-013迭代次数：795(2)反幂法：初始条件：最大迭代次数L=1000；向量u0=1，0，0运行结果：第1次迭代结果：最大特征值：1.07542e+000第2次迭代结果：最大特征值：-3.66550e+000相对误差：1.29339e+000第3次迭代结果：最大特征值：1.22709e+001相对误差：1.29871e+000第4次迭代结果：最大特征值：-1.03421e+000相对误差：1.28650e+001第5次迭代结果：最大特征值：-5.46339e-001相对误差：8.92983e-001第995次迭代结果：最大特征值：-3.24922e-001相对误差：6.61387e-003第996次迭代结果：最大特征值：-3.27176e-001相对误差：6.88964e-003第997次迭代结果：最大特征值：-3.29570e-001相对误差：7.26527e-003第998次迭代结果：最大特征值：-3.32147e-001相对误差：7.75893e-003第999次迭代结果：最大特征值：-3.34960e-001相对误差：8.39573e-003第1000次迭代结果：最大特征值：-3.38073e-001相对误差：9.20985e-003*超过最大设定迭代次数，迭代失败！(3)带双步位移的QR法：初始条件：最大迭代次数L=1000；向量u0=1，0，0运行结果：全部特征值：0.747738+i*0.000000-0.316674+i*0.025741-0.316674+i*-0.0257410.488580+i*0.1393760.488580+i*-0.1393761.045898+i*0.000000-0.630981+i*0.0000000.071427+i*0.5398270.071427+i*-0.5398271.265389+i*0.000000-1.459878+i*0.000000-1.521321+i*0.000000-1.412499+i*0.148349-1.412499+i*-0.148349-1.976376+i*0.0000001.810854+i*0.0000001.362562+i*1.0661171.362562+i*-1.0661170.169262+i*1.9094240.169262+i*-1.909424特征向量（经谱范数归一化）：实特征值0.747738对应特征向量：-0.062705 -0.022368 0.304372 0.064466 0.521833 -0.157024 0.136942 -0.218108 0.250264 -0.043064 -0.228688 -0.184632 -0.072871 0.124721 0.029070 0.102566 -0.136358 0.167727 0.085747 0.546165 实特征值1.045898对应特征向量：-0.018001 0.019652 0.273447 0.070528 0.274896 -0.144015 0.048385 0.376439 -0.583051 -0.054008 -0.168682 -0.113430 -0.034709 0.009204 0.472291 0.125664 -0.190617 0.113145 0.046278 0.059871 实特征值-0.630981对应特征向量：0.106861 0.087709 -0.024967 -0.020897 0.064302 0.034047 0.535143 0.046383 0.028832 0.003479 -0.097276 -0.383801 0.089445 -0.039560 -0.036928 -0.021330 0.014811 0.705836 -0.108904 0.082022 实特征值1.265389对应特征向量：-0.055201 0.003399 0.242191 0.102847 0.372470 -0.372826 0.113953 0.240659 -0.310401 -0.076590 -0.244632 -0.192549 -0.077259 0.263328 0.201662 0.154166 -0.407814 0.186782 0.094649 0.173302 实特征值-1.459878对应特征向量：0.427828 -0.546801 0.007822 -0.382580 0.025199 0.012788 0.033241 0.005389 -0.004065 0.043524 -0.032112 -0.044233 0.135395 -0.006564 0.001214 0.020165 0.011678 0.050001 -0.585765 0.013115 实特征值-1.521321对应特征向量：0.236032 -0.139250 -0.008143 0.638527 -0.009049 -0.002911 -0.001307 0.003054 0.006515 -0.030134 0.012712 0.011368 -0.018792 -0.001753 -0.005749 -0.014290 -0.005292 -0.014591 0.717590 0.001369 实特征值-1.976376对应特征向量：-0.227404 -0.048154 0.022615 0.297305 0.070372 0.039927 0.078503 0.015822 -0.012182 0.605334 -0.083616 -0.106270 -0.573963 -0.019907 0.003839 0.051362 0.036567 0.115613 0.332707 0.036954 实特征值1.810854对应特征向量：-0.027768 -0.051081 -0.159642 -0.054573 -0.084441 0.118378 0.029553 0.211088 0.203867 0.048627 0.075470 0.026824 -0.011103 -0.584846 -0.196009 -0.097668 0.673159 -0.036833 0.003394 0.0812442.结果分析以上三种方法中，幂法计算共进行了795次迭代才达到收敛，计算量较大，收敛性不好；反幂法计算结果未能收敛，通过进一步分析发现，这是因为反幂法迭代程序未考虑按模最小特征值为复数的情况，造成迭代失败。按双步位移的QR法迭代程序则考虑了特征值为复数的情况，并一次性给出了全部特征值。三、结论幂法、反幂法与QR法相比存在以下不足：第一，需要设定初始条件，不适合自动运算；第二，幂法 反幂法的迭代是否收敛依赖于特征值的分布情况，因此实际使用起来很不方便。QR方法的收敛速度与特征值分布无关，是较好的特征值求解数值方法。四、附录（VC+程序代码）(1)幂法：#include#include#define N 20 #define e 1e-12void main()int i,j,k,r,rt;double b,max,beta,betat,er,aNN,yN,ytN;double uN=1;printf(请输入最大迭代次数：n);scanf(%d,&L);for(i=0;iN;i+)for(j=0;jN;j+)if(i=j) aij=1.5*cos(double(i+1)+1.2*double(j+1);else aij=sin(0.5*double(i+1)+0.2*double(j+1);for(k=0;k+)max=u0;r=1;for(i=1;iN;i+)if(fabs(max)fabs(ui)max=ui;r=i;for(i=0;iN;i+)yi=ui/fabs(max);for(i=0;iN;i+)ui=0;for(j=0;j0)/从k=1开始可以计算特征值printf(第%d次迭代结果：n,k);beta=max*ytrt;printf(最大特征值：%1.5en,beta);if(k1)/从k=2开始可以计算相对误差er=fabs(beta-betat)/beta);printf(相对误差：%1.5en,er);if(ere)printf(超过最大设定迭代次数，迭代失败！n);break;for(i=0;iN;i+)yti=yi;rt=r;betat=beta;(2)反幂法：#include#include#define N 20#define e 1e-12 void LU(double cN+1,double b,double x);/LU求解线性方程组函数的声明void main()int i,j,k,L;double er,ita,beta,l,temp,aN+1N+1,yN+1;double uN+1=0,1;for(i=1;iN+1;i+)for(j=1;jN+1;j+)if(i=j) aij=1.5*cos(double(i)+1.2*double(j);else aij=sin(0.5*double(i)+0.2*double(j);printf(请输入最大迭代次数：n);scanf(%d,&L);for(k=0;k+)ita=0;for(i=1;iN+1;i+)ita=ita+pow(ui,2);ita=sqrt(ita);for(i=1;i0)/从k=1开始可以计算特征值beta=0;for(i=1;i1)/从k=2开始可以计算相对误差er=fabs(1/beta-1/temp)*beta);printf(相对误差：%1.5en,er);if(ere)printf(超过最大设定迭代次数，迭代失败！n);break;temp=beta;/LU求解线性方程组函数的定义void LU(double cN+1,double d,double x)double aN+1N+1,bN+1,yN+1,lN+1N+1,uN+1N+1,sN+1,temp,max;int i,j,k,t=0,MN+1;for(i=1;iN+1;i+)for(j=1;jN+1;j+)aij=cij;for(i=1;iN+1;i+)bi=di;for(k=1;kN+1;k+)for(i=k;iN+1;i+) si=aik; for(t=1;tk;t+) si=si-lit*utk;max=0;for(i=k;iN+1;i+)if(maxfabs(si)max=fabs(si);Mk=i;if(Mk!=k)for(t=1;tk;t+)temp=lkt;lkt=lMkt;lMkt=temp;for(t=k;tN+1;t+)temp=akt;akt=aMkt;aMkt=temp;temp=sk;sk=sMk;sMk=temp;ukk=sk;for(j=k+1;kN,jN+1;j+)ukj=akj;for(t=1;tk;t+)ukj=ukj-lkt*utj;for(i=k+1;kN,iN+1;i+)lik=si/ukk;for(k=1;kN;k+)t=Mk;temp=bk;bk=bt;bt=temp;y1=b1;for(i=2;iN+1;i+)yi=bi;for(t=1;t0;i-)xi=yi;for(t=i+1;tN+1;t+)xi=xi-uit*xt;xi=xi/uii;(3)带双步位移的QR法：#include#include#include#define N 20/稀疏矩阵大小#define L 100/最大迭代次数#define e 1e-12/精度水平void Hessenberg(double aN);/拟上三角函数的声明void QR(double bN,double aN,int m);/QR分解函数的声明void LU(double aN-1,double b,double x);/Doolittle分解函数的声明/主函数void main(void)int i,j,k,m,r,n;double det,tri,pd,sum;double vecN,dN-1,reN,imN,lamN;double aNN,bNN,bpNN,cN-1N-1;/系数矩阵for(i=0;iN;i+)for(j=0;jN;j+)if(i=j) aij=bpij=1.5*cos(double(i+1)+1.2*double(j+1);else aij=bpij=sin(0.5*double(i+1)+0.2*double(j+1);Hessenberg(a);/调用拟上三角函数printf(拟上三角化后的矩阵：n);for(i=0;iN;i+)for(j=0;jN;j+)printf( %f,aij);printf(n);k=1;m=N-1;r=0;third:if(fabs(amm-1)=e)rer=amm;imr=0;r+;m-;else goto fifth;fourth:if(m=0)rer=a00;imr=0;r+;goto eleventh;if(m0)goto third;fifth:det=am-1m-1*amm-am-1m*amm-1;tri=am-1m-1+amm;if(m=1|fabs(am-1m-2)=e)pd=pow(tri,2)-4*det;if(pd0)rer=rer+1=tri/2;imr=sqrt(-pd)/2;imr+1=-imr;r+=2;elserer=(tri+sqrt(pd)/2;rer+1=(tri-sqrt(pd)/2;imr+1=imr=0;r+=2;if(m=1)goto eleventh;if(fabs(am-1m-2)=e)m-=2;goto fourth;if(k=L)printf(超过给定最大迭代次数，迭代失败！n);goto end;elsegoto nineth;nineth:for(i=0;i=m;i+)for(j=0;j=m;j+)sum=0;for(n=0;n=m;n+)sum+=ain*anj;bij=sum-tri*aij;for(i=0;i=m;i+)bii+=det;QR(b,a,m);/调用QR分解函数k+;goto third;eleventh:printf(n全部特征值：n);for(i=0;iN;i+)printf(%f+i*%fn,rei,imi);printf(nQR分解后的矩阵：n);for(i=0;iN;i+)for(j=0;jN;j+)printf(%f ,aij);printf(rn);/求特征向量r=0;for(i=0;iN;i+)if(imi=0)lamr=rei;r+;printf(n特征向量（经谱范数归一化）：n);for(k=0;kr;k+)for(i=0;iN-1;i+)for(j=0;jN-1;j+)cij=bpij;for(i=0;iN-1;i+)cii-=lamk;di=-bpiN-1;LU(c,d,vec);vecN-1=1;sum=0;for(i=0;iN;i+)sum+=pow(veci,2);sum=sqrt(sum);printf(实特征值%f对应特征向量：n,lamk);for(i=0;iN;i+)printf(%f ,veci/sum);printf(n);end:;void Hessenberg(double aN)/拟上三角函数的定义 int i,j,r,f; double d,c,h,t,sum; double uN,pN,qN,wN,bpN; for(r=0;rN-2;r+) f=1; for(i=r+2;iN;i+) if(air!=0) f=0;break; if(f=0) sum=0; for(i=r+1;i0) c=-d; else c=d; h=c*(c-ar+1r); for(i=r+2;iN;i+) ui=air; ur+1=ar+1r-c; for(i=r+1;iN;i+) bpi=ui/h; for(j=0;jN;j+)/p sum=0; for(i=r+1;iN;i+) sum+=aij*bpi; pj=sum; for(i=0;iN;i+)/q sum=0; for(j=r+1;jN;j+) sum+=aij*bpj; qi=sum; sum=0; for(i=r+1;iN;i+) sum+=pi*bpi; t=sum;/t for(i=0;i=r;i+)/w wi=qi; for(i=r+1;iN;i+) wi=qi-t*ui; for(i=r+1;iN;i+) for(j=r+1;jN;j+) aij-=ui*pj; for(i=0;iN;i+) for(j=r+1;jN;j+) aij-=wi*uj; ar+1r=c; for(i=r+2;iN;i+) air=0; void QR(double bN,double