數(shù)值分析冪法和反冪法

冪法和反冪法求矩陣特征值隨機(jī)產(chǎn)生一對(duì)稱矩陣，對(duì)不同的原點(diǎn)位移和初值（至少取3個(gè)）分別使用冪法求計(jì)算矩陣的主特征值及主特征向量，用反冪法求計(jì)算矩陣的按模最小特征值及特征向量，并比較不同的原點(diǎn)位移和初值說明收斂。認(rèn)真讀題，了解問題的數(shù)學(xué)原形;選擇合適問題求解的數(shù)值計(jì)算方法;設(shè)計(jì)程序并進(jìn)行計(jì)算;對(duì)結(jié)果進(jìn)行解釋說明;采采用方法及結(jié)果說明對(duì)于冪法和反冪法求解矩陣特征值和特征向量的問題將從問題分析，算法設(shè)計(jì)和流程圖，理論依據(jù)，程序及結(jié)果進(jìn)行闡述該問題。一.問題的分析：求n階方陣A的特征值和特征向量，是實(shí)際計(jì)算中常常碰到的問題，如：機(jī)械、結(jié)構(gòu)或電磁振動(dòng)中的固有值問題等。對(duì)于n階矩陣A,若存在數(shù)人和n維向量x滿足TOC\o"1-5"\h\zAx=人x（1）則稱人為矩陣A的特征值，x為相應(yīng)的特征向量。由高等代數(shù)知識(shí)可知，特征值是代數(shù)方程|人I-A|=人n+a人n-i+—+a人+a=0（2）的根。從表面上看，矩陣特征值與特征向量的求解問1題似乎很簡單，只需求解方程（2）的根，就能得到特征值人，再解齊次方程組（人I-A）x=0（3）的解，就可得到相應(yīng)的特征向量。上述方法對(duì)于n很小時(shí)是可以的。但當(dāng)n稍大時(shí)，計(jì)算工作量將以驚人的速度增大，并且由于計(jì)算帶有誤差，方程（2）未必是精確的特征方程，自然就不必說求解方程（2）與（3）的困難了。冪法是一種計(jì)算矩陣主特征值（矩陣按模最大的特征值）及對(duì)應(yīng)特征向量的迭代方法，特別是用于大型稀疏矩陣。反冪法是計(jì)算海森伯格陣或三角陣的對(duì)應(yīng)一個(gè)給定近似特征值的特征向量的有效方法之一。二.算法設(shè)計(jì)及流程圖1、冪法算法（1）取初始向量u（o）（例如取u（o）=（1,1,...1）Q,置精度要求e，置k=1.（2）計(jì)算v（k）=Au（k-i）,m=max（v（k））,u（k）=v（k）/m（3）若|m廣m.1|<e,則停止計(jì)算（m.作為絕對(duì)值最大特征值人1,u俄）作為相應(yīng)的特征向量）否則置k=k+1，轉(zhuǎn)（2）2、反冪法算法（1）取初始向量U（0）（例如取U（0）^1,1,…1）T），置精度要求E，置k=1.（2）對(duì)A作LU分解，即A=LU（3）解線性方程組Ly（k）=u（k-i）,Uv（k）=y（k）（4）計(jì)算m=max（v（k））,u（k）=v（k）/m（5）若|m廣m；1|<e，則停止計(jì)算（1/m「作為絕對(duì)值最小特征值牛u（k）作為相應(yīng)的特征向量）；否則置k=k+1,轉(zhuǎn)（3）."冪法流程圖:開始m1=m;k=k+1結(jié)束反冪法流程圖輸入入；[m,u,index]=pow(A,1e-6)k=0;m1=rv=A*ur[vmax,i]=max(abs(v))1rm=v(i);u=v/m一abs(m-m1)Z1e-e'''''''''''vindex=1;break;輸出：m，u，indexabs(m-m1)<1e-6開始k=0;m1=0v=invA*um1=m;k=k+1輸入A；結(jié)束（一）冪法算法的理論依據(jù)及推導(dǎo)index=1;break;輸出：m=v(i);u=v/mm，abs(m-m1)<1e-6開始k=0;m1=0v=invA*um1=m;k=k+1輸入A；結(jié)束（一）冪法算法的理論依據(jù)及推導(dǎo)index=1;break;輸出：m=v(i);u=v/mm，u，index三、算法的理論依據(jù)及其推導(dǎo)[vmax,i]=max(abs(v))輸入A；[m,u,index]=pow_inv(A,1e-6)冪法是用來確定矩陣的主特征值的一種迭代方法，也即，絕對(duì)值最大的特征值。稍微修改該方法，也可以用來確定其他特征值。冪法的一個(gè)很有用的特性是它不僅可以生成特征值，而且可以生成相應(yīng)的特征向量。實(shí)際上，冪法經(jīng)常用來求通過其他方法確定的特征值的特征向量。1、冪法的迭代格式與收斂性質(zhì)設(shè)n階矩陣A的特征值％，人2，…,人〃是按絕對(duì)值大小編號(hào)的，x‘（iF,2，…,n）為對(duì)應(yīng)人.的特征向量，且七為單根，即IX|>|人|2,??2|人|則計(jì)算最大特征值與特征向量的迭代格式為v（k）=Au（k-i）,m=max（v（k））,u（k）=v（k）/m⑴上'其中max（v（k））表示向量v（k）絕對(duì)值的最大分量。2、對(duì)于冪法的定理按式（1）計(jì)算出mk和u（k）滿足limm=X,limu（k）=——ik—Xk—X（二）反冪法算法的理論依據(jù)及推導(dǎo)反冪法是用來計(jì)算絕對(duì)值最小的特征值忽然相應(yīng)的特征向量的方法。是對(duì)冪法的修改，可以給出更快的收斂性。1、反冪法的迭代格式與收斂性質(zhì)設(shè)A是非奇異矩陣，則零不是特征值，并設(shè)特征值為|X|n|x|n，??n|x|>|X|TOC\o"1-5"\h\z則按A-i的特征值絕對(duì)'值的大小排序，有“^|:|>|《|>W|1|nn一11對(duì)A-i實(shí)行冪法，就可得A-i的絕對(duì)值最大的特征值1/Xn和相應(yīng)的特征向量，即A的絕對(duì)值最小的特征值和相應(yīng)的特征向量。n由于用A-i代替A作冪法計(jì)算，因此該方法稱為反冪法，反冪法的迭代格式為v（k）=A-iu（k-i）,m=max（v（k））,u（k）=v（k）/m（2）2、對(duì)于反冪法的定理*按式（2）計(jì)算出的mk和u（k）滿足：1xk-*k廠'k-*max（x）nn在式（2）中，需要用到A-i，這給計(jì)算帶來很大的不方便，因此，把（2）式的第一式改為求解線性方程組Av（k）=u（k-i）（3）但由于在反冪法中，每一步迭代都需求解線性方程組（3）式，迭代做了大量的重復(fù)計(jì)算，為了節(jié)省工作量，可事先把矩陣A作LU分解，即A=LU所以線性方程組（3）改為Ly（k）=u（k-i）,Uv（k）=y（k）四、算法程序設(shè)計(jì)代碼冪法程序，在matlab中建立一個(gè)M文件并保存。%pow.mfunction[m,u,index,k]=pow(A,u,ep,it_max)ifnargin<4it_max=1000;endifnargin<3ep=1e-5;endn=length(A);index=0;k=0;m1=0;m0=0;I=eye(n);T=A-m0*I;whilek<=it_maxv=T*u;[vmax,i]=max(abs(v));m=v(i);u=v/m;ifabs(m-m1)<ep;index=1;break;endm=m+m0;m1=m;k=k+1;end在matlab輸入面板，輸入A=rand(4)；%產(chǎn)生一個(gè)4維隨機(jī)矩陣B=A+A’；u=[1111]’；％設(shè)立初始向量[m,u,index,k]=pow(B，u,ep,it_max)%最多可省略2個(gè)參數(shù)程序結(jié)束。在M文件中可以通過改變m0的值改變?cè)c(diǎn)位移，從而達(dá)到原點(diǎn)位移加速。反冪法程序設(shè)計(jì)代碼：

在matlab中建立一個(gè)M文件并保存。%pow_inv.mfunction[m,u,index,k]=pow_inv(A,u,ep,it_max)ifnargin<4it_max=1000;endifnargin<3ep=1e-5;endn=length(A);index=0;k=0;m1=0;m0=0;I=eye(n);T=A-m0*I;invT=inv(T);whilek<=it_maxv=invT*u;[vmax,i]=max(abs(v));m=v(i);u=v/m;ifabs(m-m1)<epindex=1;break;endm1=m;k=k+1;endm=1/m;m=m+m0;在matlab輸入面板，輸入A=rand（4）；%產(chǎn)生一個(gè)4維隨機(jī)矩陣B=A+A’；u=[1111]’；％設(shè)立初始向量[m,u,index,k]=pow_inv(B，u,ep,it_max）%最多可省略2個(gè)參數(shù)[m,u,index,k]=pow_inv(B，u,ep,it_max）%最多可省略2個(gè)參數(shù)程序結(jié)束。在M文件中可以通過改變m0的值改變?cè)c(diǎn)位移，從而達(dá)到原點(diǎn)位移加速。【結(jié)果顯示】%在M0=1e-4>>B=rand(4);>>A=B+B’A=0.26750.57760.63441.31300.57761.15030.76410.13670.63440.76410.02570.41931.31300.13670.41931.2248>>u=[1111]';>>[m,u,index,k]=pow(A,u)m=2.6813u=0.85760.69340.56231.0000index1=k=49修改M0=1e-3m=2.6814u=0.85760.6934

0.56231.0000二k=1001修改M0=0%此時(shí)為冪法m=2.6815u=0.85760.69350.56231.0000inTk=10修改U=[1234]修改M0=1e-4m=2.6813u=

0.85760.69340.56231.0000k=9修改M0=1e-3m=2.6805u=0.85760.69340.56221.0000inTk=修改M0=0m=2.6814u=0.85760.6934

0.56231.0000index=k=9修改U=[3567]修改M0=1e-4m=2.6819u=0.85770.69370.56241.0000ind:k=修改M0=1e-3m=2.6814=xapui0000'I死9如0ZG69'0ZZG80=n0源9?乙=in0=0/現(xiàn)割1001=習(xí)0=xapui0000'I￡乙9如0州69'09藥80

Um0muindexk[0.85760.69340.5623[1111]0.00012.68131.0000][0.58760.69340.56231490.0012.68141.0000][0.85760.69350.56230100102.68151.0000][0.85760.69340.5623110[1234]0.00012.68131.0000][0.85760.69340.5622190.0012.68051.0000][0.85760.69340.56231702.68141.0000][0.85770.69370.562419[3567]0.00012.68191.0000][0.85760.69340.5623170.0012.69141.0000][0.85770.69370.56240100102.692反冪法結(jié)果顯示：在m0為0時(shí)1.0000]17

>>A=rand(4)A=0.81800.34120.8385□.54660.66020.53410.56810.44490.34200.72710.3704□.63460.28970.30930.T027□.6213?E=A+A'B=1.63591.00141.18050.83631.00141.Ot5S21.Z9520.75421.18051.29520.MOB1.33730.83630.TC421.39F31.242G>>u=[111I]311—1111>>.JiijUjiiideK,k]=pow_inv(Bjle—4)m=0.3037u=-0.54751.CICICICI0.12S6-0.5040ind.es=1k二16M0=0.001U=[1111]1.63591.00141.1805□.83631.00141.06S21.2952□.5M21.18051.2952□.fias1.53730.83630.57421.39731.2426?u=[l11iru=1111>>[n,UjindeKjk]=povr_inv(BjLe-4)n=0.3847u=-0.8996i.oooa0.272S-0.2364index=1k=1.63591.00141.1805□.83631.00141.06821.2952□.5M21.1S051.29520.74081.53730.S3630.5M21.33^31.2426?n=[l11ir1111[iRjUjLndeXjk]=pov_mv(Bj1e-4)n=0.3847u=-0.89951.00000.2725-0.2364irides=1B=1.C3591.00141.18050.83631.C014i.oess1.29520.57421.18051.29520.74081.39730.83630.57421.39731.2426?u.=[l357]nu=1357>>[ikjv.sindeZjk]=pow_inv（瓦u）m=0.384?u=-0.89961.00000.2726-0.2354iniex-1B=1.63591.00141.1S050.03631.00141.06S21.29520.57421.10051.29520.74081.39730.83630.57421.39731.2426?u=[l35r],u一1357>>[iRj烏indeXjk]=pow_inv(BJu)m=0.384?u=-0.89951.00000.2726-0.2364index=1B=1.63591.00141.10050.S3631.00141.06821.29520.57421.1S051.29520.740S1.39730.83630.57421.39731.2426?u.=[l357]nu=1357?[iRjUjindeKjk]=pmv(Bju)in=0.3347-0.09951.00000.2F26-0.2364indeK=1B=1.63591.OOM1.18050.83631.00141.06821.29520.57421.18051.29520.74081.39730.83C30.57421.39731.2426?u=[2345],u-234-EiiLjUjinde^jk]=pov_inv(Bju)0.3847-0.89961.00000.2726-0.23C4indeK=1k二1.63591.00141.18050.S3631.00141.06321.29520.57421.18051.29E20.74081.39730.03630.57421.39731.2426?34S]‘u=2345>>[mjUjindsKjk]=pcw_inT7(BJn=0.3047u=-0.B9951.00000.2726-0.2364indeK=1k=1.63591.00141.1S050.03631.00141.06821.29520.57421.18051.29520.74081.39730.83630.574