矩阵的计算方法研究.doc

矩阵的计算方法研究.doc 山东轻工业学院2012届本科生毕业论文 摘 要 自从1946年第一台电子计算机的问世～经过半个多世纪的发展～使得科学工程 计算已成为当今世界最重要的科学进步之一。计算物理学家～诺贝尔奖金获得者 Wilson在二十世纪八十年代指出:科学计算、理论研究、科学实验并列为当今世界科 学活动三种主要方式。许多科学和工程领域如果没有科学计算不可能有一流的研究成 果。 矩阵计算是科学与工程计算的核心～可以毫不夸张地讲～大部分科学与工程问题 都要归结为一个矩阵计算问题～其中具有挑战性的问题是大规模矩阵的计算问题。 在大量的实际问题中～经常会碰到求矩阵特征值和特征向量的问题～这类问题我 们统称为特征值问题。本文首先阐述了此类问题提出的背景以及目前的研究现状～然 后针对常用的方法进行了分析与比较,最后介绍了该领域正在进行的研究工作。 第一章 引言—初步介绍矩阵特征值和矩阵特征向量～并且初步引入矩阵特征值计 算的一些算法。 第二章 详细介绍计算矩阵特征值的常用算法,幂法、Jacobi法、QR法,的基本 原理。 第三章是本文的主要部分～具体比较以上三种常用的计算矩阵特征值的方法～并 且引入具体的例子来加以证明。 第四章 介绍矩阵特征值计算方法这个领域正在进行的研究工作和基本状况。 关键词:特征值 幂法 Jacobi法 QR法 扰动 ABSTRACT More than half a century has seen the science and engineering computing to be one of the important sign of scientific progress, since the 1946, the year the first computer was invented. The science computing、theoretical study and scientific experiment, the calculating physicist Wilson, a Nobel winner, argued in the 1980s, are the major way of scientific activities worldwide. If there was not science computing, a increasing number of science and engineering fields would not have the first-class research achievements. Matrix calculation is the core of the science and engineering computing, and it is not exaggeration to say that most problems of science and engineering computing are attributed to matrix calculation and the most challenging ones are the calculation of large-scale matrix. When a great many practical problems emerge, we would have eigenvalue problems and eigenvectors problems. This study firstly shows the background of these problems as well as theirs research status nowadays; then the author would analyze and compare some commonly used methods; at last the writer would present the situation of this field’s research work that is being done nowadays. 1 山东轻工业学院2012届本科生毕业论文 Chapter1 Introduction- we show matrix eigenvalue and matrix eigenvectors problems and introduce simply some method of matrix eigenvalue calculation by giving basic knowledge of matrix. Chapter2 We would describe the basic theory of matrix eigenvalue including Power method and Jacobi method as well as QR method. Chapter3 Major part, the three commonly used methods would be compared in detail, and we would introduce some examples to exemplify. Chapter4 The study shows the research work and basic situation nowadays of matrix eigenvalue calculation. Key words:matrix eigenvalue ;Power method ; Jacobi method ; QR method ;Disturbance 第一章 引言—矩阵特征值的计算方法的预备知识 1.1前言 随着科学技术的发展和良好性能计算机的日益普及，大规模的计算问题正越来越 多的引起人们的重视。而矩阵在科学计算中起着重要作用，所以，矩阵理论与应用越 来越受到数学学者、工程技术人员和科技人员的关注。矩阵理论不仅是一门重要的数 学理论，而且在数值分析、最优化方法、数学模型等数学分支上有极其重要的应用。 2 山东轻工业学院2012届本科生毕业论文 由于利用矩阵理论与方法来处理错综复杂的工程问题时，具有表达简洁、对工程问题的实质刻画深刻等特点，因此利用矩阵理论方法来处理工程技术的各种问题越来越受到工程界人士的重视。数值代数和矩阵理论与应用已成为众多科学领域的数学工具。矩阵计算的理论和方法与方程组的求解是矩阵理论的一个重要方向，已经成为经济学、生物学、现代物理学等领域处理数学问题不可缺少的强大工具，成为数学计算的一个重要分支。许多实际问题最后常常归结为一个或一些大型系数矩阵为特殊矩阵的线性 [1]方程组的求解问题。 随着科学的发展，矩阵理论已被广泛地运用到应用数学、计算机科学、经济学、工程学、系统科学等诸多方面，成为现代科技领域处理大量有限维形式与数量关系的强有力的工具。对矩阵理论的现代研究与系统工程、优化方法及稳定理论、群论、图论等有着密切的相互关系。作为数学中的一个分支，包含了丰富的内容，成为一门最有实用价值的数学理论。特征值问题是矩阵理论的一个主要研究领域，对它的研究具有重要的理论意义和实用价值。许多科学和工程问题如结构力学中的固有频率分析以及控制系统中的稳定性问题，最终都转化为特征值问题。 在实际应用中，特征值问题有着广泛的应用背景，例如微分方程的刚性比和数值方法的稳定性;动力系统和结构系统的振动问题;电力系统的静态稳定分析;结构振动、量子化学、电路网络、动力系统稳定性、化学反应、宏观经济平衡等，实质上都是矩阵的特征值问题。例如，为了预防地震破坏，建筑物要有固有频率分析报告，美国加州法律规定，每个新设计的建筑，若高于50英尺，都要提供建筑物的固有频率分析报告，防止与地震波共振。这样，许多建筑公司自己或请人计算建筑物的固有频率也促进了矩阵特征值问题计算的研究。 一般的，代数方法只适用阶数较低的矩阵，当矩阵阶数较高时用代数方法求矩阵的特征值和特征向量是极其困难的，就必须使用效率较高的数值方法。 1.2预备知识 在这里我们将讨论这些内容的初等基础理论，它包括两个方面:第一部分是向量，矩阵及运算的定义;第二部分是由向量或矩阵的思想引出的各种概念之间的抽象关系，例如线性相关，列空间等。为了理解后面所提到的算法的描述，必须熟悉矩阵的运算，也要深入掌握矩阵的理论。 [1]nn，n定义1.2.1设，若,,c,0,,,c，使得成立，则称的特A,,,,A,c,是A ,征值，的对应的特征值的特征向量。 ,是A 将上式的形式改写为 ，，，， ,I,A,,0,,,0 这表明特征向量是齐次线性方程组 , ，， ,I,Ax,0 的非零解，由线性方程组的理论知，的特征值的充分必要条件是行列式 ,为A 3 山东轻工业学院2012届本科生毕业论文 ，，det,I-A,0 AA称为矩阵的特征方程，特征方程的根就是的特征值。多项式称，，，，f,,det,I-AAAA为的特征多项式，而矩阵称为的特征矩阵。 ，，,I-A nn[1]-1nnn，n，，，，，，f,,,-trA,，...，-1tedA性质1.2.1设，则，其中为trA,aA,c,An,1i A矩阵的迹。 [1]A性质1.2.2阶矩阵有且仅有n个特征值，其中重特征值的代数重数为m。 nm 然而5次或5次以上的多项式方程一般是没有公式求解的。所以对于阶数较大的矩阵，实际上求特征值是非常困难的，因而就要研究特征值的各种近似求法，在第二章中的幂法和QR算法正是求特征值近似值的基本方法。 [1]，，A,a性质1.2.3设,,,,...,,为的n个特征值(未必互异)，则,ij12nin，n nnAdetA,0，显然当且仅当具有零特征值。 trA,,,detA,,,,ii,1ii,1 __n，n，，，，A,a,cA,A设，用表示的元素的共轭复数为元素组成的矩阵，而ijaij T_HA称为的共轭转置矩阵。 ，，A,A 矩阵的共轭转置运算性质: T_,,HA,A1) ,,,, HHH，，A，B,A，B2) _HH，，3)kA,kA HHH，，AB,BA4) HH5) ，，A,A ,1H,1HA6)如果可逆，则 ，，，，A,A [1]A,性质1.2.4设为阶矩阵的特征值，则 n T,1)也是矩阵的特征值; A _H2)为矩阵的特征值; ,A ,1,1A3)若非奇异，为矩阵的特征值。 A, 第二章 常用的计算矩阵特征值的方法简介 本章主要是介绍常用的计算矩阵特征值的方法，包括幂法和Jacobi法和QR算法这三种方法。这三种方法是计算不同类型矩阵经常会使用的，并且效率较高的方法。 2.1幂法 4 山东轻工业学院2012届本科生毕业论文 矩阵A的特征值是它的特征多项式的根。我们知道阶数超过四次的根式多项式的根一般不能用有限次运算求得，因而矩阵特征值的计算方法本质上总是采用迭代格式。向量迭代法是指不破坏原始矩阵，而是利用进行运算产生一些迭代向量的求解方AA 法。这类方法大多数来求矩阵的部分特征值(常是较大或较小的一部分)和相应的特征向量，由于迭代过程中不破坏原始矩阵A，因而这类方法对高阶稀疏矩阵特别适用。 [2]其中，幂法是计算按模最大的一个或按模最大的几个特征值(称为优势特征值)和相应特征向量的方法。 2.1.1基本原理 [2]幂法主要用于求矩阵的按摸最大的特征值与相应的特征向量的数值方法。它是通过迭代产生序列，由此计算特征值和特征向量。假设A有n个线性无关的特征向量 ,,,,,,...,,x,x,...x，其对应特征值为，且 (2-1) ,,,,...,123n12n12n 其基本思想是:若我们求某个n阶方阵A的特征值和特征向量,先任取一个初始向量，构造如下序列: u,00 u,Au,u,Au,...u,Au... 1021kk,1 由假设可知可用向量组x,x,...x线性表出，设: u012n u,ax，ax，...，ax， 01122nn A这里假设，由于u随即选取，的可能性很小，将上式用左乘得:a,0a,0011 u,Au,a,x，a,x，...，a,x，同理可得: 10111222nnn kn,,,,,jkkkkk,,,,，，u,Au,,ax，a,，...，a,x,,ax，x,,ax，, (2-2),kk-111122nnn11j111k,,,,,,2j1,,,, kn,,,j,,,k,,时,,,其中，因为，故当，则有: ,x,kkj,,,,2j1,, uk ， lim,lim(ax，,),ax11k11k,1 ku,a,xu即当k充分大时有近似等式，从而可以看出k的值越大，越接近特征向111kk ，，uu量的一个倍数。用表示的第i个分量。则由(2-2)式有: ukk1i ,uaa，，，，，，，kk，111，1iii,, 1，，，，，，uaa，,kk11iii u，，k，1ilim,,于是有:。 1k,,，，uki 当k增大时,序列的收敛情况与绝对值最大的特征值有密切关系,分析这一序列的 ，，u，1ki极限,即可求出按模最大的特征值和特征向量。即:当k充分大时，用作为的,1，，uki 5 山东轻工业学院2012届本科生毕业论文 近似值。 2.1.2幂法的基本计算格式 k，nn，,,,kkkki,,,,,,,,1,Au,,， 由式，当k,,时，若的分Auaxaxax10,,0iii111ii,,,,,i,1i,21,,，， k,,1时，Au量会趋于无穷大，而的分量又趋于零。因此，在实际计算时需要做适10 当的规范化处理，以免发生上溢或下溢的现象。 ,,2 幂法的迭: 代 y,Au;kk1, ;m,yk,1,2,,...kk, yk;u,kmk k,,在这种情况下，当时 x1;u,kx 1, m,,k1 幂法对于模相等但特征值不相等的优势特征值不能直接使用，因为这时向量迭代 法不收敛。用中模最大的分量，现在分两种情况进行讨论: max(x)表示x ,,,,,,,,,i,3,4,...,n第一种情况: 121i kn,,,ki,,，，ax1axax，,，,1122ii,,,i,3,,1 u, kkn,,,,,ki,,,,，，maxax1axx，,，,1122i,,,,,i,31,,,, a,0,a,0,当k,,时u如果不收敛，但是 12k k，2n,,,ki,,，，ax-1axax，，,1122ii,,,i,3,,122,,Au， k1kn,,,,,ki,,,,，，maxax，-1ax，ax,1122i,,,,,i,31,,i,, 22k,,j使得u的第ju因为已规范化，则存在某个分量是1，则当时有:，，，Au,,kkk1j ,3且收敛率为。另一方面 ,1 6 山东轻工业学院2012届本科生毕业论文 kn,,,i,,,,,，，2axax，，,111ii1i,,,i,3,,1,，,Auu k1kkn,,,,,ki,,,,，，maxax，-1ax，ax,1122ii,,,,,i,31,,,, kn,,,ki,,,,,，，，，-12axa-x，,222ii1i,,,i,3,,1, Au-u,k1kkn,,,,,ki,,,,，，maxax-1axax，，,1122ii,,,,,i,31,,,, 由此可见，Au，,u和Au-,u可以分别作为,和,的近似特征向量。 k1kk1k12 __ 第二种情况:。 ,,,,,,,1322 A因为是实矩阵，所以复特征值总是成共轭对出现，它们的特征向量也可以取成 n__,i,共轭向量，即x,x，因此，记于是u,ax，ax，ax,,re11,011211ii,3i kn___，，,,,kkik,ik,,2i,，，Aureaxeaxax，如果,和,是二次方程,，b,，c,0,,111,,,011ii1r,,i,3,,，， u，bu，cu,0的根。当k充分大时，可推出，常数的确定可按最小二乘法，b和ck，2k，1k 2,,，，，，，，使u，bu，u取极小值，求出，这是因为b,c后即可求求出,,k，2k，1k1iii b12，，，，Re,-Im,,4c-b,。的实部，虚部，对应于的特,，,,-b,,,,c,,1112121122 u和u征向量可由相邻的向量求出。的特征向量的算法如下: x,kk，111 ，，iau-ykk，1,,a，i,,y,Au设，则。 ,，xu1k，1k1k, 2.2 Jacobi算法 矩阵的相似变换不会改变矩阵的特征值。根据这一原理，我们可以利用一系列的 A特殊相似变换把原矩阵化为易求特征值的特殊矩阵，然后再对其这类特俗形状的矩 A阵求解特征值问题。如任意实对称矩阵总可以通过正交相似变换化为对角型。因此 TR，，ARAR,diag,，，diag,寻找正交矩阵，使得，对角矩阵的特征值就是对角阵的ii R对角元素。的各列就是对应的特征向量。Jacobi于1946年提出了用一系列平面旋转 R来构造矩阵。在正交相似变换下，矩阵元素的平方和保持不变。因此寻找这样的正 A交相似变换，使得对称矩阵经过变换后使得矩阵的非对角线元素的平方和减少，对角线元素的平方和增大，且保持对称性不变。不断地施行这种正交相似变换，最终使非对角线元素的平方和任意接近于零，对角线元素平方和取极大值，这是Jacobi算法 7 山东轻工业学院2012届本科生毕业论文 的基本想法。 2.2.1基本原理 Jacobi法是求实对称矩阵全部特征值的一种有效方法。它不但能求出特征值，而且同时能求出特征向量，虽然和其他方法相比计算量大了一些，但它是一种可靠的方法，如果矩阵有较多的零元素，或者矩阵本身接近一个对角阵，那么利用法更加显示其优越性。 Jacobi方法的基本思想是通过一组平面旋转变换(正交相似变换)将实对称矩阵A n，nA,R化为对角矩阵，得其全部特征值。由代数学知道，若为对称矩阵，则存在一 T，，RAR,diag,,,...,,,D正交矩阵R，使，其中D的对角元素,(i,1,2,...n)就是A12ni T的特征值，的列向量v就是A的对应特征向量。于是求实对称矩阵A的特征值,Rii T问题在于寻找正交矩阵R，使为对角矩阵。 RAR,D [1]下面是Jacobi方法的步骤: T，，，，A,Rp.qARp,q令A,于是，依次构造矩阵序列，使得:，其中A,A，1sss1 ,的一个旋转。在选择平面和确定旋转角度中使用，，，，Rp,q是在p,q平面上转过角度, (s)aA的策略是简单的，通过寻求位于主对角线以上的元素，以确定最大模的项(由于pqs ,考虑到对称性，我们只需要在A的上部的元素中来寻求)。然后，选择旋转角，使得，，sa,0p,q，利用平面旋转的相似变换仅仅影响位于第行和列的元素。 ，，Rp,qpq 被修改的元素由 s，1sss，1a,acos,，asin,,aipipippi s，1sss，1a,,asin,，acos,,a,i,p,qiqipiqqi s，1s2ss2a,acos,，2acos,sin,，asin,pppppqqq s，1s2ss2a,asin,，2acos,sin,，acos,qqpppqqq s，1sss22s，1a,，，，，a,acos,sin,，acos,,sin,,apqqqpppqpq s，1a,,给出。如果把旋转角选择的能消去，那么要求选择角满足:pq sassspq,,,a,0,,-a0,,,,此外，若a,那么选择,。注意，若那么pqppqqs44,，，apq4 sa,0就不需要旋转了。(当然按目前的策略，若那么A已经是对角型了)Jacobi算法pq 产生一个趋向于确定的对角形矩阵的矩阵序列，而这个确定的对角矩阵是与初始矩阵相似的，且这一过程是稳定的，用算法最后能得到一个具有一定精度的对角矩阵，因 TTTT而有，其中R是平面旋转矩阵的乘积，由于，故的各列就是AR,RDRRAR,0 矩阵A的特征向量。 2.2.2关于特征向量的计算 8 山东轻工业学院2012届本科生毕业论文 rJacobi算法可以同时给出特征向量，且不需要增加计算量。设经过次旋转变换 TTT后迭代过程终止，则，若记。因为近似的看作对A,RR...A...RRAP,PAA,11,1rrrrrrrrr 角阵，的对角元素作为A的特征值的近似值。因此的第列就对应于的近似特i,APirr 征向量，且这些特征向量都是正交的。 TTT实际计算时不时保留全部，最后再相乘，而是在计算过程中按照R、R、...、RRi12r TP,PR的特殊结构形成的。记P,I，则，即: -10kkk kk-1，，，，p,p,j,p,q;ijij kk-1k-1，，，，，，,,p,pcos，psin; i,1,2,...,nipipiq ，，，，，，kk-1k-1p,-psin,，pcos,,iqipiq 这样每一次旋转变换只修改P的第p、q列，不必保留R。如果只求特征值，不求特k-1k 征向量，那么上式这一步可以省略了。 2.2.3Jacobi算法的变形 古典Jacobi算法的优点是收敛速度比较快，关于舍入误差有较强的稳定性，因而所求的精度一般都很高且特征向量的正交性好。其缺点是不能有效地利用矩阵的各种特性，如稀疏带状性。计算工作量大，尤其是寻找非对角元素的按模最大值时，花费的机器时间多，对绝对值较小的特征值精度差一些，因而实际计算需要对古典Jacobi算法做一些必要的修正。 [3]第一种是循环Jacobi法，它不需要寻求最大模的元素，而是按照矩阵元素的排列次序把元素消去。通常的排列次序是选择旋转，以便依次消去位于 maxa(1,2),(1,3),…(1,n),(2,3),…,(n-1,n)的元素。其优点在于省去了寻找，缺点是不论ij,ij a的大小如何，均要进行相应的变换。 ij ,,4第二种是过关Jacobi。在算法的每一步，都要设置一个阀值，并进行依次搜法 查，以便找出一个按模超过阀值的元素。然后消去该元素并重复这个过程。由于非对角元素的平方和趋近于零，因此阀值的高度是随着算法的进行而不断的降低，其缺点是阀值不易选取。 2.3 QR方法 A我们知道矩阵在一定条件下可以分解为一个下三角阵和一个上三角阵的乘积，A,LR即，若记，则 A,A1 ,A,LAL,RL;A,LR， 21111111 ,A,LAL,RL; ， A,LR32222222 ？？ A,LRA,RL; ， s,1,2,...ssss,1ss RL其中为单位下三角阵，为上三角矩阵，可以证明: ss 9 山东轻工业学院2012届本科生毕业论文 ， 。 limL,IlimR,limAsss,,,,,,sss 按上面的分解式，计算矩阵的特征值的方法就是LR方法。它是1958年由H.Rutishauser提出的，这种方法的完成总是假定每一步均有可能进行LR分解。LR算 趋向于一个上三角阵，其对角线的元素趋向于的特征值。虽然其法的重要性在于AAs1 计算量小，收敛快，但稳定性差，由于部分选主元必不可少，适用范围小。 1961年，J.G.Francis在LR算法的基础上，提出了QR算法。它把LR算法中的单位下三角阵L用正交矩阵Q代替，这样以来，QR算法就成了目前计算特征值的最有效的方法。 2.3.1基本原理 QR方法是目前计算一般矩阵(中小型矩阵)全部特征值问题的最有效方法之一。有着十分广泛的应用，它有点类似与雅克比方法，通过一系列正交相似变换为容易求出其特征值的矩阵形式，不同的是雅克比直接用正交阵对所给矩阵进行相似变换，而QR方法是通过所给矩阵的正交三角分解和逆向乘法来实现相似变换的。 为非奇异矩阵，则可分解为一正交矩阵Q与上三角矩阵R由代数可知，如果AA T的乘积，即对A进行QR分解:。于是可得到新矩阵B,QR,QAQ。 A,QR 显然，B是由A经过我相似变换得到的，因此B与A的特征值完全相同。再对B进行QR分解，又可得到一新矩阵，重复这一过程可得到矩阵序列，QR算法就是用此 ,,,5QA,QR,Q,RQ序定义的。当然，在实际应用中确定矩阵，使得列ssss，1sss,QA,R。 sss ,A,RQ,QAQ由于，所以矩阵序列A是彼此相邻的，令，进一步A,As1ssssss，1 ,,,,，，，，A,QAQ,QQAQQ,Q,...QAQQ...Q s，1sssss,1s,1s,1ss1s,112s (s)(s)Q,QQ...Q,R,RR...RQQ...QA,AQQ...Q或。定义,则有: 12sss,1112ss，1112s (s)(s)QR,QQ...QQRR...R 12s,1sss,11 ,QQ...QAR...R 12s,1ss,11 ,AQQ...QR...R 112s,1s,11 (s)(s),AQR 1 (s)(s)s最后可得:QR,A, 1 (s)(s)sRQ和R给出了A的U分解可见，,如果我们要求上三角矩阵的对角元素都是j1 正的，则这一分解是唯一的。QR算法有一个非常重要的性质是:如果是三对角矩阵，A1 A{A}那么所有的矩阵都是三对角矩阵，且如果A非奇异，则由QR算法产生的收敛sk于对角阵。 2.3.2特征向量的计算 TA,QAQ设是由正交相似变换得到的上Hessenberg矩阵，若以求得的近似AA111 10 山东轻工业学院2012届本科生毕业论文 特征值，，，利用逆幂法求出的相应的特征向量，再由,i,1,2,...,n,yAiii1 A，求出的相应于的特征向量。 x,Qy,i,1,2,...,n,xiiii 第三章 常用计算矩阵特征值的方法比较及其例证 矩阵特征值的计算问题看似简单，但要获得好的数值结果并不容易。随着科学的进步，需求和挑战促使新的方法、新理论、新想法不断出现，积累了大量的材料。本章在矩阵特征值计算问题方面，从浩瀚的材料中挑选出三种方法，从收敛性分析上做了系统性的比较，并且引入适当的例子，来加以证明。 3.1 三种方法分析和比较 3.1.1幂法的加速 k,[6],,2根据幂法的基本原理知，幂法收敛的快慢显然取决于的主要项收,a,,k2,1,, ,,,222敛于零的速度，即取决于的大小。愈接近1，幂法收敛愈慢;反之愈,,,111接近0，幂法收敛愈快。 ,,61)原点平移加 速 A设,是的特征值且按模的大小排列，即 i ,,,,...,,， 12n ,,p,,p则的特征值为，如果的排列顺序为 A,pIii ,,p,,,p,...,,p， 12n ,,,p22p则适当的选择，使比值远小于比值，从而用幂法计算矩阵的按模A,pI,,p,11 最大特征值时，收敛速度就会显著的提高。 例3.1矩阵A为 ,310，， ,, A,1,3,3,,, ,,0,34，， 用幂法迭代，发现有规律的摆动。我们取，再用幂法求的按模最大特征p,,4A,pI值，迭代六次就得到精确到小数点后三位的结果。 ,,5.1251 p这种方法称为原点平移法，其关键是的选取。若我们能用某种方法估算按模次 p大特征值时，这时的选择最好接近，若对次大特征值的信息无法得到，则应注,,22 意在迭代过程中发现次大特征值。例如在计算的前后两次迭代有规律的摆动就可以判 2pp明接近且符号相反，这时选为接近。也可以用试验的方法，随时修正,和,,121 11 山东轻工业学院2012届本科生毕业论文 值，达到加速为止。 22)加速法 ,, 此法是Aitken1937年提出的，故亦称Aitken加速。此法的基本细想是用相邻三次 [6]迭代的近似向量进行组合来提高其精度，从而使迭代过程加速。设的a,a,a是akk，1k，2 三个近似值，令 ,, ,,,，,,,,，,,,,，,k，1k，1k，2k，2kk 且 ,,k，1k，2 ,,,kk，1显然有 ,,,,--k，1k，2 ,,,,,,--kk，1即得 2,,,-kk，2k，1,,。 ,，,-2,kk，2k，1 22,,,这就是加速技巧。现就说明如何进行加速。 ,-,-12 uk设u，u，是用幂法迭代产生的三个向量，当适当大时使得x,x,...,x可kk，1k，234n 忽略不计，即有如下形式 ,x，x，，12u,, k，，maxx，,x12 kk,,,,,,,,,222,,,,,,,,,0.其中 ,,,,,,,,,111,,,,,, ku设已规范化，当充分大时的最大分量与的最大分量所在位置相同，并x,xxk121 ，，dd,1设对应于该位置的分量是，从而有 x2 x,dx21,,u,x，,,x，;k11k,1，d x,dx21,,u,x，,r,x，; k，111k，1,1，rd x,dx221u,x，,r,,x，,.k，211k，221，rd, ,2r,其中。 ,1 ,,2k，1k，2u,u,u向量具有是上述的形式，且成立。现用法得 ,,,kk，1k，2,,kk，1 12 山东轻工业学院2012届本科生毕业论文 2uu,u，，，，，，kkk，2，1iii 。 ,,,i,1,2,...,ni，，，，，，u，u,2u，2，1kkkiii 2可以证明，而，，，故更接近。 u,x，0,，，,比u,,x，0,xk，21k，211 1.01.00.5,,,,例3.2矩阵，根据幂法可以求得主特征值的结果，迭代第A,1.01.00.25,, ,,0.50.252.0,, 18次时，，迭代第19次时，，迭代第20次时，,,2.5365456,,2.53653741819,,2.5365323，那么可以求得主特征值的近似值20 2,,,,，，，，，，182019,,=2.5364516。 0，，，，，，,，,,2,182019 [6][7]3)Rayleigh商加速技术 A设是对称矩阵，则可选它的特征向量x,x,...,x满足 12n ，，x,x,,,i,j,1,2,...,n. ijij 应用Rayleigh商，则有 22k1，,,Au,u，，,iikkRu ，，,,k22k，，u,u,,,kkii 2k1，n,，，2k122，,,,,,i，,,,,1,1i,1,,i2,，，,2kn,，，2k22,,i,,,，,, ,,,11i,1,,i2,，， 2k,,,2,,,,，01,,,,,1 A易见，当对称时，用u的Rayleigh商改进近似特征值比原来精度提高一倍。 ,k1 0.2510.5,,,,例3.3 矩阵，根据幂法可以求得第10次迭代的结果，特A,10.250.25,, ,,0.50.251.25,, ，，,,1.7865914u,0.7482,0.6497,1征值，和其规范化的特征向量，由Rayleigh商可10k ,Auu，，,kk得，其主特征值的近似值。 ,,,1.7833509，，uu,kk 3.1.2幂法的收缩 A应用幂法可求得的优势特征值，现在讨论如何求次优势特征值。为叙述方便， ,,,,...,,不妨设是互不相同的。这时不能再直接应用没魔法的计算格式，而采用12n A,,,,...,,,收缩技术，即把矩阵降低一阶，使降阶后的矩阵只含有再按幂法的计算23n 13 山东轻工业学院2012届本科生毕业论文 公式求得次优特征值及相应的特征向量。 最通用的压缩方法是以相似变换为基础，我们假定的优势特征值为和相A,A,11 应的特征向量已被算出。现在设为非奇异矩阵，使得 xR11 ， Rx,ke111 1,其中k,0，那么 RARRx,,Rx11111111 T，，r,1,1,1于是，因而A,RAR,。其中是n-1阶矩阵，RARe,,eB,,211111111120B2，， n,1，因而具有同样的特征值，所以具有特征值,,,,...,,，从而可以A与ABr,R23n212对进行计算，以求得 B2 By,,y2222的特征值和相应的特征向量。 y,22 剩下的是求对应的特征向量。设对应于的特征向量，有 A和,u是Ax,122222 T，，,r1u,,u ,,2220B2，，由，可取 By,,y2222 ,，，,u, 2,,y2，， T，，其中是一个数，由,-,a，ry,0确定，最后得到 ,122 -1。 x,Ru212这就是我们所要求的特征向量。 现在我们讨论的选取问题。最简单的方法是选择为一个初等下三角阵和一个RR11 I置换阵的乘积，其中p是的分量中按模最大的位置。于是有 x1p1 y,Ix; Ly,ke1p111其中 1，， ,,y2,1,,y1 L,.,,1？？,,yn,,,？1y1，， ，，k,y,x而。 11p I的引入，有效地确定了一个主元的方法，能使的元素按模不超过1，以便保R1p1证所得出的算法是稳定的。 例3.4假定矩阵 14 山东轻工业学院2012届本科生毕业论文 232，， ,, ,1034A1,, ,,361，， T，，,,11.0,x,0.5,1.0,0.75的优势特征值为对应的特征向量。这里，因而 p,211 1，， ,,， ,,0.51L1,, ,,,0.7501，， 于是 ,1 A,LIAIL21121121 131041，，，，，， ,,,,,,,,0.513220.51,,,,,, ,,,,,,,0.75016310.7501，，，，，， 11104，， ,,,0,30,,, ,,0,4.5,2，， 所以 ,30，， ,B2,,,4.5,2，， 由此可知，其余的特征值为-3和-2。 总之，用幂法求实矩阵的按模最大的特征值和相应特征向量时，因为是事先不知道特征根的分布情况，实际使用时很不方便，特别是不适宜自动计算，收敛速度较慢，其收敛速度依赖于特征值的分布，然而幂法的结构简单，不破坏原始矩阵，利用幂法的基本思想还可以诱导出一些更有效的算法。 3.1.3 经典Jacobi法的收敛性 ，，Ak,1,2,...A按照经典Jacobi算法迭代产生的矩阵序列中，中的零元素在Akk,1k A中可能会变成非零元素。因此，一般地说不能指望通过有限次旋转变换把原矩阵化为对角矩阵，但是我们可以证明下面的收敛定理。 A{A}定理3.1设是实对称矩阵，则由经典Jacobi算法产生的矩阵序列的非对k A角元素趋于零，即趋于对角阵。 k ，，kaA证明 记的元素为，并令 ijk k，，，，A,diaga，E， kijk E其中为对角元素全为零的对称矩阵，由公式 k s，1sss，1a,acos,，asin,,a ipipippi s，1sss，1a,,asin,，acos,,a,i,p,q iqipiqqi 15 山东轻工业学院2012届本科生毕业论文 s，1s2ss2a,acos,，2acos,sin,，asin, pppppqqq s，1s2ss2a,asin,，2acos,sin,，acos, qqpppqqq s，1sss22s，1，，，，a,a,acos,sin,，acos,,sin,,a pqqqpppqpq可得 2kk，，，，a，a，，，，piqi 22，，k22，，k,122,acos，sin，acos，sin ，，，，，，，， ,,,,piqi 22，，k,1，，k,1,，，a，，，a,i,p,qpiqi 从而有 222，，k,1 ，，E,E,2a,kkpq,1FF ，，k,1a是E其中按模最大元素，因此 pqk-1 2222-E,EE kkk-1-12FFF-nn 22,,1- ,E,,k-12F-nn,, k22,,,1,E ,,02Fn,n,, 2故当。定理证毕。 k,,,E,0kF 在定理证明中的上述估计式是很保守的，但是它可以估计迭代次数。我们记 2n,n N,2这时上述的估计式就可写为 kk,1222,,nE,,E,eE1 ,,k00FFFN,, 2当时一定有 k,,NIn, 222， E,,Ek0FF ,t若取，则有迭代次数估计式 ,,2 ，，。 k,2tIn2N,1.39tN 这个迭代次数的估计式是很保守的，实际迭代次数远远小于这个估计得下界。 A{A}定理3.2 设是实对称矩阵，由Jacobi算法产生的矩阵序列的每个对角元素k ，，k，，ai,1,2,...,nAAk,,趋于的一个固定的特征值，即对于的特征值的某个排列，当ij 时，有 ，，A,diag,,,,...,,。 k12n 16 山东轻工业学院2012届本科生毕业论文 证明 首先假定，，是互异的，令 ,i,1,2,..,ni ， ,,min,,,ij,ij ,E,,对，由定理3.1知存在整数，当时，。其次，对任意，kk,kk,k0,,,k00024 ，，kA有的特征值满足，也就是在每个区间，，中只含有，，,,,a,,1,i,n,,,,,，,iiiiii ，，k，，k，1一个对角元a，现在只能让a亦属于区间，，即证明了对任意的m,1，,,,,,，,iiiiii，，k，ma也属于，，。 ,,,,,，,iiii ，，，，kka和ak，1第次迭代仅改变对角元中两个元素，故 ppqq kk，，，，，1，，a,a,,,,,,，,,i,p,q; iiiiii ，，，，K，1k，1，，，，a,,,,,,,，,,,,,,,，, 。 ppqqppqq因此只须证明 ，，k，1，，a,,,,,,，, qqpp 222由公式知 E,,Ek0FF kkkk，，，，，，，，，122a,,,asin,,2sin,cos,a，acos,,, qqppppqqqp k2，，k2，，k，， ，，,a,,sin,，a,,，,,,cos,,2asin,cos,pppqqpqqpq由此 ，，k，12k，，2k2，，k，， a,,,,,,cos,,a,,sin,,a,,cos,,a qqpqppppqqqpq 222,,,,,,,,4cos,sin,cos, ,2,cos2, k，，2a,pqtan21,2,,,因为， ,,，，，，kk，，aa,4ppqq 1所以从而有 cos2,,, 2 k，，a,,,2,。 qqp k，，，，a,,,,,,，,这就证明了。 qqpp As,s,...,s现在假定的特征值中有m个不同的特征值，，，其重数分别为。m,n12m 令 ,,min,,,， ij,ij ,,kk,k根据定理3.1，对任意的正数,，存在整数，有E,对任意的成立。,00k282 ，，s,,,,,，,由圆盘定理知在每个区间内恰有个对角元，应用上面的证明方法亦可iii 证明这些对角元素经过迭代仍属原区间，由于的任意性定理得证。 , 17 山东轻工业学院2012届本科生毕业论文 进一步分析可知，Jacobi算法的收敛性较高，并且对舍入误差有较强的数值稳定 性，因而求得的解精度较高，且所求的解的正交性较好。 3.1.4 QR方法的收敛性分析 定义 若当k,,时，收敛于分块上三角阵(对角块为一阶或二阶子块)，则Ak 称QR算法是收敛的。若当k,,时，趋于分块上三角形(对角块为一阶或二阶子Ak 快，主对角块上方元素可能不收敛)，则称QR算法是基本收敛的。 定理3.3 假定 ,1，，1)，其中 A,X,X,,diag,;i1 ,,,,...,,,02) 12n ,1X,Y,LR3) yy 则QR算法基本收敛。 证明 我们知道当令，则 A,A1 TA,QR,A,QAQ,RQ;111211111 TA,QR,A,QAQ,RQ;222321222 ？？ A,QR,A,RQ k,1,2,...kkkk，1kk由上式可见，只需分析{A}Q的极限情形，这是因为Q决定序列的收敛性质。先假kkk ____k，1A,QRQ,QQ...Q,R,RR...R设A,QR，则，其中，可知k，1k，112k，112k，11k，1k，1k，1k，1k，1 _kkQ是分解的第一个因子，因而需要分析的极限情况。 AAk11 注意到 kk,1 A,X,X 1 k,X,LRyy k,kk，，，，,X,L,,Ryy k,1，，,L,,I，EE,e令。因而的元素 ykkkij k,,,i,,e,l,i,j. ijij,,,j,, 由定理2)知，因而 limE,0kk,, kk,1A,X,X 1 k,X,LRyy k,kk，，，，,X,L,,Ryy可改写成 18 山东轻工业学院2012届本科生毕业论文 kk，，A,XI，E，，,R， 1ky由于X的列线性无关，故可作正交分解，即，且约定的对角元素为正。则X,QRRxxx有 kk，，A,QRI，E，，,R 1xxky k,1,Q，，I，RERR,R， xxkxxy ,1因为当k充分大时，I，RER为非奇异阵，故有唯一的分解式 xkx ,1，，，，kkI，RER,QR， xkx ，，，，kk，，kk,,其中的对角元为正。当时，。因为 Q,I,R,IR kkkk，，，，，，，，A,QQRR,R， 1xsy k妨碍上式称为的QR分解仅仅是第二个因子(上三角阵)的对角元可能非正。为此A1 我们引进两个对角正交阵: ,,,,,n12,,D,diag,,...,;1,,,,,12n,, ,,，，，，RRyy,,nn11D,diag,...,2,,，，，，RRyy11nn,, ，，RR其中为的对角元素。于是 yyii kkkkkk，，,1,，，，，，，A,QQDDDDRR,R xxy12112 ，，k,1,kkkQQDD是的唯一QR分解，且是的正交因子。由QR性质,矩阵序列{A}是AAx21k11 正交相似的，即 __,1,,,111A,QAQ,Q...QAQQ...Q,QAQ kkkkkkkk，111121 _ Q,QQ...Q其中，可以得到 12kk k，，__1，，kk,,k A,DDQQAQQDD112122k，xx k，，_,,11k,k，，,,k,,,, DDQRRQDD1221xx,,,, k，，_k，，,1,1,1k,,A因为时，(上三角阵)，所以的对角线DQR,RQD,R,R,Rkxxxx22 A下方元素收敛于零。故基本收敛于R，定理证完。 k 总之QR算法是目前计算特征值的最有效的方法。 3.2 实例分析 前面我们介绍了几种常用的求特征值的方法，并就其基本理论与收敛稳定性性进 行了分析。为了研究各种方法的优缺点与它们之间的差异，下面我们给出一个具体的 19 山东轻工业学院2012届本科生毕业论文 案例: 422,,,,例3.5:已知矩阵，分别用Jacobi算法，幂法，QR方法计算A的A,251,, ,,216,, 特征值。 解: 1)Jacobi算法 422,,,,，，0解:记 ,251A,, ,,216,, (1)首先取 ，，0 i,1,j,2,a,212 00，，，，a,a11122cot2 a,,,,,，，042a12 2 b,tan,,,a，1,a,,0.780776 1 c,cos,,,0.78820621，b d,sin,,,0.615412 ,,cossin00cd，，，， ,,,,0，，,,,sincos0,,0 ，，VV,,,dc12,,,, ,,,,001001，，，， 0422,0cdcd，，,,，，,,T,,,,，，，，，，，，1000,,,02510 AVAVdcdc,,,,,,,,,,,,001216001，，,,，， 2.43844800.961,,,,，，1,06.5615522.020190 A,, ,,0.9612.0201906,, ，，1i,2,j,3,a,2.020190(2)取 23 11，，，，a,a2233 a,,0.138985 ，，12a23 2 b,tan,,a，1,a,0.88033183 c,0.754208 d,0.656635 与(1)类似，可以得到: 20 山东轻工业学院2012届本科生毕业论文 2.4384480.6310260.724794,,,,，，2 ,0.6310268.3203860A,, ,,00.2096144.241166,,重复上面的过程，如此继续下去，可得 2.1831850.5951920,,,,，，3 ,0.5915928.3203860.209614A,, ,,00.2096144.496424,, 2.1259960,0.020048,,,,，，4,08.3775760.0208653A,, ,,,0.0200480.2086534.496424,, 2.125995,0.001073,0.020019,,,,，，5 ,,0.0010738.3887610A,, ,,,0.02001904.485401,, 2.125825,0.0010720,,,,，，6,,0,0010728.3887610 A,, ,,00.0000094.485401,, 2.12582500,,,,，，7,08.3887610.000009 A,, ,,00.0000094.485401,, ,,2.125825,,,8.388761,,,4.485401故得到A的特征值: 1232)幂法 Mathematica语句 A,{{4,2,2},{2,5,1},{2,1,6}}; EigenvaluesNA,,,, x,{1,1,1}; ,,,,Doy,A.x;Printk,"",y"",x; yx,,,,,MaxAbsy//N,{k,1,20} 运行结果 4 2 2 2 5 1 2 1 6 {8.38762, 4.48646, 2.12592} 1{8, 8, 9}{1, 1, 1} 2{7.33333, 7.22222, 8.66667}{0.888889, 0.888889, 1.} 21 山东轻工业学院2012届本科生毕业论文 3{7.05128, 6.85897, 8.52564}{0.846154, 0.833333, 1.} 4{6.91729, 6.67669, 8.45865}{0.827068, 0.804511, 1.} 5{6.84978, 6.58222, 8.42489}{0.817778, 0.789333, 1.} 6{6.81473, 6.5325, 8.40736}{0.813041, 0.781283, 1.} 7{6.79626, 6.50612, 8.39813}{0.810567, 0.776997, 1.} 8{6.78646, 6.49207, 8.39323}{0.809259, 0.77471, 1.} 9{6.78123, 6.48457, 8.39062}{0.808563, 0.773489, 1.} 10{6.77844, 6.48057, 8.38922}{0.808192, 0.772836, 1.} 11{6.77695, 6.47842, 8.38848}{0.807994, 0.772487, 1.} 12{6.77615, 6.47728, 8.38808}{0.807888, 0.772301, 1.} 13{6.77573, 6.47667, 8.38786}{0.807832, 0.772201, 1.} 14{6.7755, 6.47634, 8.38775}{0.807801, 0.772147, 1.} 15{6.77538, 6.47616, 8.38769}{0.807785, 0.772119, 1.} 16{6.77531, 6.47607, 8.38766}{0.807777, 0.772104, 1.} 17{6.77528, 6.47602, 8.38764}{0.807772, 0.772095, 1.} 18{6.77526, 6.47599, 8.38763}{0.807769, 0.772091, 1.} 19{6.77525, 6.47598, 8.38762}{0.807768, 0.772089, 1.} 20{6.77524, 6.47597, 8.38762}{0.807767, 0.772087, 1.} 由以上程序运行结果可以看出，运用幂法求解此矩阵的特征值时，经过19次循环 后，结果趋于稳定，得到矩阵A的最大特征值为8.38762,与准确值吻合。 3)QR方法: Mathematica语句: Clear A,{{4,2,2},{2,5,1},{2,1,6}} MatrixFrom%,, EigensystemA//MatrixForm//N//Chop,, {q,j},JordanDecompositionA//N,, q//MatrixForm//Chop j//MatrixForm//Chop 运行结果: 4 2 2 2 5 1 2 1 6 8.38762 2.12592 4.48646 {0.807767, 0.772086, 1.} {-2.70384, 1.5336, 1.} {-0.228931, -1.05568, 1.} 22 山东轻工业学院2012届本科生毕业论文 -0.228931 -2.70384 0.807767 1.05568 1.5336 0.772086 1. 1. 1. 4.48646 0 0 0 2.12592 0 0 0 8.38762 在运用QR方法求特征值时，省去了迭代或循环的冗余，能够快捷的计算出矩阵所以的特征值，并且具有较高的精确度。 总结 初级经济法重点总结下载党员个人总结TXt高中句型全总结.doc高中句型全总结.doc理论力学知识点总结pdf :幂法是用于求矩阵的按摸最大的实特征值及相应的特征向量的一种有效的方法，它最大的优点是方法简单，特别适用于大型稀疏矩阵。为了防止“溢出”，通常采用改进的幂法，但有时收敛速度较慢，幂法的收敛速度由第一特征值和第二特征 迭代的收敛速度就越慢。Jacobi算法值比值的绝对值来确定。这个比值越接近于1， 是求实对称矩阵的全部特征值及特征向量的方法，是通过一系列正交相似变换(平面旋转矩阵)把实对称矩阵化为对角阵，从而求出矩阵的全部特征值与特征向量的近似值。在矩阵的特征值问题的数值解法中，QR方法是至今最有效的方法，它利用矩阵正交相似变换把矩阵化为上三角阵，它可求出一般矩阵的全部特征值及特征向量。然后对矩阵进行QR分解。QR方法的特点在于收敛快，精度高。 第四章 矩阵特征值算法的研究现状 在进行特征值计算时，由于初始数据的误差，计算机字长的限制带来的舍入误差是不可避免的，这些误差对矩阵的特征值有什么影响，矩阵特征值对矩阵元素的扰动是否敏感，这些问题是需要加以考虑的，这是一个很复杂的问题，我们称之为矩阵扰动分析。矩阵扰动分析主要是研究矩阵元素的变化对于矩阵问题的解的影响，它不仅仅与矩阵论和算子理论密切相关，而且对于矩阵计算同样是有重要的意义。 矩阵扰动中的矩阵特征值问题不仅可直接解决数学中诸如非线性规划、优化、常微分方程 ,以及各类数学计算问题 ,而且在结构力学、工程设计、计算物理和量子力学中具有重要作用 ,目前矩阵特征值问题的应用大多来自于解数学物理方程、差分方程、Markov过程等。正因为它具有重要意义和广泛的应用 ,所以矩阵特征值扰动问题是具有深度理论意义和广泛应用背景的研究任务之一。 矩阵特征值扰动理论在上个世纪后半叶得到充分的发展，国外的发展体系比较完善，建立了矩阵特征值扰动理论的基本框架，国内在上世纪80年代中期以后，一批致力于基础理论研究的数学工作者在这一领域取得了长足的发展，使矩阵特征值扰动理论的分析方法、研究范围都有突破性的进展，为其在其他学科上的应用起到了导向和借鉴作用。 23 山东轻工业学院2012届本科生毕业论文 参考文献 [1] 朱元国，饶玲，严涛等.矩阵分析与计算[M].国防工业出版社，2010:90-100. [2] 美G.H.戈卢布C.F.范洛恩.矩阵计算[M].科学出版社，2001:359-534 [3] 曹志浩，张玉德，李瑞遐.矩阵计算和方程求根[M].高等教育出版社，1979:108-196 [4] 彭文华.对称矩阵的两特征值问题[J] .大学数学，2004，(03) :59-60. [5] 徐仲，张凯院，陆全，冷国伟.矩阵论简明教程[M].科学出版社，2001:120-137. [6] 蒋尔雄.矩阵计算[M].科学出版社,2005:155-192. [7] 蒋尔雄.关于Rayleigh商迭代的理论基础[J].高等学校计算数学学报，1990，12(1): 78-96 [8] 刑志栋，曹见荣，矩阵数值分析[M].陕西科学技术出版社，2005:172-182. [9] 斯图尔特著.矩阵计算引论[M].上海科技出版社，国防工业出版社,1980:46-66. [10] Reinsch C H.A stable rational QR algorithm for the computation of the eigenvalues of an hermitian tridiagonal matrix[J].Numer Math,1971,25:591-597. 24 山东轻工业学院2012届本科生毕业论文 致 谢 首先感谢山东轻工业学院对我的教育和培养，在这种良好的学术氛围中学到了丰富的知识，锻炼和培养了良好的学术素养。本论文的完成，要感谢我的导师田霞副教授，是她带我走进了学术研究的殿堂!本文从选提到完成都是在田老师的悉心指导下进行的。田老师还提供了大量详实的资料供我参阅，每一篇论文她都认真阅读和批改，我所取得的每一点小小的成绩都凝聚了她无数的心血和汗水!她对学术工作的一丝不苟、锲而不舍的精神，是我们学习的榜样，深深的感染了我。在此，我向她表示深深的感谢! 在此，还要深深的感谢田霞老师学习上给予的关怀和支持!同时，我还要感谢我的父母和家人，理学院的领导和老师在我他们的关爱和支持始终激励着我勤奋刻苦地学习、进取!感谢本科生期间所有给过我帮助的老师、同学和朋友! 最后，感谢百忙之中为评阅论文而辛勤劳动的各位专家教授们! 二零一二年五月于济南 25