关于幂法与反幂法的研究

东北大学秦皇岛分校数值计算课程设计报告幂法及反幂法学院数学与统计学院专业信息与计算科学学号******姓名***指导教师******成绩教师评语：指导教师签字：2014年07月07日数学与统计学院课程设计报告第1页1绪论1.1课题的背景矩阵特征值的数值算法，在科学和工程技术中很多问题在数学上都归结为矩阵的特征值问题。例如，结构的振动波形和频率可分别由适当矩阵的特征向量和特征值来决定，结构的稳定性由特征值决定；又如机械和机件的振动问题，无线电工及光学系统第电磁振荡问题和物理学中各种临界值都牵涉到特征值计算。所以说研究利用数学软件解决求特征值的问题是非常必要的。求矩阵特征值的一种方法是从原始矩阵出发，求出其特征多项式及其根，即得到矩阵的特征值。但高次多项式求根问题尚有困难，而且重根的计算精度较低。另外，原始矩阵求特征多项式系数的过程，对舍入误差非常敏感，对最终结果影响很大。所以，从数值计算的观点来看，这种求矩阵特征值的方法不够好。实际问题中，有时需要的并不是所有的特征根，而是最大最小的实特征根。称模最大的特征根为主特征值。解决特征值计算的算法有很多种，古老的雅可比方法、兰乔斯方法以及较为常用的幂法、QR方法。QR方法是一种变换法，可求全部的特征值；幂法和反幂法是迭代法，只求模最大与模最小的特征值及特征向量。下面主要来研究一下幂法、反幂法，利用MATLAB解决矩阵特征值问题。幂法是一种计算矩阵主特征值（矩阵按模最大的特征值）及对应特征向量的迭代方法，特别适用于大型稀疏矩阵。反幂法是计算海森伯格阵或三对角阵的对应一个给定近似特征值的特征向量的有效方法之一。1.2概念的认识对于n阶矩阵A，若存在数和n维向量x满足：xAx，则称为矩阵A的特征值，x为相应的特征向量。病态矩阵：求解方程组时对数据的小扰动很敏感的矩阵。例如希尔伯特矩阵就是一类著名的病态矩阵。本次课题不对病态矩阵做深入研究。非亏损矩阵：矩阵存在n个线性无关的特征向量，即有一个完全的特征向量组。2MATLAB特征值计算工具简介查阅MATLABHELP可以知道，利用eig函数可以快速求解矩阵的特征值和特征向量。可利用该函数对以下所做的幂法及反幂法程序进行检验。数学与统计学院课程设计报告第2页调用格式为：,()VDeigA,,nobalanceAeigDVAeig说明：(1)输入参量A必须是方阵；(2)输出参量V是一个矩阵，它的各列是方阵A的特征向量；(3)输出方阵D是一个对角阵，其元素是方阵A的特征值，V与D同列向量相对应；(4)当不写输出格式DV,时，只输出由A得特征值为元素的列阵；(5)如果A中含有小到跟截断误差相当的元素时，加写输入参数“nobalance”,它可以提高小元素的作用，通常可以省略该参数，以免使结果误差变大。例题检验：求方阵410131012A的特征值和特征向量。MATLAB实现过程如下：A=[210;131;014];[V,D]=eig(A)V=0.7887-0.57740.2113-0.5774-0.57740.57740.21130.57740.7887D=1.26790003.00000004.73213幂法3.1幂法算法的理论基础及推导设实矩阵nnijaA)(有一个完备的特征向量组（矩阵A有n个线性无关的特征向量），其特征值为n,,21，相应的特征向量为nxxx,,21。已知A的主特征值是实根，且满足条件n321。现在讨论求1及1x的方法：数学与统计学院课程设计报告第3页幂法的基本思想是任取一个非零的初始向量0v，由矩阵A构造一向量序列011021201vAAvvvAAvvAvvkkk称为迭代向量。由假设，0v可表示为nnxxxv22110（设01）（3.1）于是nknnkkkkkxxxvAAvv22211101)(11121111kkinikiikxxx（3.2）其中inikiikx21。由假设),,2,1(1/1nii，故lim0kk,从而111limkkkvx。这说明序列1kkv越来越接近A的对应于1的特征向量，或者说当k充分大时111kkvx，即迭代向量kv为1的特征值的近似向量。下面再考虑主特征值1的计算，用kiv表示kv的第i个分量，则1111111kkiiikikiivxvx，故11limkikkivv。也就是说两相邻的迭代向量分量的比值收敛到主特征值。通过以上推论可以得出结论，设nnRA有n个线性无关的特征向量（即非亏损的），主特征值1满足n321，则对任何非零初始向量v01，构造的向量序列0vAvkk收敛到主特征向量1x；ikikvv1收敛到主特征值1。（定理一）幂法只能对非亏损矩阵求实的主特征值，且常用于实对称矩阵。3.2幂法算法的迭代向量规范化应用幂法计算A的主特征值1及对应的特征向量时，如果11（或11），迭代向量kv的各个不等于零的分量将随k而趋向于无穷（或趋向于零），这样在计算机实现时就可能“溢出”。为了克服这个缺点，就需要将迭代向量加以规范化。设有一向量0v，将其规范化得到向量vvumax，其中vmax表示向量v的绝对值数学与统计学院课程设计报告第4页最大的分量，即如果有iniivv1max0，则0maxivv，且0i为所有绝对值最大的分量中的最小下标。任取一初始向量001v，构造向量序列vmax：0001002022220021200111001maxmaxmaxmaxmaxmaxmaxvAvAuvAvAvvAvAvvuAvvAAuvAvAvvvuAvAuvkkkkkk由（3.1）式有niikiikinikiikxxxvA2111110，niikiikniikiikkkkxxxxvAvAu211112111100maxmax121112111maxmaxxxxxxxiniikiiniikiik同理，可得到niikiikniikiikkxxxxv211111121111max12111121111maxmaxmaxniikiiniikiikxxxxvk结论：设nnRA有n个线性无关的特征向量，主特征值1满足n321，则对任意非零初始向量0100uv，按下述方法构造的向量数学与统计学院课程设计报告第5页序列ku，kv：kkkkkkkvuvAuvuvmax0100,2,1k则有11maxlim1xxukk；1lim2kk。3.3幂法算法1、算法步骤(1)取初始向量0u（例如取Tu1,,1,10）,置精度要求，置1k。(2)计算1kkAuv,kkvmmax,kkkmvu(3)若1kkmm，则停止计算（km作为绝对值最大特征值1，ku作为相应的特征向量）否则置1kk，转(2)。2、实现过程在MATLAB中没有提供现成的函数幂法，通过自定义编写函数实现幂法。其代码如下：function[l,v,s]=mifa(A,x0,eps)%A为已知矩阵%x0为迭代初始向量%eps为迭代精度%l为求得的矩阵主特征值%v为求得的矩阵主特征向量%s为迭代步数ifnargin==2eps=1.0e-6;endv=x0;%v为主特征向量M=5000;%对迭代步数限制m=0;数学与统计学院课程设计报告第6页l=0;fork=1:My=A*v;m=max(y);%m为模的最大分量v=y/m;if(abs(m-l)eps)l=m;%到所需精度，退出，l为主特征值s=k;%s为迭代步数returnelseifk==mdisp('迭代步数太多，收敛速度太慢');l=m;s=M;elsel=m;endendend【注意】：1、当矩阵为病态矩阵时，幂法收敛速度非常慢。2、初始迭代向量不能取零向量。3、例题检验：用幂法计算131716251A的主特征值和相应的特征向量。(1)先利用MATLAB自带函数进行运算：A=[152;6-17;131];[V,D]=eig(A)V=-0.6047-0.76040.4778-0.6981-0.0472-0.8372-0.38340.64770.2659D=数学与统计学院课程设计报告第7页8.0407000-0.3929000-6.6478(2)再使用幂法的MATLAB程序代码实现如下：clearA=[152;6-17;131];x0=[111]';[l,v,s]=mifa(A,x0)l=8.0407v=0.86611.00000.5491s=84(3)将两种结果对比，很明显最大特征值一样，而特征向量不一样，这是由于矩阵的特征向量有无数个，选不同的基底，结果不一样，但是两种向量对应比值相同，即5491.03834.00000.16981.08661.06047.0，因此可以判断这种算法是可行的。3.4幂法扩展由上述讨论已知，幂法是一种求主特征值的方法，试想一下，能不能把幂法改进一下求出所有的特征值？接下来继续讨论：设矩阵A有n个特征值按模从大到小排序为：n321。其相应的n个线性无关的特征向量为nxxx,,21。在计算A的最大特征值1及相应特征向量1x后进行运算，从而得到一个首行全为0的矩阵B，此时A的阶数就被降了一阶，使矩阵得到收缩。在降阶后，继续用乘幂法计算2及其相应的特征向量2x。具体操作如下：收缩法算法：(1)用幂法求出A的主特征值1和主特征向量1x。数学与统计学院课程设计报告第8页(2)令nnnnnnnnnnnTaxaaxaaxaaxaaxaaxaAxAB111212111121212212211212111000其中1`A代表矩阵A的第一行组成的列向量，即nTaaaA112111,,,；为了不失一般性设Tnxxx1211.,,1；(3)去掉1`A的第1行和第1列的1n阶矩阵记为：nnnnnnnnnnnnaxaaxaaxaaxaaxaaxaaxaaxaaxaB1113131212131313313312313212121321231221221则1B与有A除1项外的其余1n相同个特征值n321，可以用幂法计算2及其对应的特征向量2x，如此经过n次收缩，就可以把A的所有特征值求出来。收缩法可用来求矩阵的所有的特征值。它是基于幂法的求解特征值的方法。在MATLAB中实现这个算法，代码如下：functionT=shousuo(A,eps)%A是已知矩阵%eps是精度%T是对角阵，对角元素为矩阵A的特征值ifnargin==1eps=1.0e-6;endN=size(A);