一求解一般右端項(xiàng)的Toeplitz方程組及Toeplitz矩陣的逆-Read

一求解一般右端項(xiàng)的103「1讓方程組及Toeplitz矩陣的逆1實(shí)驗(yàn)?zāi)康?熟悉求解Yule-Walker方程組，一般右端項(xiàng)的Toeplit方程組及Toeplit矩陣的逆的求解步驟；?通過實(shí)驗(yàn)體會方程組的性質(zhì)對求解的影響；2實(shí)驗(yàn)原理1)求解一般右端項(xiàng)的Toeplit方程組：Tnx^其中匕(1, J是已知向量.通過f/ 、 ,一~ 、 xEy求解Tx(, )t,k1,2,,n)來計(jì)算*可得；xkkkk,其kk1k k1 k1k中；yk是k階Yule-Walker方程組的解， Ek是k階單位陣，k由k(k1；^/「/(1確定，這里,%,‘)丁.這樣，便可通過弋出發(fā)遞推地求得方程組的解。ToeplitZ郵車的逆以n5為例)通過求解一個(gè)5階的Yule-Walke方程組得到工，再利用： 1/(1R1T11),v Enj”1求出、1的最后一列和最后一行元素，然后再利用、1的廣義對稱性求出相應(yīng)元素；利用ij ……(vivjvnvnj)/及(1得到的結(jié)果，然后再利用、1的廣義對稱性，可得到丁/第二層的全部元素；再利用ij …ni(vivjvnvnj)/,可求得T51最中間的元素，則就求得了T51的全部元素。3數(shù)值算例①求解Yule-Walked程組：(求解四階Yule-Walker方程組R1[1]保存方程組的階數(shù))pleaseinputtheflag1to3:1pleaseinputthestringR1[N]:4-138R1[0]=5.000000R1[1]=4.000000R1[2]=-1.000000R1[3]=3.000000R1[4]=8.000000結(jié)果y1[0]=5.000000y1[1]=-0.786753y1[2]=1.675522y1[3]=2.357687y1[4]=-2.539950②求解Toeplit矩陣的逆（5階）pleaseinputtheflag1to3:2pleaseinputthestringR1[N]:0.41.23.41.2R1[0]=5.000000R1[1]=0.400000R1[2]=-1.200000R1[3]=3.400000R1[4]=1.20000得到的結(jié)果為T[0][0]=-0.089410T[0][1]=-0.015410T[0][2]=0.119910T[0][3]=0.319410T[0][4]=-0.111810T[1][0]=-0.015410T[1][1]=1.137010T[1][2]=-0.205710T[1][3]=-1.434810T[1][4]=0.319410T[2][0]=0.119910T[2][1]=-0.205710T[2][2]=-1.732610T[2][3]=0.205710T[2][4]=0.119910T[3][0]=0.319410T[3][1]=-1.434810T[3][2]=-0.205710T[3][3]=1.137010T[3][4]=0.015410T[4][0]=-0.111810T[4][1]=0.319410 T[4][2]=0.119910T[4][3]=0.015410T[4][4]=-0.089410③一般右端項(xiàng)的Toeplit方程組：（4階右端項(xiàng)的Toeplit方程組R1[1]保存方程組的階數(shù)）pleaseinputtheflag1to3:3pleaseinputthestringR矩陣系數(shù)”3416R1[N]R1[0]=5.000000R1[1]=3.000000R1[2]=4.000000R1[3]=1.000000R1[4]=6.000000pleaseinputthestringb[N]:6724b[N]右端項(xiàng)”b[0]=5.000000b[1]=6.000000b[2]=7.000000b[3]=2.000000b[4]=4.000000結(jié)果：x[0]=5.000000x[1]=6.000000x[2]=1.375000x[3]=4.512821x[4]=-0.100273將上面結(jié)果與直接用乂@埒@匕計(jì)算的結(jié)果比較誤差較大。4心得體會本次試驗(yàn)我是用c語言完成的，實(shí)驗(yàn)過程中，，正是因?yàn)門oeplit方程組本身具有的良好的性質(zhì):廣義對稱性.使得在編寫程序的過程中節(jié)省了很多步驟，整過過程求解的過程即用到向量乘法，這只要兩輪循環(huán)就可以實(shí)現(xiàn)，用c語言很容易實(shí)現(xiàn)的。，特別是在求逆的過程中,因?yàn)樗膶ΨQ性，我們可以僅僅去求它的一個(gè)倒三角矩陣,這樣在一個(gè)循環(huán)中用三條賦值語句就可以實(shí)現(xiàn)。因此，方程組或者矩陣的某些方面的性質(zhì)會對我們的計(jì)算結(jié)果，計(jì)算方法產(chǎn)生很大的作用。然而，由于自己對算法的不熟悉和c與語言掌握得不精通。在向量運(yùn)算的過程中對數(shù)組中元素與腳標(biāo)的關(guān)系的處理過程中遇到很多麻煩，對結(jié)果的保存處理得不當(dāng)。我本想用malloc函數(shù)開辟內(nèi)存區(qū)以保存結(jié)果，但該函數(shù)在wintC-191下無法實(shí)現(xiàn)，最后用全局變量保存結(jié)果，這使得函數(shù)的通用性降低了，而且將得了程序的清晰性。我會在以后學(xué)習(xí)的過程中多請教老師多和同學(xué)交流，提高自己的編程技巧和運(yùn)用知識的能力。二實(shí)用共軛梯度法的求解線性方程組的算法實(shí)現(xiàn)及其分析1實(shí)驗(yàn)?zāi)康模菏煜?shí)用共軛梯度法的基本原理及思想方法；通過與一般方法的比較，體會實(shí)用共輾梯度法的優(yōu)劣；掌握稀疏矩陣存儲方式；2實(shí)驗(yàn)原理：實(shí)用共輾梯度法是對系數(shù)矩陣為對稱正定的一類線性方程的求解的一種非常有效的算法。在實(shí)際使用時(shí)，由于誤差的出現(xiàn)，吏得負(fù)梯度方向r的正交性很快損失，以致于其有限k步終止性不再成立.因此，實(shí)際應(yīng)用時(shí)，只是將共軛梯度法作為一種迭代的方法使用，其基本步驟是：(1)首先輸入二次發(fā)國:2XTAxbTX中的A,b及初始點(diǎn)X。，注意這里A應(yīng)為對稱正定的稀疏矩陣；(2)令xx,rAx, rrr,k1；(3)如果k1^^pr；否則：k1/k2，Prp；wAp,a /pTw,xxap,rraw, rrr(5)如果k0，則輸出x；否則：kk1，轉(zhuǎn)步到(2),實(shí)際上，每迭代一次，只需作依次矩陣和向量乘法運(yùn)算和5n次數(shù)量乘法運(yùn)算，在存儲上僅需增加4n個(gè)存儲單元即可。3數(shù)值算例①A為一般的小型對稱正定矩陣3002050090A051000001290021寫成稀疏矩陣的形式如下：112234451523345539251121選取常量如下0000x00e0.001006b543計(jì)算結(jié)果如下：x=CG_4_1(B,b,x0,e);-0.266670.82609x0.869572.26670.86667結(jié)果分析：(wucha=A*x-b)wucha=1.0e-013*0.3109-0.0444-0.04440.05330.1599從上述的計(jì)算結(jié)果來看，效果還比較理想。②A為Hilberj^^令H=hilb(5；H=1.00000.50000.33330.25000.20000.50000.33330.25000.20000.16670.33330.25000.20000.16670.14290.25000.20000.16670.14290.12500.20000.16670.14290.12500.1111將其轉(zhuǎn)化為稀疏矩陣的表現(xiàn)形式n=length(H);m=n*(n+1)/2;B=zeros(m,3);k=1;whilek<=mfori=1:nforj=i:nB(k,2)=j;B(k,1)=i;B(k,3)=H(i,j);k=k+1;endendendB=1.00001.00001.00001.00002.00000.50001.00003.00000.33331.00004.00000.25001.00005.00000.20002.00002.00000.33332.00003.00000.25002.00004.00000.20002.00005.00000.16673.00003.00000.20003.00004.00000.16673.00005.00000.14294.00004.00000.14294.00005.00000.12505.00005.00000.11110.0001，選取右端常數(shù)為b[1,2,3,4,5]，初始值x0[0,0,0,0,0]，以及誤差控制量0.0001，得到如下的結(jié)果：x=1.0e+004*0.0125-0.2880 1.4490-2.4640 1.3230為了進(jìn)行比較，用乂@埒@匕直接算出了(這里乂乂inv(H)*b)xx=1.0e+004*0.0125-0.2880 1.4490-2.4640 1.3230計(jì)算二者的殘差，得到 (這里。。xxx)cc=1.0e-004*0.04530.1679 0.0977 0.0697 0.05650.0453從結(jié)果看來已經(jīng)達(dá)到了所要求的精度要求了。4心得體會共軛梯度法是解大型稀疏現(xiàn)行方程組的有效方法。通過學(xué)習(xí)和上機(jī)試驗(yàn)，我更深刻共軛梯度法的優(yōu)點(diǎn)：由于它所處理的對象是稀疏矩陣，存儲時(shí)只要存儲不為零的元素即可；不需要預(yù)先估計(jì)別的參數(shù)就可以進(jìn)行計(jì)算；每次迭代所需的計(jì)算，主要是向量之間的運(yùn)算，便于并行化。這些優(yōu)點(diǎn)使得它的算法快速，并且其精確度相對也比較高。.同時(shí)，也發(fā)現(xiàn)即使是對Hilber矩陣這樣典型的病態(tài)方程組，運(yùn)用共軛梯度法還是可以得到較精確的解，這也是它優(yōu)于其他迭代法的地方。這是我第一次用乂@埒@匕編寫程序，在編寫程序的過程中，我充分感受到乂@埒@匕實(shí)現(xiàn)矩陣與向量運(yùn)算的優(yōu)越性。我先通過一般共軛梯度法的Matlab程序?qū)崿F(xiàn)小型矩陣的計(jì)算，在此過程中并未考慮存儲矩陣的方式，在實(shí)現(xiàn)實(shí)用共軛梯度法時(shí)對矩陣存儲的方式及轉(zhuǎn)化的過程遇到很多麻煩，其間我在圖書館查閱有關(guān)矩陣計(jì)算的乂@垣@匕程序?qū)嵗⑴c同學(xué)交流經(jīng)驗(yàn),特別是戴仕全同學(xué)幾次幫我修改程序并提出很多意見，讓我能在較短的時(shí)間內(nèi)運(yùn)用Matlab編寫程序。三雙重隱式QR迭代算法實(shí)現(xiàn)及其分析1實(shí)驗(yàn)?zāi)康模菏煜R算法求解矩陣特征值基本思想；編程實(shí)現(xiàn)雙重隱式QR迭代算法并比較其優(yōu)越性；熟悉卜。尤6卜。^6變換，Givens變換，和Hessenber吩解的基本思想及其應(yīng)用;2實(shí)驗(yàn)原理：QR迭代基本思想是:給定ARnn與AA,構(gòu)造迭代：%QkRk,k1,2,,1 ARQk1kk其中Qk是正交矩陣，Rk是上三角矩陣。對一般的方陣進(jìn)行一次QR迭代所需的運(yùn)算量是O(n3),這顯然不是我們所希望，因此需要在進(jìn)行迭代前對初始矩陣進(jìn)行一下處理，使之具有較多的零元素，這樣就可以使運(yùn)算量大大減少。這一步就稱之為上呢$$6廿6嗎化。主要是通過Householder變換實(shí)現(xiàn)這一過程。HQtAQ同時(shí)，用AI代替A來加快QR迭代的收斂速度，迭代為:Hk同時(shí)，用AI代替A來加快QR迭代的收斂速度，迭代為:HRQIk1kkH為了不引進(jìn)復(fù)運(yùn)算，用，連續(xù)作兩次位移，即進(jìn)行口H為了不引進(jìn)復(fù)運(yùn)算，用，連續(xù)作兩次位移，即進(jìn)行口11 2 H1H2IURRUI11 1IUR2 22R2U2 2I這里HHk，這樣就找到了一種實(shí)現(xiàn)由H到H2的變換方法，既減少了運(yùn)算量又防止了復(fù)運(yùn)算的出現(xiàn)，這就是雙重步位移的隱式QR算法的基本思想。3數(shù)值算例①用基本QR算法求解A矩陣的全部特征值:A=3234567111234562891234-4291113158-1-2-3-1-1-1-1323413158-2-2-3-4-5-3-3>>[iterD]=qralg(A)iter=537D=Columns1through318.4123 11.1805 1.7100-4.2522iColumns4through61.7100+4.2522i4.4982 -2.2327Column7-0.2783②用雙重隱式QR算法求解A矩陣的全部特征值：iter=9D=Columns1through318.4123 11.1805 1.7099-4.2522iColumns4through61.7099+4.2522i4.4983 -2.2327Column7-0.27834心得體會：QR算法是目前計(jì)算中小型稠密矩陣全部特征值和特征向量的最有效的方法?；綫R算法每次迭代都需要明顯地作矩陣的QR分解，其收斂是線性的，這種算法實(shí)在實(shí)數(shù)域內(nèi)進(jìn)行運(yùn)算，當(dāng)矩陣有復(fù)數(shù)特征值時(shí)，即使收斂也是很慢的，雙重隱式QR迭代算法其迭代過程都是實(shí)數(shù)運(yùn)算，原點(diǎn)位移加速了收斂性，隱式方法可減少計(jì)算工作量并節(jié)省存儲空間。QR算法包括著householde變換，上Hessenberg分解，QR變換，QR迭代收斂性判別等一系列問題，與前面兩個(gè)實(shí)驗(yàn)一比，難度當(dāng)然加大了，我先通過編寫基本QR算法程序，體會QR算法的原理和householde變換，化Hessenberg型方法，在編寫雙重隱式QR迭代算法對矩陣作初始正交相似變換時(shí)是最棘手的。隱式QR算法是通過劃分和收縮以及原點(diǎn)平移對QR算法計(jì)算工作量和收斂性加以改善。這也給了我們研究和改進(jìn)其他計(jì)算方法的一條思路，我會認(rèn)真體會它的思想，學(xué)以致用。四矩陣的奇異值分解1實(shí)驗(yàn)?zāi)康模豪斫饩仃嚻娈愔捣纸獾幕舅枷?；熟悉奇異值分解的?nèi)在規(guī)律；熟練矩陣的Schur分解和Householde逡換;2實(shí)驗(yàn)原理：設(shè)ARmn,則A的奇異值分解可通過可通過實(shí)對稱矩陣CAtA的Schur分解導(dǎo)出,因此可先利用對稱QR方法來實(shí)現(xiàn)C的Schur分解，然后借助C的Schur分解來實(shí)現(xiàn)A的奇異值分解.但是由于計(jì)算CAtA的過程中會引入較大的誤差，并且計(jì)算量也挺大，因此，通過隱含的將對稱QR方法應(yīng)用于AtA上來解決這個(gè)問題，主要分為兩步:(1)將A二對角化，這一步主要是通過Householder變換來解決；(2)對TBTB進(jìn)行帶聽出$。田立移的對稱QR迭代,當(dāng)然這一步還是希望不明確的將T計(jì)算出來就可進(jìn)行下去.這里又可分為三步：2r2r2第一步：選取位移，取門1n1nirn第二步:確定Givens變換GJ12,)；rn1n的靠近2r2的特征值作為2r2 nnnn第三步:確定正交矩陣Q使、丁GjTG'Q是對稱三對角陣目Qe]eiO這樣就得到了奇異值分解的最基本的迭代算法。3數(shù)值算例任做100*70的隨機(jī)矩陣，誤差設(shè)為1e-16>A=rand(100,70);>u=1e-16;[U,V,B,B1]=svd1(A,u);迭代次數(shù)：tt=1409計(jì)算誤差：>norm(U'*U-eye(100))ans=1.2424e-013>norm(V'*V-eye(70))ans=101.2569e-013>A2=U*A*V;>A1=A2(1:70,1:70);>norm(A1-B)ans=1.0442e-013由程序生成的B的對角線元素：>B1B1=Columns1through1642.0305 5.2185 4.9596 4.7958 4.7379 4.5545 4.4544 4.40184.3203 4.2068 4.1756 4.0837 3.9697 3.8604 3.8464 3.8079Columns17through323.6789 3.6041 3.5366 3.4952 3.4367 3.3254 3.3076 3.26773.2126 3.2078 3.0966 3.0310 2.9434 2.9119 2.8335 2.7346Columns33through482.7135 2.6312 2.5903 2.5415 2.4359 2.3976 2.3398 2.29792.2760 2.2214 2.2017 2.1386 2.0803 2.0066 1.9668 1.8415Columns49through641.8296 1.7528 1.7078 1.6529 1.5861 1.5548 1.5247 1.46351.3945 1.2653 1.2365 1.1640 1.0889 1.0669 1.0171 0.9268Columns65through700.8697 0.8462 0.6894 0.6695 0.5948 0.474011

五利用Lanczos算法求矩陣A的前10個(gè)和后10個(gè)特征值1實(shí)驗(yàn)?zāi)康?理解Lanczos算法的基本思想；了解Lanczos算法在計(jì)算大型稀疏對稱矩陣特征值問題上的優(yōu)勢;能編程實(shí)現(xiàn)Lanczos算法并用它計(jì)算大型稀疏對稱矩陣特征值；對大型稀疏對稱矩陣矩陣的存儲方法有更深入的理解和研究；2實(shí)驗(yàn)原理:Lanczos方法是適用于求解大型稀疏對稱矩陣特征值問題的比較理想的辦法。設(shè)ASRnn,則通常情況下，我們都將A三對角化來求其的特征值和特征向量，即求一個(gè)正交矩陣Q,使得丁QtAQ,因?yàn)橐浞掷肁的稀疏性以減少內(nèi)存的占用量，所以直接對T和Q進(jìn)行三對角化，這是Lanczos迭代的第一步，即產(chǎn)生一列對稱的三對角矩陣T^j1,2,,m；第二步：對某一個(gè)km,計(jì)算T的部分特征值或全部特征值；第三步：選擇這些特征值k中的某些個(gè)作為A的近似特征值；第四步：如果對應(yīng)的特征向量亦需要，則對第三步選定的每個(gè)特征值，求對應(yīng)的Ritz向量z,然后以2作為A對應(yīng)于的近似特征向量。3數(shù)值算例4 114 1,14 1144 114 1,14 114100100IHI計(jì)算：A （其中HIHIIH4040I為100階單位矩陣）的前10個(gè)及后10個(gè)特征值。取k=300,然后計(jì)算它的前十個(gè)特征值作為A的前十個(gè)特征值的近似.取口的第一個(gè)分量為1,則結(jié)果為：前10個(gè)特征值:4.9744545.09842814.84842 24.91109 34.9744545.09842854.72038 65.0367 74.6556112

4.590359后4.590359后10個(gè)特征值為:3.1515813.2796253.409659103.025552.901583.025552.901582 3 462.9633 73.34439 82.840371013六總結(jié)通過本學(xué)期對矩陣計(jì)算這門課程的學(xué)習(xí)我對數(shù)值代數(shù)方面的知識有了一個(gè)系統(tǒng)的回顧和總結(jié)，對求解線性方程組的直接法和迭代法有更深入的理解和探索，對求解矩陣特征值和特征向量的分二治之法，1@^20$算法，QR方法，以及奇異值分解的基本思想有深刻地認(rèn)識和研究，對矩陣的常用分解和變換如Schur分解，Householder變換，Givens變換掌握得比較扎實(shí)，對稀疏矩陣的存儲方法也有一定的掌握，對誤差分析和數(shù)值計(jì)算方法的穩(wěn)定性和收斂性分析有一定的理解。在編程過程中我對C語言的應(yīng)用更加熟練,特別是掌握了Matlab編程的基本技巧，雖然還不能應(yīng)用得如魚得水，這對我以后的學(xué)習(xí)和科研有很大的幫助。我本科畢業(yè)于湖南科技大學(xué)，在本科期間我在我班是名列前茅的學(xué)生，但剛?cè)氪筮B理工大學(xué)攻讀研究生學(xué)位時(shí)以前的那份“榮耀”已不復(fù)在，相反，很多同學(xué)都是重點(diǎn)大學(xué)的保送生，與之相比我卻成了“落后分子”，而我一向心比天高，剛開學(xué)時(shí)很難接受這段落差，一方面想抓緊時(shí)間，迎頭趕上。一方面又感覺很緊張很自卑，自暴自棄，有很長時(shí)間我都很消沉，然而在學(xué)習(xí)中特別是在本次實(shí)驗(yàn)報(bào)告編寫程序的過程中，我經(jīng)常請教同學(xué)和同學(xué)交流，我認(rèn)識到每個(gè)人都有自己長處和不足，其他同學(xué)也不是我想象中的那么強(qiáng)，這讓我有更大的信心投入學(xué)習(xí)，然而我自己的學(xué)習(xí)自覺性還是不夠一碰到難點(diǎn)就忍不住去請教別人，特別是最后兩個(gè)程序我請教了很多同學(xué)，這也是我以后學(xué)習(xí)和科研的大忌，與有的同學(xué)比我覺得我看書還是很細(xì)致，思維也比較嚴(yán)謹(jǐn)，我會在以后學(xué)習(xí)中取長補(bǔ)短，不斷進(jìn)步。本次實(shí)驗(yàn)報(bào)告的順利完成離不開程明松老師的耐心指導(dǎo)，程老師在上課時(shí)不僅對算法的基本原理，算法的計(jì)算量，不同算法的優(yōu)劣性，算法的穩(wěn)定性和收斂性進(jìn)行了細(xì)致的分析，還每個(gè)算法都在黑板上認(rèn)真板書，這便于我們更深刻的理解算法，也教育了我們做學(xué)問要有一絲不茍，精益求精的態(tài)度。在此表示衷心感謝，另外，我還要感謝我的同寢室好友鮑鏡如，沒有她的電腦和我遇到困難時(shí)給的鼓勵(lì)我的實(shí)驗(yàn)報(bào)告也不能如期完成。14七附錄部分程序代碼1求解一般右端項(xiàng)的10?口1設(shè)2方程組及Toeplitz矩陣的逆①#defineN5%主函數(shù)部分#include"stdio.h"floatx[N];floatY[N];floatT[N][N];floaty[N];main(){voidfunction1(floatR[N]);voidfunction2(floatR[N],floatb[N]);voidfunction3(floatR[N]);floatR1[N];floatb[N];inti,j,flag;R1[0]=N;b[0]=N;printf("pleaseinputtheflag1to3:\n");scanf("%d”,&flag);switch(flag){printf("pleaseinputthestringR1[N]:\n");for(i=1;i<N;i++)scanf("%f”,&R1[i]);printf(〃\n〃)；printf("輸出向量R1[N]\n");for(i=0;i<N;i++)printf("R1[%d]=%f”,i,R1[i]);printf(〃\n〃)；function1(R1);printf(〃輸出方程組解\n〃)；for(i=0;i<N;i++)printf("y1[%d]=%f”,i,Y[i]);printf(〃\n〃)；getch();break;printf("求矩陣的逆：\n");15printf("pleaseinputthestringR1[N]:\n");for(i=1;i<N;i++)scanf("%f”,&R1[i]);printf(〃\n〃)；printf("輸出向量R1[N]\n");for(i=0;i<N;i++)printf("R1[%d]=%f”,i,R1[i]);printf(〃\n〃)；function3(R1);getch();break;printf("pleaseinputthestringR1[N]:\n");for(i=1;i<N;i++)scanf("%f”,&R1[i]);printf(〃\n〃)；printf("輸出向量R1[N]\n");for(i=0;i<N;i++)printf("R1[%d]=%f”,i,R1[i]);printf(〃\n〃)；printf("pleaseinputthestringb[N]:\n");for(j=1;j<N;j++)scanf("%f”,&b[j]);printf(“輸出向量b[N]\n");for(i=0;i<N;i++)printf("b[%d]=%f”,i,b[i]);function2(R1,b);printf(〃輸出方程組222解\n〃)；for(i=0;i<N;i++)printf("x[%d]=%f〃，i,x[i]);getch();break;))②voidfunction1(floatR[N])%求解Yule-Walker5■程組{floatb=1,a,x;intk,i,j;Y[0]=R[0];Y[1]=-1*R[1];a=-1*R[1];for(k=1;k<N-1;k++){b=(1-a*a)*b;for(j=1;j<=k;j++)16{x=0.0;x=x+R[k+1-j]*Y[j];)a=-1*(R[k+1]+x)/b;for(i=1;i<=k;i++)Y[i]=Y[i]+a*Y[k+1-i];Y[k+1]=a;))③voidfunction2(floatR[N],floatb[N])%求解一般右端項(xiàng)的Toeplitz方程組{floatx1;floaty1;floata,d,u;inti,j,k;x[0]=b[0];y[0]=R[0];x[1]=b[1];y[1]=-1*R[1];a=-1*R[1];d=1;for(k=1;k<N-1;k++){d=(1-a*a)*d;for(i=1;i<=k;i++){x1=0;y1=0;x1=x1+R[i]*x[k+1-i];y1=y1+R[i]*y[k+1-i];)u=(b[k+1]-x1)/d;a=-1*(R[k+1]+y1)/d;for(j=1;j<=k;j++){x[i]=x[i]+u*x[k+1-i];y[i]=y[i]+u*y[k+1-i];)x[k+1]=u;y[k+1]=a;))④voidfunction3(floatR[N])%求解Toeplitz矩陣的逆 {17floata;floatx;inti,j,k;for(k=1;k<N;k++){x=0;x=x+R[k]*Y[k];)a=1/(1+x);T[0][0]=a;for(j=1;j<N;j++){T[0][j]=Y[j];T[j][0]=T[0][j];T[N-1][N-1-j]=T[j][0];T[N-1-j][N-1]=T[0][j];)for(i=1;i<(N-1)/2+1;i++)for(j=i;j<N-i+1;j++){T[i][j]=T[i-1][j-1]+(Y[i]*Y[j]-Y[N-i]*Y[N-j])/a;T[j][i]=T[i][j];T[N-j-1][N-i-1]=T[i][j];T[N-i-1][N-j-1]=T[j][i];)for(i=0;i<N;i++){for(j=0;j<N;j++)printf("T[%d][%d]=%f",i,j,T[i][j]);printf(〃\n〃)；)getch();)2實(shí)用共軛梯度法①實(shí)用的共軛梯度法functionx=CG_4_1(A,b,x0,e)x=x0;r=b-MV_4_1(A,x0);rr(1)=r'*r;k=1;whilerr(k)>rr(1)*eifk==1p=r;elseB=rr(k)/rr(k-1);18p=r+B*p;endW=MV_4_1(A,p);a=rr(k)/(p'*W);x=x+a*p;r=r-a*W;rr(k+1)=r'*r;k=k+1;endx;②稀疏矩陣與向量的乘積functionv=MV_4_1(B,b)m=length(B);v=zeros(length(b),1);fornum=1:mifB(num,1)==B(num,2)v(B(num,2))=v(B(num,2))+B(num,3)*b(B(num,1));elsev(B(num,1))=v(B(num,1))+B(num,3)*b(B(num,2));v(B(num,2))=v(B(num,2))+B(num,3)*b(B(num,1));endend3QR方法：①functionH=househ(x)構(gòu)造Householder變換矩陣Hn=length(x);I=diag(ones(1,n));sn=sign(x(1));ifsn==0sn=1;endz=x(2:n);ifnorm(z,inf)==0H=I;break;endsigma=-sn*norm(x);u=x;u(1)=u(1)-sigma;lo=sigma*(sigma-x(1));H=I-u*u'/lo;functionH=househ(x)n=length(x);I=diag(ones(1,n));sn=sign(x(1));ifsn==0sn=1;endz=x(2:n);ifnorm(z,inf)==0H=I;break;endsigma=-sn*norm(x);u=x;u(1)=u(1)-sigma;lo=sigma*(sigma-x(1));H=I-u*u'/lo;②functionA=hessen(A)用)Householder變換矩陣A化肥$$6廿6呼型[n,n]=size(A);19fork=1:n-2x=A(k+1:n,k);H=househ(x);A(k+1:n,k:n)=H*A(k+1:n,k:n);A(1:n,k+1:n)=A(1:n,k+1:n)*H;end③functionA=qrtran(m,A先用Givens變換作QR分解，再作相似變換Q=diag(ones(1,m));fori=1:m-1xi=A(i,i);xk=A(i+1,i);ifxk~=0d=sqrt(xi"2+xk"2);c=xi/d;s=xk/d;J=[cs;-sc];A(i:i+1,i:m)=J*A(i:i+1,i:m);Q(1:m,i:i+1)=Q(1:m,i:i+1)*J';endendA(1:m,1:m)=A(1:m,1:m)*Q;?function[iter,D]=qralg用基本QR算法求解矩陣的全部特征值ep=0.5*10"(-4);[n,n]=size(A);D=zeros(1,n);i=n;m=n;iter=0;A=hessen(A);while(n>0)ifm<=2la=eig(A(1:m,1:m));D(1:m)=la';break;enditer=iter+1;A=qrtran(m,A);fork=m-1:-1:1ifabs(A(k+1,k))<epifm-k<=2la=eig(A(k+1:m,k+1:m));j=i-m+k+1;D(j:i)=la';i=j-1;m=k;break;endendendend⑤雙重步位移QR迭代function[P,HH]=QR_2_2(H)20n=length(H);m=n-1;s=H(m,m)+H(n,n);t=H(m,m)*H(n,n)-H(m,n)*H(n,m);x=H(1,1廠2+H(1,2)*H(2,1)-s*H(1,1)+t;y=H(2,1)*(H(1,1)+H(2,2)-s);z=H(2,1)*H(3,2);P=eye(n);fork=0:n-3X=[x,y,z]';[v,beta]=house(X);p=eye(3)-beta*v*v';%Y=[eye(k),p,eye(n-k-3);]ifk<1q=1;elseq=k;endH(k+1:k+3,q:n)=p*H(