清華大學(xué)數(shù)值分析實(shí)驗(yàn)報(bào)告

1、數(shù)值分析實(shí)驗(yàn)報(bào)告 一、 實(shí)驗(yàn)3.1題目：考慮線性方程組，編制一個(gè)能自動(dòng)選取主元，又能手動(dòng)選取主元的求解線性代數(shù)方程組的Gauss消去過(guò)程。（1）取矩陣，則方程有解。取計(jì)算矩陣的條件數(shù)。分別用順序Gauss消元、列主元Gauss消元和完全選主元Gauss消元方法求解，結(jié)果如何？（2）現(xiàn)選擇程序中手動(dòng)選取主元的功能，每步消去過(guò)程都選取模最小或按模盡可能小的元素作為主元進(jìn)行消元，觀察并記錄計(jì)算結(jié)果，若每步消去過(guò)程總選取按模最大的元素作為主元，結(jié)果又如何？分析實(shí)驗(yàn)的結(jié)果。（3）取矩陣階數(shù)n=20或者更大，重復(fù)上述實(shí)驗(yàn)過(guò)程，觀察記錄并分析不同的問(wèn)題及消去過(guò)程中選擇不同的主元時(shí)計(jì)算結(jié)果的差異，說(shuō)明主元素的

2、選取在消去過(guò)程中的作用。（4）選取其他你感興趣的問(wèn)題或者隨機(jī)生成的矩陣，計(jì)算其條件數(shù)，重復(fù)上述實(shí)驗(yàn)，觀察記錄并分析實(shí)驗(yàn)的結(jié)果。1. 算法介紹首先，分析各種算法消去過(guò)程的計(jì)算公式，順序高斯消去法：第k步消去中，設(shè)增廣矩陣中的元素（若等于零則可以判定系數(shù)矩陣為奇異矩陣，停止計(jì)算），則對(duì)k行以下各行計(jì)算，分別用乘以增廣矩陣的第行并加到第行，則可將增廣矩陣中第列中以下的元素消為零；重復(fù)此方法，從第1步進(jìn)行到第n-1步，則可以得到最終的增廣矩陣，即；列主元高斯消去法：第k步消去中，在增廣矩陣中的子方陣中，選取使得，當(dāng)時(shí)，對(duì)中第行與第行交換，然后按照和順序消去法相同的步驟進(jìn)行。重復(fù)此方法，從第1步進(jìn)行第n

3、-1步，就可以得到最終的增廣矩陣，即；完全主元高斯消去法：第k步消去中，在增廣矩陣中對(duì)應(yīng)的子方陣中，選取使得，若或，則對(duì)中第行與第行、第列與第列交換，然后按照和順序消去法相同的步驟進(jìn)行即可。重復(fù)此方法，從第1步進(jìn)行到第n-1步，就可以得到最終的增廣矩陣，即；接下來(lái)，分析回代過(guò)程求解的公式，容易看出，對(duì)上述任一種消元法，均有以下計(jì)算公式：2. 實(shí)驗(yàn)程序的設(shè)計(jì)一、輸入實(shí)驗(yàn)要求及初始條件；二、計(jì)算系數(shù)矩陣A的條件數(shù)及方程組的理論解；三、對(duì)各不同方法編程計(jì)算，并輸出最終計(jì)算結(jié)果。3. 計(jì)算結(jié)果及分析（1）先計(jì)算系數(shù)矩陣的條件數(shù)，結(jié)果如下，可知系數(shù)矩陣的條件數(shù)較大，故此問(wèn)題屬于病態(tài)問(wèn)題， b或A的擾動(dòng)都

4、可能引起解的較大誤差；采用順序高斯消去法，計(jì)算結(jié)果為：最終解為x=(1.0000, 1.0000, 1.0000, 1.0001, 0.9998, 1.0004, 0.9993, 1.0012, 0.9979, 1.0028)T使用無(wú)窮范數(shù)衡量誤差，得到=2.0401e-14，可以發(fā)現(xiàn)，采用順序高斯消元法求得的解與精確解之間誤差較小。通過(guò)進(jìn)一步觀察，可以發(fā)現(xiàn)，按照順序高斯消去法計(jì)算時(shí)，其選取的主元值和矩陣中其他元素大小相近，因此順序高斯消去法方式并沒(méi)有對(duì)結(jié)果造成特別大的影響。 若采用列主元高斯消元法，則結(jié)果為：最終解為x=(1.0000, 1.0000, 1.0000, 1.0000, 1.0

5、000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000)T同樣使用無(wú)窮范數(shù)衡量誤差，有=0； 若使用完全主元高斯消元法，則結(jié)果為最終解x=(1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000, 1.0000)T同樣使用無(wú)窮范數(shù)衡量誤差，有=0；（2）若每步都選取模最小或盡可能小的元素為主元，則計(jì)算結(jié)果為最終解x=(1.0000 1.0000 1.0000 1.0001 0.9998 1.0004 0.9993 1.0012 0.9979 1.0028)T使用無(wú)窮范數(shù)衡量誤差

6、，有為2.0401e-14；而完全主元消去法的誤差為=0。從（1）和（2）的實(shí)驗(yàn)結(jié)果可以發(fā)現(xiàn)，列主元消去法和完全主元消去法都得到了精確解，而順序高斯消去法和以模盡量小的元素為主元的消去法沒(méi)有得到精確解。在后兩種消去法中，由于程序計(jì)算時(shí)的舍入誤差，對(duì)最終結(jié)果產(chǎn)生了一定的影響，但由于方程組的維度較低，并且元素之間相差不大，所以誤差仍比較小。為進(jìn)一步分析，計(jì)算上述4種方法每步選取的主元數(shù)值，并列表進(jìn)行比較，結(jié)果如下：第n次消元順序列主元完全主元模最小16.0000886.000024.6667884.666734.4286884.428644.3333884.333354.2258884.22586

7、4.6032884.603274.6063884.606384.4902884.490294.4853884.4853104.30990.54690.54694.3099從上表可以發(fā)現(xiàn)，對(duì)這個(gè)方程組而言，順序高斯消去選取的主元恰好事模盡量小的元素，而由于列主元和完全主元選取的元素為8，與4在數(shù)量級(jí)上差別小，所以計(jì)算過(guò)程中的累積誤差也較小，最終4種方法的輸出結(jié)果均較為精確。在這里，具體解釋一下順序法與模最小法的計(jì)算結(jié)果完全一致的原因。該矩陣在消元過(guò)程中，每次選取主元的一列只有兩個(gè)非零元素，對(duì)角線上的元素為4左右，而其正下方的元素為8，該列其余位置的元素均為0。在這樣的情況下，默認(rèn)的主元也就是該列

8、最小的主元，因此兩種方法所得到的計(jì)算結(jié)果是一致的。理論上說(shuō)，完全高斯消去法的誤差最小，其次是列主元高斯消去法，而選取模最小的元素作為主元時(shí)的誤差最大，但是由于方程組的特殊性（元素相差不大并且維度不高），這個(gè)理論現(xiàn)象在這里并沒(méi)有充分體現(xiàn)出來(lái)。（3）時(shí)，重復(fù)上述實(shí)驗(yàn)過(guò)程，各種方法的計(jì)算結(jié)果如下所示，在這里，仍采用無(wú)窮范數(shù)衡量絕對(duì)誤差。順序高斯消去法列主元高斯消去完全主元高斯消去選取模最小或盡可能小元素作為主元消去X1.0000 1.0000 1.0000 1.0001 0.9998 1.0004 0.9993 1.0014 0.9972 1.0057 0.98861.0227 0.9547 1.0

9、902 0.8209 1.3524 0.3179 1.2732 0.8173 1.910211111111111111111111111111111111111111111.0000 1.0000 1.0000 1.0001 0.9998 1.0004 0.9993 1.0014 0.9972 1.0057 0.98861.0227 0.9547 1.0902 0.8209 1.3524 0.3179 1.2732 0.8173 1.91022.9430e-11002.9430e-11 可以看出，此時(shí)列主元和完全主元的計(jì)算結(jié)果仍為精確值，而順序高斯消去和模盡可能小方法仍然產(chǎn)生了一定的誤差，并且

10、兩者的誤差一致。與n=10時(shí)候的誤差比相比，n=20時(shí)的誤差增長(zhǎng)了大約1000倍，這是由于計(jì)算過(guò)程中舍入誤差的不斷累積所致。所以，如果進(jìn)一步增加矩陣的維數(shù)，應(yīng)該可以看出更明顯的現(xiàn)象。（4）不同矩陣維度下的誤差如下，在這里，為方便起見(jiàn)，選取2-條件數(shù)對(duì)不同維度的系數(shù)矩陣進(jìn)行比較。維度條件數(shù)順序消去列主元完全主元模盡量小1.7e+32.84e-14002.84e-141.8e+62.91e-11002.91e-115.7e+79.31e-10009.31e-101.8e+92.98e-08002.98e-081.9e+123.05e-05003.05e-053.8e+163.28e+043.88e

11、-123.88e-123.28e+048.5e+163.52e+134.2e-34.2e-33.52e+13 從上表可以看出，隨著維度的增加，不同方法對(duì)計(jì)算誤差的影響逐漸體現(xiàn)，并且增長(zhǎng)較快，這是由于舍入誤差逐步累計(jì)而造成的。不過(guò)，方法二與方法三在維度小于40的情況下都得到了精確解，這兩種方法的累計(jì)誤差遠(yuǎn)比方法一和方法四慢；同樣地，出于與前面相同的原因，方法一與方法四的計(jì)算結(jié)果保持一致，方法二與方法三的計(jì)算結(jié)果保持一致。4. 結(jié)論本文矩陣中的元素差別不大，模最大和模最小的元素并沒(méi)有數(shù)量級(jí)上的差異，因此，不同的主元選取方式對(duì)計(jì)算結(jié)果的影響在維度較低的情況下并不明顯，四種方法都足夠精確。對(duì)比四種方法

12、，可以發(fā)現(xiàn)采用列主元高斯消去或者完全主元高斯消去法，可以盡量抑制誤差，算法最為精確。不過(guò)，對(duì)于低階的矩陣來(lái)說(shuō)，四種方法求解出來(lái)的結(jié)果誤差均較小。另外，由于完全選主元方法在選主元的過(guò)程中計(jì)算量較大，而且可以發(fā)現(xiàn)列主元法已經(jīng)可以達(dá)到很高的精確程度，因而在實(shí)際計(jì)算中可以選用列主元法進(jìn)行計(jì)算。附錄：程序代碼clearclc;format long;%方法選擇n=input(矩陣A階數(shù):n=);disp(選取求解方式);disp(1 順序Gauss消元法，2 列主元Gauss消元法，3 完全選主元Gauss消元法，4 模最小或近可能小的元素作為主元);a=input(求解方式序號(hào)：);%賦值A(chǔ)和bA=z

13、eros(n,n);b=zeros(n,1);for i=1:n A(i,i)=6; if i1 A(i,i-1)=8; end if i1，所以這是一個(gè)病態(tài)問(wèn)題。由上表可以看出，四種方法的求解都存在一定的誤差。下面分析誤差的來(lái)源：LU分解方法的誤差存在主要是由于Hilbert矩陣各元素由分?jǐn)?shù)形式轉(zhuǎn)換為小數(shù)形式時(shí)，不能除盡情況下會(huì)出現(xiàn)舍入誤差，在進(jìn)行LU分解時(shí)也存在這個(gè)問(wèn)題，所以最后得到的結(jié)果不是方程的精確解 ，但結(jié)果顯示該方法的誤差非常??；Jacobi迭代矩陣的譜半徑為4.30853，故此迭代法不收斂；GS迭代法在迭代次數(shù)為541次時(shí)得到了方程的近似解，其誤差約為0.05 ，比較大。GS迭代

14、矩陣的譜半徑為0.，很接近1，所以GS迭代法收斂速度較慢；SOR迭代法在迭代次數(shù)為100次時(shí)誤差約為0.08，誤差較大。SOR迭代矩陣的譜半徑為0.，也很接近1，所以時(shí)SOR迭代法收斂速度不是很快，但是相比于GS法，在迭代速度方面已經(jīng)有了明顯的提高；另外，對(duì)不同的，SOR方法的迭代速度會(huì)相應(yīng)有變化，如果選用最佳松弛因子，可以實(shí)現(xiàn)更快的收斂；（2）考慮不同維度的情況，時(shí)，算法Gauss消去J法GS法SOR法（w=0.5）計(jì)算結(jié)果0.6269 1.9060 0.2676 1.8103 0.4116 1.2164 0.7125 1.4747-0.53491.68390.60151.50361.257

15、61.66120.53640.55361.56340.51710.30811.51541.39741.55140.84690.7655迭代次數(shù)-356100譜半徑-6.0421310.8776-時(shí)，算法Gauss消去法Jacobi法GS法SOR法（w=0.5）計(jì)算結(jié)果0.1197 1.6354 0.3700 1.8631 0.4329 1.7805 0.9370 1.1556 0.5812 1.9020 0.5643-0.67511.84090.67390.32370.13531.65911.66081.85331.94220.62500.06971.65720.44830.79100.272

16、71.35371.64681.50691.11920.24240.41420.1922迭代次數(shù)-1019100譜半徑-8.6496410.9966-時(shí)算法Gauss消去法Jacobi法GS法SOR法（w=0.5）計(jì)算結(jié)果0.3799 1.4896 0.9769 0.7706 1.9297 0.5194 2.5890 -2.2128 7.1018 -7.9081 7.0484-1.5142 0.7065 1.1601 0.3362-不收斂1.15901.75540.77290.90530.84571.79031.87641.39441.97681.60561.51760.20410.80030.

17、40650.3220迭代次數(shù)-262100譜半徑-6.0421311.00008.9082-0.6780分析以上結(jié)果可以發(fā)現(xiàn)，隨著n值的增加，Gauss消去法誤差逐漸增大，而且誤差增大的速度很快，在維數(shù)小于等于10情況下，Gauss消去法得到的結(jié)果誤差較?。坏?dāng)維數(shù)達(dá)到15時(shí)，計(jì)算結(jié)果誤差已經(jīng)達(dá)到精確解的很多倍；J法迭代不收斂，無(wú)論n如何取值，其譜半徑始終大于1，因而J法不收斂，所以J迭代法不能用于Hilbert矩陣的求解；對(duì)于GS迭代法和SOR迭代法，兩種方法均收斂，GS迭代法是SOR迭代法松弛因子取值為1的特例，SOR方法受到取值的影響，會(huì)有不同的收斂情況?？梢缘贸鯣S迭代矩陣的譜半徑小于

18、1但是很接近1，收斂速度很慢。雖然隨著維數(shù)的增大，所需迭代的次數(shù)逐漸減少，但是當(dāng)維數(shù)達(dá)到15的時(shí)候，GS法已經(jīng)不再收斂。因此可以得出結(jié)論，GS迭代方法在Hilbert矩陣維數(shù)較低時(shí)，能夠在一定程度上滿足迭代求解的需求，不過(guò)迭代的速度很慢。另外，隨著矩陣維數(shù)的增加， SOR法的誤差水平基本穩(wěn)定，而且誤差在可以接受的范圍之內(nèi)。經(jīng)過(guò)比較可以得出結(jié)論，如果求解較低維度的Hibert矩陣問(wèn)題，Gauss消去法、GS迭代法和SOR迭代法均可使用，且Gauss消去法的結(jié)果精確度較高；如果需要求解較高維度的Hibert矩陣問(wèn)題，只有采用SOR迭代法。（3）系數(shù)矩陣的條件數(shù)較大時(shí)，為病態(tài)方程。由實(shí)驗(yàn)可知，Gau

19、ss法在解上述方程時(shí)，結(jié)果存在很大的誤差。而對(duì)于收斂的迭代法，可以通過(guò)選取最優(yōu)松弛因子的方法來(lái)求解，雖然迭代次數(shù)相對(duì)較多，但是結(jié)果較為精確?？傮w來(lái)看，對(duì)于一般病態(tài)方程組的求解，可以采用以下方式：1. 低維度下采用Gauss消去法直接求解是可行的；Jacobi迭代方法不適宜于求解病態(tài)問(wèn)題；GS迭代方法可以解決維數(shù)較低的病態(tài)問(wèn)題，但其譜半徑非常趨近于1，導(dǎo)致迭代算法收斂速度很慢，維數(shù)較大的時(shí)候，GS法也不再收斂；SOR方法較適合于求解病態(tài)問(wèn)題，特別是矩陣維數(shù)較高的時(shí)候，其優(yōu)勢(shì)更為明顯。2. 采用高精度的運(yùn)算，如選用雙倍或更多倍字長(zhǎng)的運(yùn)算，可以提高收斂速度；3. 可以對(duì)原方程組作某些預(yù)處理，從而有效

20、降低系數(shù)矩陣的條件數(shù)。4. 實(shí)驗(yàn)結(jié)論（1）對(duì)Hibert矩陣問(wèn)題，其條件數(shù)會(huì)隨著維度的增加迅速增加，病態(tài)性會(huì)越來(lái)越明顯；在維度較低的時(shí)候，Gauss消去法、GS迭代法和SOR迭代法均可使用，且可以優(yōu)先使用Gauss消去法；如果需要求解較高維度的Hibert矩陣問(wèn)題，只有SOR迭代法能夠求解。（2）SOR方法比較適合于求解病態(tài)問(wèn)題，特別是矩陣維數(shù)較高的時(shí)候，其優(yōu)點(diǎn)更為明顯。從本次實(shí)驗(yàn)可以看出，隨著矩陣維數(shù)的增大，SOR方法所需的迭代次數(shù)減少，而且誤差基本穩(wěn)定，是解決病態(tài)問(wèn)題的適宜方法。附錄：程序代碼clear allclc;format long;%矩陣賦值n=input(矩陣H的階數(shù):n=);

21、for i=1:n for j=1:n H(i,j)=1/(i+j-1); endendb=H*ones(n,1);disp(H矩陣為：);Hdisp(向量b：);b%方法選擇disp(選取求解方式);disp(1 Gauss消去法，2 J迭代法，3 GS迭代法，4 SOR迭代法);a=input(求解方式序號(hào)：);%Gauss消去法if a=1; H1=H;b1=b; for k=1:n if H1(k,k)=0 disp(主元為零，Gauss消去法無(wú)法進(jìn)行); break end fprintf(第%d次消元所選取的主元是：%gn,k,H1(k,k) for p=k+1:n m5=-H1(

22、p,k)/H1(k,k); H1(p,k:n)=H1(p,k:n)+m5*H1(k,k:n); b1(p)=b1(p)+m5*b1(k); end endx1(n)=b1(n)/H1(n,n);for k=n-1:-1:1 for v=k+1:n b1(k)=b1(k)-H1(k,v)*x1(v); end x1(k)=b1(k)/H1(k,k);end disp(Gauss消去法解為:); disp(x1); disp(解與精確解之差的無(wú)窮范數(shù)); norm(x1-a),inf)endD=diag(diag(H);L=-tril(H,-1);U=-triu(H,1); %J迭代法if a=2

23、; %給定初始x0 ini=input(初始值設(shè)定:x0=); x0(:,1)=ini*diag(ones(n); disp(初始解向量為:); x0 xj(:,1)=x0(:,1); B=(D(-1)*(L+U); f=(D(-1)*b; fprintf(J法B矩陣譜半徑為：%gn,vrho(B); if vrho(B)1; for m2=1:5000 xj(:,m2+1)=B*xj(:,m2)+fj; if norm(xj(:,m2+1)-xj(:,m2),inf)0.0001 break end end disp(J法計(jì)算結(jié)果為:); xj(:,m2+1) disp(解與精確解之差的無(wú)窮范數(shù)); norm(xj(:,m2+1)-diag(ones(n),inf) disp(J迭代法迭代次數(shù):); m2 else disp(由于B矩陣譜半徑大于1，因而J法不收斂); endend %GS迭代法if a=3; %給定初始x0 ini=input(初始值設(shè)定:x0=); x0(:,1)=ini*diag(ones(n); disp(初始解向量為:); x0 xG(:,1)=x0(:,1); G=inv(D-L)*U; fG=inv(D-L)*b; fprintf(GS法G矩陣譜半徑為：%gn,vrho(G); if vr