1、矩阵的特征值与特征向量的计算摘 要矩阵是高等代数学中的常见工具，也常见于统计分析等应用数学学科中.在物理学中，矩阵于电路学、力学、光学和量子物理中都有应用；计算机科学中，三维动画制作也需要用到矩阵. 矩阵的运算是数值分析领域的重要问题.将矩阵分解为简单矩阵的组合可以在理论和实际应用上简化矩阵的运算.在本论文中，我们主要讨论矩阵的特征值和特征向量的计算，我们知道，有很多现实中的问题都可以用到矩阵特征值与特征向量计算的知识，比如，在物理、力学和工程技术方面有很多的应用，并且发挥着极其重要的作用.因为这些问题都可归结为求矩阵特征值的问题，具体到一些具体问题，如振动问题，物理中某些临界值的确定问题以及
2、一些理论物理中的问题.在本论文中，我们主要介绍求矩阵的特征值与特征向量的一些原理和方法，原理涉及高得代数中矩阵的相关定理，方法主要介绍冥法及反冥法，Jacobi 方法和 QR 算法，并利用 MATLAB，VC 等软件编写相关算法的程序来求解相关问题，加以验证.关键词： 矩阵；特征值；特征向量；冥法；反冥法；Jacobi 方法；QR 算法；VC软件；MATLAB 软件THE CALCULATIONS OF EIGENVALUE AND EIGENVECTOR OF MATRIXABSTRACTThe matrix is an usual tool in Advanced Algebra, whi
3、ch also used by applied mathematics such as Statistics Analysis. In Physics, we can see the important usage of matrix including Electric Circuits, Mechanics, Optics and Quantum Physics. Making three dimension needs matrix in Computer. The arithmetic of matrix is a very important part in Numerical An
4、alysis. It can simplify the calculation of matrix that we decompose the matrix into several simple parts.In this thesis, we mainly talk about the calculation of eigenvalue and eigenvector of matrix. As we all know, there are lots of realistic problems which need the knowledge of the thesis to solve.
5、 We can see the important usage of matrix including Electric Circuits, Mechanics, Optics and Quantum Physics. It play an important role in these problems inferred above. Because these problems can regarded as the calculation of eigenvalue and eigenvector of matrix, like vibrating problems and critic
6、al value problems and so on.We primarily introduce the principle and approach of the calculation of eigenvalue and eigenvector of matrix that infer the relevant principle in Advanced Algebra. We mainly talk about iteration methods, Jacobi method and QR method by using MATLAB.Key words: Matrix；Eigenv
7、alue；Eigenvector；Iteration methods; Jacobi method;QR method;MATLAB目 录1 引言12 相关定理。13 符号说明24 冥法及反冥法24.1冥法.34.2反冥法.85 QR算法.14参考文献.18附录191 引言在本论文中，我们主要讨论矩阵的特征值和特征向量的计算，我们知道，有很多现实中的问题都可以用到矩阵特征值与特征向量计算的知识，比如，在物理、力学和工程技术方面有很多的应用，并且发挥着极其重要的作用.因为这些问题都可归结为求矩阵特征值的问题，具体到一些具体问题，如振动问题，物理中某些临界值的确定问题以及一些理论物理中的问题.在本
8、论文中，我们主要介绍求矩阵的特征值与特征向量的一些原理和方法，原理涉及高得代数中矩阵的相关定理，方法主要介绍冥法及反冥法，Jacobi 方法和 QR 算法，并利用 MATLAB，VC 等软件编写相关算法的程序来求解相关问题，加以验证.2 相关定理定理2.1 如果 是矩阵A的特征值，则有i),.21(ntranii1o.det22nAo定理2.2 设A与B为相似矩阵 ，则)(1ATBA与B有相同的特征值；o1若 是 的一个特征向量，则 是A的特征向量2xx定理2.3 设 ，则A的每一个特征值必属于下述某个圆盘之中：nija)().,21(1niaijji定义2.1 设A是n阶是对称矩阵，对于任意
9、非零向量x，称 为对应于向量xAxR,)()的Rayleigh商.定理2.4 设 为对称矩阵（其特征值次序记作 ，对应的特征向量nRn21组成规范化正交组，即 ），则nx,21 ijjx),(（对于任何非零向量x）；1),(Ano;),(ma21oxRno.),(i30AxRnno3 符号说明A:n阶矩阵B:n阶矩阵I：n阶单位阵:矩阵特征值i),.21(x:实数域上的n维向量：实数域上的n维向量),0(iv：实属上的规范化向量1ku4 冥法及反冥法4.1 冥法幂法是一种计算矩阵 的主特征值的一种迭代法，它最大优点是方法简单，AnR适合于计算大型稀疏矩阵的主特征值.设 ，其特征值为 ，对应特征
10、向量为 即nRaijA)( i),1(nixiix),1且 线性无关.设 特征值满足：（即 为强占优）,ni A1（4.1.1）|21n幂法的基本思想，是任取一个非零初始向量 ，由矩阵 的乘幂构造一向量nRv0A序列(4.1.2)0112201vAvkk称 为迭代向量.下面来分折 .关 系与及 1kvx由设 为 中一个基本，于是， 有展开式 ,1n R0v niixav10（且设 ）01且有ikniKkxvAv101)()(1211 nkkkx1kka由假设（4.1.1）式,则 ),2(0)1niimk即 lik(4.1.3 )且收敛速度由比值 确定.且有|12r1limxvk这说明，当 充分
11、大时，有 ，或 越来越接近特征向量 .1/xvkkv1/1x下面考虑主特征值 的计算.1用 表示 的第 个分量，考虑相邻迭代向量的分量的比值.ikv)(ki)0(,)()(11 ikikiik vx设从而是(4.1.5)11)(limikkv说明相邻迭代向量分量的比值收敛到主特征 ,且收敛速度由比值 来度量, 越1|12rr小收敛越快,但 越小收敛越快,但 ,而接近于 1 时,收敛可能很慢.r|12r定理 4.1 （1）设 n 个线性无关的特征向量：RA（2）设 特征值满足 |21n（3）幂法： )0(10且v）kA,2则 （1） ；11)(limxvik（2） 1)(liik如果 主特征值为
12、实的重根，即有A| 121 nr又设 A 有 个线性无关的特征向量， 其中n ,21nx（41.4）),1(),1( nrixArixAii 对于任意初始向量 ),(10不 全 为 零且 riniixv则由幂法有 rinriikikKk xxvA110 )()(1kiik且有（设 不全为零）,lim1irikxvr,1)0(kk当由此，当 充分大时， 接近于与 对应的特征向量的某个线性组合.kv1/1应用幂法计算 的主特征值 及对应的特征向量时，如果 ） ，迭代A11(或向量的各个不等于零的分量将随 而趋于无究（或趋于零） ，这样电算时就可能溢k出.为此，就南非要将迭代向量加以规范化.设有非零
13、向量 )()max2等或归 范 化 vuvuv其中 表示向量 绝对值最大的元素，即如果有草药 则)ax( ,)(max)(0iviv0iv其中 为所有绝对值最大的分量中最小指标.0i显然有下面性性质：设 ，则nrvt,为 实 数)max()(t在定理 4.1 条件下幂法可改进为：n1任取初始向量 .)0(10且vu迭代： 规范化：， 01Av)max(11v)(0Av,)max(0212u)()(0222uk ,)ax(011vAuvkk )max()(00vAvukkk 于是，由上式产生迭代向量序列 及规范化向量kku且改进幂法计算公式为：设 )0(10且vu对于 ,2k（4.1.7） :规
14、 范 化迭 代 kkkkvAu/)max(1下面考查 与计算 的关系.,u1x及由 niixv10且有 （4.1.8） )(110kkniikA其中 )(0)(1xiknik当(1)考查规范化向量序列:由(4.1.7)及(4.1.8)式,则有 )(max1kkku（4.16）)()max)ax(11 kk当(2)考查迭代向量序列: )(ax)ax( 1101 kkkKkvAv)m(111k)(,ax()1kvkkk 当定理 (改进幂法)(1)设 有 个线性无关特征向量;nRA(2)设 特征值满足n21且 ),1(ixAi (3) 由改进幂法得到(4.1.7)式),则有,kvu(a) )max(
15、li1k(b) 1lili kkkv且收敛速度由比值 确定.|12r实现幂法，每迭代一次主要是计算一次矩阵乘向量 ，可编一个子程序.)(Au例 1.用 MATLAB 编写冥法程序求矩阵主特征值及近似主特征量用幂法计算下列矩阵的主特征值和对应的特征向量的近似向量，精度 .并510把输出的结果真实结果进行比较. 6312B解 输入 MATLAB 程序 B=1 2 3;2 1 3;3 3 6; V0=1,1,1; k,lambda,Vk,Wc=mifa(B,V0,0.00001,100), V,D = eig (B), 于是， Dzd=max(diag(D), wuD= abs(Dzd- lambd
16、a), wuV=V(:,3)./Vk,运行后屏幕显示结果请注意：迭代次数 k,主特征值的近似值 lambda,主特征向量的近似向量Vk,相邻两次迭代的误差 Wc 如下：k = lambda = Wc = Dzd = wuD =3 9 0 9 0Vk = wuV =0.50000000000000 0.816496580927730.50000000000000 0.816496580927731.00000000000000 0.81649658092773V =0.70710678118655 0.57735026918963 0.40824829046386-0.7071067811865
17、5 0.57735026918963 0.408248290463860 -0.57735026918963 0.816496580927734.2 反冥法及位移反冥法(1)反幂法可用来计算矩阵按模最小的特征值及对应的特征向量.设 为非厅异矩阵， 特征值满足nRAA021n对应特征向量 为线性无关，则 特征求值为x,1 1|21n特征向量为 .,x因此计算 的按模最小的特征值 的部题就是计算 按模最大的特征值部题.An1A对于 应用幂法迭代（称为反幂法） ，可求矩阵 的主特征值 .1 n/反幂法迭代公式：任取初始向量 ，)0(0nuv且1，2，k（4.2.1）kkkkuvA/)max(1其中迭
18、代向量 可通过解方程组求得:kv如果 个线性无关特征向量且 特征值满足:nRA有A011n则由反幂法(2.11)构造的向量序列 满足,kvunknkbxua1lim)()a(且收敛速度由比值 确定.|1nr（2）应用反幂法求一个的似特征值对应的特征向量.设已知 的特征值 的一个近似值 （通常是用其它方法得到） ，现要求对nRAjj应的特征向量 （近似） ，在反幂法中也可用原点平移法来加速收敛.jx如果 存在，显然，特征值为1)(pIpn,21对应的特征向量 .x,21现取 （但不能取 ） ，且设 与其它特征值是分离的，即pj jj|j)(|,ii1kuv即 |1pj)(,jii说明 是 的主特
19、征值.j 1)(IA现对 应用幂法得到反幂法计算公式：1p取初始向量 ),0(0jvu且,21k（4.2.2）kkkuvpIA/)max(1与定理 8 证明类似,可得下述结果.定理 10 （1）设 有 个线性无关特征向量即 .nRA ),1(nixAi（2）取 （为 特征值 一个近似值） ，设 存在且pjj1)(pI|j)(|ii则由反幂法迭代公式（2，12）构造向量序列 满足：,kvu)(1)max()( kpjkvbjku当当或 jkp)(当且收敛速度由比值|minprijj确定.由定理可知，反幂法计算公式（4.2.2）可用计算特征向量 .选择 是 的一个jxpj近似且 的特征值分离情况较
20、好，一般 很小，所以迭代过程收敛较快，同时改进特征Ar值.反幂法迭代公式中 是以通过解方程组kv1)(kupIA求得.为了节省计算量，可先将 进行三角分解.)(pIALUIP)(其中 为置换阵，于是每次迭代求 相当于求解两个三角形方程组pkvkyvu1可按下述方法取 ，即选 使0v0uTPuLUv)1,(01回代求解即求得 .反幂法计算公式：1分解计算，且保存 及 信息LUpIAP)(,P2反幂法迭代（1） 求Tv)1,(1v)max(/1u（2） ,k1） 求1kPLyky求Uv2） )max(kk3） vu/对于计算对称三对角阵，或计算 Hessenberg 阵对应于一个给定的近似特征值的
21、特征向量，反幂 法是一个有效方法.例 2 .用反幂法计算 对应于近似特征值 （精确特征值为A26.13）的特征向量2679413.33410A解 取 ，用部分选主元分解法实现 ，26.p LUpIAP)(其中 ,12876.04.1L,104821U01p（1）求解 TUv)1,(1)291.0,5486.,720.Tu)(7.13.,1（2）求解 2PLUvT)(7034.,89.,763.153Tu)21.0,.,（3） 求解特征向量（真解）是3Tx)(32,1,06794,58.,由此， 相当好的近似 .是2u3x 267943.12.133/例3.用原点位移反幂法的迭代公式，根据给定的
22、下列矩阵的特征值 的初始值 ，nn计算与 对应的特征向量 的近似向量,精确到0.000 1.nnX,2104.0解 输入 MATLAB 程序 A=1 -1 0;-2 4 -2;0 -1 2;V0=1,1,1;k,lambda,Vk,Wc=ydwyfmf(A,V0,0.2,0.0001,10000)运行后屏幕显示结果请注意：因为 A-aE 的各阶主子式都不等于零，所以 A-aE 能进行 LU 分解.A-aE 的秩 R(A-aE)和各阶顺序主子式值 hl、迭代次数 k,按模最小特征值的近似值lambda,特征向量的近似向量 Vk,相邻两次迭代的误差 Wc 如下：k = lambda = Wc =
23、hl =3 0.2384 1.0213e-007 0.8000 1.0400 0.2720Vk = V = D =1.0000 -0.2424 -1.0000 -0.5707 5.1249 0 00.7616 1.0000 -0.7616 0.3633 0 0.2384 00.4323 -0.3200 -0.4323 1.0000 0 0 1.6367例 4. 用原点位移反幂法的迭代公式（5.28） ，计算 的分别对应于126475A特征值 的特征向量 的近似向量，相邻迭代误差为 0.001.将计算结果与 1.01X精确特征向量比较.解 计算特征值 对应的特征向量 的近似向量.输入 MATLA
24、B 程序 .01 1X A=0 11 -5;-2 17 -7;-4 26 -10;V0=1,1,1;k,lambda,Vk,Wc= ydwyfmf(A,V0,1.001, 0.001,100),V,D=eig(A);Dzd=min(diag(D), wuD= abs(Dzd- lambda),VD=V(:,1),wuV=V(:,1)./Vk,运行后屏幕显示结果请注意：因为 A-aE 的各阶主子式都不等于零，所以 A-aE 能进行 LU 分解.A-aE 的秩 R(A-aE)和各阶顺序主子式值 hl、迭代次数 k,按模最小特征值的近似值lambda,特征向量的近似向量 Vk,相邻两次迭代的误差 W
25、c 如下：hl =-1.00100000000000 5.98500100000000 -0.00299600100000k = lambda = RA1 = 5 1.00200000000000 -0.00299600100000Vk = VD = wuV =-0.50000000000000 -0.40824829046386 0.81649658092773-0.50000000000000 -0.40824829046386 0.81649658092773-1.00000000000000 -0.81649658092773 0.81649658092773Wc = Dzd = wu
26、D =1.378794763695562e-009 1.00000000000000 0.00200000000000 从输出的结果可见，迭代 5 次，特征向量 的近似向量 的相邻两次迭代的误1X1差 Wc 1.379 e-009，由 wuV 可以看出， = Vk 与 VD 的对应分量的比值相等.特征值的近似值 lambda 1.002 与初始值 1.001 的绝对误差为 0.001，而与 的绝对误11 1差为 0.002，其中,1X T)01.0 ,0.5- ,00.5( .5 QR 方法用最末元位移QR方法求下列实对称矩阵的全部近似特征值，并将计算结果与 全部真A特征值比较.其中精度为,2
27、 1 3 4- 5A510解 首先保存用最末元位移 QR 方法求实对称矩阵 全部特征值的 MATLAB 主程序A为 M 文件，取名为 qr4.m.在 MATLAB 工作窗口输入程序 A=5 2 2 1;2 -4 1 1;2 1 3 1;1 1 1 2; tzg=qr4(A,5,100)运行后屏幕显示结果请注意：下面的 i 表示求第 i 个特征值，k 是迭代次数，sk 是原点位移量，Bk=Ak-sk*eye(n),Qk 和 Rk 是 Bk 的 QR 分解，At=Rk*Qk+sk*eye(n)迭代矩阵：i =3tzgk =-4.54918876933872k =5sk =-4.5491887662
28、1637B =7.16946633623142tzgk =7.16946633623142请注意：n 阶实对称矩阵 A 的全部真特征值 lamoda 和至少含 t 个有效数字的近似特征值 tzg 如下：lamoda =-4.549188769343661.379722432931842.000000000000007.16946633641181tzg =-4.549188769338721.379722433107302.000000000000007.16946633623142用求根位移 QR 方法求实对称矩阵 全部特征值，精度为 .并将计算结果与A410全部真特征值比较.其中A2135
29、解 首先把用求根位移 QR 方法求实对称矩阵 全部特征值的 MATLAB 主程序保存为AM 文件，命名为 qr8.m.然后在工作窗口输入 MATLAB 程序 A=5 2 2 1;2 -3 1 1;2 1 3 1;1 1 1 2; tzg=qr8(A,0.0001,100)运行后屏幕显示结果如下：请注意：下面的 i 表示求第 i 个特征值，如果迭代矩阵 Ak 的阶数2,且 m 阶矩阵 Ak的 m 行第 m-1 列的元近似等于零 .则原 n 阶矩阵 A 的第 j 个特征值 j=skj，j=1,2,.,n-2;下面的矩阵 Ak 降一阶.i =3tzgk =2.0004Ak =4.8235 -2.02
30、82-2.0282 -5.2085如果迭代矩阵 Ak 的阶数=2,则原 n 阶矩阵 A 的最后两个特征值 j=k+xj，k=n-2,j=1,2.x1 =5.2180x2 =-5.6030tzg1 =7.2183tzg2 =-3.6027请注意：n 阶实对称矩阵 A 的全部真特征值 lamoda 和精度为 jd 的近似特征值 tzg 如下：lamoda =-3.60271.38432.00007.2183tzg =-3.60271.38432.00047.2183参 考 文 献1 姜启源，谢金星，叶俊编数学模型（第三版）M北京：高等教育出版社，2005：1-202.2 王建卫，曲中水 凌滨编著.
31、 MATLAB 7.X 程序设计M . 北京：中国水利水电出版社，2007：55-80.3 李庆扬，王能超，易大义编著.数值分析（第四版）M. 武汉：华中科技大学出版社，2006：219-245.4 王萼芳编著 .高等代数M. 北京：高等教育出版社,2009.161-210.5 张军编著 .数值计算M .北京：清华大学出版社，2008.7.1-200.6 莫勒编著 .MATLAB 数值计算 M.北京：机械工业出版社，2006.6.1-150.附 录程序1： 用幂法计算矩阵 的主特征值和对应的特征向量的MATLAB主程序Afunction k,lambda,Vk,Wc=mifa(A,V0,jd,
32、max1)lambda=0;k=1;Wc =1; ,jd=jd*0.1;state=1; V=V0;while(kjd)state=1;endk=k+1;Wc=Wc;endif(Wcjd)state=1;endk=k+1;%Vk=Vk2,mk=mk1,endif(Wcjd)state=1;endk=k+1;endif(Wc1)b1=abs(Ak(n,n-1); b2=abs(Ak(n,n); b3=abs(Ak(n-1,n-1);b4=min(b2, b3); jd=10(-t); jd1=jd*b4;if(b1=jd1)sk=Ak(n,n); Bk=Ak-sk*eye(n); Qk,Rk=q
33、r(Bk);At=Rk*Qk+sk*eye(n); k=k+1;tzgk=Ak(n,n);disp(请注意：下面的i表示求第i 个特征值，k是迭代次数， sk是原点位移量，)disp( Bk=Ak-sk*eye(n),Qk和Rk是Bk 的QR分解，At=Rk*Qk+sk*eye(n)迭代矩阵： )i,state=1;k,sk,Bk,Qk,Rk,At,Ak=At;elsedisp(请注意：i表示求第i 个特征值，tzgk是矩阵A的特征值的近似值，k是迭代次数， )disp( 下面的矩阵B是m 阶矩阵At 的（m-1）阶主子矩阵，即At 降一阶.)i,tzgk=Ak(n,n),tzg(n,1)=t
34、zgk;k=k,sk,Ak;B=Ak(1:n-1,1:n-1),Ak=B;n=n-1;state=1; i=i+1;endendendtzg(1,1)=Ak;tzg=sort(tzg(:,1);tzgk=Akdisp(请注意：n阶实对称矩阵A 的全部真特征值lamoda和至少含t 个有效数字的近似特征值tzg如下： )lamoda=sort(eig(A)程序 5.用求根位移 QR 方法求实对称矩阵 全部特征值的 MATLAB 主程序Afunction tzg=qr8(A,jd,max1)n,n=size(A); Ak=A; k=0; tzg=zeros(n); state=1;i=1;s0=0
35、;while(k2)bn=abs(Ak(n,n-1);if(bn=jd)S=eig(Ak(n-1:n,n-1:n);sk=s0; j,t=min(abs(Ak(n,n)*1,1-S);t=t;sk=S(t); Bk= Ak-sk*eye(n); Qk,Rk=qr(Bk); Ak=Rk*Qk;k=k+1;tzgk=sk+s0;s0=tzgk;disp(请注意：下面的i表示求第i 个特征值，k是迭代次数， sk是原点位移量，Ak迭代矩阵： )i,state=1;k,sk,tzgk;Ak,elsedisp(请注意：下面的i表示求第i 个特征值，如果迭代矩阵 Ak的阶数2, 且m 阶矩阵Ak的m 行第
36、m-1列的元近似等于零. 则原n 阶矩阵A 的第j个特征值j=skj，j=1,2,.,n-2; 下面的矩阵Ak 降一阶. )i=i+1, k;sk;tzgk,tzg(n,1)=tzgk;Ak;B=Ak(1:n-1,1:n-1); Ak=B,n=n-1;state=1;endenddisp(如果迭代矩阵Ak的阶数 =2,则原n阶矩阵A的最后两个特征值j=k+xj，k=n-2,j=1,2.)for n=2:2D=eig(Ak);x1=D(1),x2=D(2),tzg1=tzgk+x1,tzg2=tzgk+x2,endtzg(1,1)= tzg1; tzg(2,1)= tzg2;tzg=sort(tzg(:,1);disp(请注意：n阶实对称矩阵A 的全部真特征值lamoda和精度为jd 的近似特征值tzg如下： )lamoda=sort(eig(A)
