冪法和反冪法求矩陣特征值課程

20/20題目?jī)绶ê头磧绶ㄇ缶仃囂卣髦稻唧w內(nèi)容隨機(jī)產(chǎn)生一對(duì)稱矩陣,對(duì)不同的原點(diǎn)位移和初值(至少取3個(gè)）分別使用冪法求計(jì)算矩陣的主特征值及主特征向量，用反冪法求計(jì)算矩陣的按模最小特征值及特征向量,并比較不同的原點(diǎn)位移和初值說(shuō)明收斂。要求1。認(rèn)真讀題，了解問(wèn)題的數(shù)學(xué)原形；２．選擇合適問(wèn)題求解的數(shù)值計(jì)算方法；３.設(shè)計(jì)程序并進(jìn)行計(jì)算；4．對(duì)結(jié)果進(jìn)行解釋說(shuō)明;采用方法及結(jié)果說(shuō)明對(duì)于冪法和反冪法求解矩陣特征值和特征向量的問(wèn)題將從問(wèn)題分析，算法設(shè)計(jì)和流程圖，理論依據(jù)，程序及結(jié)果進(jìn)行闡述該問(wèn)題.一．問(wèn)題的分析：求ｎ階方陣A的特征值和特征向量，是實(shí)際計(jì)算中常常碰到的問(wèn)題,如：機(jī)械、結(jié)構(gòu)或電磁振動(dòng)中的固有值問(wèn)題等。對(duì)于ｎ階矩陣A,若存在數(shù)和n維向量x滿足Ax=x(1）則稱為矩陣A的特征值,x為相應(yīng)的特征向量。由高等代數(shù)知識(shí)可知，特征值是代數(shù)方程｜Ｉ-A｜=+a+…＋a+a=０（２）的根。從表面上看,矩陣特征值與特征向量的求解問(wèn)題似乎很簡(jiǎn)單，只需求解方程(2)的根,就能得到特征值，再解齊次方程組（I-Ａ）x=0（3）的解，就可得到相應(yīng)的特征向量。上述方法對(duì)于ｎ很小時(shí)是可以的.但當(dāng)n稍大時(shí)，計(jì)算工作量將以驚人的速度增大，并且由于計(jì)算帶有誤差,方程(2）未必是精確的特征方程，自然就不必說(shuō)求解方程（2)與（3）的困難了.冪法是一種計(jì)算矩陣主特征值（矩陣按模最大的特征值）及對(duì)應(yīng)特征向量的迭代方法，特別是用于大型稀疏矩陣。反冪法是計(jì)算海森伯格陣或三角陣的對(duì)應(yīng)一個(gè)給定近似特征值的特征向量的有效方法之一。二。算法設(shè)計(jì)及流程圖1、冪法算法（1)取初始向量u（例如?。?(1,1,…1))，置精度要求，置k=1.（2）計(jì)算ｖ=Au,m=max（v)，ｕ=ｖ/ｍ(3)若｜m=m|〈，則停止計(jì)算(m作為絕對(duì)值最大特征值,ｕ作為相應(yīng)的特征向量)否則置k=k+１,轉(zhuǎn)(２）2、反冪法算法(1）取初始向量ｕ(例如取u＝（1，1，…1）)，置精度要求，置k=1。（2）對(duì)A作LU分解，即A=LU(３）解線性方程組Ly=ｕ，Uv＝y(tǒng)(４）計(jì)算ｍ=maｘ（v)，u=v/m（5）若|ｍ=ｍ|〈,則停止計(jì)算（1/m作為絕對(duì)值最小特征值，u作為相應(yīng)的特征向量）；否則置k＝k+１，轉(zhuǎn)(3).冪法流程圖：開(kāi)始開(kāi)始輸入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結(jié)束m1=m;k=k+1反冪法流程圖開(kāi)始開(kāi)始輸入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結(jié)束m1=m;k=k+1輸入A；[m,u,index]=pow(A,1e-6)三、算法的理論依據(jù)及其推導(dǎo)（一）冪法算法的理論依據(jù)及推導(dǎo)冪法是用來(lái)確定矩陣的主特征值的一種迭代方法，也即,絕對(duì)值最大的特征值。稍微修改該方法，也可以用來(lái)確定其他特征值。冪法的一個(gè)很有用的特性是它不僅可以生成特征值，而且可以生成相應(yīng)的特征向量。實(shí)際上，冪法經(jīng)常用來(lái)求通過(guò)其他方法確定的特征值的特征向量。1、冪法的迭代格式與收斂性質(zhì)設(shè)n階矩陣A的特征值,，…,是按絕對(duì)值大小編號(hào)的,ｘ（i＝1，2，…,n）為對(duì)應(yīng)的特征向量，且為單根，即|｜>｜｜≥…≥｜｜則計(jì)算最大特征值與特征向量的迭代格式為v=Ａu,m=mａx（v），ｕ=v/m（1）其中ｍax(ｖ）表示向量v絕對(duì)值的最大分量.2、對(duì)于冪法的定理按式（1）計(jì)算出ｍ和u滿足m=，u＝（二）反冪法算法的理論依據(jù)及推導(dǎo)反冪法是用來(lái)計(jì)算絕對(duì)值最小的特征值忽然相應(yīng)的特征向量的方法。是對(duì)冪法的修改，可以給出更快的收斂性.1、反冪法的迭代格式與收斂性質(zhì)設(shè)A是非奇異矩陣，則零不是特征值,并設(shè)特征值為｜|≥｜|≥…≥｜｜>|｜則按Ａ的特征值絕對(duì)值的大小排序，有｜|＞||≥…≥|｜對(duì)A實(shí)行冪法，就可得A的絕對(duì)值最大的特征值1／和相應(yīng)的特征向量，即A的絕對(duì)值最小的特征值和相應(yīng)的特征向量。由于用A代替A作冪法計(jì)算,因此該方法稱為反冪法,反冪法的迭代格式為v=Aｕ,ｍ＝maｘ(ｖ),u=ｖ/m(２）2、對(duì)于反冪法的定理按式（2）計(jì)算出的m和ｕ滿足:m＝，u=在式(2）中,需要用到A,這給計(jì)算帶來(lái)很大的不方便，因此，把（２)式的第一式改為求解線性方程組Av=u(3）但由于在反冪法中，每一步迭代都需求解線性方程組（3)式,迭代做了大量的重復(fù)計(jì)算,為了節(jié)省工作量，可事先把矩陣A作LＵ分解,即A＝LU所以線性方程組(3）改為Ｌy=ｕ，Uv=y四、算法程序設(shè)計(jì)代碼冪法程序,在ｍatlａb中建立一個(gè)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ａt(yī)lａｂ輸入面板，輸入A=rａnｄ(4）；%產(chǎn)生一個(gè)4維隨機(jī)矩陣B=A+A’;u＝［111１］’;％設(shè)立初始向量［m，u，indeｘ，k］=pow（B，ｕ,eｐ，iｔ＿ｍax)%最多可省略2個(gè)參數(shù)程序結(jié)束。在M文件中可以通過(guò)改變m0的值改變?cè)c(diǎn)位移，從而達(dá)到原點(diǎn)位移加速.反冪法程序設(shè)計(jì)代碼:在matlab中建立一個(gè)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）；％產(chǎn)生一個(gè)4維隨機(jī)矩陣B=A+A’；u=［１１11］’；％設(shè)立初始向量［m，u,ｉndｅx,k]=pｏw＿inv（B，u，ep,ｉt＿maｘ）％最多可省略2個(gè)參數(shù)程序結(jié)束.在Ｍ文件中可以通過(guò)改變m０的值改變?cè)c(diǎn)位移，從而達(dá)到原點(diǎn)位移加速。【結(jié)果顯示】%在M０=1e—4>〉B=ranｄ（4）；＞〉A(chǔ)=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修改Ｍ０=０％此時(shí)為冪法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總結(jié)以上，冪法如下：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反冪法結(jié)果顯示：在ｍ０為0時(shí)Ｍ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]綜上，反冪法結(jié)果如下：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