精编资料推荐幂法,反幂法求解矩阵最大最小特征值及其对应的特征向量

1、专业资料推荐数值计算解矩阵的按模最大最小特征值及对应的特征向量一 .幂法1.幕法简介：当矩阵A满足一定条件时，在工程中可用幕法计算其主特征值(按模最大)及其特征向量。矩阵A需要满足的条件为：(1)|2卜门|- 0, i为A的特征值存在n个线性无关的特征向量，设为X2，Xn1.1计算过程：n对任意向量x(0)，有x(0) -7 rUi,i不全为0，则有i 4x(k Ax(k)Ak1x(0)nnk 1k 1八 A am = a A 5(八2/ n、k 1()3nUn1i =1i T1u1可见，当丨厘丨越小时，收敛越快；且当k充分大时，有x(k 1) k1:-x(kx(k)二 k ： U1 1U1x

2、(k)；(k)1，对应的特征向量即是x102算法实现(1).输入矩阵A，初始向量X，误差限最大迭代次数N(k).k 1,0; y(k)max(abs(x(k).计算x Ay, max(x);若-：；,输出,y,否则,转(5)(5)若k N ,置k- k 1/ ：-,转3,否则输出失败信息，停机.3 matlab程序代码function t,y=lpowerA,xO,eps,N)% t为所求特征值，y是对应特征向量k=1;z=0;y=x0./max(abs(x0); x=A*y;% z相当于%规范化初始向量%迭代格式b=max(x);% b相当于：if abs(z-b)eps & k lx, y

3、=lpower (A, xO, eps, X)ans =2. 53650. 74820. 6497L 00000. 01661. 48012. 5365 A#y-t*yans -1. 0e-004 +-0.1603-0, 16840结果正确，表明算法和代码正确，然后利用此程序计算15阶Hilb矩阵，与eig(A)的得到结果比较，再计算 A*y-t*y ，验证y是否是对应的特征向量。设置初始向量为x0=ones(15,1),结果显示如下 A二hilb(15); x0=ones (15,1):】；eps-le6: N=30; Lt, y -1 power (A, xO, eps, Xt 二L 84

4、59y =V - eig(A) A*y-T*y1. 0000ans =ans =0. 62280. 47060. 3838-0. 00DO0. 0000L 0e-007 *a0. 32640. 0000-0“ 55160. 0000F 64 3G0. 28510. 0000-0“ 64650. 25370. 0000池 62450. 22900. 0000- 59530. 20890. 0000-0.56500.19220. 00007 53590. 17800, 0000-0.5Q阴0. 16590, 0004-0” 48340, 0056-0.46030. 1554a 0572-0. 43

5、900. 1462a. 4266-0.41950. 13811. 8459-0.401bjn. rn rs可见，结果正确。得到了 15阶Hilb矩阵的按模最大特征值和对应的特征向 量。二反幂法1反幕法简介及其理论在工程计算中，可以利用反幕法计算矩阵按模最小特征值及其对应特征向 量。其基本理论如下，与幕法基本相同：1由Ax二X= x = A(X),则Axx，可知，A和A-1的特征值互为倒数，Z求A按模最小特征值即求 A-1的按模最大特征值，取倒数即为A的按模最小特征 值所以算法基本相同，区别就是在计算 x(k 1时，不是令x(k1)=Ay(k),而是x(kA-1y(k)具体计算时，变换为人严1)

6、y(k)；对a做LU分解，来计算x(k 1)2算法实现(1) .输入矩阵A,初始向量x,误差限；,最大迭代次数N,(2) 置k -1,， “0, y,max(abs( x)(3) .作三角分解A二LU(4) .解方程组 LUx = y( Lz 二 y,Ux 二 z),(5) .： max(x), - :一 -(6)若| :-咒0卜；,输出,y,停机,否则转(7),转(4)；max(abs(x)(7).若k ： N,置k k 1,，屮否则输出失败信息，停机.3 matlab程序代码fun ctio n s,y=i nvpower(A,xO,eps, n) % s为按模最小特征值，y是对应特征向量

7、 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按模最大的倒数.if abs(u-r)eps & k s, y = invpower (A? eps n)-5.9673e-0180. 0000-0.00000.0000-0.00060.0050-0. 02450. 0699-0. 0968-0. 04210. 4497-0. 9084L 000065270.2379-0.0375 eig(A)ans =-0. 00000. 00000. 00000. 00000. 0

8、0000+ 00000+ 00000, 00000” 00000. 00000, 00040. 00560. 05720. 42661 8459 A*vs*y已nw =L 0e-017 *0. 3903-0. 56380. 0000-0.5208-0.47416 35400. 4537-0. 2746CL 9724-0. 6556-0.62886 66180. 06590. 1419可见，结果真确。得到了 量。15阶Hilb矩阵的按模最大特征值和对应的特征向三计算条件数矩阵A的条件数等于A的范数与A的逆的范数的乘积，即 cond(A)= IIAA(-1) II,对应矩阵的3种范数，可以定义3种

9、条件数。函数cond(A,1)、cond(A)或cond(A inf)是判断矩阵病态与否的一种度量，条件数越 大表明矩阵的病态程度越大这里我们选择矩阵的2范数，即cond(A)分别为 矩阵At A的最大和最小特征值而如果A为对称矩阵，如Hilb矩阵，AtA的最大最小特征值，分别为 A的 最大最小特征值的平方。所以 cond(A)为A的最大最小特征值得比值。对于本 例中的15阶Hilb矩阵来说，利用上面计算结果得其条件数(选择第二种条件数)为：3.0934e+017；这与直接利用cond(A)得到的结果：2.5083e+017在同一数量级，再次表明了上述算得得最大最小特征值的正确性，同时又表明

10、阵是病态矩阵。Hilb 矩四.Aitken商加速法1.简介与原理若：ak收敛与a,且lim 一二c = 0;即：ak ：线性收敛， k护 ak - a当k充分大时，有玉口玉心ak - aakd1 - a、，、，(Xn 雄xn*)yn 2 =XnXn42 一 2Xn* + Xn一 r r(ak*ak)2“二 a ak=-= Okak 2 一 2ak -1 ak用ak逼近a,这种方法称为Aitken加速法.同幕法和反幕法计算最大和最小特征值类似，如果计算最大特征值,为x(k1) = Ay(k)；计算最小特征值时，迭代格式为则迭代格式x(k1 A-1y(k),即y(k Ax(k1)2.算法实现计算按

11、模最大特征值算法如下：(1).输入A- (aj),初始向量x,误差限；,最大迭代次数N,(2).置7 十00厂 max(x).计算 x=Ay,置 max(x)二:2( -Of ).计算O。-丄，1 *0:;,输出 , y停机，否则转(6),(5) .若卩-町(6) .若 k : N,置-：1 =-； 0, - 2 =1，儿=八 0,k * “mxxr 弹否则,输出失败信息，停机.类似幕法和反幕法可以写出按模最小特征值算法，此处不再赘述。专业资料推荐3.matlab程序代码function r,y=aitken(A,xO,eps,n) % r按模最大特征值$为对应特征向量 k=1;a0=0;%

12、a相当于0a仁1; % a1相当于:1r0=1;%相当于2中的y=x0./max(abs(x0);%规范化初始向量x=A*y;a2=max(abs(x); % a2相当于：2r=a0-(a1-a0)A2/(a2-2*a1+a0); % 相当于if (a2-2*a1+a0)=0%若上式中分母为0,则迭代失败，返回disp 初始向量迭代失败return;endif abs(r-r0)eps & k0 |aa(i ndex)=0r=r;elser=-r;endendend9专业资料推荐类似可得按模最小特征值和特征向量的代码如下：与上面类似，所不同的只 是迭代格式不同fun cti on r,y=i

13、nv aitke n( A,xO,eps, n)k=1;a0=0;a1=1;r0=1;y=x0./max(abs(x0);L,U=lu(A);%迭代格式的不同z=Ly;x=Uz;a2=max(abs(x);r=a0-(a1-a0)A2/(a2-2*a1+a0);if (a2-2*a1+a0)=0disp 初始向量迭代失败return;endif abs(r-r0)eps & k0 |aa(i ndex)=0 r=1/r;elser=-1/r;endend 4.计算Hilb矩阵特征值ones(15,1)此处不再举例，而是直接应用于15阶Hilb矩阵，初始向量选为结果如下，并将结果与幕法和反幕法得

14、到结果比较14 A二hilb (15);xO=ones (1& I); eps=le-6; n=30;y_ =aitken(AJ xOj eps, n)1. 8459L 00000. 62290. 47060, 38380. 32640. 28510, 25370. 22900.20890. 19220. 17800. 16590. 15540. 14620. 1381同理，最小这与幕法得到的特征值和特征向量一致， 表明算法和代码正确; 特征值结果如下： r, y = invaitken & x0f epSj n)5. 9673e-018专业资料推荐0+ 0000 -0, 00000. 000

15、0 -0, 00060. 0050 -0. 02450. 0699-0. 0968-Q. 04210. 44977+ 9084L 0000- 6527Q. 2379 -0. 0375这与反幕法得到的结果一致，表明结果正确五，对称矩阵的Rayleigh商加速法1. 简介与原理A为对称矩阵，设x = 0,则称R(x)二xtAxTx x为关于A的Rayleigh商x15原理如下:设Xj =0为的特征向量，即AXj二iXi用Xi左乘上式有:(k)y()x(k 1)即k)_x(k)max(x(k)二 Ay(k)_(k)T/(k)、(y() (y()这称为Rayleigh商加速法(k)其中 R(y(k)1

16、?y(k) max(x(k)专业资料推荐2. 算法实现(1).输入矩阵A,初始向量X,误差限；,最大迭代次数N,(2).置匕 1, r0,yxmax(abs(x)20(3) x A* y, rxmax(abs(x)(4) .若 |r _r卜：；,输出r,y,停机,否则转(5),(5) .若 k : N,置kk 1, r r, y否则输出失败信息，停机.3. Matlab程序代码function r,y=rayleigh(A,xO,eps,n)% r 是特征值，y是特征向量k=1; r0=0;y=xO./max(abs(xO);x=A*y; %迭代格式计算新的xr=dot(y,x)/dot(y,

17、y); % Reyleigh 商if abs(r-r0)eps & knk=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)eps & k I); eps=le6;n30 _r, y =rayleigh (A x0P eps n)8459v 二L 00000. 62290. 47060. 38390, 32650. 2852