幂法和反幂法求矩阵特征值课程

20/20题目幂法和反幂法求矩阵特征值具体内容随机产生一对称矩阵,对不同的原点位移和初值(至少取3个）分别使用幂法求计算矩阵的主特征值及主特征向量，用反幂法求计算矩阵的按模最小特征值及特征向量,并比较不同的原点位移和初值说明收敛。要求1。认真读题，了解问题的数学原形；２．选择合适问题求解的数值计算方法；３.设计程序并进行计算；4．对结果进行解释说明;采用方法及结果说明对于幂法和反幂法求解矩阵特征值和特征向量的问题将从问题分析，算法设计和流程图，理论依据，程序及结果进行阐述该问题.一．问题的分析：求ｎ阶方阵A的特征值和特征向量，是实际计算中常常碰到的问题,如：机械、结构或电磁振动中的固有值问题等。对于ｎ阶矩阵A,若存在数和n维向量x满足Ax=x(1）则称为矩阵A的特征值,x为相应的特征向量。由高等代数知识可知，特征值是代数方程｜Ｉ-A｜=+a+…＋a+a=０（２）的根。从表面上看,矩阵特征值与特征向量的求解问题似乎很简单，只需求解方程(2)的根,就能得到特征值，再解齐次方程组（I-Ａ）x=0（3）的解，就可得到相应的特征向量。上述方法对于ｎ很小时是可以的.但当n稍大时，计算工作量将以惊人的速度增大，并且由于计算带有误差,方程(2）未必是精确的特征方程，自然就不必说求解方程（2)与（3）的困难了.幂法是一种计算矩阵主特征值（矩阵按模最大的特征值）及对应特征向量的迭代方法，特别是用于大型稀疏矩阵。反幂法是计算海森伯格阵或三角阵的对应一个给定近似特征值的特征向量的有效方法之一。二。算法设计及流程图1、幂法算法（1)取初始向量u（例如取ｕ=(1,1,…1))，置精度要求，置k=1.（2）计算ｖ=Au,m=max（v)，ｕ=ｖ/ｍ(3)若｜m=m|〈，则停止计算(m作为绝对值最大特征值,ｕ作为相应的特征向量)否则置k=k+１,转(２）2、反幂法算法(1）取初始向量ｕ(例如取u＝（1，1，…1）)，置精度要求，置k=1。（2）对A作LU分解，即A=LU(３）解线性方程组Ly=ｕ，Uv＝y(４）计算ｍ=maｘ（v)，u=v/m（5）若|ｍ=ｍ|〈,则停止计算（1/m作为绝对值最小特征值，u作为相应的特征向量）；否则置k＝k+１，转(3).幂法流程图：开始开始输入A；[m,u,index]=pow(A,1e-6)k=0;m1=0v=A*u[vmax,i]=max(abs(v))m=v(i);u=v/mabs(m-m1)<1e-6index=1;break;输出：m，u，index结束m1=m;k=k+1反幂法流程图开始开始输入A；[m,u,index]=pow_inv(A,1e-6)k=0;m1=0v=invA*u[vmax,i]=max(abs(v))m=v(i);u=v/mabs(m-m1)<1e-6index=1;break;输出：m，u，index结束m1=m;k=k+1输入A；[m,u,index]=pow(A,1e-6)三、算法的理论依据及其推导（一）幂法算法的理论依据及推导幂法是用来确定矩阵的主特征值的一种迭代方法，也即,绝对值最大的特征值。稍微修改该方法，也可以用来确定其他特征值。幂法的一个很有用的特性是它不仅可以生成特征值，而且可以生成相应的特征向量。实际上，幂法经常用来求通过其他方法确定的特征值的特征向量。1、幂法的迭代格式与收敛性质设n阶矩阵A的特征值,，…,是按绝对值大小编号的,ｘ（i＝1，2，…,n）为对应的特征向量，且为单根，即|｜>｜｜≥…≥｜｜则计算最大特征值与特征向量的迭代格式为v=Ａu,m=mａx（v），ｕ=v/m（1）其中ｍax(ｖ）表示向量v绝对值的最大分量.2、对于幂法的定理按式（1）计算出ｍ和u满足m=，u＝（二）反幂法算法的理论依据及推导反幂法是用来计算绝对值最小的特征值忽然相应的特征向量的方法。是对幂法的修改，可以给出更快的收敛性.1、反幂法的迭代格式与收敛性质设A是非奇异矩阵，则零不是特征值,并设特征值为｜|≥｜|≥…≥｜｜>|｜则按Ａ的特征值绝对值的大小排序，有｜|＞||≥…≥|｜对A实行幂法，就可得A的绝对值最大的特征值1／和相应的特征向量，即A的绝对值最小的特征值和相应的特征向量。由于用A代替A作幂法计算,因此该方法称为反幂法,反幂法的迭代格式为v=Aｕ,ｍ＝maｘ(ｖ),u=ｖ/m(２）2、对于反幂法的定理按式（2）计算出的m和ｕ满足:m＝，u=在式(2）中,需要用到A,这给计算带来很大的不方便，因此，把（２)式的第一式改为求解线性方程组Av=u(3）但由于在反幂法中，每一步迭代都需求解线性方程组（3)式,迭代做了大量的重复计算,为了节省工作量，可事先把矩阵A作LＵ分解,即A＝LU所以线性方程组(3）改为Ｌy=ｕ，Uv=y四、算法程序设计代码幂法程序,在ｍatlａb中建立一个M文件并保存。%pｏｗ。mfuｎctiｏn[m，u,indeｘ，k]=pow(A，u，ep，it_ｍax)ifｎａrgｉn<4iｔ＿mａx=1000；enｄiｆnargiｎ<3ep=１ｅ—5；endn=lengｔh（Ａ）；iｎｄｅx=0;k=0；m1=0;ｍ０=０;I=eyｅ（ｎ）;Ｔ＝A—m０＊I;whilｅｋ<=it_ｍaｘv=Ｔ*u；[vmax,i］＝ｍax(ａbs(v)）；m＝ｖ(i);u＝v/m；ifabs(m—m１）<ep;ｉndex=１；bｒeak；ｅndm＝ｍ+m0；m1＝m；ｋ=k+1；enｄ在mａtlａｂ输入面板，输入A=rａnｄ(4）；%产生一个4维随机矩阵B=A+A’;u＝［111１］’;％设立初始向量［m，u，indeｘ，k］=pow（B，ｕ,eｐ，iｔ＿ｍax)%最多可省略2个参数程序结束。在M文件中可以通过改变m0的值改变原点位移，从而达到原点位移加速.反幂法程序设计代码:在matlab中建立一个M文件并保存。%ｐｏw_inｖ.mｆunction［m，ｕ，ｉnｄex,k］=pow_ｉnv（A，u,eｐ,iｔ_mａx）ｉfnargｉn＜4it＿ｍax=1000；endifnarｇｉｎ<3ep＝1e-5；enｄn=ｌｅngtｈ（A)；indｅｘ＝0；ｋ＝0;m1＝０；m０=０;I＝ｅye(ｎ);T=A-m０＊I；iｎｖＴ=iｎｖ(T)；whilek<=it_ｍaxｖ=invT＊u；［ｖｍaｘ,i］=maｘ(aｂs（v））;m=v(i）；ｕ=v/ｍ;ifaｂs（ｍ－m1)＜epindex=1;ｂｒeａk;ｅndm１=m；k=ｋ+1；eｎdm=1/ｍ;ｍ=m＋ｍ0；在mａｔlab输入面板,输入Ａ=rａnd（4）；％产生一个4维随机矩阵B=A+A’；u=［１１11］’；％设立初始向量［m，u,ｉndｅx,k]=pｏw＿inv（B，u，ep,ｉt＿maｘ）％最多可省略2个参数程序结束.在Ｍ文件中可以通过改变m０的值改变原点位移，从而达到原点位移加速。【结果显示】%在M０=1e—4>〉B=ranｄ（4）；＞〉A=B＋B’A=0。２６750.57760.６34４1.３１３00.５7761。15030.7６410.1３670．６3440．764１0．02５70.４１931．３1300。13６7０.41９31.２248＞>ｕ=［1111］＇；>〉［m,u，inｄｅx,ｋ］=pow(A,u)m＝２.6813ｕ=０。857６0.6９340.５6231.0０00inｄex=1k＝４９修改M０＝１e-3m=2．６8１４ｕ=0．8576０.6９３40．5623１.０000indｅｘ＝０k=１0０1修改Ｍ０=０％此时为幂法m=2.６815ｕ＝0．8５７60．69350。562３1．0000ｉｎdｅx=１k=10修改U=[1２34］修改M0=1e-4m＝２。6８13u=0。85760。69340.５6231．0０0０ｉnｄex＝１ｋ=9修改Ｍ０=１e-3m＝2．6805u=0。85７60。69３40.56221.00００ｉndｅx=１k=7修改M０=0ｍ=2．６８14u＝0．8５７６0.6９３40。5623１。０000indｅx=1k=9修改U＝［35６7］修改M0=1e-4m＝2.6819u=0．857７0。６9370.56241．０000ｉnｄex＝1k＝7修改M0＝１e—3m=2.６814ｕ=０。85760．69３40.562３１．000０indｅx=0k＝１0０1修改M0=0ｍ=２.6８２０u=０。8５770.６937０.56241.0000ｉｎdex=1k=7总结以上，幂法如下：Um０ｍuindeｘk[１111]0．00012.681３［0.８576０。6９３4０．56２31。0000］14９０.0012。6814[0.５8760。6９3４0。５6２3１．0000］０1０0１０２。68１5［0.８5760。6９３50。56231．０0０0]１10［1234］0．00012。６81３［0。857６0.６9３40．5623１.0０00]１90．0012．680５［0。85760。６９340。56221。０0０0］1702．6814[０。８５760。69340.56２31。0000]１9［3567］0。０0０１2.68１9［0。85770．6937０.５６２41.0000]1７０。0012．6９１４[０。85760．69340.５６231．０000］01０0102．692[0.8５７70.６９3７０．５6２41。０0０0]17反幂法结果显示：在ｍ０为0时Ｍ0=0.0０１Ｕ=[11１1］Ｍ0=0。１u＝［1111］M0=0u=［1357]M0=0．1u＝[１35７］Ｍ0=0.5u＝［1357］Ｍ０＝0u=［２345]M０＝0。１u=[2345］M０=0．7ｕ=［2345]综上，反幂法结果如下：um0mｕiｎdｅxk［１111］0．10.3847[-0．8９９６1．00000.27２6－0。２364］1150．0010．３847[-0。8９9６1．00000.２7２６—0.2364］1１60０.3８47［-０。８９9６1.000０0。272６-０.2364]１１6［13５7］0．50.3847［—0．899５１.00００0.2726—0.2364］1270。1０。38４７［—0.89９61。０0０00.2７26—0。23６4］117０0。3847［—0.８９961.00000。2７2６－０。2364］120［23４5］0。70．７0９1[－０。6962-0.４４９7０。2