迭代解法的matlab實(shí)現(xiàn)

1、逐次求出的A為低階稠解線性方程組AX =b逐次求出的A為低階稠適合于算中常會(huì)遇至吠型稀近似解逼近精適合于算中常會(huì)遇至吠型稀密矩陣（指n不大且元疏矩陣（指n很大且零兀較多）的方程組，迭代法在計(jì)算和存貯兩方面都適合后一種情 況.因?yàn)榈ㄊ峭ㄟ^逐次迭代來逼近方程組的解，所以收斂性和收斂速度是構(gòu)造迭代法 時(shí)應(yīng)該注意的問題.另外，因?yàn)椴煌南禂?shù)矩陣具有不同的性態(tài)，所以大多數(shù)迭代方法都 具有一定的適用范圍.有時(shí)，某種方法對(duì)于一類方程組迭代收斂，而對(duì)另一類方程組迭代 時(shí)就發(fā)散.因此，我們應(yīng)該學(xué)會(huì)針對(duì)具有不同性質(zhì)的線性方程組構(gòu)造不同的迭代4.1迭代法和斂散性及其MATLAB程序4.1.2迭代法斂散性的判別及

2、其MATLAB程序根據(jù)定理4.1和譜半徑定義，現(xiàn)提供一個(gè)名為pddpb.m的M文件，用于判別迭代公 式（4.7）產(chǎn)生的迭代序列的斂散性.用譜半徑判別迭代法產(chǎn)生的迭代序列的斂散性的MATLAB主程序輸入的量：線性方程組AX = b的迭代公式（4.7）中的B ；輸出的量：矩陣B的所有特征值和譜半徑mH = P （B）及其有關(guān)迭代法產(chǎn)生的迭代序 列的斂散性的相關(guān)信息.function H=ddpbj（B）H=eig（B）;mH=norm（H,inf）;if mH=1disp（請(qǐng)注意：因?yàn)樽V半徑不小于1,所以迭代序列發(fā)散,譜半徑mH和B的所 有的特征值H如下：）elsedisp（請(qǐng)注意：因?yàn)樽V半徑小于

3、1，所以迭代序列收斂,譜半徑mH和B的所有 的特征值H如下：）endmH4.1.3與迭代法有關(guān)的MATLAB命令（一）提?。óa(chǎn)生）對(duì)角矩陣和特征值可以用表4 - 1的MATLAB命令提取對(duì)角矩陣和特征值.表4-1提?。óa(chǎn)生）對(duì)角矩陣和特征值MATLAB命令功能DX=diag(X)若輸入向量X，則輸出DX是以X為對(duì)角元的對(duì)角矩陣； 若輸入矩陣X，則輸出DX是以X的對(duì)角元構(gòu)成的向量；DX=diag(diag(X)輸入矩陣X，輸出DX是以X的對(duì)角元構(gòu)成的對(duì)角矩陣， 可用于迭代法中從A中提取D.lm=eig(A)輸入矩陣A，輸出lm是A的所有特征值.（二）提?。óa(chǎn)生）上（下）三角形矩陣可以用表4-2的

4、MATLAB命令提取矩陣的上三角形矩陣和下三角形矩陣.表4-2 提取矩陣的上三角形矩陣和下三角形矩陣MATLAB命令功能U=triu(A)輸入矩陣A，輸出U是A的上三角形矩陣.L=tril(A)輸入矩陣A，輸出L是A的下三角形矩陣.U=triu(A,1)輸入矩陣A，輸出U是A的上三角形矩陣，但對(duì)角元為0，可 用于迭代法中從A中提取U .L=tril(A,-1)輸入矩陣A，輸出L是A的下三角形矩陣，但對(duì)角元為0，可 用于迭代法中從A中提取L.（三）稀疏矩陣的處理對(duì)稀疏矩陣在存貯和運(yùn)算上的特殊處理，是MATLAB進(jìn)行大規(guī)?？茖W(xué)計(jì)算時(shí)的特點(diǎn)和優(yōu)勢之一.可以用表4-3的MATLAB命令，輸入稀疏矩陣的

5、非零元（零元不必輸入），即可進(jìn)行運(yùn)算.表4-3 稀疏矩陣的MATLAB命令MATLAB命令功能ZA=sparse(r,c,v,m,n)表示在第r仃、第c列輸入數(shù)值力 矩陣共m仃n 列，輸出ZA，給出Sc）及V為一稀疏矩陣.MA=full(ZA)輸入稀疏矩陣Z4,輸出為滿矩陣M4 （包含零元）4.2雅可比（Jacobi）迭代及其MATLAB程序4.2.2雅可比迭代的收斂性及其MATLAB程序根據(jù)定理4.3和公式（4.14），現(xiàn)提供一個(gè)名為jspb.m的M文件如下:判別雅可比迭代收斂性的MATLAB主程序輸入的量：線性方程組AX=b的系數(shù)矩陣A；輸出的量：系數(shù)矩陣A =(.)萬nxnj輸出的量：系

6、數(shù)矩陣A =(.)萬nxnj=ii豐kkjkk(k = 1, 2,n）的值和有關(guān)雅可比迭代收斂性的相關(guān)信息.function a=jspb(A)n m=size(A);or j=1:ma(j)=sum(abs(A(:,j)-2*(abs(A(j,j);end for i=1:nif a（i）=0disp（請(qǐng)注意：系數(shù)矩陣A不是嚴(yán)格對(duì)角占優(yōu)的，此雅可比迭代不一定收斂）returnendendif a（i） A=10 -1 -2;-1 10 -2;-1 -1 5;a=jspb（A）運(yùn)行后輸出結(jié)果請(qǐng)注意：系數(shù)矩陣A是嚴(yán)格對(duì)角占優(yōu)的，此方程組有唯一解，且雅可比迭代收 斂a =-8-8-1（2）在MAT

7、LAB工作窗口輸入程序 A=10 -1 -2;-1 10 -2;-1 -1 0.5;a=jspb（A）運(yùn)行后輸出結(jié)果請(qǐng)注意：系數(shù)矩陣A不是嚴(yán)格對(duì)角占優(yōu)的，此雅可比迭代不一定收斂a =-8.0000e+000 -8.0000e+000 3.5000e+0004.2.3雅可比迭代的兩種MATLAB程序利用MATLAB程序和雅可比迭代解線性方程組AX = b的常用的方法有兩種，一種方法是根據(jù)雅可比迭代公式(4.11)、(4.12)式、定理4.3和公式(4.14)編寫一個(gè)名為jacdd.m的M文件并保存，然后在MATLAB工作窗口輸入對(duì)應(yīng)的命令，執(zhí)行此文件；另一 種方法是根據(jù)雅可比迭代的定義，利用提取

8、對(duì)角矩陣和上、下三角矩陣的MATLAB程序解 線性方程組AX = b .下面我們分別介紹這兩種方法.(一)雅可比迭代公式的MATLAB程序用雅可比迭代解線性方程組AX = b的MATLAB主程序輸入的量：線性方程組AX = b的系數(shù)矩陣A和b,初始向量X0范數(shù)的名稱P= 1,才aakk(k = 1, 2, , n)的值和有關(guān)雅萬nxn2, inf或fro.，近似解X的誤差(精度)wucha和迭代的最大次數(shù)maxi； (k = 1, 2, , n)的值和有關(guān)雅萬nxnj=ii豐k可比迭代收斂性的相關(guān)信息及其AX = b的精確解jX和近似解X.根據(jù)雅可比迭代(4.11)、(4.12)式、定理4.3

9、和公式(4.14)，現(xiàn)提供一個(gè)名為jacdd.m 的M文件如下：function X=jacdd(A,b,X0,P,wucha,max1) n m=size(A);or j=1:ma(j)=sum(abs(A(:,j)-2*(abs(A(j,j);end for i=1:nif a(i)=0disp(請(qǐng)注意：系數(shù)矩陣A不是嚴(yán)格對(duì)角占優(yōu)的，此雅可比迭代不一定收斂)returnendendif a(i)0disp(請(qǐng)注意：系數(shù)矩陣A是嚴(yán)格對(duì)角占優(yōu)的，此方程組有唯一解，且雅 可比迭代收斂)endfor k=1:max1kfor j=1:mX(j)=(b(j)-A(j,1:j-1,j+1:m)*X0(

10、1: j-1,j+1:m)/A(j,j);endX,djwcX=norm(X-X0,P); xdwcX=djwcX/(norm(X,P)+eps); X0=X;X1=Ab;if (djwcXwucha)&(xdwcXwucha)&(xdwcXwucha)disp(請(qǐng)注意：雅可比迭代次數(shù)已經(jīng)超過最大迭代次數(shù)max1 )enda,X=X;jX=X1,例4.2.3用8范數(shù)和判別雅可比迭代的MATLAB主程序解例4.2.2中的方程組，解 的精度為0.001，分別取最大迭代次數(shù)max1=100, 5,初始向量X0=(0 0 0)T，并比較 它們的收斂速度.解 (1)取最大迭代次數(shù)max1=100時(shí).首先

11、保存名為jacdd.m的M文件，然后在MATLAB工作窗口輸入程序 A=10 -1 -2;-1 10 -2;-1 -1 5; b=7.2;8.3;4.2;X0=0 0 0; X=jacdd(A,b,X0,inf,0.001,100)運(yùn)行后輸出結(jié)果請(qǐng)注意：系數(shù)矩陣A是嚴(yán)格對(duì)角占優(yōu)的，此方程組有唯一解，且雅可比迭代收斂請(qǐng)注意：雅可比迭代收斂，此方程組的精確解jX和近似解X如下： a = -8-8-1jX =1.10001.20001.3000X =1.09941.19941.2993在MATLAB工作窗口輸入程序 A=10 -1 -2;-1 10 -2;-1 -1 0.5； b=7.2;8.3；4

12、.2; X0=0 0 0;X=jacdd（A,b,X0,inf, 0.001,100）運(yùn)行后輸出結(jié)果請(qǐng)注意：系數(shù)矩陣A不是嚴(yán)格對(duì)角占優(yōu)的，此雅可比迭代不一定收斂請(qǐng)注意：雅可比迭代收斂，此方程組的精確解jX和近似解X如下： a =-8.0000-8.00003.5000jX =24.500024.6000 106.6000X = 24.073824.1738 104.7974（2）取最大迭代次數(shù)max1=5時(shí)，在MATLAB工作窗口輸入程序 A=10 -1 -2;-1 10 -2;-1 -1 5;b=7.2;8.3;4.2; X0=0 0 0; X=jacdd（A,b,X0,inf,0.001,

13、5） 運(yùn)行后輸出結(jié)果請(qǐng)注意：系數(shù)矩陣A是嚴(yán)格對(duì)角占優(yōu)的，此方程組有唯一解，雅可比迭代收斂 請(qǐng)注意：雅可比迭代次數(shù)已經(jīng)超過最大迭代次數(shù)max1 a =-8-8-1jX = 1.10001.20001.3000X =1.09511.19511.2941在MATLAB工作窗口輸入程序 A=10 -1 -2;-1 10 -2;-1 -1 0.5; b=7.2;8.3;4.2;X0=0 0 0; X=jacdd（A,b,X0,inf, 0.001,5） 運(yùn)行后輸出結(jié)果請(qǐng)注意：系數(shù)矩陣A不是嚴(yán)格對(duì)角占優(yōu)的，此雅可比迭代不一定收斂 請(qǐng)注意：雅可比迭代次數(shù)已經(jīng)超過最大迭代次數(shù)max1a =-8.0000-8.

14、00003.5000jX =24.500024.6000106.6000X =5.54905.649027.6553a = fj=1 i壬a = fj=1 i壬k-|akJ （k = 1,2, , n）的值越小，雅可比迭代收斂的速度越快.（二）利用雅可比迭代定義編寫的解線性方程組的MATLAB程序利用雅可比迭代定義編寫解線性方程組（4.5）的MATLAB程序的一般步驟步驟1將線性方程組（4.5）的系數(shù)矩陣A分解為A = D - L-U，其中(0)(0 a a )121na 00:21,U =: I. I. an-1,na a -a 00 A=a11 a12 -a1n; a21 a22 a2n;

15、an1 D=diag(A) U=triu(A,1), L=tril(A,-1)運(yùn)行后即可輸出A的D,L, U ；步驟2若對(duì)角矩陣D非奇異(即a,。0, i = 1,n)，則(4.5)化為X = D -1(L + U) X + D-ib .若記B = D -1 (L + U), f = D -1b .則方程組的迭代形式可寫作11X (k+1) = B X (k) + f(k = 0,1,2 )可以利用下面的MATLAB程序完成1dD=det(D);ifdD=0disp(請(qǐng)注意：因?yàn)閷?duì)角矩陣D奇異，所以此方程組無解.)elsedisp(請(qǐng)注意：因?yàn)閷?duì)角矩陣D非奇異，所以此方程組有解.)iD=inv

16、(D); B1=iD*(L+U);f1=iD*b;for k=1:max1X= B1*X0+ f1; X0=X; djwcX=norm(X-X0,P); xdwcX=djwcX/(norm(X,P)+eps); X1=Ab;if (djwcXwucha)&(xdwcXwucha)|(xdwcXwucha)disp(請(qǐng)注意：雅可比迭代次數(shù)已經(jīng)超過最大迭代次數(shù)max1 )endenda,X=X;jX=X1,4.3高斯-塞德爾(Gauss-Seidel)迭代及其MATLAB程序4.3.3高斯-塞德爾迭代兩種MATLAB程序利用MATLAB程序和高斯-塞德爾迭代解線性方程組AX = b的常用方法有兩種

17、，一 種方法是根據(jù)高斯-塞德爾迭代公式(4.16)、(4.17)、定理4.3和公式(4.14)編寫一 個(gè)名為gsdd.m的M文件并保存，然后在MATLAB工作窗口輸入對(duì)應(yīng)的命令，執(zhí)行此文件; 另一種方法是根據(jù)高斯-塞德爾迭代的定義，利用提取對(duì)角矩陣和上、下三角矩陣的 MATLAB程序解線性方程組AX = b.下面我們分別介紹這兩種方法.(一)高斯-塞德爾迭代定義的MATLAB程序1根據(jù)高斯-塞德爾迭代定義，現(xiàn)提供名為gsdddy.m的M文件如下：用高斯-塞德爾迭代定義解線性方程組AX = b的MATLAB主程序1輸入的量:線性方程組AX = b的系數(shù)矩陣A和b,初始向量X0,范數(shù)的名稱P =

18、1, 2,inf,或fro,近似解X的誤差(精呷 wucha和迭代的最大次數(shù)maxi.輸出的量：以系數(shù)矩陣A =L. 的對(duì)角元構(gòu)成的對(duì)角矩陣D. A的上三角形矩陣 nxnU,但對(duì)角元為0、A的下三角形矩陣L,但對(duì)角元為0和有關(guān)高斯-塞德爾迭代收斂性的 相關(guān)信息及其AX = b的精確解jX和近似解X.function X=gsdddy(A,b,X0,P,wucha,max1)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1); dD=det(D);if dD=0disp(請(qǐng)注意：因?yàn)閷?duì)角矩陣D奇異，所以此方程組無解.)elsedisp(請(qǐng)注意：因?yàn)閷?duì)角矩陣D非奇異

19、，所以此方程組有解.) iD=inv(D-L); B2=iD*U;f2=iD*b;jX=Ab; X=X0; n m=size(A);for k=1:max1X1= B2*X+f2; djwcX=norm(X1-X,P);xdwcX=djwcX/(norm(X,P)+eps);if (djwcXwucha)|(xdwcXwucha)returnelsek,X1,k=k+1;X=X1;endendif (djwcXwucha)|(xdwcX A=10 -1 -2;-1 10 -2;-1 -1 0.5; b=7.2;8.3;4.2; X0=0 0 0;X=gsdddy(A,b,X0,inf, 0.0

20、01,100)運(yùn)行后輸出結(jié)果請(qǐng)注意：因?yàn)閷?duì)角矩陣D非奇異，所以此方程組有解.請(qǐng)注意：高斯-塞德爾迭代收斂，此A的分解矩陣D,U,L和方程組的精確解jX和 近似解X如下：D =L =10.000000000010.00000100000.5000110U =jX =01224.500024.6000 106.6000002X =00024.499624.5996 106.5984此近似解與例4.2.3中的（1）中的解X = （24.0738, 24.1738, 104.7974） t比 較，在相同的條件下，高斯-塞德爾迭代比雅可比迭代得到的近似解的精度更高.（2）在MATLAB工作窗口輸入程序

21、A=3 4 -5 7;2 -8 3 -2;4 51 -13 16；7 -2 21 3；b=5；2;-1;21;X0=0 0 0 0;X=gsdddy（A,b,X0,inf,0.001,100）運(yùn)行后輸出結(jié)果請(qǐng)注意：因?yàn)閷?duì)角矩陣D非奇異，所以此方程組有解.請(qǐng)注意：高斯-塞德爾迭代發(fā)散，并且迭代次數(shù)已經(jīng)超過最大迭代次數(shù)max1, 方程組的精確解jX和迭代向量X如下：jX =0.1821-0.25710.72861.3036X = 1.0e+142 *0.28830.10620.3622-3.1374（二）高斯-塞德爾迭代公式的MATLAB程序2根據(jù)高斯-塞德爾迭代公式（4.16）、（4.17）、定

22、理4.3和公式（4.14），現(xiàn)提供名為gsdd.m的M文件如下用高斯-塞德爾迭代解線性方程組AX = b的MATLAB主程序2輸入的量:線性方程組AX = b的系數(shù)矩陣A和b,初始向量X0,范數(shù)的名稱P = 1,2, inf,或fro.，近似解X的誤差（精度）wucha和迭代的最大次數(shù)maxi.才aakjkkj=ii豐k（k = 1, 2, , n）j=ii豐k（k = 1, 2, , n）的值和有關(guān)高萬nxn斯-塞德爾迭代收斂性的相關(guān)信息及其AX = b的精確解/X和近似解X.function X=gsdd(A,b,X0,P,wucha,max1)n m=size(A);for j=1:ma

23、(j)=sum(abs(A(:,j)-2*(abs(A(j,j);endfor i=1:nif a(i)=0disp(請(qǐng)注意：系數(shù)矩陣A不是嚴(yán)格對(duì)角占優(yōu)的，此高斯塞德爾迭代不定收斂)return endend if a(i)0disp(請(qǐng)注意：系數(shù)矩陣A是嚴(yán)格對(duì)角占優(yōu)的，此方程組有唯一解，且高 斯-塞德爾迭代收斂)endfor k=1:max1for j=1:mif j=1X(1)=(b (1)-A(1,2:m)*X0(2:m)/A(1,1)endif j=mX(m)=(b(m)-A(m,1:M1)*X(1:M1)/A(m,m);endfor j=2:M1X(j)=(b(j)-A(j,1:j-

24、1)*X(1:j-1) -A(j,j+1:m)*X(j+1:m)/A(j,j);endenddjwcX=norm(X-X0,P);xdwcX=djwcX/(norm(X,P)+eps); X0=X;X1=Ab;if (djwcXwucha)|(xdwcXwucha)&(xdwcXwucha)disp(請(qǐng)注意：高斯-塞德爾迭代次數(shù)已經(jīng)超過最大迭代次數(shù)maxi )enda,X=X;jX=X1,4.4解方程組的超松弛迭代法及其MATLAB程序用雅可比迭代法和高斯-塞德爾迭代法解線性方程組時(shí)，有時(shí)收斂速度很慢，為了提 高收斂速度，我們介紹超松弛迭代法，它是對(duì)高斯-塞德爾迭代的一種加速方法，適用于 大型

25、稀疏矩陣方程組的求解.4.4.2超松弛迭代法收斂性及其MATLAB程序根據(jù)定理4.5和譜半徑定義，現(xiàn)提供名為ddpbj.m的M文件，用于判別超松弛迭代 公式(4.23)產(chǎn)生的迭代序列的斂散性.用譜半徑判別超松弛迭代法產(chǎn)生的迭代序列的斂散性的MATLAB主程序輸入的量：線性方程組AX = A的系數(shù)矩陣A和松弛因子om；輸出的量：矩陣% = (-誠)-1U + (1 -)D 的所有特征值和譜半徑 mH = P (B“)及其有關(guān)超松弛迭代法產(chǎn)生的迭代序列的斂散性的相關(guān)信息function H=ddpbj(A,om)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);

26、iD=inv(D-om*L);B2=iD*(om*U+(1-om)*D);H=eig(B2);mH=norm(H,inf);if mH=idisp(請(qǐng)注意：因?yàn)樽V半徑不小于1,所以超松弛迭代序列發(fā)散,譜半徑mH 和B的所有的特征值H如下：)elsedisp(請(qǐng)注意：因?yàn)樽V半徑小于1,所以超松弛迭代序列收斂,譜半徑mH 和B的所有的特征值H如下：)end mH例4.4.1當(dāng)取=1.15, 5時(shí)，判別用超松弛迭代法解下列方程組產(chǎn)生的迭代序列是 否收斂？5 % + x 2 - x 3 一 2 x 4 = 42 x + 8 x + x + 3 x = 1 A=5 1 -1 -2;2 8 1 3;1 -

27、2 -4 -1;-1 3 2 7; H=ddpbj(A,1.15)運(yùn)行后輸出結(jié)果請(qǐng)注意：因?yàn)樽V半徑小于1，所以超松弛迭代序列收斂,譜半徑mH和B的所有 的特征值H如下：mH =0.1596H =0.1049 + 0.1203i 0.1049 - 0.1203i -0.1295 + 0.0556i -0.1295 - 0.0556i(2)當(dāng)取 =5時(shí)，然后在MATLAB工作窗口輸入程序 H=ddpbj(A, 5)運(yùn)行后輸出結(jié)果請(qǐng)注意：因?yàn)樽V半徑不小于1，所以超松弛迭代序列發(fā)散,譜半徑mH和B的所有的特征值H如下：mH =14.1082H =-14.1082-2.5107 0.5996 + 2.6

28、206i 0.5996 - 2.6206i4.4.3超松弛迭代法的MATLAB程序根據(jù)超松弛迭代公式(4.23)和定理4.5,現(xiàn)提供名為cscdd.m的M文件如下：用超松弛迭代法解線性方程組AX = b的MATLAB主程序輸入的量：線性方程組AX = b的系數(shù)矩陣A和b,初始向量X,范數(shù)的名稱P = 1, 2, inf,或fro.，松弛因子om,近似解X的誤差(精度)wucha和迭代的最大次數(shù)maxi.輸出的量：譜半徑mH，以系數(shù)矩陣A的對(duì)角元構(gòu)成的對(duì)角矩陣D、A的上三角形矩 陣，但對(duì)角元為0、A的下三角形矩陣乙，但對(duì)角元為0,迭代次數(shù)z，有關(guān)超松弛迭代收 斂性的相關(guān)信息及其AX = b的精確

29、解jX和近似解X.function X=cscdd (A,b,X,om,wucha,max1)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1); jX=Ab;n m=size(A);iD=inv(D-om*L); B2=iD*(om*U+(1-om)*D);H=eig(B2);mH=norm(H,inf);for k=1:max1iD=inv(D-om*L); B2=iD*(om*U+(1-om)*D);f2= om*iD*b; X1= B2*X+f2;X=X1; djwcX=norm(X1-jX,inf); xdwcX=djwcX/(norm(X,inf)+

30、eps);if (djwcXwucha)|(xdwcX max1disp(迭代次數(shù)已經(jīng)超過最大迭代次數(shù)max1，譜半徑mH，方程組的精確 解jX,迭代次數(shù)i如下：)mH,D,U,L,jX=jX, i=k-1,endendendif mH=1disp(請(qǐng)注意：因?yàn)樽V半徑不小于1,所以超松弛迭代序列發(fā)散.)disp(譜半徑mH，A的分解矩陣D,U,L和方程組的精確解jX,迭代次數(shù)i和迭代序列X如下：)i=k-1,mH,D,U,L,jX,elsedisp(因?yàn)樽V半徑小于1,所以超松弛迭代序列收斂，近似解乂如下：)end或function X=cscdd1 (A,b,X,om,wucha,max1)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1); jX=Ab;n m=size(A