數值分析-第三章 線性方程組迭代法

1、1（五）（五） 迭代法的收斂條件迭代法的收斂條件（一）（一） 迭代法的一般形式迭代法的一般形式主要內容：主要內容：（二）（二） 雅可比（雅可比（JacobiJacobi）迭代法）迭代法（三）（三） 高斯高斯- -賽德爾賽德爾(Gsuss-Seidel)(Gsuss-Seidel)迭代法迭代法 （四）（四） 松弛法松弛法( (六六) ) 小結小結線性方程組的迭代解法 第三章2（一）一般迭代形式（一）一般迭代形式 x=Mx+g x=Mx+g （3-2）TnnnijbbbaA),(,1nnnnRxRgM)0(.,這里代入迭代公式x x（k+1k+1）=Mx=Mx(k)(k)+g+g (k=0,1,2

2、,-) (k=0,1,2,-) （3-3）1、對線性方程組 Ax=b （3-1）構造同解方程組3產生向量序列 ，當k充分大時，以 作為方程組（3-1）的近似解-一般迭代法。 kxkx主要解決的問題：(1).在多種方法中，構造那種迭代格式的好？(2).他們的收斂性如何？(3).在收斂的條件下，最佳計算運行格式？ 計算近似值的迭代次數.(4).多種迭代格式誤差的類比。42、向量序列、矩陣序列的收斂性、向量序列、矩陣序列的收斂性,如果定義定義3、1 設 為 中向量序列，( )kxnRnxR( )lim0kkxx(3-4)其中 為范數，則稱 收斂于x，記為 ( )kx( )limkkxx( )lim0

3、kkAA其中 為矩陣范數，則稱 收斂于A，記為 ( )kA()limkkAA(3-5)定義定義3、2 設 為n階方陣序列，A為n階方陣, 如果 () kA5定理定理3.1nR中的向量序列( )kx 收斂于 中的xnR當且僅當( )( )( )( )( )1212lim,(1,2, ),kiikTkkkknTnxxinwherexxxxxx xx充要條件定理定理3.2 設下面( )( ),kkijijAaAa( )kA 收斂于矩陣A的均為n階方陣，則矩陣序列( )lim,( ,1,2, )kijijkaai jn6）（因1322112222212111212111nnnnnnnnnnbxaxax

4、abxaxaxabxaxaxa系數矩陣A是非奇異的，不妨設將方程組（3-1）變?yōu)?, 2 , 1( , 0niaiibi)迭代法二、雅可比(Jaco7）（23112211223231212112121nnnnnnnnnnngxbxbxbxgxbxbxbxgxbxbx)., 2 , 1( ,), 2 , 1,( ,niabgnjijiaabiiiiiiijij其中gBxxggggbbbbbbBnnnnn,0002121221112記：)0(x取初始向量 代入（3-2）右邊有新向量。 ) 0 (x8 如此下去，產生一個序列 ，滿足為中向量序列 （3-33-3）上述過程給出的迭代法稱雅可比迭代法（簡

5、單迭代法）。 ), 2 , 1 , 0( ,)()1(kgBxxkk kxnR), 2 , 1 , 0(),()(1)(1)()1(kbDxADIgBxxkkk其矩陣表示為),(),(1111111nnnnaadiagDaadiagD其中（3-4）9i從1到n循環(huán)輸入：, , ,A b n置初值：0;0(1, )ikxin輸出k,xi(i=1,-,n)結束e1111iniiijjijjjj iiiyba xa xa 1;0kkei從1到n循環(huán)max(, ),iiiiexye xy輸入：, , ,A b n置初值：0;0(1, )ikxinyesNoJacobi迭代法的算法框圖10位有效數字。精

6、確到迭代法解方程組用3 ,12232121xxxxJacobi例例1解解計算結果如下：取迭代格式為, 0, 0, 2 , 1 , 0,2/ )1 (3/ )2()0(2)0(1)(1)1(2)(2)1(1xxkxxxxJacobikkkk 110.194450.250000.166670.5000000.583330.611110.500000.66667043210k)(1kx)(2kx567890.601850.597220.600310.599540.60050.208330.199080.201390.199850.20023200. 0,600. 0*2*1xx因而12例例2 2.

7、.用雅可比迭代法求解方程組4258321072210321321321xxxxxxxxx10111 102 ,115A 解解：100001000010D111110000100005D131BID A111110001012101001 102100511502 .02 .02 .001 .02 .01 .0014代入迭代式取初值Tx)0,0,0()0(TgBxx)4 . 8 , 3 . 8 , 2 . 7()0()1(TgBxx)50.11,70.10,17. 9 () 1 ()2(TgBxx)9992.12,9994.11,9994.10()1()9(，如此下去， 易知，方程組的精確解為

8、Tx)13,12,11(迭代結果（當迭代次數增大時）越來越接近精確解。15clear;fprintf(Jacobi迭代法)x1(1)=0;x2(1)=0;x3(1)=0;for i=1:10 x1(i+1)=7.2+0.1*x2(i)+0.2*x3(i); x2(i+1)=8.3+0.1*x1(i)+0.2*x3(i); x3(i+1)=8.4+0.2*x1(i)+0.2*x2(i);endx=x1,x2,x3x = 0 0 0 7.2000 8.3000 8.4000 9.7100 10.7000 11.5000 10.5700 11.5710 12.4820 10.8535 11.8534

9、 12.8282 10.9510 11.9510 12.9414 10.9834 11.9834 12.9804 10.9944 11.9944 12.9933 10.9981 11.9981 12.9978 10.9994 11.9994 12.9992 10.9998 11.9998 12.9997161231231238322041133631236xxxxxxxxx3,2,1Tx 例：用Jacobi迭代法求解方程組，精確解為17(1)( )( )123(1)( )( )213(1)( )( )3121(3220)81( 433)111( 6336)12kkkkkkkkkxxxxxxxx

10、x(0)0,0,0Tx(10)3.0000321.9998380.999813x解：按迭代過程取初始向量，迭代10次得實際計算結果表明Jacobi迭代法收斂與精確解。18迭代法 Seidel)-賽德爾(Gauss-三、高斯 研究研究JacobiJacobi迭代法就會發(fā)現在逐個求的分量迭代法就會發(fā)現在逐個求的分量時時, ,當計算時都已求出當計算時都已求出, ,但沒有被利用但沒有被利用. .直觀上看直觀上看, ,新算出的分量可能比舊的分量要準確新算出的分量可能比舊的分量要準確. . 因此因此, ,設想一旦當新分量已求出設想一旦當新分量已求出, ,馬上馬上就用它來就用它來代替代替,也就是在也就是在J

11、acobi迭代法中求迭代法中求(1)(1)(1)(1)()()12112,kkkkkkiixxxxxx時用代替()1,kixG aussSeidel這 就 是迭 代 法 。19GaussSeidel迭代格式如下：) 43 (), 2 , 1, 0(,)(1)(1)(1) 1(11) 1(22) 1(11) 1(2)(2)(323) 1(12122) 1(21)(1)(313)(21211) 1(1niabxaxaxaaxbxaxaxaaxbxaxaxaaxiinknnnknknnnknknnkkkknnkkk公式（3-4）表示為 20),2, 1 ,0( ,)()1()1(kgxUxLxkkk

12、000,00011122121nnnnnbbbUbbbL其中)(,11111LDDLDDDLIUDULDL若從而得迭代式), 2 , 1 , 0( ,)()(1)(1)1(kbLDUxLDxkk上式中矩陣 為Gauss-Seidel迭代矩陣。ULDM1)(21yesNo輸出k,xi(i=1,-,n)結束e1111m ax,iniiijjijjjjiiiiiybaxaxaexSexS1;0kkei從1到n循環(huán)輸入：, , ,A b n置初值：0;0(1, )ikxinGauss-Seidel迭代法的算法框圖22例例3.位有效數字。精確到迭代法解方程組用3 12232121xxxxSeidelGa

13、uss(1)( )12(1)(1)21(2)/3,0,1,2,(1)/2kkkkGauss Siedelxxkxx迭代格式為 解解(0)(0)120,0,xx取計算結果如下：230.199850.6003140.199750.199080.194450.1666700.600050.601850.611110.66667053210k)(1kx)(2kx200. 0,600. 0*2*1xx 因而較快。收斂速度迭代法迭代次數少，即法比迭代可得，與例比較例JacobiSeidelGauss 31 24例例4 4 用用Gauss-SeidelGauss-Seidel迭代法求解方程組迭代法求解方程組

14、4258321072210321321321xxxxxxxxx取初值 代入迭代式有 Tx) 0 , 0 , 0 () 0 (6440.11)42(510200. 9)832(1012000. 772101)722(101) 1 (2) 1 (1) 1 (3)0(3) 1 (1) 1 (2)0(3)0(2) 1 (1xxxxxxxxx解解：25(6)(5)xBxb易知，方程組的精確解為 Tx)13,12,11(迭代結果（當迭代次數為迭代結果（當迭代次數為6 6次時）就趕上雅可比迭次時）就趕上雅可比迭代代9 9次的結果。（其它情況以后再討論）次的結果。（其它情況以后再討論）如此下去，(10.999

15、9,11.9999,13.000)T26clear;fprintf(Gauss-Seidel迭代法)x1(1)=0;x2(1)=0;x3(1)=0;for i=1:7 x1(i+1)=7.2+0.1*x2(i)+0.2*x3(i); x2(i+1)=8.3+0.1*x1(i+1)+0.2*x3(i); x3(i+1)=8.4+0.2*x1(i+1)+0.2*x2(i+1);endx=x1,x2,x3Gauss-Seidel迭代法迭代法x=0007.20009.020011.644010.430811.671912.820510.931311.957212.977710.991311.99471

16、2.997210.998911.999312.999610.999911.999913.000011.000012.000013.000027clear;fprintf(SOR(超松弛)迭代法)x1(1)=0;x2(1)=0;x3(1)=0;for i=1:7 x1(i+1)=(1-1.42)*x1(i)+1.42/10*(7.2+0.1*x2(i)+0.2*x3(i); x2(i+1)=(1-1.42)*x1(i)+1.42/10*(8.3+0.1*x1(i+1)+0.2*x3(i); x3(i+1)=(1-1.42)*x1(i)+1.42/5*(8.4+0.2*x1(i+1)+0.2*x2

17、(i+1);endx=x1,x2,x3SOR(超松弛)迭代法x = 0 0 0 1.0224 1.1931 2.5114 0.6813 0.8302 2.0420 0.8061 0.9619 2.1999 0.7600 0.9133 2.1421 0.7770 0.9313 2.1634 0.7707 0.9246 2.1556 0.7730 0.9271 2.1585281231231238322041133631236xxxxxxxxx3,2,1Tx 例：用Gauss-Seidel迭代法求解方程組，精確解為。29(1)( )( )123(1)(1)( )213(1)(1)(1)3121(3

18、220)81( 433)111( 6336)12kkkkkkkkkxxxxxxxxx(0)0,0,0Tx(5)2.9998432.0000721.000061x解：其迭代過程取初始向量，迭代5次得與Jacobi迭代10次的精度相同。30注：1.Jacobi迭代法收斂，迭代法不一定收斂；2.編程計算時，Gauss-Seidel迭代法用一 套存儲單元，而Jacobi迭代法用兩套存儲 單元。31四四. .松弛法松弛法為加速迭代過程的收斂，引入參數 在上得到一種新算法,)()1(xxxkkikikixxx)()1(), 2 , 1(),()1 ()(1)1(11)(nixaxabaxkjnijijkj

19、ijijiiiki（3-5）按（按（3-5）計算方程組（）計算方程組（3-1）的近似解序列的方法為松弛法。）的近似解序列的方法為松弛法。 為低松弛，為低松弛， 為為Gauss-SeidelGauss-Seidel迭代法，迭代法， 為超松弛法為超松弛法(SOR)(SOR)。, 1, 1, 132例例5. 用超松弛法求解方程組8 . 1202123232121xxxxxxx(0 )1.4,(1,1,1)Tx其 中有迭代公式（3-5）有(1)( )( )113(1)( )(1)( )2213(1)( )(1)3320.40.7(1)0.40.7(),0.40.7(1.8)(0,1,2,)kkkkkk

20、kkkkxxxxxxxxxxk ( 0 )(1,1,1)Tx將代 入 上 式解：解：3356.1)18 .1 (7 .04 .01)11 (7 .04 .01)11 (7 .04 .0)1(3)1(2)1(1xxx繼續(xù)下去，方程組有解 6001.1,3996.1,2000.1)9(3)9(2)9(1xxx與精確解 6 . 1, 4 . 1, 2 . 1321xxx比較，誤差為 3)9(10210004.0 xx34clear;fprintf(SOR(超松弛)迭代法)x1(1)=1;x2(1)=1;x3(1)=1;for i=1:9 x1(i+1)=(1-1.4)*x1(i)+1.4/2*(1+

21、x2(i); x2(i+1)=(1-1.4)*x2(i)+1.4/2*(x1(i+1)+x3(i); x3(i+1)=(1-1.4)*x3(i)+1.4/2*(1.8+x2(i+1);endx=x1,x2,x3SOR(超松弛)迭代法x =1.0000 1.0000 1.0000 1.0000 1.0000 1.5600 1.0000 1.3920 1.6104 1.2744 1.4626 1.6396 1.2140 1.4125 1.5929 1.2032 1.3922 1.5974 1.1933 1.3966 1.5987 1.2003 1.4006 1.6010 1.2003 1.4007

22、 1.600111.11.21.31.411.21.41.61.811.21.41.61.8x1SOR-V919x2x3351.41231231234202422233xxxxxxxxx (1)( )71,2,3max10kkiiixx(0)1,1,1Tx例：分別用Gauss-Seidel迭代法和SOR迭代法()求解方程組當取初始向量時，停止迭代36結論：v1。Gauss-Seidel迭代法要迭代72次得v2SOR迭代法( )，只須迭代72次得(72)0.9999520.9999451.999995x1.41.4(25)0.9999940.9999982.000000 x適當選擇，SOR迭代法

23、具有明顯的加速收斂效果。373GaussSeidel、迭代法2Jacobi、迭代法小結 4.松弛迭代法松弛迭代法(SOR)1、迭代法的一般形式38112233131132axbaxbaxb11 11221112221 12222(0)a xa xba aa xa xb122111221a aa a思考與練習思考與練習 1若用Jacobi迭代法求解方程組迭代收斂的充要條件是討論實數a與收斂性的關系。2若用Jacobi迭代法求解方程組39五五. .迭代法的收斂條件迭代法的收斂條件 5 51 1 矩陣的譜半徑矩陣的譜半徑 稱的特征值為階矩陣設), 2 , 1(nininnRA 定義定義6（譜半徑（譜

24、半徑)ini1max）（A的為矩陣A譜半徑譜半徑。之間的關系。下面討論譜半徑與范數值，則由，的特征的為對應于的任意特征值，為設AXAAXX 再由相容性條件有可得,AXX 40XAAXX故有為非零向量，因為, 0XXA的任何特征值，于是有是由于AAA )(的任何一種范數。的譜半徑不超過即矩陣AA矩陣矩陣A A的譜半徑與范數有如下關系：的譜半徑與范數有如下關系： )()(AAA410limkkA定理定理1 1 ： 設A為n階方陣，則的充要條件為 1)(A0limkkA0limkkA由定義kkkAAA)()(0而由極限存在準則，有 0)(limkkA1)(A所以5 5、2 2迭代法的收斂條件迭代法的

25、收斂條件證明：證明：若若必要性必要性,42若若 1)(A12)(1A取矩陣A的譜半徑與范數的關系有12)(1)(AAA0limlimkkkkAA于是kkAA而0limkkA所以充分性充分性。43 定理定理2 2 對任意初始向量 和右邊g，由迭代 格式 產生向量序列 收斂條件是)0(x), 2 , 1 , 0(,)()1(kgMxxkk1)(M)(kx設存在n維向量x*，使得*)(limxxkk則x*滿足 gMxx*由迭代公式 (0)*()kMxx證明：證明：()*kxx(1)*kM xgM xg(1)*()kM xx2(2)*()kMxx440)(lim)(lim*)(*)0(xxxxMkkk

26、k于是因x0為任意n維向量，上式成立必須0limkkM。由定理知1)(M)(kx推論推論1 1、在定理、在定理2 2條件下，若條件下，若 ，則，則 收斂。收斂。推論推論2 2、松弛法收斂的必要條件為、松弛法收斂的必要條件為 。1M2045給定線性方程組例例5。1、21132121xx是否收斂？迭代法求解迭代法和問用JacobiSeidelGauss 矩陣為所給線性方程組的系數解：1321迭代法的迭代矩陣為）Jacobi10320)(1ULDJ46的特征方程為因而其迭代矩陣J0326621，解得 16)(J 于是迭代法發(fā)散。因此Jacobi迭代法的迭代矩陣為SeidelGauss )26020)

27、(1ULDG47的特征方程為因而其迭代矩陣G06026021，解得 16)(A 于是迭代法發(fā)散。因此SeidelGauss48 解方程組3222122321321321xxxxxxxxx討論雅可比迭代法與Gauss-Seidel迭代法的收斂性。 （1）由定理2知迭代法是收斂等價迭代矩陣 的譜半徑1)(M,100010001,122111221DA例例5。2、解解：4900 01 0 0 ,22 0L 1JacobiBID A特征方程為.02211223BI1)(JB因此雅可比迭代法收斂。 022001 ,000U022101220 50（2）由Gauss-Seidel迭代法， 111221DL

28、1()GAUSSMDLU11()11021DL 100022110010210 02223251Gauss-Seidel迭代法發(fā)散。12)(GAUSSM1230,2,22023002IM特征方程為2(2)0. 52注注：例也說明確實 只是松弛法收斂的必要條件，而非充分的，因為Gauss-Seidel迭代法（ ）。 20153（1-1）若矩陣為嚴格對角占優(yōu),即對所有的方陣滿足為方便給出以下判據：為方便給出以下判據： （1-2）若矩陣為弱對角占優(yōu)，若至少由一個值，使下式成立。 1,(1,2, )niiijjj iaain1,(1,2, )niiijjj iaain（1）：541112220AAAA

29、（）： （2-1）若矩陣為可約,若矩陣能通過行、 列的互換成為 1122,AA其中為方陣.（-2）若矩陣為不可約,若矩陣不能通過行、列的互換成為 1112220AAAA55矩陣為對稱正定的，則松弛法收斂矩陣為對稱正定的，則松弛法收斂20如：例中系數矩陣如：例中系數矩陣 51121012110A嚴格對角占優(yōu)，嚴格對角占優(yōu)，Jacobi迭代法迭代法和和Gauss-Seidel迭代法迭代法均收斂。均收斂。設有線性方程組，下例結論成立：設有線性方程組，下例結論成立：（）若矩陣為（）若矩陣為嚴格對角占優(yōu)嚴格對角占優(yōu)或為或為不可約弱對角占優(yōu)不可約弱對角占優(yōu)，則則Jacobi迭代法迭代法和和Gauss-Se

30、idel迭代法迭代法均收斂。均收斂。20（）若矩陣為對稱正定的，（）若矩陣為對稱正定的，則松弛法收斂則松弛法收斂（）若矩陣為嚴格對角占優(yōu)，（）若矩陣為嚴格對角占優(yōu)，則松弛法收斂。，則松弛法收斂。1056為非嚴格對角占優(yōu)陣，但為對稱正定矩陣，210121012A4 .1例中例中系數矩陣松弛法收斂。試討論用三種迭代法求解的收斂性。例、設有方程組，其中15 .05 .05 .015 .05 .05 .01A57因為對稱正定矩陣（各階主子式大于零），由判別條件知Gauss-Seidel迭代法收斂， 松弛法收斂。但不是弱對角占優(yōu)陣，則不能用判別條件斷定迭代法收斂性。只能通過計算迭代矩陣為2005 . 0

31、5 . 05 . 005 . 05 . 05 . 001ADIB解解：581,213211)(JB由定理知雅可比迭代法不收斂。. 0) 1()2()4/3(42222222113111111BI特征方程為59、誤差估計、誤差估計 收斂于收斂于，則有誤差估計，則有誤差估計gMxxkk)()1()63(1)0()1(*)(xxMMxxkk定理定理、設有迭代格式、設有迭代格式)(, 1kxM 若注注由誤差估計式（由誤差估計式（3 36 6），據給定的），據給定的精度精度，可估計出迭代次數可估計出迭代次數： ）（73ln)1(ln)0()1(MxxMEk60事實上事實上，1IM*(),IM xg*1(

32、)(*1)xIMg( )*(1)*()()kkxxMxgMxg( )*1(0)(1)()kkxxMIMxx( )*(0)(1)1(*3)1kkxxMxxM, 1xMxg(0)(1)()(*2)fxx61定理、公式的應用定理、公式的應用( )*(0)(1)11kkxxMxxM(1)(0)(1)lnlnMxxkM(1)( )kkxMxg62如在例中，若要求各分量的誤差絕對值不超過，則由 4 .83 .82 .7.000,02 .02 .02 .001 .02 .01 .00)1()0(xxM4.8,4.0)0()1(xxM有代入（3-7）得932.124 . 0ln4 . 8)4 . 01 (10

33、ln4k所需迭代1次才能保證各分量的誤差絕對值不超過。 63采用采用事后誤差估計事后誤差估計方法方法-用用相鄰兩次迭代值之差相鄰兩次迭代值之差作為達到精度標準。作為達到精度標準。 收斂于，則有誤差估計gMxxkk)()1()(, 1kxM若。)1()(*)(1kkkxxMMxx更好的結論更好的結論：設有迭代格式64等價的極值問題等價的極值問題-求解方程組兩種方法求解方程組兩種方法1、最速下降法；、最速下降法；2、共軛梯度法、共軛梯度法設設x*方程組方程組Ax=b的精確解，的精確解，Ax*=b，若，若A是正定矩陣是正定矩陣僅當僅當x=x*，二次函數，二次函數F0(x)=(x-x*)A(x-x*)

34、=xAx-2bx+(x*)Ax*達到極小值達到極小值F0(x*)=0。這里這里F0(x)與二次函數與二次函數F(x)=xAx-2bx僅差一個常數僅差一個常數(x*)Ax*，它們的極小值點是相同的，它們的極小值點是相同的，所以解所以解Ax*=b等價于求解二次函數等價于求解二次函數F(x)的極小值點的極小值點x*。65例、用共軛梯度法求解對稱正定方程組例、用共軛梯度法求解對稱正定方程組Txxxxx)0 , 0(,023)0(2121A=3 -1 ;-1,1 ;b=2 0;x,k=getd(A,b)d = 0.2222 0.6667d = 1.0e-015 * 0.0000 -0.1110d = 1

35、.0e-015 * -0.1110 -0.1110 x = 1.0000 1.0000k = 366 %A=2-1-1;-120;-101;b=010;x,k=getd(A,b)%d=0.50000.25000d=0.22220.11110.3333d=1.0e-016*0.00000.0000-0.5551d=1.0e-016*-0.55510.0000-0.5551x=111k=4例、用共軛梯度法求解對稱正定方程組例、用共軛梯度法求解對稱正定方程組1231231232020100 xxxxxxxxx67function x,k=getd(A,b,x0,ep,Nmax)%用共軛梯度法求解正定

36、系數矩陣線性方程組AX=b%A為線性方程組的系數矩陣，正定對稱，b為方程組的右端向量%x為解向量，k為迭代次數，x0為迭代初值%ep為精度，Nmax為迭代次數上限以防發(fā)散（默認值為500）n=length(A);k=0;if nargin5 Nmax=500;endif nargin4 ep=1e-10;endif narginep&kNmax k=k+1;x0=x; alpha=(r*r)/(d*A*d);r1=r; s=alpha*d;x=x+s;r=r-A*s; beta=(r*r)/(r1*r1);d=r+beta*dendif k=Nmax warning(已迭代上限次數);end 68三、小結斂（2）迭代法收的判定.陣斂譜徑嚴對（1）向量-矩序列的極限、迭代法收、 半、格角占優(yōu);69112233131132axbaxbaxb11 112211