數(shù)值分析希爾伯特病態(tài)線性方程組

1、病態(tài)線性代數(shù)方程組的求解理論的分析表明，求解病態(tài)的線性代數(shù)方程組是困難的?？紤]方程組Hx=1.一b的求解，其中H為Hilbert矩陣，H=(hj)nM，hj=,i,j=1,2,.,nij-11 .估計Hilbert矩陣2-條件數(shù)與階數(shù)的關(guān)系；2 .選擇問題的不同維數(shù)，分別用Gaus醉I去法，Jacob迭代，GS迭代和SOR迭代求解，比較結(jié)果；3 .討論病態(tài)問題求解的算法。解：1、取Hilbert矩陣階數(shù)最高分別為n=20和n=100。采用Hilbert矩陣的2-條件數(shù)作為實驗的比較對象，畫出的曲線如下圖所示：lg(cond(Hn)nlg(cond(Hn)n關(guān)系圖20204,nHdnorTMWo

2、818661210n從圖中可以看出，在n013之刖，圖像接近直線，在n>13之后，圖像趨于平緩，在一定的范圍內(nèi)上下波動。為了比較圖像的線性部分，作出一條直線與已知曲線進行比較。比較直線的關(guān)系式為：lg(cond(Hn)=1.519n-1.833,結(jié)果下圖所示lg(cond(Hn)-n關(guān)系圖30.,.IIIIIIIBI25一_20._H15.一g10-5-0.-51111111102468101214161820n從圖2中可以看出，當n較小時，lg(cond(Hn)n之間近似滿足線性關(guān)系。當n繼續(xù)增大到100時，lg(cond(Hn)n關(guān)系圖下圖所示：lg(cond(Hn)-n關(guān)系圖25,

3、11,11,20-.二d-.RJi”n15.-Hog10-5一r01111I11110102030405060708090100n從圖中可以看出，圖像的走勢符合在n=20時的猜想，在n大于一定的值之后，圖像趨于平緩，且在一定范圍內(nèi)震蕩，同時又有一定上升趨勢，但上升速度很慢。2、選擇不同的階數(shù)n,設(shè)方程組的精確解為xz=(1,1；1)T進行計算，用四種方法解x_Guass1x_Jacobi1x_GS1、x_SOR1寸比表如下表所小。四種方法解x_Guassx_Jacobi、x_GSx_SOR對比表nx_Guassx_Jacobix_GSx_SOR31-8.22E+3070.9999911-1.5

4、52e+30810.999991-Inf0.999971411.0172e+30810.999991Inf0.999651.00011Inf1.00080.999741Inf0.999481.000251-6.07E+3070.998791.00011-1.2685e+3081.02230.999041-1.6592e+3080.905111.0041-Inf1.14210.994111-Inf0.930951.002861-6.64E+3070.9994911-1.4322e+3081.00590.999541-Inf0.995181.00111-Inf0.951491.00091-Inf1

5、.10210.996211-Inf0.945361.002371-9.20E+3071.000711-Inf0.981921.00021-Inf1.10220.99791-Inf0.808241.00621-Inf1.0620.994641-Inf1.14480.998981-Inf0.899561.0021811.4378e+3081.00190.999971Inf0.958751.0011Inf1.19620.994081Inf0.730241.01141Inf0.954490.994641Inf1.18240.998061Inf1.14660.996681Inf0.828041.0042

6、915.09E+3071.00310.9999211.1719e+3080.939461.001911.6249e+3081.24890.990291Inf0.754461.01540.99999Inf0.838240.995491Inf1.09850.998760.99998Inf1.22470.994051Inf1.11480.999691Inf0.775611.0046101-7.82E+3071.00360.999941-Inf0.93761.00111-Inf1.21420.995551-Inf0.889821.00450.99992-Inf0.779991.00111.0002-I

7、nf0.958251.00030.99964-Inf1.16090.997481.0003-Inf1.22070.997760.99982-Inf1.08170.999271-Inf0.750191.0031117.97E+3071.00290.999951Inf0.958571.00070.99998Inf1.09130.998091.0002Inf1.09090.99980.99894Inf0.802911.0031.0037Inf0.828491.0010.99187Inf1.02360.999361.0111Inf1.18860.997760.99081Inf1.21280.99799

8、1.0042Inf1.06010.999610.99916Inf0.736721.00281218.11E+3071.00120.999971Inf0.995921.00030.99993Inf0.919041.00021.0009Inf1.29740.996210.99316Inf0.887531.0041.0298Inf0.743091.00140.91765Inf0.874511.0011.1481Inf1.08370.998340.82734Inf1.22150.997821.1258Inf1.21480.998030.9479Inf1.04310.999931.0094Inf0.71

9、5151.00291314.90E+3070.998980.999990.999971.1863e+3081.03880.99991.00121.6981e+3080.74141.00220.98043Inf1.46940.993251.1771Inf1.0031.00450.033234Inf0.706161.00154.3909Inf0.746491.0023-6.8988Inf0.953940.9991313.35Inf1.16010.99821-11.81Inf1.26350.997399.4534Inf1.21940.99831-2.2127Inf1.01981.00021.5352

10、Inf0.675951.00321415.94E+3070.9967710.999991.4525e+3081.07990.999511.0004Inf0.58371.00410.99479Inf1.59150.990681.0373Inf1.12841.00460.85719Inf0.711571.00151.2465Inf0.650911.00351.176Inf0.824241-0.85802Inf1.06130.99895.3287Inf1.24210.99733-4.4541Inf1.3020.997485.0263Inf1.21660.99841-0.64129Inf0.98559

11、1.00061.2863Inf0.621961.0035211Inf0.9894910.99978Inf1.16240.999291.0084Inf0.559161.00170.86586Inf0.938321.00142.1174Inf1.53570.9961-4.2693Inf1.44540.9985715.159Inf1.05070.99966-18.28Inf0.698661.00176.2524Inf0.526541.002113.831Inf0.538331.002110.008Inf0.681851.0014-49.838Inf0.892471.000445.871Inf1.11

12、260.99942-27.1Inf1.29840.9986154.555Inf1.41970.99808-34.976Inf1.45850.99793-31.788Inf1.40570.9981824.491Inf1.25930.9988541.386Inf1.02150.99993-43.821Inf0.697761.001413.526Inf0.295371.0032Gauss消去法求解：選擇問題白階數(shù)為38時，用Gauss消去法求得的解與精確解一致，當階數(shù)為914時，解開始出現(xiàn)偏差，而且n越大，偏差越大。用迭代法求解：取初始向量為1.2（1,1；1）t.無論n為多少階，用Jacobi迭代

13、方法迭代出現(xiàn)發(fā)散的不穩(wěn)定現(xiàn)象，無法求解；用GS迭代方法迭代不發(fā)散，能求得解，但收斂非常緩慢，當?shù)螖?shù)取得相當大（20000次）時解仍在精確解附近波動；取w=1.5,用SORg代方法迭代不發(fā)散，能求得解，收斂速度較GS迭代快一些，但仍非常緩慢。選擇問題白階數(shù)為21時：用Gauss消去法求得的解與精確解相差很大，相差10的量級。用迭代法求解時，取初始向量為1.2（1,1；1）t用Jacobi迭代方法迭代發(fā)散，無法求解；用GS迭代方法迭代不發(fā)散，能求得解，但收斂非常緩慢，迭代20000次后，算得的值與精確值1相差很大。用SOR4代方法迭代不發(fā)散，能求得解，收斂速度還行。從上面的結(jié)果可以看出當病態(tài)問

14、題的階數(shù)升高時作為直接法的Gauss消去法又能求解變成不能求解。用Jacobi迭代方法迭代發(fā)散，無法求解；而GS和SOR迭代法在階數(shù)升高時仍能求解，選擇合適的w,可使SOR迭代法收斂速度加快。但在階數(shù)較低時直接法能求得精確解而迭代法卻總存在一定的誤差?？梢娭苯臃ㄅc迭代法各有各的優(yōu)勢與不足。3、關(guān)于病態(tài)問題的求解，主要的方法是對原方程作某些預(yù)處理，以降低系數(shù)矩陣的條件數(shù)。例如選擇非奇異矩陣P,QR,使方程組Ax=b化為等價方程組(PAQ)y=Pb,原方程的解x=Qy.原則上應(yīng)使矩陣PAQ的條件數(shù)比A有所改善。一般P和Q可選擇為三角形矩陣或?qū)蔷仃?。理論上最好選擇對角陣P和Q,滿足：cond(PA

15、Q)=mincond(PAQ)。以Hilbert矩陣為例。lg(cond(hn)/cond(Hn)n關(guān)系圖2XJnwnuanocrnb001I5012503003504005045011上圖中的hn=DnHnDn(Dn為Hn的對角線兀素開萬構(gòu)成的對角矩陣)，n為希爾伯特矩陣的階數(shù)。我們觀察到，隨希爾伯特矩陣階數(shù)的增大，函數(shù)lgcond(，)/cond(Hn)在卜3,1.5區(qū)間波動，主要集中在-1.5,0.5區(qū)間。我們知道當cond(hn)Wcond(Hn)時，有l(wèi)gcond(hn)/cond(Hn)W0,在上圖中，我們可以很容易的觀察到，對于大部分n,函數(shù)值都是小于或者等于零的，這說明Hn經(jīng)過

16、預(yù)處理后的hn的條件數(shù)較小，在一定程度上改善了原希爾伯特矩陣的性質(zhì)。由于希爾伯特矩陣病態(tài)的性質(zhì)，對于對原希爾伯特矩陣進行預(yù)處理后的新希爾伯特矩陣的條件數(shù)在一定的區(qū)間呈波動的變化規(guī)律。從整體上觀察，對于大多數(shù)的n,進行預(yù)處理后的希爾伯特矩陣的條件數(shù)有明顯的降低，這就說明預(yù)處理改善了大多數(shù)階數(shù)的希爾伯特矩陣的性質(zhì)。t1=20i1=1:t1forn=1:t1H=hilb(n)K1(n)=cond(H)K(n)=log10(K1(n)endl1=1.519*i1-1.833plot(i1,K(i1),i1,l1)xlabel('n','fontsize',20)ylab

17、el('lg(cond(Hn)','fontsize',20)title('lg(cond(Hn)n關(guān)系圖','fontsize',24)set(gca,'fontsize',16)figuret2=100i2=1:t2forn=1:t2H=hilb(n)K1(n)=cond(H)K(n)=log10(K1(n)endplot(i2,K(i2)xlabel('n','fontsize',20)ylabel('lg(cond(Hn)','fontsize'

18、;,20)title('lg(cond(Hn)n關(guān)系圖','fontsize',24)set(gca,'fontsize',16)clearform=3:25n=m%Hilbert矩陣的階數(shù)n1=n+1H=hilb(n)K=cond(H)u=1.2x0=u*ones(n,1)xz=ones(n,1)b=H*xza=HD=diag(diag(a)U=-triu(a,1)L=-tril(a,-1)w=1.5d=1.0e-6max=1000c=a,b%Gaus目肖去法fork=1:nbmax=0;%選主元，并存放在變量bmax中fori=k:nifbm

19、ax-abs(c(i,k)<0bmax=abs(c(i,k);l=i;%記下絕對值最大者所在位置endendifbmax<1.0e-20'stop'pause;%如果主元絕對值小于1.0e-20,則A奇異，計算終止end%進行行交換ifl=k%如果l=k,則不需要交換forj=k:n1t=c(l,j);c(l,j)=c(k,j);c(k,j)=t;endend%進行消元t=1/c(k,k);c(k,j)=c(k,j)*t;c(k+1:n,j)=c(k+1:n,j)-c(k+1:n,k)*c(k,j);forj=k+1:n1endendC的第n+1歹U中%回代法求解上

20、三角形方程組，解存放在矩陣fork=2:ni=n1-k;forj=i+1:nc(i,n1)=c(i,n1)-c(i,j)*c(j,n1);endendxGuass=c(1:n,n1)%Jacobi迭代法BJ=D(L+U)fJ=Dbx2=x0xJacobi=BJ*x2+fJkJacobi=1while(norm(xJacobi-x2)>=d)&(kJacobi<max)x2=xJacobixJacobi=BJ*x2+fJkJacobi=kJacobi+1end%GS迭代法BG=(D-L)UfG=(D-L)bx3=x0xGS=BJ*x3+fGkGS=1while(norm(xGS-x3)>=d)&(kGS<max)x3=xGSxGS=BG*x3+fGkGS=