1、-矩阵特征值和特征向量Eigenvalues and Eigenvectors-问题的提出问题的提出矩阵特征值计算非常重要，在很多方面应用矩阵特征值计算非常重要，在很多方面应用数值分析中，和矩阵有关的迭代序列的收敛数值分析中，和矩阵有关的迭代序列的收敛取决于迭代矩阵的特征值大小取决于迭代矩阵的特征值大小动态系统中，特征值标志着系统是否是稳定动态系统中，特征值标志着系统是否是稳定的的振动系统中，微分方程的特征值或者有限元振动系统中，微分方程的特征值或者有限元模型的矩阵系数和系统的固有频率直接相关模型的矩阵系数和系统的固有频率直接相关数学中方阵的对角化、微分方程组的解等等数学中方阵的对角化、微分方
2、程组的解等等-6.1 基本概念回顾基本概念回顾DEF6.1 设设A是是n阶方阵，如果数阶方阵，如果数和一维非零向量和一维非零向量使关系式使关系式A=成立，则称数成立，则称数为方阵为方阵A的的特征值特征值，非零向量非零向量称为称为A的属于特征值的属于特征值的的特征向量特征向量.推论：推论：如果如果是矩阵是矩阵A的属于特征值的属于特征值0的特征向量，的特征向量，那么那么的任何一个非零倍数的任何一个非零倍数k也是也是A的属于的属于的特征向的特征向量。这是因为量。这是因为A=0所以所以A(k)=0(k)，这说明属这说明属于同一个特征值的特征向量不是唯一的，但一个特征于同一个特征值的特征向量不是唯一的，
3、但一个特征向量只能属于一个特征值。向量只能属于一个特征值。xAx 可以写成齐次线性方程组可以写成齐次线性方程组0 xEA )(方程组有解方程组有解0EA 即即0aaaaaaaaannn2n12n22121n2111 上式是以上式是以为未知量的一元为未知量的一元n n次方程，称为方阵次方程，称为方阵A A的的特征方程特征方程，EA 是是的的n n次多项式，记为次多项式，记为)(f称为方阵称为方阵A A的的特征多项式特征多项式。显然，方阵显然，方阵A的特征值就是其特征方程的解。特征的特征值就是其特征方程的解。特征方程在复数范围内恒有解，其解的个数为方程的方程在复数范围内恒有解，其解的个数为方程的次
4、数（重跟按重数计算），因此次数（重跟按重数计算），因此n阶方阵有阶方阵有n个特个特征值。显然，征值。显然，n阶单位矩阵阶单位矩阵E的特征值都是的特征值都是1。设设n n阶方阵阶方阵)(ijaA 的特征值为的特征值为n21,则有则有（1 1）;nn2211n21aaa （2 2）.An21 如果如果i 是方阵是方阵A A的一个特征值，的一个特征值，求得非零解求得非零解,ipx 则则ip就是就是A A的对应于特征值的对应于特征值i的特征向量。的特征向量。由以上分析知：由以上分析知：求方阵的特征值和特征向量实际上就是求行列式和求方阵的特征值和特征向量实际上就是求行列式和方程组的解。方程组的解。,)(
5、0 xEAi 程组程组由线性方由线性方例例6.1求矩阵求矩阵 6321A的特征值与特征向量。的特征值与特征向量。解解A A的特征多项式为的特征多项式为),()(76616321 故故A A的特征值为的特征值为.,7021 当当01 时时,由由0 xEA1 )(即方程组即方程组,00 xx632121解得基础解系为解得基础解系为.12p11p就是就是A A的一个属于特征值的一个属于特征值01 的特征向量，的特征向量，A A的属于特征值的属于特征值01 的所有特征向量为的所有特征向量为).0k(kp1为任意常数为任意常数 当当时，时，72 由由0 xEA2 )(即方程组即方程组,00 xx1326
6、21解得基础解系解得基础解系.31p2A A的属于特征值的属于特征值72 的所有特征向量为的所有特征向量为).0k(kp2为任意常数为任意常数 1p就是就是A A的一个属于特征值的一个属于特征值01 的特征向量，的特征向量，对于一阶矩阵对于一阶矩阵A A，如果，如果0是是A A的的k k重特征根，重特征根，个数不大于个数不大于k k，所含向量的个数不大于所含向量的个数不大于k.k.定理定理的线性无关特征向量的的线性无关特征向量的则则A A对应于对应于0 x)EA(0 的基础解系的基础解系也就是说，也就是说，0定理定理 属于不同特征值的特征向量是线性无关的。属于不同特征值的特征向量是线性无关的。
7、事实事实 方阵在复数域内总有特征根，但不一定有实方阵在复数域内总有特征根，但不一定有实特征根。特征根。例例矩阵矩阵 0110A的特征值。的特征值。A A的特征多项式为的特征多项式为.)(111f2 其有复特征根其有复特征根.,ii21 -方程一般形式方程一般形式 0 xIxAxxA 00IAxIA注意：上面用定义阐述了如何求解矩阵注意：上面用定义阐述了如何求解矩阵A A的特征值的特征值和特征向量和特征向量X X。但众所周知，高次多项式求根是。但众所周知，高次多项式求根是相当困难的，而且重根的计算精度较低。同时，相当困难的，而且重根的计算精度较低。同时，矩阵矩阵A A求特征多项式系数的过程对舍入
8、误差十分敏求特征多项式系数的过程对舍入误差十分敏感，这对最后计算结果影响很大。因此，从数值感，这对最后计算结果影响很大。因此，从数值计算角度来看，上述方法缺乏实用价值。计算角度来看，上述方法缺乏实用价值。问题的解决：问题的解决：目前，求矩阵特征值问题实际采用目前，求矩阵特征值问题实际采用的是迭代法和变换法。的是迭代法和变换法。-6.2 幂法（幂法（Power Method）的特征值的特征值计算矩阵计算矩阵例例 1110A2.7618034.02/5161803.12/5101111AI1212 ）（）（：方法方法解解 377233X,233144X53321110XAX32211110XAX2
9、1111110XAX,11X2(12)(11)(2)(3)(1)(2)(0)(1)(0)取取：方方法法在很多问题中，矩阵的按模最大特征值往往起重要在很多问题中，矩阵的按模最大特征值往往起重要的作用。例如矩阵的谱半径即按模最大特征值，决的作用。例如矩阵的谱半径即按模最大特征值，决定了迭代矩阵是否收敛。因此矩阵的按模最大的特定了迭代矩阵是否收敛。因此矩阵的按模最大的特征值比其余特征值更重要。征值比其余特征值更重要。幂法是计算按模最大特征值及相应的特征向量的数幂法是计算按模最大特征值及相应的特征向量的数值方法值方法。简单地说，任取初始向量。简单地说，任取初始向量X(0),迭代计算迭代计算X(k+1)
10、=A X(k)得到迭代序列得到迭代序列X(k+1)，k0,1,；再分析；再分析X(k+1)与与X(k)之间的关系，就可得到之间的关系，就可得到A的按模最大特征值及的按模最大特征值及特征向量的近似解特征向量的近似解-幂法分析幂法分析.n,2,1ixAxxn,;n,2,1i,Aiiiin321i ，即，即特征向量特征向量个线性无关的个线性无关的并有并有其中其中有特征值有特征值在幂法中，假设矩阵在幂法中，假设矩阵nn22110i00 xxxXxnAXX的线性组合，即的线性组合，即的特征向量的特征向量个线性无关个线性无关的的可表示成可表示成，任取初始向量任取初始向量）（）（）（模最大的特征值。模最大的
11、特征值。的变化趋势计算矩阵按的变化趋势计算矩阵按分布有关，幂法根据分布有关，幂法根据的变化趋势与特征值的的变化趋势与特征值的一般地有一般地有）（那么，那么，）（）（）（）（）（）（kknknn2k221k111kknnn222111nn2211nn221101XXxxxAXXxxxxAxAxAxxxAAXX以下考虑两种简单情况。以下考虑两种简单情况。一个一个按模最大的特征值只有按模最大的特征值只有）（）（）（）（故故，有有对充分大的对充分大的，由于，由于若若，由上式得到，由上式得到设设k1111k11k11k1kk1i1i1nk1knn2k1k2211k1nknn2k221k11kn321Xx
12、XxXn32i0kn32i10 xxxxxxX ）（）（）（）（）（）（，得到特征向量近似为，得到特征向量近似为由由，征值征值于是得到按模最大的特于是得到按模最大的特kk1k1kki1ki1XXAXXn21iXX 的大小的大小的速度取决于比值的速度取决于比值收敛于收敛于显然显然）（）（121ki1kiXX和特征向量和特征向量的按模最大的特征值的按模最大的特征值计算矩阵计算矩阵例例 1120A3.7），（，相相应应的的特特征征向向量量值值准准确确近近按按模模最最大大的的特特征征值值的的继继续续算算下下去去，越越来来越越接接），（向向量量，相相应应的的特特征征得得到到按按模模最最大大的的特特征征值
13、值11x254625460 x00073.21111 为为反反号号的的实实根根按按模模最最大大的的特特征征值值是是互互 nk1knnnk1k3322k11k1kn321211xxx)1(xX0，有有，即即，且且设设）（一个一个按模最大的特征值只有按模最大的特征值只有 221k11kk11k111k1k11k)k(i)2k(i1)k(i)2k(i21kk21222k112k12k221k111k11k22k11k1kx)1(2XXx2XXn,2,1iX/XX/XXXx)1(xXx)1(xXx)1(xXk）（）（）（）（）（）（）（）（）（又有又有，呈现有规律的摆动呈现有规律的摆动）（）（）（充分
14、大时充分大时当当 ）（）（）（）（得到相应的特征向量得到相应的特征向量k11k2k11k1XXxXXx从上述过程可得出计算矩阵从上述过程可得出计算矩阵A的按模最大特征值的方的按模最大特征值的方法法,具体步骤如下：具体步骤如下：任取一非零向量任取一非零向量X0,一般可取一般可取 X0=(1,1,1)T X(k+1)=A X(k)当当k足够大时足够大时,即可得到：即可得到：1 X(k+1)/X(k)-6.3 反幂法反幂法(Inverse Power Method)-6.4 规范化幂法规范化幂法 若按若按6.2中计算过程，有一严重缺点，当中计算过程，有一严重缺点，当|1|1时，时，X(k)中不为零的
15、分量将随中不为零的分量将随K的增大的增大而无限增大，计算机就可能出现上溢（或随而无限增大，计算机就可能出现上溢（或随K的增大而很快出现下溢），因此，在实际计算的增大而很快出现下溢），因此，在实际计算时，须按规范法计算，每步先对向量时，须按规范法计算，每步先对向量 进行进行“规规范化范化”，即用，即用X(k)中绝对值最大的一个分量记中绝对值最大的一个分量记作作max|xik|，用，用max|xik|遍除遍除X(k)的所有的所有分量，得到规范化向量分量，得到规范化向量Y(k)，并令，并令X(k+1)=A Y(k)实际计算公式实际计算公式Y(k)X(k)/|X(k)|X(k+1)=A Y(k)的特征
16、值和特征向量的特征值和特征向量的按模最大的按模最大用规范法计算矩阵用规范法计算矩阵例例 1439A4.7-反幂法的规范算法反幂法的规范算法实际计算公式实际计算公式Y(k)X(k)/|X(k)|AX(k+1)=Y(k)-6.5 幂法的加速和降阶幂法的加速和降阶幂法的收敛速率依赖于次大和最大特征值之比，幂法的收敛速率依赖于次大和最大特征值之比，当比值很小时，收敛快当比值很小时，收敛快先对矩阵进行变换，使得有很大的特征值先对矩阵进行变换，使得有很大的特征值原点移位法原点移位法：用：用A0I来代替来代替A进行迭代进行迭代原点移位法原点移位法：A0I和和A的特征值的特征值0，相应的特征向量不变相应的特征
17、向量不变)xxx()(X)IA(Xnk010nn22k010211k010k0k ）（）（为了加速收敛，适当选取为了加速收敛，适当选取0，使得，使得120102 从理论上讲，幂法可以采取从理论上讲，幂法可以采取降阶降阶的方法求出矩阵的方法求出矩阵A A的全部特征值。当求出的全部特征值。当求出1和对应的特征向量和对应的特征向量x1后，后，按同样的思想可以依次求出按同样的思想可以依次求出2,3,n以及相应以及相应的特征向量的特征向量x2,x3,xn 。在幂法中，求出矩阵。在幂法中，求出矩阵A A的主特征值的主特征值1及对应的特征向量及对应的特征向量x1后，可用后，可用压缩压缩方法方法求出求出n-1
18、n-1阶矩阵阶矩阵B B使它的特征值为使它的特征值为2,从而把求从而把求A次特征值次特征值2的问题转化为求的问题转化为求B B的主特征值，等等。的主特征值，等等。幂法小结：幂法小结：幂法适用范围幂法适用范围 为求矩阵的按模最大特征值及相应为求矩阵的按模最大特征值及相应的特征向量，其优点是算法简单，容易编写程序的特征向量，其优点是算法简单，容易编写程序在计算机上实现，缺点是收敛速度慢，其有效性在计算机上实现，缺点是收敛速度慢，其有效性依赖与矩阵特征值的分布情况依赖与矩阵特征值的分布情况反幂法的适用范围是求矩阵反幂法的适用范围是求矩阵A A的按模最小特征值及的按模最小特征值及对应的特征向量。对应的
19、特征向量。-6.6 其它方法其它方法 平行迭代法：可求出前几个较大的特征值和特征平行迭代法：可求出前几个较大的特征值和特征测量，适用于高阶对称稀疏矩阵测量，适用于高阶对称稀疏矩阵 QR算法：基于任何实非奇异矩阵都可分解为正算法：基于任何实非奇异矩阵都可分解为正交矩阵和上三角矩阵的乘积，适用于任意实非奇交矩阵和上三角矩阵的乘积，适用于任意实非奇异矩阵的全部特征值异矩阵的全部特征值 Jacobi法：用平面旋转矩阵构成的正交相似变法：用平面旋转矩阵构成的正交相似变换将对称矩阵化为对角形，适用于实对称矩阵换将对称矩阵化为对角形，适用于实对称矩阵-Power MethodThe basic comput
20、ation of the power method is summarized as lim and 1kk1-k1-kkuAuAuuThe equation can be written as:1-k1-k11-k11-kuAuuAu-Power MethodThe basic computation of the power method is summarized as lim and 1kk1-k1-kkuAuAuuThe equation can be written as:1-k1-k11-k11-kuAuuAu-Shift methodIt is possible to obta
21、in another eigenvalue from the set equations by using a technique known as shifting the matrix.xxASubtract the a vector from each side,thereby changing the maximum eigenvalue xsxIsxA-Shift methodThe eigenvalue,s,is the maximum value of the matrix A.The matrix is rewritten in a form.IABmaxUse the Pow
22、er method to obtain the largest eigenvalue of B.-Inverse Power MethodThe inverse method is similar to the power method,except that it finds the smallest eigenvalue.Using the following technique.xxA xAxAA11 xAx11 xBx-Inverse Power MethodThe algorithm is the same as the Power method and the“eigenvecto
23、r”is not the eigenvector for the smallest eigenvalue.To obtain the smallest eigenvalue from the power method.1 1-Accelerated Power MethodThe Power method can be accelerated by using the Rayleigh Quotient instead of the largest wk value.The Rayeigh Quotient is defined as:11zAzzwz1-Accelerated Power M
24、ethodThe values of the next z term is defined as:The Power method is adapted to use the new value.12wz-QR factorizationAnother form of factorization A=Q*R Produces an orthogonal matrix(“Q”)and a right upper triangular matrix(“R”)Orthogonal matrix-inverse is transposeT1QQ -Why do we care?We can use Q
25、 and R to find eigenvalues1.Get Q and R (A=Q*R)2.Let A=R*Q3.Diagonal elements of A are eigenvalue approximations 4.Iterate until convergedQR factorizationNote:QR eigenvalue method gives all eigenvalues simultaneously,not just the dominant-Householder MatrixHouseholder matrix reduces zk+1,zn to zeroT
26、o achieve the above operation,v must be a linear combination of x and ek 00yyyHxy xxxxxxvv w;ww2IHk21n1kk212 n1kk1k21kTkxxxxxxexv001000e,.,.,-对象对象双曲型方程双曲型方程:0 xuatuLu (5.1)x),x()0,x(u给给定定初初始始条条件件-建立差分格式建立差分格式将将xt平面分割成矩形网格平面分割成矩形网格,2,1,0j,jttt,2,1,0k,khxx0jk 用用(k,j)表示网格节点表示网格节点(xk,tj)，网格节点上的函数，网格节点上的
27、函数值为值为u(k,j)用差商表示导数用差商表示导数)j,x(u2hh)j,k(u)j,1k(uxu)t,k(u2)j,k(u)1j,k(utuxj,ktj,k )j,x(u6hh2)j,1k(u)j,1k(uxu)t,k(u62)1j,k(u)1j,k(utux2j,kt2j,k 方程方程(5.1)式变为式变为称为误差项称为误差项其中其中)h()j,x(u2ah)t,k(u2),h(R0),h(Rh)j,k(u)j,1k(ua)j,k(u)1j,k(uxt11 0huuauuj,kj,1kj,k1j,k (5.2)略去误差项，得到差分方程略去误差项，得到差分方程加上初始条件，构成差分格式加上
28、初始条件，构成差分格式k0,kj,kj,1kj,k1j,ku)uu(aruu -差分格式的收敛性和稳定性差分格式的收敛性和稳定性差分格式的依赖区域差分格式的依赖区域库朗条件：差分格式收敛的必要条件是差分格式的依库朗条件：差分格式收敛的必要条件是差分格式的依赖区域应包含微分方程的依赖区域赖区域应包含微分方程的依赖区域稳定性稳定性-对象对象抛物型方程抛物型方程:Tt0,0b0 xubtuLu22 (5.3)x(g)t,1(u),x(g)t,0(u1x0),x()0,x(u2x),x()0,x(u121 边界条件边界条件）初边值混合问题）初边值混合问题（）初值问题）初值问题（定解条件有两类定解条件有
29、两类-建立差分格式建立差分格式将将xt平面分割成矩形网格平面分割成矩形网格,2,1,0j,jttt,2,1,0k,khxx0jk 用用(k,j)表示网格节点表示网格节点(xk,tj)，网格节点上的函数，网格节点上的函数值为值为u(k,j)用差商表示导数用差商表示导数)j,x(u12hh)j,1k(u)j,k(u2)j,1k(uxu)t,k(u2)j,k(u)1j,k(utux)4(22j,k22tj,k 方程方程(5.3)式变为式变为称为误差项称为误差项其中其中)h()j,x(u12bh)t,k(u2),h(R0),h(Rh)j,1k(u)j,k(u2)j,1k(ub)j,k(u)1j,k(u
30、2)4(x2t112 )uu2u(bsuuj,1kj,kj,1kj,k1j,k (5.4)略去误差项，并令略去误差项，并令s/h2 得到差分方程得到差分方程边界条件差分化（第二、三类边界条件）边界条件差分化（第二、三类边界条件）)t,x(u2hh)t,x(u)t,x(uxu)t,x(u2hh)t,0(u)t,h(uxux1NNt,1xt,0 -常用的差分格式显式格式显式格式 ,2,1,0j),j(gu),j(gu1N,2,1k)kh(u,2,1,0j,1N,2,1k)uu2u(bsuu2j,N1j,00,kj,1kj,kj,1kj,k1j,k隐式格式隐式格式 ,2,1,0j),j(gu),j(
31、gu1N,2,1k)kh(u,2,1,0j,1N,2,1k)uu2u(bsuu2j,N1j,00,kj,1kj,kj,1k1j,kj,kRichardson格式格式 ,2,1,0j),j(gu),j(gu1N,2,1k)kh(u,2,1,0j,1N,2,1k)uu2u(bs2uu2j,N1j,00,kj,1kj,kj,1k1j,k1j,k菱形格式菱形格式 ,2,1,0j),j(gu),j(gu1N,2,1k)kh(u,2,1,0j,1N,2,1k)uuuu(bs2uu2j,N1j,00,kj,1k1j,k1j,kj,1k1j,k1j,k六点格式六点格式 ,2,1,0j),j(gu),j(gu1
32、N,2,1k)kh(u,2,1,0j,1N,2,1k)uu2u(2bs)uu2u(2bsuu2j,N1j,00,kj,1kj,kj,1k1j,1k1j,k1j,1kj,k1j,k-对象对象椭圆型方程椭圆型方程:)y,x(fyuxuuPoisson0yuxuuLaplace22222222 方方程程方方程程(5.5)给给定定的的边边界界条条件件题题，即即在在边边界界上上满满足足主主要要定定解解条条件件是是边边值值问问 )y,x()y,x(u)y,x()y,x(fyuxuu2222-建立差分格式建立差分格式将将xy平面分割成矩形网格平面分割成矩形网格,2,1,0j,jyy,2,1,0k,khxxj
33、k 用用(k,j)表示网格节点表示网格节点(xk,yj)，网格节点上的函数，网格节点上的函数值为值为u(k,j)用差商表示导数用差商表示导数)y,k(u12h)1j,k(u)j,k(u2)1j,k(uyu)j,x(u12hh)j,1k(u)j,k(u2)j,1k(uxux)4(22j,k22x)4(22j,k22 方程方程(5.5)式变为式变为称为误差项称为误差项其中其中)h()j,x(u12bh)y,k(u12b),h(Rf),h(R)1j,k(u)j,k(u2)1j,k(uh)j,1k(u)j,k(u2)j,1k(u22)4(x2)4(y21j,k122 j,h1j,kj,k1j,k2j,1kj,kj,1k2f)uu2u(s1)uu2u(h1 (5.6)略去误差项，得到差分方程略去误差项，得到差分方程边界条件处理边界条件处理直接转移直接转移线性插值线性插值
