幂法,反幂法求解矩阵最大最小特征值及其对应地特征向量

1、实用标准文案数值计算解矩阵的按模最大最小特征值及对应的特征向量一. 幂法1. 幂法简介：当矩阵 A 满足一定条件时，在工程中可用幂法计算其主特征值( 按模最大 )及其特征向量。矩阵A 需要满足的条件为：(1) | 1 | | 2 | .| n | 0 ,i 为 A的特征值(2) 存在 n 个线性无关的特征向量，设为1.1 计算过程 ：对任意向量 x(0 )，有 x( 0)ni ui ,i 1x(k 1)Ax(k).Ak 1 x( 0 )nk 1nAk 1u uii ii ii 1i1x1, x2 ,., xni 不全为 0，则有k1u ( 2 )k 1a u( n )k 1a un1112 2

2、n11k1u111可见，当 |2 | 越小时，收敛越快；且当 k 充分大时，有1x(k1 )k1x(k 1 )11u1，对应的特征向量即是 x(k 1 )。x(k )k1u1x(k )112 算法实现输入矩阵，初始向量 ，误差限,最大迭代次数N(1).Ax(2).k1,0; y( k)x( k )max(abs( x( k ) )计算xAy,max(x);(3).若| ,输出, y,否则转(5)(4).,若kN ,置kk1,转3 ,否则输出失败信息,停机.(5).3 matlab程序代码精彩文档实用标准文案functiont,y=lpowerA,x0,eps,N) % t为所求特征值， y是对

3、应特征向量k=1;z=0;% z相当于y=x0./max(abs(x0);%规范化初始向量x=A*y;%迭代格式b=max(x);% b相当于if abs(z-b)<eps % 判断第一次迭代后是否满足要求 t=max(x);return ;endwhile abs(z-b)>eps && k<Nk=k+1;z=b;y=x./max(abs(x);x=A*y;b=max(x);endm,index=max(abs(x);%这两步保证取出来的按模最大特征值t=x(index);%是原值，而非其绝对值。end4 举例验证选取一个矩阵 A，代入程序，得到结果，并与e

4、ig(A) 的得到结果比较，再计算 A*y-t*y ，验证 y 是否是对应的特征向量。结果如下：精彩文档实用标准文案结果正确，表明算法和代码正确， 然后利用此程序计算15 阶 Hilb 矩阵，与eig(A) 的得到结果比较，再计算A*y-t*y ，验证 y 是否是对应的特征向量。设置初始向量为 x0=ones(15,1) ，结果显示如下精彩文档实用标准文案可见，结果正确。得到了15阶 Hilb 矩阵的按模最大特征值和对应的特征向量。二. 反幂法1. 反幂法简介及其理论在工程计算中，可以利用反幂法计算矩阵按模最小特征值及其对应特征向量。其基本理论如下，与幂法基本相同：由Axx x A1，则A11

5、-1的特征值互为倒数，( x)xx ，可知，A 和 A-1求 A 按模最小特征值即求 A 的按模最大特征值 , 取倒数即为 A 的按模最小特征值所以算法基本相同，区别就是在计算x( k 1)时，不是令 x( k 1)Ay(k ) , 而是 x( k 1)A-1 y(k )具体计算时，变换为Ax( k 1)y( k ) ; 对 A做LU 分解，来计算 x( k 1)2. 算法实现精彩文档实用标准文案(1).输入矩阵 A,初始向量 x,误差限, 最大迭代次数 N ,(2).置 k1, 00, yx,max( abs( x)(3).作三角分解 ALU(4).解方程组 LUxy( Lzy,Uxz),(

6、5).max( x),(6).若 |0 |,输出 1, y,停机 ,否则转 (7),(7).若 kN , 置 kk1,0, yx, 转 (4);max(abs( x)否则输出失败信息,停机 .3 matlab程序代码functions,y=invpower(A,x0,eps,n) % s为按模最小特征值， y是对应特征向量k=1; r=0;% r相当于 0y=x0./max(abs(x0); %规范化初始向量L,U=lu(A);z=Ly;x=Uz;u=max(x);s=1/u;%按模最小为 A-1 按模最大的倒数 .ifabs(u-r)<eps % 判断第一次迭代后是否满足终止条件 re

7、turnendwhile abs(u-r)>eps && k<n%终止条件 .k=k+1;r=u;y=x./max(abs(x);z=Ly;x=Uz;u=max(x);endm,index=max(abs(x);%这两步保证取出来的按模最大特征值s=1/x(index);%是原值，而非其绝对值。end4 举例验证精彩文档实用标准文案同幂法一样，选取一个矩阵A，代入程序，得到结果，并与eig(A) 的得到结果比较，再计算A*y-t*y ，验证 y 是否是对应的特征向量。可见结果正确， 然后利用此程序计算15 阶 Hilb 矩阵，eig(A) 的得到结果比较，再计算A*

8、y-s*y ，验证 y 是否是对应的特征向量。设置初始向量为x0=ones(15,1) ，结果显示如下精彩文档实用标准文案可见，结果真确。得到了15阶 Hilb 矩阵的按模最大特征值和对应的特征向量。三.计算条件数矩阵 A 的条件数等于 A 的范数与 A 的逆的范数的乘积，即 cond(A)= A· A(-1) ，对应矩阵的 3种范数，可以定义 3种条件数。 函数 cond(A,1) 、cond(A)或 cond(A inf) 是判断矩阵病态与否的一种度量，条件数越大表明矩阵的病态程度越大 .这里我们选择矩阵的2范数，即 cond ( A)1, 1, 2分别为2,矩阵 AT A的最大

9、和最小特征值而如果 A 为对称矩阵，如Hilb 矩阵， AT A 的最大最小特征值，分别为A 的最大最小特征值的平方。所以 cond(A) 为 A 的最大最小特征值得比值。对于本例中的 15阶 Hilb 矩阵来说，利用上面计算结果得其条件数 ( 选择第二种条件数）为：3.0934e+017；这与直接利用 cond(A) 得到的结果： 2.5083e+017 在同一精彩文档实用标准文案数量级，再次表明了上述算得得最大最小特征值的正确性，同时又表明 Hilb 矩阵是病态矩阵。四 .Aitken 商加速法1. 简介与原理若 a收敛与 a,且 lim ak 1ac 0;即 a线性收敛 ,kkakak当

10、 k充分大时 ,有ak1aak2aakaaka1yn 2( xn2xn 1 )2xn 22xn 1xnxn 2( a1ak)2aakk:a?k2ak 1ak 2ak用 a?k逼近 a, 这种方法称为 Aitken 加速法 .同幂法和反幂法计算最大和最小特征值类似，如果计算最大特征值， 则迭代格 式为 x(k 1)Ay( k)；计算最小特征值时，迭代格式为x( k 1)A 1 y( k ) ,即 y( k )Ax (k 1) 。2. 算法实现计算按模最大特征值算法如下：(1).输入 A(aij ), 初始向量 x,误差限,最大迭代次数 N ,( 2).置k1,0 0,0,0 1.0,x,1yma

11、x( x)(3).计算 xAy, 置 max( x)2 ,( 4).计算( 10 )202,210(5).若0, 输出, y停机 , 否则转 (6),( 6).若kN,置 10 ,21,0 ,k 1k,xy转(3),max( x)否则 , 输出失败信息 , 停机 .类似幂法和反幂法可以写出按模最小特征值算法，此处不再赘述。精彩文档实用标准文案3.matlab程序代码function r,y=aitken(A,x0,eps,n) % r按模最大特征值 ,y 为对应特征向量k=1;a0=0;% a相当于0a1=1; % a1相当于1r0=1; %相当于 2中的0y=x0./max(abs(x0);

12、 %规范化初始向量x=A*y;a2=max(abs(x); % a2相当于2r=a0-(a1-a0)2/(a2-2*a1+a0); %相当于if (a2-2*a1+a0)=0%若上式中分母为 0，则迭代失败，返回disp" 初始向量迭代失败 "return;endif abs(r-r0)<eps% 判断第一次迭代后是否满足要求，如满足，则返回结果returnendwhile abs(r-r0)>eps && k<n %终止条件k=k+1;a0=a1;a1=a2;r0=r;y=x./max(abs(x);x=A*y;%迭代格式a2=max(a

13、bs(x);if (a2-2*a1+a0)=0 %若分母为 0，则迭代失败，返回return;endr=a0-(a1-a0)2/(a2-2*a1+a0);m,index=max(abs(eig(A); %以下代码保证取出来的按模最大特征值aa=eig(A);%是原值，而非其绝对值。if aa(index)>0 |aa(index)=0r=r;elser=-r;endend精彩文档实用标准文案end类似可得按模最小特征值和特征向量的代码如下：与上面类似，所不同的只是迭代格式不同 .functionr,y=invaitken(A,x0,eps,n)k=1;a0=0;a1=1;r0=1;y=x

14、0./max(abs(x0);L,U=lu(A);%迭代格式的不同z=Ly;x=Uz;a2=max(abs(x);r=a0-(a1-a0)2/(a2-2*a1+a0);i f (a2-2*a1+a0)=0disp" 初始向量迭代失败 "return;endifabs(r-r0)<eps% 判断第一次迭代后是否满足要求，如满足，则返回结果returnendwhile abs(r-r0)>eps && k<nk=k+1;a=b;b=c;r0=r;y=x./max(abs(x);z=Ly;x=Uz;a2=max(abs(x);if (a2-2*a

15、1+a0)=0return;endr=a0-(a1-a0)2/(a2-2*a1+a0);endm,index=min(abs(eig(A);%以下代码保证取出来的按模最大特征值aa=eig(A);%是原值，而非其绝对值。if aa(index)>0 |aa(index)=0r=1/r;elser=-1/r;endend精彩文档实用标准文案4. 计算 Hilb 矩阵特征值此处不再举例，而是直接应用于15阶Hilb 矩阵，初始向量选为 ones(15,1) ，结果如下，并将结果与幂法和反幂法得到结果比较这与幂法得到的特征值和特征向量一致， 表明算法和代码正确； 同理，最小特征值结果如下：精彩

16、文档实用标准文案这与反幂法得到的结果一致，表明结果正确。五，对称矩阵的 Rayleigh 商加速法1. 简介与原理为对称矩阵，设x0,则称 ()xT Ax 为关于 A的Rayleigh商AR xxT x原理如下：设Xi0为 i的特征向量，即 AXii Xi用Xi 左乘上式有：y( k)x(k )max( x( k) )x( k 1)Ay(k )这称为商加速法。RayleighR( y(k ) )( y(k ) )T x( k 1)( y(k ) )T ( y(k ) )其中 R(y(k ) )1, y(k )x(k )max(x( k) )精彩文档实用标准文案2. 算法实现输入矩阵A,初始向量

17、x,误差限,最大迭代次数N ,(1).置k1, r00, yx,( 2).max(abs(x)(3).xA * y, ryT xyT y若| rr0 |,输出r , y,停机否则转(5),( 4).,若kN,置kk1, r0r, yx,转(3);(5).max(abs( x)否则输出失败信息 , 停机 .3. Matlab 程序代码function r,y=rayleigh(A,x0,eps,n) % r 是特征值， y是特征向量 k=1; r0=0;y=x0./max(abs(x0);x=A*y;%迭代格式计算新的 xr=dot(y,x)/dot(y,y);% Reyleigh商if abs(r-r0)<epsreturnendwhile abs(r-r0)>eps && k<nk=k+1;r0=r;y=x./max(abs(x);x=A*y;r=dot(y,x)/dot(y,y);endend类似得计算按模最小特征值的Rayleigh 商加速法，如下：function r,y=invrayleigh(A,x0,eps,n)k=1;r0=0;y=x0./max(abs(x0);L,U=lu(A); %迭代格式不同z=Ly;x=Uz;r=dot(y,x)/dot(y,y);if abs(r-r0)<epsreturnendwhile