數(shù)值分析實(shí)驗(yàn)報(bào)告
一、實(shí)驗(yàn)3.1
題目：
考慮線性方程組，，編制一個(gè)能自動選取主元，又能手動選取主元的求解線性代數(shù)方程組的Gauss消去過程。
（1）取矩陣，則方程有解。取計(jì)算矩陣的條件數(shù)。分別用順序Gauss消元、列主元Gauss消元和完全選主元Gauss消元方法求解，結(jié)果如何？
（2）現(xiàn)選擇程序中手動選取主元的功能，每步消去過程都選取模最小或按模盡可能小的元素作為主元進(jìn)行消元，觀察并記錄計(jì)算結(jié)果，若每步消去過程總選取按模最大的元素作為主元，結(jié)果又如何？分析實(shí)驗(yàn)的結(jié)果。
（3）取矩陣階數(shù)n=20或者更大，重復(fù)上述實(shí)驗(yàn)過程，觀察記錄并分析不同的問題及消去過程中選擇不同的主元時(shí)計(jì)算結(jié)果的差異，說明主元素的選取在消去過程中的作用。
（4）選取其他你感興趣的問題或者隨機(jī)生成的矩陣，計(jì)算其條件數(shù)，重復(fù)上述實(shí)驗(yàn)，觀察記錄并分析實(shí)驗(yàn)的結(jié)果。
1.算法介紹
首先，分析各種算法消去過程的計(jì)算公式，順序高斯消去法：
第k步消去中，設(shè)增廣矩陣中的元素（若等于零則可以判定系數(shù)矩陣為奇異矩陣，停止計(jì)算），則對k行以下各行計(jì)算，分別用乘以增廣矩陣的第行并加到第行，則可將增廣矩陣中第列中以下的元素消為零；重復(fù)此方法，從第1步進(jìn)行到第n-1步，則可以得到最終的增廣矩陣，即；
列主元高斯消去法：
第k步消去中，在增廣矩陣中的子方陣中，選取使得，當(dāng)時(shí)，對中第行與第行交換，然后按照和順序消去法相同的步驟進(jìn)行。重復(fù)此方法，從第1步進(jìn)行第n-1步，就可以得到最終的增廣矩陣，即；
完全主元高斯消去法：
第k步消去中，在增廣矩陣中對應(yīng)的子方陣中，選取使得，若或，則對中第行與第行、第列與第列交換，然后按照和順序消去法相同的步驟進(jìn)行即可。重復(fù)此方法，從第1步進(jìn)行到第n-1步，就可以得到最終的增廣矩陣，即；
接下來，分析回代過程求解的公式，容易看出，對上述任一種消元法，均有以下計(jì)算公式：
2.實(shí)驗(yàn)程序的設(shè)計(jì)
一、輸入實(shí)驗(yàn)要求及初始條件；
二、計(jì)算系數(shù)矩陣A的條件數(shù)及方程組的理論解；
三、對各不同方法編程計(jì)算，并輸出最終計(jì)算結(jié)果。
3.計(jì)算結(jié)果及分析
（1）
先計(jì)算系數(shù)矩陣的條件數(shù)，結(jié)果如下，可知系數(shù)矩陣的條件數(shù)較大，故此問題屬于病態(tài)問題，b或A的擾動都可能引起解的較大誤差；
采用順序高斯消去法，計(jì)算結(jié)果為：
最終解為x=(1.***,1.***,1.***,1.***,0.***,1.***,0.***,1.***,0.***,1.***)T
使用無窮范數(shù)衡量誤差，得到=2.842***1e-14，可以發(fā)現(xiàn)，采用順序高斯消元法求得的解與精確解之間誤差較小。通過進(jìn)一步觀察，可以發(fā)現(xiàn)，按照順序高斯消去法計(jì)算時(shí)，其選取的主元值和矩陣中其他元素大小相近，因此順序高斯消去法方式并沒有對結(jié)果造成特別大的影響。
若采用列主元高斯消元法，則結(jié)果為：
最終解為x=(1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***)T
同樣使用無窮范數(shù)衡量誤差，有=0；
若使用完全主元高斯消元法，則結(jié)果為
最終解x=(1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***,1.***)T
同樣使用無窮范數(shù)衡量誤差，有=0；
（2）
若每步都選取模最小或盡可能小的元素為主元，則計(jì)算結(jié)果為
最終解x=(1.***
1.***
1.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***)T
使用無窮范數(shù)衡量誤差，有為2.842***1e-14；而完全主元消去法的誤差為=0。
從（1）和（2）的實(shí)驗(yàn)結(jié)果可以發(fā)現(xiàn)，列主元消去法和完全主元消去法都得到了精確解，而順序高斯消去法和以模盡量小的元素為主元的消去法沒有得到精確解。在后兩種消去法中，由于程序計(jì)算時(shí)的舍入誤差，對最終結(jié)果產(chǎn)生了一定的影響，但由于方程組的維度較低，并且元素之間相差不大，所以誤差仍比較小。
為進(jìn)一步分析，計(jì)算上述4種方法每步選取的主元數(shù)值，并列表進(jìn)行比較，結(jié)果如下：
第n次消元
順序
列主元
完全主元
模最小
6.***
6.***
4.***
4.***
4.***
4.***
4.***3333
4.***3333
4.***
4.***
4.***
4.***
4.0***063
4.0***063
4.***
4.***
4.0039***
4.0039***
4.***
0.0***469
0.0***469
4.***
從上表可以發(fā)現(xiàn)，對這個(gè)方程組而言，順序高斯消去選取的主元恰好事模盡量小的元素，而由于列主元和完全主元選取的元素為8，與4在數(shù)量級上差別小，所以計(jì)算過程中的累積誤差也較小，最終4種方法的輸出結(jié)果均較為精確。
在這里，具體解釋一下順序法與模最小法的計(jì)算結(jié)果完全一致的原因。該矩陣在消元過程中，每次選取主元的一列只有兩個(gè)非零元素，對角線上的元素為4左右，而其正下方的元素為8，該列其余位置的元素均為0。在這樣的情況下，默認(rèn)的主元也就是該列最小的主元，因此兩種方法所得到的計(jì)算結(jié)果是一致的。
理論上說，完全高斯消去法的誤差最小，其次是列主元高斯消去法，而選取模最小的元素作為主元時(shí)的誤差最大，但是由于方程組的特殊性（元素相差不大并且維度不高），這個(gè)理論現(xiàn)象在這里并沒有充分體現(xiàn)出來。
（3）
時(shí)，重復(fù)上述實(shí)驗(yàn)過程，各種方法的計(jì)算結(jié)果如下所示，在這里，仍采用無窮范數(shù)衡量絕對誤差。
順序高斯消去法
列主元高斯消去
完全主元高斯消去
選取模最小或盡可能小元素作為主元消去
X
1.***
1.***
1.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
1.***
1.***
1.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
2.***e-11
0
0
2.***e-11
可以看出，此時(shí)列主元和完全主元的計(jì)算結(jié)果仍為精確值，而順序高斯消去和模盡可能小方法仍然產(chǎn)生了一定的誤差，并且兩者的誤差一致。與n=10時(shí)候的誤差比相比，n=20時(shí)的誤差增長了大約1000倍，這是由于計(jì)算過程中舍入誤差的不斷累積所致。所以，如果進(jìn)一步增加矩陣的維數(shù)，應(yīng)該可以看出更明顯的現(xiàn)象。
（4）
不同矩陣維度下的誤差如下，在這里，為方便起見，選取2-條件數(shù)對不同維度的系數(shù)矩陣進(jìn)行比較。
維度
條件數(shù)
順序消去
列主元
完全主元
模盡量小
1.7e+3
2.84e-14
0
0
2.84e-14
1.8e+6
2.91e-11
0
0
2.91e-11
5.7e+7
9.31e-10
0
0
9.31e-10
1.8e+9
2.98e-08
0
0
2.98e-08
1.9e+12
3.05e-05
0
0
3.05e-05
3.8e+16
3.28e+04
3.88e-12
3.88e-12
3.28e+04
8.5e+16
3.52e+13
4.2e-3
4.2e-3
3.52e+13
從上表可以看出，隨著維度的增加，不同方法對計(jì)算誤差的影響逐漸體現(xiàn)，并且增長較快，這是由于舍入誤差逐步累計(jì)而造成的。不過，方法二與方法三在維度小于40的情況下都得到了精確解，這兩種方法的累計(jì)誤差遠(yuǎn)比方法一和方法四慢；同樣地，出于與前面相同的原因，方法一與方法四的計(jì)算結(jié)果保持一致，方法二與方法三的計(jì)算結(jié)果保持一致。
4.結(jié)論
本文矩陣中的元素差別不大，模最大和模最小的元素并沒有數(shù)量級上的差異，因此，不同的主元選取方式對計(jì)算結(jié)果的影響在維度較低的情況下并不明顯，四種方法都足夠精確。
對比四種方法，可以發(fā)現(xiàn)采用列主元高斯消去或者完全主元高斯消去法，可以盡量抑制誤差，算法最為精確。不過，對于低階的矩陣來說，四種方法求解出來的結(jié)果誤差均較小。
另外，由于完全選主元方法在選主元的過程中計(jì)算量較大，而且可以發(fā)現(xiàn)列主元法已經(jīng)可以達(dá)到很高的精確程度，因而在實(shí)際計(jì)算中可以選用列主元法進(jìn)行計(jì)算。
附錄：程序代碼
clear
clc;
format
long;
%方法選擇
n=input(＇矩陣A階數(shù):n=＇);
disp(＇選取求解方式＇);
disp(＇1
順序Gauss消元法，2
列主元Gauss消元法，3
完全選主元Gauss消元法，4
模最小或近可能小的元素作為主元＇);
a=input(＇求解方式序號：＇);
%賦值A(chǔ)和b
A=zeros(n,n);
b=zeros(n,1);
for
i=1:n
A(i,i)=6;
if
i>1
A(i,i-1)=8;
end
if
i A(i,i+1)=1; end end for i=1:n for j=1:n b(i)=b(i)+A(i,j); end end disp(＇給定系數(shù)矩陣為：＇); A disp(＇右端向量為：＇); b %求條件數(shù)及理論解 disp(＇線性方程組的精確解：＇); X=(A＼b)＇ fprintf(＇矩陣A的1-條件數(shù): %f ＼n＇,cond(A,1)); fprintf(＇矩陣A的2-條件數(shù): %f ＼n＇,cond(A)); fprintf(＇矩陣A的無窮-條件數(shù): %f ＼n＇,cond(A,inf)); %順序Gauss消元法 if a==1 A1=A;b1=b; for k=1:n if A1(k,k)==0 disp(＇主元為零，順序Gauss消元法無法進(jìn)行＇); break end fprintf(＇第%d次消元所選取的主元：%g＼n＇,k,A1(k,k)) %disp(＇此次消元后系數(shù)矩陣為：＇); %A1 for p=k+1:n l=A1(p,k)/A1(k,k); A1(p,k:n)=A1(p,k:n)-l*A1(k,k:n); b1(p)=b1(p)-l*b1(k); end end x1(n)=b1(n)/A1(n,n); for k=n-1:-1:1 for w=k+1:n b1(k)=b1(k)-A1(k,w)*x1(w); end x1(k)=b1(k)/A1(k,k); end disp(＇順序Gauss消元法解為:＇); disp(x1); disp(＇所求解與精確解之差的無窮-范數(shù)為＇); norm(x1-X,inf) end %列主元Gauss消元法 if a==2 A2=A;b2=b; for k=1:n [max_i,max_j]=find(A2(:,k)==max(abs(A2(k:n,k)))); if max_i~=k A2_change=A2(k,:); A2(k,:)=A2(max_i,:); A2(max_i,:)=A2_change; b2_change=b2(k); b2(k)=b2(max_i); b2(max_i)=b2_change; end if A2(k,k)==0 disp(＇主元為零，列主元Gauss消元法無法進(jìn)行＇); break end fprintf(＇第%d次消元所選取的主元：%g＼n＇,k,A2(k,k)) %disp(＇此次消元后系數(shù)矩陣為：＇); %A2 for p=k+1:n l=A2(p,k)/A2(k,k); A2(p,k:n)=A2(p,k:n)-l*A2(k,k:n); b2(p)=b2(p)-l*b2(k); end end x2(n)=b2(n)/A2(n,n); for k=n-1:-1:1 for w=k+1:n b2(k)=b2(k)-A2(k,w)*x2(w); end x2(k)=b2(k)/A2(k,k); end disp(＇列主元Gauss消元法解為:＇); disp(x2); disp(＇所求解與精確解之差的無窮-范數(shù)為＇); norm(x2-X,inf) end %完全選主元Gauss消元法 if a==3 A3=A;b3=b; for k=1:n VV=eye(n); [max_i,max_j]=find(A3(k:n,k:n)==max(max(abs(A3(k:n,k:n))))); if numel(max_i)==0 [max_i,max_j]=find(A3(k:n,k:n)==-max(max(abs(A3(k:n,k:n))))); end W=eye(n); W(max_i(1)+k-1,max_i(1)+k-1)=0; W(k,k)=0; W(max_i(1)+k-1,k)=1; W(k,max_i(1)+k-1)=1; V=eye(n); V(k,k)=0; V(max_j(1)+k-1,max_j(1)+k-1)=0; V(k,max_j(1)+k-1)=1; V(max_j(1)+k-1,k)=1; A3=W*A3*V; b3=W*b3; VV=VV*V; if A3(k,k)==0 disp(＇主元為零，完全選主元Gauss消元法無法進(jìn)行＇); break end fprintf(＇第%d次消元所選取的主元：%g＼n＇,k,A3(k,k)) %disp(＇此次消元后系數(shù)矩陣為：＇); %A3 for p=k+1:n l=A3(p,k)/A3(k,k); A3(p,k:n)=A3(p,k:n)-l*A3(k,k:n); b3(p)=b3(p)-l*b3(k); end end x3(n)=b3(n)/A3(n,n); for k=n-1:-1:1 for w=k+1:n b3(k)=b3(k)-A3(k,w)*x3(w); end x3(k)=b3(k)/A3(k,k); end disp(＇完全選主元Gauss消元法解為:＇); disp(x3); disp(＇所求解與精確解之差的無窮-范數(shù)為＇); norm(x3-X,inf) end %模最小或近可能小的元素作為主元 if a==4 A4=A;b4=b; for k=1:n AA=A4; AA(AA==0)=NaN; [min_i,j]=find(AA(k:n,k)==min(abs(AA(k:n,k)))); if numel(min_i)==0 [min_i,j]=find(AA(k:n,k)==-min(abs(AA(k:n,k:n)))); end W=eye(n); W(min_i(1)+k-1,min_i(1)+k-1)=0; W(k,k)=0; W(min_i(1)+k-1,k)=1; W(k,min_i(1)+k-1)=1; A4=W*A4; b4=W*b4; if A4(k,k)==0 disp(＇主元為零，模最小Gauss消元法無法進(jìn)行＇); break end fprintf(＇第%d次消元所選取的主元：%g＼n＇,k,A4(k,k)) %A4 for p=k+1:n l=A4(p,k)/A4(k,k); A4(p,k:n)=A4(p,k:n)-l*A4(k,k:n); b4(p)=b4(p)-l*b4(k); end end x4(n)=b4(n)/A4(n,n); for k=n-1:-1:1 for w=k+1:n b4(k)=b4(k)-A4(k,w)*x4(w); end x4(k)=b4(k)/A4(k,k); end disp(＇模最小Gauss消元法解為:＇); disp(x4); disp(＇所求解與精確解之差的無窮-范數(shù)為＇); norm(x4-X,inf) end 二、實(shí)驗(yàn)3.3 題目： 考慮方程組的解，其中系數(shù)矩陣H為Hilbert矩陣： 這是一個(gè)著名的病態(tài)問題。通過首先給定解（例如取為各個(gè)分量均為1）再計(jì)算出右端的辦法給出確定的問題。 （1）選擇問題的維數(shù)為6，分別用Gauss消去法（即LU分解）、J迭代法、GS迭代法和SOR迭代法求解方程組，其各自的結(jié)果如何？將計(jì)算結(jié)果與問題的解比較，結(jié)論如何。 （2）逐步增大問題的維數(shù)，仍用上述的方法來解它們，計(jì)算的結(jié)果如何？計(jì)算的結(jié)果說明的什么？ （3）討論病態(tài)問題求解的算法。 1.算法設(shè)計(jì) 對任意線性方程組，分析各種方法的計(jì)算公式如下，（1）Gauss消去法： 首先對系數(shù)矩陣進(jìn)行LU分解，有，則原方程轉(zhuǎn)化為，令，則原方程可以分為兩步回代求解： 具體方法這里不再贅述。 （2）J迭代法： 首先分解，再構(gòu)造迭代矩陣，其中，進(jìn)行迭代計(jì)算，直到誤差滿足要求。 （3）GS迭代法： 首先分解，再構(gòu)造迭代矩陣，其中，進(jìn)行迭代計(jì)算，直到誤差滿足要求。 （4）SOR迭代法： 首先分解，再構(gòu)造迭代矩陣，其中，進(jìn)行迭代計(jì)算，直到誤差滿足要求。 2.實(shí)驗(yàn)過程 一、根據(jù)維度n確定矩陣H的各個(gè)元素和b的各個(gè)分量值； 二、選擇計(jì)算方法（Gauss消去法，J迭代法，GS迭代法，SOR迭代法），對迭代法設(shè)定初值，此外SOR方法還需要設(shè)定松弛因子； 三、進(jìn)行計(jì)算，直至滿足誤差要求（對迭代法，設(shè)定相鄰兩次迭代結(jié)果之差的無窮范數(shù)小于0.0001；
3.計(jì)算結(jié)果及分析
（1）時(shí)，問題可以具體定義為
計(jì)算結(jié)果如下，Gauss消去法
第1次消元所選取的主元是：1
第2次消元所選取的主元是：0.0833333
第3次消元所選取的主元是：0.00555556
第4次消元所選取的主元是：0.000357143
第5次消元所選取的主元是：2.26757e-05
第6次消元所選取的主元是：1.43155e-06
解得X=(0.***
1.***
0.***
1.***
0.***
1.***)T
使用無窮范數(shù)衡量誤差，可得=4.***e-10；
J迭代法
設(shè)定迭代初值為零，計(jì)算得到
J法的迭代矩陣B的譜半徑為4.30853＞1，所以J法不收斂；
GS迭代法
設(shè)定迭代初值為零，計(jì)算得到GS法的迭代矩陣G的譜半徑為：0.999998＜1，故GS法收斂，經(jīng)過541次迭代計(jì)算后，結(jié)果為X=(1.001***6
0.***
0.***
1.***
1.***
0.***)T
使用無窮范數(shù)衡量誤差，有=0.***；
SOR迭代法
設(shè)定迭代初值為零向量，并設(shè)定，計(jì)算得到SOR法迭代矩陣譜半徑為0.***，經(jīng)過100次迭代后的計(jì)算結(jié)果為
X=(1.***
0.***
1.03***59
1.06***81
1.***
0.9***527)T；
使用無窮范數(shù)衡量誤差，有=0.***；
對SOR方法，可變，改變值，計(jì)算結(jié)果可以列表如下
迭代次數(shù)
迭代矩陣的譜半徑
0.***
0.***
0.***
0.***
X
1.***
0.***
1.01***40
1.***
1.0***681
0.***
1.***
0.***
1.***
1.***
1.***
0.***
1.***
0.***
1.***
1.***
0.***
0.***
1.05***66
0.***
1.***
0.***
1.***
0.***
0.***
0.***
0.***
0.***
可以發(fā)現(xiàn)，松弛因子的取值對迭代速度造成了不同的影響，上述四種方法中，松弛因子=0.5時(shí)，收斂相對較快。
綜上，四種算法的結(jié)果列表如下：
算法
Gauss消去法
Jacobi法
GS法
SOR法（取）
迭代次數(shù)
--
不收斂
541
迭代矩陣的譜半徑
--
4.30853
0.999998
0.***
X
0.***
1.***
0.***
1.***
0.***
1.***
--
1.001***6
0.***
0.***
1.***
1.***
0.***
1.***
0.***
1.03***59
1.06***81
1.***
0.9***527
4.***e-10
--
0.***
0.***
計(jì)算可得，矩陣H的條件數(shù)為>>1，所以這是一個(gè)病態(tài)問題。由上表可以看出，四種方法的求解都存在一定的誤差。下面分析誤差的來源：
LU分解方法的誤差存在主要是由于Hilbert矩陣各元素由分?jǐn)?shù)形式轉(zhuǎn)換為小數(shù)形式時(shí)，不能除盡情況下會出現(xiàn)舍入誤差，在進(jìn)行LU分解時(shí)也存在這個(gè)問題，所以最后得到的結(jié)果不是方程的精確解，但結(jié)果顯示該方法的誤差非常??；
Jacobi迭代矩陣的譜半徑為4.30853，故此迭代法不收斂；
GS迭代法在迭代次數(shù)為541次時(shí)得到了方程的近似解，其誤差約為0.05，比較大。GS迭代矩陣的譜半徑為0.999998，很接近1，所以GS迭代法收斂速度較慢；
SOR迭代法在迭代次數(shù)為100次時(shí)誤差約為0.08，誤差較大。SOR迭代矩陣的譜半徑為0.999999，也很接近1，所以時(shí)SOR迭代法收斂速度不是很快，但是相比于GS法，在迭代速度方面已經(jīng)有了明顯的提高；另外，對不同的，SOR方法的迭代速度會相應(yīng)有變化，如果選用最佳松弛因子，可以實(shí)現(xiàn)更快的收斂；
（2）
考慮不同維度的情況，時(shí)，算法
Gauss消去
J法
GS法
SOR法（w=0.5）
計(jì)算結(jié)果
0.***
1.***
0.***
1.***
0.***
1.***
0.***
1.***
--
0.***
1.***
0.***
1.***
1.***
1.***
0.9968***
0.***
1.***
0.9397***
0.***
1.***
1.***
1.***
0.***
0.***
迭代次數(shù)
--
--
356
譜半徑
--
6.04213
0.***
--
時(shí)，算法
Gauss消去法
Jacobi法
GS法
SOR法（w=0.5）
計(jì)算結(jié)果
0.***
1.***
0.***
1.000***1
0.***
1.***
0.***
1.***
0.***
1.***
0.***
--
0.***
1.***
0.***
0.***
0.***
1.02***91
1.***
1.***
1.***
0.***
0.947***7
1.0***572
0.***
0.***
0.***
1.***
1.***
1.***
1.***
0.***
0.***
0.***
迭代次數(shù)
--
--
1019
譜半徑
--
8.64964
0.***
--
時(shí)
算法
Gauss消去法
Jacobi法
GS法
SOR法（w=0.5）
計(jì)算結(jié)果
0.***
1.***
0.***
0.***
1.***
0.***
2.***
-2.***
7.***
-7.***
7.***
-1.***
0.***
1.***
0.***
--
不收斂
1.***
1.***
0.907***9
0.***
0.***
1.***
1.09***64
1.***
1.***
1.***
1.0385***
0.***
0.942***3
0.***
0.***
迭代次數(shù)
--
--
262
譜半徑
--
6.04213
＞1
1.***
8.***
--
--
0.***
分析以上結(jié)果可以發(fā)現(xiàn)，隨著n值的增加，Gauss消去法誤差逐漸增大，而且誤差增大的速度很快，在維數(shù)小于等于10情況下，Gauss消去法得到的結(jié)果誤差較小；但當(dāng)維數(shù)達(dá)到15時(shí)，計(jì)算結(jié)果誤差已經(jīng)達(dá)到精確解的很多倍；
J法迭代不收斂，無論n如何取值，其譜半徑始終大于1，因而J法不收斂，所以J迭代法不能用于Hilbert矩陣的求解；
對于GS迭代法和SOR迭代法，兩種方法均收斂，GS迭代法是SOR迭代法松弛因子取值為1的特例，SOR方法受到取值的影響，會有不同的收斂情況?？梢缘贸鯣S迭代矩陣的譜半徑小于1但是很接近1，收斂速度很慢。雖然隨著維數(shù)的增大，所需迭代的次數(shù)逐漸減少，但是當(dāng)維數(shù)達(dá)到15的時(shí)候，GS法已經(jīng)不再收斂。因此可以得出結(jié)論，GS迭代方法在Hilbert矩陣維數(shù)較低時(shí)，能夠在一定程度上滿足迭代求解的需求，不過迭代的速度很慢。另外，隨著矩陣維數(shù)的增加，SOR法的誤差水平基本穩(wěn)定，而且誤差在可以接受的范圍之內(nèi)。
經(jīng)過比較可以得出結(jié)論，如果求解較低維度的Hibert矩陣問題，Gauss消去法、GS迭代法和SOR迭代法均可使用，且Gauss消去法的結(jié)果精確度較高；如果需要求解較高維度的Hibert矩陣問題，只有采用SOR迭代法。
（3）
系數(shù)矩陣的條件數(shù)較大時(shí)，為病態(tài)方程。由實(shí)驗(yàn)可知，Gauss法在解上述方程時(shí)，結(jié)果存在很大的誤差。而對于收斂的迭代法，可以通過選取最優(yōu)松弛因子的方法來求解，雖然迭代次數(shù)相對較多，但是結(jié)果較為精確。
總體來看，對于一般病態(tài)方程組的求解，可以采用以下方式：
1.低維度下采用Gauss消去法直接求解是可行的；
Jacobi迭代方法不適宜于求解病態(tài)問題；
GS迭代方法可以解決維數(shù)較低的病態(tài)問題，但其譜半徑非常趨近于1，導(dǎo)致迭代算法收斂速度很慢，維數(shù)較大的時(shí)候，GS法也不再收斂；
SOR方法較適合于求解病態(tài)問題，特別是矩陣維數(shù)較高的時(shí)候，其優(yōu)勢更為明顯。
2.采用高精度的運(yùn)算，如選用雙倍或更多倍字長的運(yùn)算，可以提高收斂速度；
3.可以對原方程組作某些預(yù)處理，從而有效降低系數(shù)矩陣的條件數(shù)。
4.實(shí)驗(yàn)結(jié)論
（1）對Hibert矩陣問題，其條件數(shù)會隨著維度的增加迅速增加，病態(tài)性會越來越明顯；在維度較低的時(shí)候，Gauss消去法、GS迭代法和SOR迭代法均可使用，且可以優(yōu)先使用Gauss消去法；如果需要求解較高維度的Hibert矩陣問題，只有SOR迭代法能夠求解。
（2）SOR方法比較適合于求解病態(tài)問題，特別是矩陣維數(shù)較高的時(shí)候，其優(yōu)點(diǎn)更為明顯。從本次實(shí)驗(yàn)可以看出，隨著矩陣維數(shù)的增大，SOR方法所需的迭代次數(shù)減少，而且誤差基本穩(wěn)定，是解決病態(tài)問題的適宜方法。
附錄：程序代碼
clear
all
clc;
format
long;
%矩陣賦值
n=input(＇矩陣H的階數(shù):n=＇);
for
i=1:n
for
j=1:n
H(i,j)=1/(i+j-1);
end
end
b=H*ones(n,1);
disp(＇H矩陣為：＇);
H
disp(＇向量b：＇);
b
%方法選擇
disp(＇選取求解方式＇);
disp(＇1
Gauss消去法，2
J迭代法，3
GS迭代法，4
SOR迭代法＇);
a=input(＇求解方式序號：＇);
%Gauss消去法
if
a==1;
H1=H;b1=b;
for
k=1:n
if
H1(k,k)==0
disp(＇主元為零，Gauss消去法無法進(jìn)行＇);
break
end
fprintf(＇第%d次消元所選取的主元是：%g＼n＇,k,H1(k,k))
for
p=k+1:n
m5=-H1(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
end
x1(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(＇解與精確解之差的無窮范數(shù)＇);
norm((x1-a),inf)
end
D=diag(diag(H));
L=-tril(H,-1);
U=-triu(H,1);
%J迭代法
if
a==2;
%給定初始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矩陣譜半徑為：%g＼n＇,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(＇解與精確解之差的無窮范數(shù)＇);
norm((xj(:,m2+1)-diag(ones(n))),inf)
disp(＇J迭代法迭代次數(shù):＇);
m2
else
disp(＇由于B矩陣譜半徑大于1，因而J法不收斂＇);
end
end
%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矩陣譜半徑為：%g＼n＇,vrho(G));
if
vrho(G)<1
for
m3=1:5000
xG(:,m3+1)=G*xG(:,m3)+fG;
if
norm((xG(:,m3+1)-xG(:,m3)),inf)<0.0001
break;
end
end
disp(＇GS迭代法計(jì)算結(jié)果:＇);
xG(:,m3+1)
disp(＇解與精確解之差的無窮范數(shù)＇);
norm((xG(:,m3+1)-diag(ones(n))),inf)
disp(＇GS迭代法迭代次數(shù):＇);
m3
else
disp(＇由于G矩陣譜半徑大于1，因而GS法不收斂＇);
end
end
%SOR迭代法
if
a==4;
%給定初始x0
ini=input(＇初始值設(shè)定:x0=＇);
x0(:,1)=ini*diag(ones(n));
disp(＇初始解向量為:＇);
x0
A=H;
for
i=1:n
b(i)=sum(A(i,:));
end
x_star=ones(n,1);
format
long
w=input(＇松弛因子:w=＇);
Lw=inv(D-w*L)*((1-w)*D+w*U);
f=w*inv(D-w*L)*b;
disp(＇迭代矩陣的譜半徑：＇)
p=vrho(Lw)
time_max=100;%迭代次數(shù)
x=zeros(n,1);%迭代初值
for
i=1:time_max
x=Lw*x+f;
end
disp(＇SOR迭代法得到的解為＇);
x
disp(＇解與精確解之差的無窮范數(shù)＇);
norm((x_star-x),inf)
end
pause
三、實(shí)驗(yàn)4.1
題目：
對牛頓法和擬牛頓法。進(jìn)行非線性方程組的數(shù)值求解
（1）用上述兩種方法，分別計(jì)算下面的兩個(gè)例子。在達(dá)到精度相同的前提下，比較其迭代次數(shù)、CPU時(shí)間等。
（2）取其他初值，結(jié)果又如何？反復(fù)選取不同的初值，比較其結(jié)果。
（3）總結(jié)歸納你的實(shí)驗(yàn)結(jié)果，試說明各種方法適用的問題。
1.算法設(shè)計(jì)
對需要求解的非線性方程組而言，牛頓法和擬牛頓法的迭代公式如下，（1）牛頓法：
牛頓法為單步迭代法，需要取一個(gè)初值。
（2）擬牛頓法：（Broyden秩1法）
其中，擬牛頓法不需要求解的導(dǎo)數(shù)，因此節(jié)省了大量的運(yùn)算時(shí)間，但需要給定矩陣的初值，取為。
2.實(shí)驗(yàn)過程
一、輸入初值；
二、根據(jù)誤差要求，按公式進(jìn)行迭代計(jì)算；
三、輸出數(shù)據(jù)；
3.計(jì)算結(jié)果及分析
（1）首先求解方程組（1），在這里，設(shè)定精度要求為，方法
牛頓法
擬牛頓法
初始值
計(jì)算結(jié)果X
x1
0.***
0.***
x2
1.***
1.0852***
x3
0.***
0.***
迭代次數(shù)
CPU計(jì)算時(shí)間/s
3.777815
2.739349
可以看出，在初始值相同情況下，牛頓法和擬牛頓法在達(dá)到同樣計(jì)算精度情況下得到的結(jié)果基本相同，但牛頓法的迭代次數(shù)明顯要少一些，但是，由于每次迭代都需要求解矩陣的逆，所以牛頓法每次迭代的CPU計(jì)算時(shí)間更長。
之后求解方程組（2），同樣設(shè)定精度要求為
方法
牛頓法
擬牛頓法
初始值
計(jì)算結(jié)果X
x1
0.***
0.***
x2
0.***
0.***
x3
-0.***
-0.***
迭代次數(shù)
CPU計(jì)算時(shí)間/s
2.722437
3.920195
同樣地，可以看出，在初始值相同情況下，牛頓法和擬牛頓法在達(dá)到同樣計(jì)算精度情況下得到的結(jié)果是基本相同的，但牛頓法的迭代次數(shù)明顯要少，但同樣的，由于每次迭代中有求解矩陣的逆的運(yùn)算，牛頓法每次迭代的CPU計(jì)算時(shí)間較長。
（2）對方程組（1），取其他初值，計(jì)算結(jié)果列表如下，同樣設(shè)定精度要求為
初始值
方法
牛頓法
擬牛頓法
計(jì)算結(jié)果
0.***
1.***
0.***
9.21***94
-5.***
18.1***205
迭代次數(shù)
CPU計(jì)算時(shí)間/s
3.907164
4.818019
計(jì)算結(jié)果
0.***
1.***
0.***
9.21***91
-5.***
18.1***807
迭代次數(shù)
2735
CPU計(jì)算時(shí)間/s
8.127286
5.626023
計(jì)算結(jié)果
0.***
1.***
0.***
0.***
1.0852***
0.***
迭代次數(shù)
CPU計(jì)算時(shí)間/s
3.777815
2.739349
計(jì)算結(jié)果
0.***
1.***
0.***
0.***
1.***
0.***
迭代次數(shù)
188
CPU計(jì)算時(shí)間/s
3.835697
2.879070
計(jì)算結(jié)果
9.21***22
-5.***
18.1***605
Matlab警告矩陣接近奇異值，程序進(jìn)入長期循環(huán)計(jì)算中
迭代次數(shù)
--
CPU計(jì)算時(shí)間/s
4.033868
--
計(jì)算結(jié)果
0.***
1.***
0.***
Matlab警告矩陣接近奇異值，程序進(jìn)入長期循環(huán)計(jì)算中
迭代次數(shù)
--
CPU計(jì)算時(shí)間/s
12.243263
--
從上表可以發(fā)現(xiàn)，方程組（1）存在另一個(gè)在(9.2,-5.6,18.1)T附近的不動點(diǎn)，初值的選取會直接影響到牛頓法和擬牛頓法最后的收斂點(diǎn)。
總的來說，設(shè)定的初值離不動點(diǎn)越遠(yuǎn)，需要的迭代次數(shù)越多，因而初始值的選取非常重要，合適的初值可以更快地收斂，如果初始值偏離精確解較遠(yuǎn)，會出現(xiàn)迭代次數(shù)增加直至無法收斂的情況；
由于擬牛頓法是一種近似方法，擬牛頓法需要的的迭代次數(shù)明顯更多，而且收斂情況不如牛頓法好（初值不夠接近時(shí)，甚至?xí)霈F(xiàn)奇異矩陣的情況），但由于牛頓法的求解比較復(fù)雜，計(jì)算時(shí)間較長；
同樣的，對方程組（2），取其他初值，計(jì)算結(jié)果列表如下，同樣設(shè)定精度要求為
初始值
方法
牛頓法
擬牛頓法
計(jì)算結(jié)果
0.***
0.***
-0.***
0.***
0.***
-0.***
迭代次數(shù)
CPU計(jì)算時(shí)間/s
2.722437
3.920195
計(jì)算結(jié)果
0.***
0.***
-0.***
0.***
-0.***
76.***
迭代次數(shù)
CPU計(jì)算時(shí)間/s
5.047111
5.619752
計(jì)算結(jié)果
0.***
0.***
-0.***
1.0e+02
*
-0.***
-0.000***6
1.754***3
迭代次數(shù)
CPU計(jì)算時(shí)間/s
3.540668
3.387829
計(jì)算結(jié)果
0.***
0.***
-0.***
1.0e+04
*
0.***
-0.***
1.***
迭代次數(shù)
CPU計(jì)算時(shí)間/s
2.200571
2.640901
計(jì)算結(jié)果
0.***
0.***
-0.***
矩陣為奇異值，無法輸出準(zhǔn)確結(jié)果
迭代次數(shù)
--
CPU計(jì)算時(shí)間/s
1.719072
--
計(jì)算結(jié)果
0.***
0.***
-0.***
矩陣為奇異值，無法輸出準(zhǔn)確結(jié)果
迭代次數(shù)
149
--
CPU計(jì)算時(shí)間/s
2.797116
--
計(jì)算結(jié)果
矩陣為奇異值，無法輸出準(zhǔn)確結(jié)果
矩陣為奇異值，無法輸出準(zhǔn)確結(jié)果
迭代次數(shù)
--
--
CPU計(jì)算時(shí)間/s
--
--
在這里，與前文類似的發(fā)現(xiàn)不再贅述。
從這里看出，牛頓法可以在更大的區(qū)間上實(shí)現(xiàn)壓縮映射原理，可以在更大的范圍上選取初值并最終收斂到精確解附近；
在初始值較接近于不動點(diǎn)時(shí)，牛頓法和擬牛頓法計(jì)算所得到的結(jié)果是基本相同的，雖然迭代次數(shù)有所差別，但計(jì)算總的所需時(shí)間相近。
（3）
牛頓法在迭代過程中用到了矩陣的求逆，其迭代收斂的充分條件是迭代滿足區(qū)間上的映內(nèi)性，對于矩陣的求逆過程比較簡單，所以在較大區(qū)間內(nèi)滿足映內(nèi)性的問題適合應(yīng)用牛頓法進(jìn)行計(jì)算。一般而言，對于函數(shù)單調(diào)或者具有單值特性的函數(shù)適合應(yīng)用牛頓法，其對初始值敏感程度較低，算法具有很好的收斂性。
另外，需要說明的是，每次計(jì)算給出的CPU時(shí)間與計(jì)算機(jī)當(dāng)時(shí)的運(yùn)行狀態(tài)有關(guān)，同時(shí)，不同代碼的運(yùn)行時(shí)間也不一定一致，所以這個(gè)數(shù)據(jù)并不具有很大的參考價(jià)值。
4.實(shí)驗(yàn)結(jié)論
對牛頓法和擬牛頓法，都存在初始值越接近精確解，所需的迭代次數(shù)越小的現(xiàn)象；
在應(yīng)用上，牛頓法和擬牛頓法各有優(yōu)勢。就迭代次數(shù)來說，牛頓法由于更加精確，所需的迭代次數(shù)更少；但就單次迭代來說，牛頓法由于計(jì)算步驟更多，且計(jì)算更加復(fù)雜，因而每次迭代所需的時(shí)間更長，而擬牛頓法由于采用了簡化的近似公式，其每次迭代更加迅速。當(dāng)非線性方程組求逆過程比較簡單時(shí)，如方程組1的情況時(shí)，擬牛頓法不具有明顯的優(yōu)勢；而當(dāng)非線性方程組求逆過程比較復(fù)雜時(shí)，如方程組2的情況，擬牛頓法就可以體現(xiàn)出優(yōu)勢，雖然循環(huán)次數(shù)有所增加，但是CPU耗時(shí)反而更少。
另外，就方程組壓縮映射區(qū)間來說，一般而言，對于在區(qū)間內(nèi)函數(shù)呈現(xiàn)單調(diào)或者具有單值特性的函數(shù)適合應(yīng)用牛頓法，其對初始值敏感程度較低，使算法具有很好的收斂性；而擬牛頓法由于不需要在迭代過程中對矩陣求逆，而是利用差商替代了對矩陣的求導(dǎo)，所以即使初始誤差較大時(shí)，其倒數(shù)矩陣與差商偏差也較小，所以對初始值的敏感程度較小。
附錄：程序代碼
%方程1，牛頓法
tic;
format
long;
%%初值
disp(＇請輸入初值＇);
a=input(＇第1個(gè)分量為：＇);
b=input(＇第2個(gè)分量為：＇);
c=input(＇第3個(gè)分量為：＇);
disp(＇所選定初值為＇);
x=[a;b;c]
%%誤差要求
E=0.0001;
%%迭代
i=0;
e=2*E;
while
e>E
F=[12*x(1)-x(2)^2-4*x(3)-7;x(1)^2+10*x(2)-x(3)-11;x(2)^3+10*x(3)-8];
f=[12,-2*x(2),-4;2*x(1),10,-1;0,3*x(2)^2,10];
det_x=((f)^(-1))*(-F);
x=x+det_x;
e=max(norm(det_x));
i=i+1;
end
disp(＇迭代次數(shù)＇);
i
disp(＇迭代次數(shù)＇);
x
toc;
%方程1，擬牛頓法
tic;
format
long;
%%初值
%%初值
disp(＇請輸入初值＇);
a=input(＇第1個(gè)分量為：＇);
b=input(＇第2個(gè)分量為：＇);
c=input(＇第3個(gè)分量為：＇);
disp(＇所選定初值為＇);
x0=[a;b;c]
%%誤差要求
E=0.0001;
%%迭代
i=0;
e=2*E;
A0=eye(3);
while
e>E
F0=[12*x0(1)-x0(2)^2-4*x0(3)-7;x0(1)^2+10*x0(2)-x0(3)-11;x0(2)^3+10*x0(3)-8];
x1=x0-A0^(-1)*F0;
s=x1-x0;
F1=[12*x1(1)-x1(2)^2-4*x1(3)-7;x1(1)^2+10*x1(2)-x1(3)-11;x1(2)^3+10*x1(3)-8];
y=F1-F0;
A1=A0+(y-A0*s)*s＇/(s＇*s);
x0=x1;
A0=A1;
e=max(norm(s));
i=i+1;
end
disp(＇迭代次數(shù)＇);
i
disp(＇迭代次數(shù)＇);
x0
toc;
%方程2，牛頓法
tic;
format
long;
%%初值
disp(＇請輸入初值＇);
a=input(＇第1個(gè)分量為：＇);
b=input(＇第2個(gè)分量為：＇);
c=input(＇第3個(gè)分量為：＇);
disp(＇所選定初值為＇);
x=[a;b;c]
%%誤差要求
E=0.0001;
%%迭代
i=0;
e=2*E;
while
e>E
F=[3*x(1)-cos(x(2)*x(3))-0.5;x(1)^2-81*(x(2)+0.1)^2+sin(x(3))+1.06;exp(1)^(-x(1)*x(2))+20*x(3)+(10*pi-3)/3];
f=[3,x(3)*sin(x(2)*x(3)),x(2)*sin(x(2)*x(3));2*x(1),-162*x(2)-81/5,cos(x(3));-x(2)*exp(1)^(-x(1)*x(2)),-x(1)*exp(1)^(-x(1)*x(2)),20];
det_x=((f)^(-1))*(-F);
x=x+det_x;
e=max(norm(det_x));
i=i+1;
end
disp(＇迭代次數(shù)＇);
i
disp(＇迭代次數(shù)＇);
x
toc;
%方程2，擬牛頓法
tic;
format
long;
%%初值
%%初值
disp(＇請輸入初值＇);
a=input(＇第1個(gè)分量為：＇);
b=input(＇第2個(gè)分量為：＇);
c=input(＇第3個(gè)分量為：＇);
disp(＇所選定初值為＇);
x0=[a;b;c]
%%誤差要求
E=0.0001;
%%迭代
i=0;
e=2*E;
A0=eye(3);
while
e>E
F0=[3*x0(1)-cos(x0(2)*x0(3))-0.5;x0(1)^2-81*(x0(2)+0.1)^2+sin(x0(3))+1.06;exp(1)^(-x0(1)*x0(2))+20*x0(3)+(10*pi-3)/3];
x1=x0-A0^(-1)*F0;
s=x1-x0;
F1=[3*x1(1)-cos(x1(2)*x1(3))-0.5;x1(1)^2-81*(x1(2)+0.1)^2+sin(x1(3))+1.06;exp(1)^(-x1(1)*x1(2))+20*x1(3)+(10*pi-3)/3];
y=F1-F0;
A1=A0+(y-A0*s)*s＇/(s＇*s);
x0=x1;
A0=A1;
e=max(norm(s));
i=i+1;
end
disp(＇迭代次數(shù)＇);
i
disp(＇迭代次數(shù)＇);
x0
toc;
點(diǎn)擊咨詢