幂法反幂法求解矩阵最大最小特征值及其对应的特征向量.docx
- 文档编号:11931916
- 上传时间:2023-06-03
- 格式:DOCX
- 页数:21
- 大小:120.59KB
幂法反幂法求解矩阵最大最小特征值及其对应的特征向量.docx
《幂法反幂法求解矩阵最大最小特征值及其对应的特征向量.docx》由会员分享,可在线阅读,更多相关《幂法反幂法求解矩阵最大最小特征值及其对应的特征向量.docx(21页珍藏版)》请在冰点文库上搜索。
幂法反幂法求解矩阵最大最小特征值及其对应的特征向量
数值计算解矩阵的按模最大最小特征值及对应的特征向量
—一.幂法
1.幕法简介:
当矩阵A满足一定条件时,在工程中可用幕法计算其主特征值(按模最大)及其特征向量。
矩阵A需要满足的条件为:
⑴I1II2|n|-0,i为A的特征值
(2)存在n个线性无关的特征向量,设为Xi,X2,…,Xn
1.1计算过程:
n
对任意向量x(0),有x(0)八:
-M—不全为0,则有
i4
X(k岀)=Ax(k)==Ak岀乂。
)
nn
Ak1aqa扌15
i=1i=1
■k
1
2
可见,当1—1越小时,收敛越快;且当k充分大时,有?
"1
x(k1)
>(k+1)
1,对应的特征向量即是x
2算法实现
⑶.计算xAy,…max(x);
⑷若|•一十:
;,输出-,y,否则,转(5)
(5)若N,置k「k1^-,转3,否则输出失败信息,停机.
3matlab程序代码
k=1;
z=0;
y=x0./max(abs(x0));x=A*y;
%z相当于■
%规范化初始向量
%迭代格式
b=max(x);
%b相当于:
ifabs(z-b) %判断第一次迭代后是否满足要求 end whileabs(z-b)>eps&&k z=b; y=x./max(abs(x)); x=A*y; b=max(x); end [m,index]=max(abs(x));%这两步保证取出来的按模最大特征值 t=x(index);end %是原值,而非其绝对值。 4举例验证 选取一个矩阵A,代入程序,得到结果,并与eig(A)的得到结果比较,再计算A*y-t*y,验证y是否是对应的特征向量。 结果如下: »A=[l10.5;110.25;0.5,0.25,2] A= 1.0000 E0000 0.5000 1.0000 1.0000 0,2500 0.5000 0.2500 2.0000 »x0=[l;l;ll;eps=le-5;220;>>y]=lpower(A,xO,eps,X) ans= 2.5365 0.7482 0.6497 L0000 ^0.0166 1.4801 2.5365 »A#y-t*y ans- 1.0e-004+ -0.1603 -0,1684 0 结果正确,表明算法和代码正确,然后利用此程序计算15阶Hilb矩阵,与 eig(A)的得到结果比较,再计算A*y-t*y,验证y是否是对应的特征向量。 设 置初始向量为x0=ones(15,1),结果显示如下 »A二hilb(15); »x0=ones(15,1): 】;eps-le^6: »N=30; >>Lt,y]-1power(A,xO,eps,X} t二 L8459 y- »eig(A) 2>A*y-T*y 1.0000 ans= ans= 0.6228 0.4706 0.3838 -0.00DO 0.0000 1.Ge-00? * 0 0.3264 0.0000 -0“5516 0.0000 F6436 0.2851 0.0000 -0“6465 0,2537 0.0000 P・624S 0.2290 0.0000 -0.5953 0.2089 0.0000 -0*5650 O1922 0.0000 -0.5359 0.1780 0,0000 -Q.5086 0.1659 0,0004 -0”4834 60056 -0.4603 0.1554 0.0572 -0.4390 0.1462 0.4266 -0.4195 61381 1.8459 -0.4015 Cflf 可见,结果正确。 得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量。 二•反幂法 1■反幕法简介及其理论 在工程计算中,可以利用反幕法计算矩阵按模最小特征值及其对应特征向量。 其基本理论如下,与幕法基本相同: 1 由AxVx=x=A‘(’X),则A’xx,可知,A和A-1的特征值互为倒数, Z 求A按模最小特征值即求A-1的按模最大特征值,取倒数即为A的按模最小特征值所以算法基本相同,区别就是在计算x(k1)时,不是令X(k1)=Ay(k),而是X(kJA-1y(k)具体计算时,变换为Ax(k-1^y(k);对A做LU分解,来计算x(k1) 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),-■ 1 (6)若-’0卜: ;,输出—,y,停机,否则转(7), ma莎面,转⑷; (7).若k: : N,置k・k1,•,y・ 否则输出失败信息,停机. 3matlab程序代码 function[s,y]=invpower(A,xO,eps,n)%s为按模最小特征值,y是对应特征向量k=1;r=0;%r相当于-o y=x0./max(abs(x0));%规范化初始向量 [L,U]=lu(A); z=L\y; x=U\z; u=max(x); s=1/u;%按模最小为A-1按模最大的倒数. ifabs(u-r) return end whileabs(u-r)>eps&&k k=k+1; r=u; y=x./max(abs(x)); z=L\y; x=U\z; u=max(x); end [m,index]=max(abs(x));%这两步保证取出来的按模最大特征值 s=1/x(index);%是原值,而非其绝对值。 4举例验证 同幕法一样,选取一个矩阵 end A,代入程序,得到结果,并与eig(A)的得到结果比 较,再计算A*y-t*y,验证y是否是对应的特征向量 y二 L0000 -0.9517 -0.1300 1.0e-016宕 -0.1388 -0.0694 -0.1908 可见结果正确,然后利用此程序计算15阶Hilb矩阵,eig(A)的得到结果比 较,再计算A*y-s*y,验证y是否是对应的特征向量。 设置初始向量为 x0=ones(15,1),结果显示如下 »A=liilb(15): x0=ones(1571);eps=le~6;n=30; >>[s,y]=invpower(A? eps^n) -5.9673e-018 可见,结果真确。 得到了15阶Hilb矩阵的按模最大特征值和对应的特征向量。 三•计算条件数 矩阵A的条件数等于A的范数与A的逆的范数的乘积,即cond(A)=IIAA(-1)II,对应矩阵的3种范数,可以定义3种条件数。 函数 cond(A,1)、cond(A)或cond(Ainf)是判断矩阵病态与否的一种度量,条件数越大表明矩阵的病态程度越大• 这里我们选择矩阵的2范数,即cond(A)j分别为 矩阵AtA的最大和最小特征值 而如果A为对称矩阵,如Hilb矩阵,AtA的最大最小特征值,分别为A的最大最小特征值的平方。 所以cond(A)为A的最大最小特征值得比值。 对于本例中的15阶Hilb矩阵来说,利用上面计算结果得其条件数(选择第二种条件数)为: 3.0934e+017;这与直接利用cond(A)得到的结果: 2.5083e+017在同一 数量级,再次表明了上述算得得最大最小特征值的正确性,同时又表明Hilb矩 阵是病态矩阵。 四.Aitken商加速法 1.简介与原理 若: ak收敛与a,且lim一二c=0;即: ak: 线性收敛,k护ak_a 当k充分大时,有亟口玉心 ak-aakd1-a 、,、,(Xn雄一xn出) yn=Xn2- Xn七一2Xn*+Xn (ak1一ak)2Oak ak2-2ak-1'ak 用ak逼近a,这种方法称为Aitken加速法. 同幕法和反幕法计算最大和最小特征值类似,如果计算最大特征值,则迭代格式 为x(k°=Ay(k);计算最小特征值时,迭代格式为 x(k1^A-1y(k),即y(k)=Ax(" 2.算法实现 计算按模最大特征值算法如下: (1).输入A-©),初始向量x,误差限;,最大迭代次数N, ⑵.置k=1, X 0=0,-1=0,°=1.0,y= max(x) ⑶.计算x=Ay,置max(x)=匚2, (4).计算丸= 2 (G-Of) °(10丿0) °2_2口1乜0 (5).若丸一町 <名,输出&,y停机,否则转(6), (6).若k: : N,置-: : 1=0,二2=■■1,八=八0, 否则,输出失败信息,停机. 类似幕法和反幕法可以写出按模最小特征值算法,此处不再赘述 3.matlab程序代码function[r,y]=aitken(A,xO,eps,n)%r按模最大特征值$为对应特征向量 k=1; a0=0;%a相当于>0 a仁1;%a1相当于: '1 r0=1;%相当于2中的‘0 y=x0./max(abs(x0));%规范化初始向量 x=A*y; a2=max(abs(x));%a2相当于: 2 r=a0-(a1-a0)A2/(a2-2*a1+a0);%相当于' if(a2-2*a1+a0)==0%若上式中分母为0,则迭代失败,返回 disp"初始向量迭代失败" return; end ifabs(r-r0) end whileabs(r-r0)>eps&&k k=k+1; a0=a1; a1=a2; r0=r; y=x./max(abs(x)); x=A*y;%迭代格式 a2=max(abs(x)); if(a2-2*a1+a0)==0%若分母为0,则迭代失败,返回return; end r=a0-(a1-a0)A2/(a2-2*a1+a0); [m,index]=max(abs(eig(A)));%以下代码保证取出来的按模最大特征值aa=eig(A);%是原值,而非其绝对值。 ifaa(index)>0||aa(index)==0 r=r; else r=-r; end end end 类似可得按模最小特征值和特征向量的代码如下: 与上面类似,所不同的只是迭代格式不同• function[r,y]=invaitken(A,xO,eps,n) k=1; a0=0; a1=1; r0=1; y=x0./max(abs(x0)); [L,U]=lu(A);%迭代格式的不同 z=L\y; x=U\z; a2=max(abs(x)); r=a0-(a1-a0)A2/(a2-2*a1+a0); if(a2-2*a1+a0)==0 disp"初始向量迭代失败" return; end ifabs(r-r0) return end whileabs(r-r0)>eps&&k k=k+1; a=b; b=c; r0=r; y=x./max(abs(x)); z=L\y; x=U\z; a2=max(abs(x)); if(a2-2*a1+a0)==0 return; end r=a0-(a1-a0)A2/(a2-2*a1+a0); end [m,index]=min(abs(eig(A)));%以下代码保证取出来的按模最大特征值 aa=eig(A);%是原值,而非其绝对值。 ifaa(index)>0||aa(index)==0 r=1/r; else r=-1/r; end end4.计算Hilb矩阵特征值 此处不再举例,而是直接应用于15阶Hilb矩阵,初始向量选为 ones(15,1) 结果如下,并将结果与幕法和反幕法得到结果比较 »A=hilb(15); xO=ones(1&I); »eps=le-6; »n=30; _Yjy_=aitken(AJxO,eps,n) 1.S459 y二 L0000 0.6229 0,4706 0,3838 0.3264 0.2851 0,2537 0.2290 0.2089 0.1922 0.1780 0.1659 0,1554 0.1462 0.1381 这与幕法得到的特征值和特征向量一致,表明算法和代码正确;特征值结果如下: 同理,最小 >>[r,y]=invaitken&xO,epSjn) 5.9673e-018 0,0000 -0,0000 0.0000 -0,0006 0.0050 -0.0245 0.0699 -0.0968 -0.0421 0.4497 7+9084 L0000 -0+6527 Q.2379 -0.0375 这与反幕法得到的结果一致,表明结果正确 五,对称矩阵的Rayleigh商加速法 1.简介与原理 xTAx A为对称矩阵,设x=0,则称R(x)二耳x为关于A的Rayleigh商 xx 原理如下: 2.算法实现 (1).输入矩阵A,初始向量X,误差限;,最大迭代次数N, (4).若|r-r0卜;,输出r,y,停机,否则转(5), 转(3); max(abs(x)) (5).若kN,置k「k1,r0•r,y「 否则输出失败信息,停机. 3.Matlab程序代码 function[r,y]=rayleigh(A,xO,eps,n)%r是特征值,y是特征向量 k=1;r0=0; y=x0./max(abs(x0)); x=A*y;%迭代格式计算新的x r=dot(y,x)/dot(y,y);%Reyleigh商 ifabs(r-r0) return end whileabs(r-r0)>eps&&k k=k+1; r0=r; y=x./max(abs(x)); x=A*y; r=dot(y,x)/dot(y,y); end end 类似得计算按模最小特征值的Rayleigh商加速法,如下: function[r,y]=invrayleigh(A,x0,eps,n) k=1;r0=0; y=x0./max(abs(x0)); [L,U]=lu(A);%迭代格式不同 z=L\y; x=U\z; r=dot(y,x)/dot(y,y); ifabs(r-r0) return end whileabs(r-r0)>eps&&k k=k+1; r0=r; y=x./max(abs(x)); z=L\y; x=U\z; r=dot(y,x)/dot(y,y); end r=1/r; end 4.计算Hilb矩阵特征值 ones(15,1) 此处不再举例,而是直接应用于15阶Hilb矩阵,初始向量选为 结果如下,并将结果与幕法和反幕法得到结果比较 »A=hilb(15): xO=ones(15>I);eps=le~6; n—30■ >>_r,y]=rayleigh(A>x0Peps^n) 8459 v二 L0000 0.6229 0.4706 0.3839 0,3265 0.2852 0.2538 0.2290 0.2089 0.1922 0.1781 0.1660 0.1555 0.1463 这与幕法得到结果一致,表明算法和代码正确 0.1381 同理,最小特征值如下: »[rfy]=invray1eigh(A,xO,eps,n) -5.9673e-018 0.0000-0.0000 0.0000 -0.0006 0.0050 -0.0245 0.0699 -0.0968 -0.0421 0.4497 -0.9084 1.0000 -0.6527 0.2379 与反幕法得到结果一致,表明代码和算法正确 70375
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 幂法反幂法 求解 矩阵 最大 最小 特征值 及其 对应 特征向量