1、數值計算第一次大作業實驗目的 以Hilbert矩陣為例，研究處理病態問題可能遇到的困難。內容 Hilbert矩陣的定義是它是一個對稱正定矩陣，而且隨著n的增加迅速增加，其逆矩陣，這里1) 畫出之間的曲線（可以用任何的一種范數）。你能猜出之間有何種關系嗎？提出你的猜想并想法驗證。用行范數for n=1:50for i=1:nfor j=1:nA(i,j)=1/(i+j-1);B(i,j)=factorial(n+i-1)*factorial(n+j-1)/(i+j-1)*(factorial(i-1)*factorial(j-1)2*factorial(n-i)*factorial(n-j);e

2、ndendresult1=0;for j=1:nresult1=result1+A(1,j);endresult1=log(result1);result2=0;for i=1:nfor j=1:nresult2=B(i,j)+result2;endresult(i)=log(result2);endm=max(result);x(n)=result1+m;endplot(1:50,x)對于更大的n值，由于Hilbert逆矩陣中的元素過大，溢出，故在此取50以內的n。圖1 關系曲線圖猜想之間存在線性關系驗證：設在以上程序基礎上，再添加>>>> y=x'>&

3、gt; l=1:40;>> k=l'>> p=polyfit(k,y,1)%一次多項式擬合p =3.5446 -3.0931% P=polyfit(k,y,2) %二次多項式擬合p = -0.0008 3.5778 -3.3253% P=polyfit(k,y,3) %三次多項式擬合 0.0000 -0.0033 3.6198 -3.4777% P=polyfit(k,y,4) %四次多項式擬合-0.0000 0.0002 -0.0082 3.6654 -3.5815% P=polyfit(k,y,5) %五次多項式擬合p = 0.0000 -0.0000 0.

4、0007 -0.0156 3.7107 -3.6542從上式可以看出，高次項系數相對于一次項和常數項系數要小很多，所以取2）設是的對角線元素開方構成的矩陣。，不難看出依然是對稱矩陣，而且對角線元素都是1。把變成的技術稱為預處理。畫出之間的曲線（可以用任何一種范數）。你能對于預處理得出什么印象？本小題用2范數>> clear>> n=500;>> c=;>> for k=2:nH=hilb(k);D=diag(sqrt(diag(H);D1=inv(D);H1=D1*H*D1;c=c,cond(H1)/cond(H);endC=log(c);k=2

5、:n;plot(k,C,'r-')圖 2 隨n的變化曲線圖從圖中給出了函數的變化曲線。我們觀察到隨著Hilbert矩陣階數的增大，函數值在-6,4區間波動，主要集中在-3,1區間。我們知道在時，有，在上圖中，我們可以容易觀察到，對于大部分，函數值都是小于或者等于零的，這說明經過預處理后的地條件數較小，由于條件數愈大，方程組的病態愈嚴重，也就愈難得到方程組比較準確地解，所以預處理在一定程度上改善了原Hilbert矩陣的特性。3）對于，給定不同的右端項。分別用以及，求解，比較計算結果。取n=4，b=1;2;3;4>> b=1,2,3,4'>> H=h

6、ilb(4);>> H1=inv(H);>> x1=H1*bx1 = 1.0e+003 * -0.0640 0.9000 -2.5200 1.8200>> D=diag(sqrt(diag(H);>> D1=inv(D);>> H2=D1*H*D1;>> H3=inv(H2);>> x2=D1*H3*D1*bx2= 1.0e+003 *-0.0640 0.9000 -2.5200 1.8200x3=Hbx3 = 1.0e+003 * -0.0640 0.9000 -2.5200 1.8200同理，當n=5時，b

7、=1;2;3;4;5x1 = 1.0e+004 * 0.0125 -0.2880 1.4490 -2.4640 1.3230x2 = 1.0e+004 * 0.0125 -0.2880 1.4490 -2.4640 1.3230x3 = 1.0e+004 * 0.0125 -0.2880 1.4490 -2.4640 1.3230當n=6時，b=1;2;3;4;5;6x1 = 1.0e+005 * -0.0022 0.0735 -0.5712 1.6632 -2.0160 0.8593x2 = 1.0e+005 * -0.0022 0.0735 -0.5712 1.6632 -2.0160 0

8、.8593x3 = 1.0e+005 * -0.0022 0.0735 -0.5712 1.6632 -2.0160 0.8593當n=7時，b=1;2;3;4;5;6;7x1 = 1.0e+006 * 0.0003 -0.0161 0.1777 -0.7728 1.5593 -1.4636 0.5165x2 = 1.0e+006 * 0.0003 -0.0161 0.1777 -0.7728 1.5593 -1.4636 0.5165x3 = 1.0e+006 * 0.0003 -0.0161 0.1777 -0.7728 1.5593 -1.4636 0.5165當n=8時，b=1;2;3

9、;4;5;6;7;8x1 = 1.0e+007 * -0.0001 0.0032 -0.0469 0.2818 -0.8316 1.2757 -0.9754 0.2934x2 = 1.0e+007 * -0.0001 0.0032 -0.0469 0.2818 -0.8316 1.2757 -0.9754 0.2934x3 = 1.0e+007 * -0.0001 0.0032 -0.0469 0.2818 -0.8316 1.2757 -0.9754 0.2934當n=9時，b=1;2;3;4;5;6;7;8;9x1 = 1.0e+007 * 0.0001 -0.0058 0.1095 -0

10、.8649 3.4685 -7.6684 9.4594 -6.0952 1.5972x2 = 1.0e+007 * 0.0001 -0.0058 0.1095 -0.8649 3.4685 -7.6684 9.4594 -6.0952 1.5972x3 = 1.0e+007 * 0.0001 -0.0058 0.1095 -0.8649 3.4685 -7.6685 9.4595 -6.0952 1.5972當n=10時，b=1;2;3;4;5;6;7;8;9;10x1 = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -1.2107 3.5942 -6.3

11、224 6.5105 -3.6228 0.8406x2 = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -1.2107 3.5943 -6.3227 6.5108 -3.6229 0.8406x3 = 1.0e+008 * -0.0000 0.0010 -0.0233 0.2330 -1.2108 3.5947 -6.3233 6.5114 -3.6232 0.8407當n=11時，b=1;2;3;4;5;6;7;8;9;10;11x1 = 1.0e+009 * 0.0000 -0.0002 0.0046 -0.0564 0.3674 -1.3991 3.

12、2758 -4.7725 4.2140 -2.0629 0.4294x2 = 1.0e+009 * 0.0000 -0.0002 0.0046 -0.0567 0.3687 -1.4038 3.2861 -4.7867 4.2259 -2.0685 0.4305x3 = 1.0e+009 * 0.0000 -0.0002 0.0046 -0.0567 0.3687 -1.4037 3.2858 -4.7862 4.2255 -2.0682 0.4305當n=12時，b=1;2;3;4;5;6;7;8;9;10;11;12Warning: Matrix is close to singular

13、or badly scaled. Results may be inaccurate. RCOND = 2.692153e-017.x1 = 1.0e+010 * -0.0000 0.0000 -0.0007 0.0110 -0.0886 0.4225 -1.2686 2.4585 -3.0691 2.3821 -1.0452 0.1980x2 = 1.0e+010 * -0.0000 0.0000 -0.0009 0.0133 -0.1055 0.4975 -1.4788 2.8408 -3.5190 2.7127 -1.1831 0.2229x3 = 1.0e+010 * -0.0000

14、0.0000 -0.0008 0.0123 -0.0978 0.4630 -1.3817 2.6636 -3.3099 2.5586 -1.1187 0.2113由此可見，當n較小時，三種方法得出的結果基本相同，隨著n的增大，三種方法得出的結果的偏差也越來越大。4）取不同的并以的第一列為右端向量，用高斯-塞德爾迭代法求解，觀察其收斂性。最后你能對于有關Hilbert矩陣的計算得出哪些結論。輸入：A為方程組的系數矩陣，b為方程組右端的列向量，X為迭代初值構成的列向量，max為最大迭代次數，w為誤差精度輸出：x為求得的方程組的解構成的列向量，n為迭代次數取n=4，則>> A=hilb(

15、4);>> b=A(:,1)；>> X=0;0;0;0;>> max=50;>> w=10-6;>> gaussseidel(A,b,X,max,w)迭代次數為n = 2方程組的解為x = 1.0000 -0.0000 0.0000 -0.0000取n=8，則>> X=0;0;0;0;0;0;0;0;>> max=50;>> w=10-6;>> gaussseidel(A,b,X,max,w)迭代次數為n = 2方程組的解為x = 1.0000 -0.0000 0.0000 -0.000

16、0 0.0000 -0.0000 -0.0000 -0.0000>> A=hilb(1000);>> b=A(:,1);>> X=zeros(1000,1);>> max=50;>> w=10-6;>> gaussseidel(A,b,X,max,w)迭代次數為n = 2迭代結果為第一行為1，其余為0的向量。取元素全為一的向量作為起始迭代向量時，有：>> A=hilb(4);>> b=A(:,1);>> X=ones(4,1);>> max=50;>> w=10-

17、6;>> gaussseidel(A,b,X,max,w)在最大迭代次數內不收斂!最大迭代次數后的結果為x = 1.0206 -0.0457 -0.0884 0.1310>> A=hilb(8);>> b=A(:,1);>> X=ones(8,1);>> max=50;>> w=10-6;>> gaussseidel(A,b,X,max,w)在最大迭代次數內不收斂!最大迭代次數后的結果為x = 0.9510 0.2777 -0.1542 -0.2055 -0.1283 -0.0147 0.1016 0.2091

18、>> A=hilb(200);>> b=A(:,1);>> X=ones(200,1);>> max=50;>> w=10-6;>> gaussseidel(A,b,X,max,w)在最大迭代次數內不收斂!最大迭代次數后的結果為x = 0.6054 1.3942 0.2032 -0.4150（由于元素太多故不一一列出）通過推導，不難發現，起始的迭代向量為零時，對于不同的n值，均只需迭代兩次就可以得出答案，且結果均是第一個元素為1，其余元素為0。這是因為右端向量b取的是Hilbert 矩陣的第一列。如果b取其它向量，則可以知

19、道用迭代法求解時的求解誤差比較大，不收斂。當起始的迭代向量不為零時，可以得到不同的答案，如上題中的取元素全為一的向量時，結果是發散。由此可以知道對于不同的初始條件，迭代結果也不同。所以用高斯-賽德爾迭代法解此方程組并不是對任意的初始向量和右端向量都收斂。Hilbert 矩陣是病態的，對于對原Hilbert 矩陣進行預處理后的新Hilbert 矩陣的條件數在一定范圍內呈波動的變化規律。從整體上觀察，對于大多數的n，進行預處理后的Hilbert 矩陣地條件數有明顯的降低，這就說明預處理改善了大多數Hilbert 矩陣的性質。用迭代法求解方程時，如果系數矩陣為Hilbert 矩陣，則求解結果的誤差較大，根據迭代法基本定理可知用高斯-賽德爾迭代法解系