数值方法课程设计幂法反幂法计算矩阵特征值和特征向量附Matlab程序Word格式文档下载.docx

1、在本论文中，我们主要讨论矩阵的特征值和特征向量的计算，我们知道，有很多现 实中的问题都可以用到矩阵特征值与特征向量计算的知识，比如，在物理、力学和工程 技术方面有很多的应用，并且发挥着极其重要的作用 . 因为这些问题都可归结为求矩阵 特征值的问题，具体到一些具体问题，如振动问题，物理中某些临界值的确定问题以及 一些理论物理中的问题 .在本论文中，我们主要介绍求矩阵的特征值与特征向量的一些原理和方法，原理涉及高得代数中矩阵的相关定理，方法主要介绍冥法及反冥法并利用 MATLAB算法的程序来求解相关问题，加以验证 .2相关定理定理2.1 如果 i （i 1,2,., n） 是矩阵A的特征值，则有n

2、n1 i aii trAi 1 i 12o detA 1 2 n.定理2.2 设A与B为相似矩阵 （B T 1AT），则1o A 与 B有相同的特征值；2o若 x是 B 的一个特征向量，则 Tx是A的特征向量定理2.3 设A （aij）n n ，则A的每一个特征值必属于下述某个圆盘之中：aij aij （i 1,2,.,n）.j1ji定义2.1 设A是n阶是对称矩阵，对于任意非零向量 x，称 R（x） （Ax,x） 为对应于向量 x x,x的Rayleigh 商 .定理2.4 设A Rn n为对称矩阵（其特征值次序记作 1 2 n ，对应的特征向量x1,x2, xn组成规范化正交组，即 （x,

3、xj ） ij ），则3符号说明A:n阶矩阵B:I ： n阶单位阵i （i 1,2,.,n） : 矩阵特征值x: 实数域上的 n维向量vi（i 0,1, n ） ：实数域上的 n维向量uk（k 0,1, ,n, ） ：实属上的规范化向量4冥法及反冥法4.1冥法幂法是一种计算矩阵 A Rn n 的主特征值的一种迭代法，它最大优点是方法简单，适合于计算大型稀疏矩阵的主特征值设 A （aij ） Rn n ，其特征值为 i ，对应特征向量为 xi （i 1, ,n）,即Axi i xi （i 1, ,n）4.1.1 ）且xi, ,xn 线性无关.设A特征值满足：（即 1为强占优）| 1 | | 2

4、| | n |称 vk 为迭代向量 .面来分折 1及 x1与vk关系 .且设 1 0 ）且有（4.1.3 ）nvk Avk 1 A v0 i i xii1vk 1k（ 1x1 2（ 2）k x2 n（ n）kxn）1 （a1x1 k ）由假设（ 4.1.1 ）式, 则kim（ i ） 0（i 2, ,n）lkim k 0且收敛速度由比值 r | |确定. 且有1这说明，当 k充分大时，有 vk / 1k 1x1，或 vk / 1k 越来越接近特征向量 1x1. 下面考虑主特征值 1的计算 .用（vk）i表示vk的第i 个分量，考虑相邻迭代向量的分量的比值 .从而是说明相邻迭代向量分量的比值收敛

5、到主特征 1,且收敛速度由比值 r | 2 |来度量, r 越小收敛越快 ,但r 越小收敛越快 , 但r | 2 | 1, 而接近于 1时, 收敛可能很慢 .定理 4.1 （1）设 A Rn nn 个线性无关的特征向量： （2）设 A 特征值满足| 1| | 2| | n|（3）幂法： v0 0（且 1 0）vk Avk 1 （k 1,2, ） 则 （1） lim （vk 1）i 1x1 ；k （vk ）i 1 1又设 A有 n个线性无关的特征向量， x1,x2, ,xn,其中Axi 1xi（i 1, ,r ）, Axi i xi（i r 1, ,n）对于任意初始向量v0 ixi（且 i, ,

6、 r不全为零 ）则由幂法有rnvk AKv0 1k ixi i（ i）k xii1 i r 1 ir1k（ i xi k）lim vkk i xi , （设 1, r 不全为零）k 1k i 1 i i 1 r（ k 0当 k ）由此，当 k 充分大时， vk / 1k 接近于与 1 对应的特征向量的某个线性组合应用幂法计算 A 的主特征值 1及对应的特征向量时， 如果 1 1（或 1 1），迭代向量的各个不等于零的分量将随 k 而趋于无究（或趋于零） ，这样电算时就可能溢出为此，就南非要将迭代向量加以规范化 .设有非零向量其 中 m axv（） 表 示 向 量 v 绝 对 值 最 大 的 元

7、 素 ， 即 如 果 有 草 药 （v）i0 m ax（v）i , 则max（ v） （v）i 01 i n其中 i0 为所有绝对值最大的分量中最小指标显然有下面性性质：设 t 为实数 ,v r n ，则max（tv） t max（v）在定理 4.1 条件下幂法可改进为：任取初始向量 u0 v0 0（且 1 0） .迭代： 规范化：v1 Av0 v1 Au0 ， u1max（v1） max（ Av0 ）A2v0 v2 A2v0v2 Au1 ,u2 2 kmax（ Av0 ） max（v2 ） max（ Ak v0 ） （4.16）面考查 uk,vk 与计算 1及x1的关系.由 v0 i xi（

8、2） 考查迭代向量序列 :v AK v0 k1 （ 1x1 k ）k max（ Ak 1v0） max（ 1k 1（ 1x1 k 1）1x1 kmax（ 1x1 k 1）max（ 1x1 k ）于是， k max（ vk ） 1 1,（当k ）max（ 1x1 k 1定理 （ 改进幂法 ）（1） 设 A Rn n有n 个线性无关特征向量 ;（2） 设 A 特征值满足且 Axi i xi （i 1,2, ,n）（3） uk, vk由改进幂法得到 （4.1.7） 式）, 则有且收敛速度由比值 r | 2 | 确定.个子程序求矩阵按模最实现幂法，每迭代一次主要是计算一次矩阵乘向量 （Au） ，可编

9、大特征值如下： %这个函数用于使用幂法求矩阵特征向量和特征值 %A-矩阵， v- 初始向量， e- 精度 function t,p=pm（A,v,e）u=v./max（abs（v）; %old = 0; %记录上一次迭代得到的特征值while （1） v=A*u;if （abs（max（v）-old）e）break ;endold = max（v）;p = u;t = max（v）;例 1. 为检验以上代码的正确性， 我们使用以上代码计算以下矩阵的最大特征值和特征向 量1,3,3A 2,1,54,3,1结果为：例 2. 利用你所编制的子程序求如下矩阵（从 60 到 70 阶）211 2 1 A

10、1 2 112按模最大、按模最小的特征值及对应特征向量。解：代码见附录，运行得到的结果如下：以上仅给出特征值的计算结果。特征向量见附录，这里给出 70 阶的特征向量：0.58-0.941.00 -0.81 0.54 -0.290.13-0.050.01-0.000.000.00 -0.00 0.00 0.000.00 0.00 0.00 0.000.00 0.00 0.000.00 0.000.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 -0.00 0.00 -0.000.00 -0.00 0.01 -0.05 0.13 -0.29 0.54 -0.8

11、1 1.00 -0.94 0.58 4.2 反冥法（1） 反幂法可用来计算矩阵按模最小的特征值及对应的特征向量设 A Rn n 为非厅异矩阵， A特征值满足2对应特征向量 x1,x1, ,xn 为线性无关，则 A 1特征求值为特征向量为 x1,x2, xn.因此计算 A的按模最小的特征值 n 的部题就是计算 A 1 按模最大的特征值部题 对于 A 1应用幂法迭代（称为反幂法） ，可求矩阵 A 1的主特征值 1/ n. 反幂法迭代公式：文案大全任取初始向量 v0 u0 0（且 n 0） ，k 1， 2，4.2.1 ）vk A uk 1 k max（vk ）uk vk /uk其中迭代向量 vk 可

12、通过解方程组求得Av k uk 1如果 A Rn n有n 个线性无关特征向量且 A 特征值满足 :1 n 1 n 0则由反幂法 （2.11） 构造的向量序列 uk,vk 满足，现要求对应| j p| | i p|,（ j i）| j p| i p,（i j）说明 是（A pI） 1的主特征值 . jp现对（A pI） 1应用幂法得到反幂法计算公式： 取初始向量 u0 v0 0（且 0）,jk 1,2, ,vk （A pI） uk 1k max（vk ） （ 4.2.2 ）与定理 8证明类似,可得下述结果 .定理 10 （1）设 A Rn n有n个线性无关特征向量即 Axi ixi（i 1, ,

13、n） .（2）取 p j （为A特征值 j一个近似值），设（A pI） 1存在且| j p| | i p|（j i）则由反幂法迭代公式（ 2，12）构造向量序列 uk,vk 满足：xj（a）uk （当k ）max（ xj ）（b） max（v ） （当k ）k k p反幂法迭代公式中 vk 是以通过解方程组（A pI）vk uk 1求得.为了节省计算量，可先将 （A pI ）进行三角分解 .P（A pI）LU其中 p为置换阵，于是每次迭代求 vk 相当于求解两个三角形方程组Ly k pu k 1Uv k yk可按下述方法取 v0 u0，即选 u0 使Uv1 L 1Pu0 （1, ,1）T回代求

14、解即求得 v1.反幂法计算公式：1分解计算P（A pI）LU ，且保存 L,U 及 P信息2反幂法迭代（1）Uv1 （1, ,1）T求 v1u1 v1 / 1, 1 max（ v1）（2） k 2,1 ）Lyk Puk 1求 ykUvk yk 求vk2 ） k max（ vk ）3 ） uk vk / k对于计算对称三对角阵，或计算 Hessenberg 阵对应于一个给定的近似特征值的特 征向量，反幂法是一个有效方法 .使用 Matlab 编写一个使用反幂法求矩阵最小特征值和特征向量的程序如下： function s,y=fpm（A,x0,eps） % s 为按模最小特征值， y 是对应特征向

15、量 文案大全k=1;r=0; % r 相当于 0?y=x0./max（abs（x0）; % 规范化初始向量L,U=lu（A）;z=Ly; x=Uz; u=max（x）;s=1/u; % 按模最小为 A-1 按模最大的倒数 .if abs（u-r）eps k=k+1;r=u; y=x./max（abs（x）; z=Ly;x=Uz;u=max（x）;end m,index=max（abs（x）; s=1/x（index）;同样，取一个矩阵进行测试：计算结果为：A12 按模最小的特征值及对应特征向量。 代码见附录，程序结果如下图：同样只给出 70 阶时的特征值，具体结果见附录0.04 0.09 0.

16、13 0.18 0.22 0.26 0.30 0.35 0.39 0.43 0.47 0.51 0.54 0.580.620.650.680.720.750.770.800.830.850.870.890.910.930.950.960.970.980.991.000.580.540.510.470.430.390.350.300.260.220.180.090.04 参考文献1 姜启源，谢金星，叶俊编数学模型（第三版）M北京：高等教育出版社， 2005： 1-202.2 王建卫，曲中水 凌滨编著 . MATLAB 7.X 程序设计 M. 北京：中国水利水电出版 社，2007：55-80.3

17、李庆扬，王能超，易大义编著 . 数值分析（第四版） M. 武汉：华中科技大学出版社，2006：219-245.附录%这个函数用来生成老师要求记算的那个矩阵， n 是指定阶数function A=createMatrix（n）A = zeros（n）; %先全部初始化为 0for i=1:for j=1:if （i=j）A（i,j）=2; %设置主对角线上的值为 2else if （i=j-1 | i=j+1） %设置主对角线傍边的两条斜线上的的值为 -1A（i,j）=-1;%这个函数用于使用幂法求矩阵特征向量和特征值%A-矩阵， v- 初始向量， e- 精度 function t,p=pm（A

18、,v,e） u=v./max（abs（v）; % old = 0; u=v./max（abs（v）;%这个程序用于求 60-60 阶矩阵的特征值和特征向量clccleare = 0.01;for i=60:70A = createMatrix（i）; %生成要计算的矩阵v = ones（i,1）; %生成初始微量v（1） = 1;t,p = pm（A,v,e）; %计算fprintf（ %d 阶 特征值： %fn ,i,t）; %输出特征值 %以下三句代码为输出特征值和特征微量% fprintf（%d 阶 :%f ,i,t）;%.2f ,p）;n）;end % 使用反幂法求矩阵按模最小特征值function s,y=fpm（A,x0,eps） % s 为按模最小特征值， y 是对应特征向量 k=1;eps% 终止条件 .k=k+1;y=x./max（abs（x）;60-60 阶矩阵的特征值和特征向量m,index=max（abs（x）;s=1/x（index）;%这个程序用于使用反幂法求t,p = fpm（A,v,e）;%d 阶:fprintf（ ,p）; ）;使用幂法求矩阵最大特征值和特征向量结果：60阶:3.754011 0.58 -0.94 1.00 -0.81 0.54 -0.29 0.13 -0.05 0.01 -0.00