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

......zz數(shù)值分析實(shí)驗(yàn)報(bào)告一、實(shí)驗(yàn)3.1題目：考慮線性程組AxbARnnbRn，編制一個(gè)能自動(dòng)選取主元，又能手動(dòng)選取主元的求解線性代數(shù)程組的Gauss消去過程。6 1 7 1取矩陣A b則程有解x*,T取n10 8 6 1 8 6

計(jì)算矩陣的條件數(shù)。分別用順序Gauss消元、列主元Gauss消元和完全選主元Gauss消元法求解，結(jié)果如？現(xiàn)選擇程序中手動(dòng)選取主元的功能，每步消去過程都選取模最小或按模盡按模最大的元素作為主元，結(jié)果又如？分析實(shí)驗(yàn)的結(jié)果。取矩陣階數(shù)n=20過程中的作用。實(shí)驗(yàn)，觀察記錄并分析實(shí)驗(yàn)的結(jié)果。算法介紹首先，分析各種算法消去過程的計(jì)算公式，順序高斯消去法：第k步消去中，設(shè)增廣矩陣B中的元素ak0（若等于零則可以判定系數(shù)kkk行以下各行計(jì)算lik

ak ikakkk

,ik1,k2, ,n，分別用lik

B的第k行并加到第k1,k

,n行，則可將增廣矩陣B中第k列中ak以下的元素消為零；重復(fù)此法，從第1步進(jìn)行到第n-1步，則kkBnAn,bn；

ak

akkk第k步消去中，在增廣矩陣B中的子陣

kn中，選取ak使得

ikkak akknk nnakmaxa(k)

，當(dāng)i

k時(shí)，對(duì)B中第k行與第

行交換，然后按照和順序消去ikk

kin ik k k法相同的步驟進(jìn)行。重復(fù)此法，從第1步進(jìn)行第n-1步，就可以得到最終的增廣矩陣，即Bn1An1,bn1；

ak

ak第k步消去中，在增廣矩陣

kkB中對(duì)應(yīng)的子陣

kn

ak使

ijkkak akkknk nn得ak

maxa(k)

k或j

kB中第k行與第

行、第k列與第j列kkij kin ij k k k kkkkjn交換，然后按照和順序消去法相同的步驟進(jìn)行即可。重復(fù)此法，從第1步進(jìn)行到第n-1Bn2An2,bn2； 接下來，分析回代過程求解的公式，容易看出，對(duì)上述任一種消元法，均有以下計(jì)算公式：bnx n ;n annnkbnkx k

njkankk

anxkj

;kn1,n2, ,1實(shí)驗(yàn)程序的設(shè)計(jì)一、輸入實(shí)驗(yàn)要求及初始條件；二、計(jì)算系數(shù)矩陣A的條件數(shù)及程組的理論解；三、對(duì)各不同法編程計(jì)算，并輸出最終計(jì)算結(jié)果。計(jì)算結(jié)果及分析（1）先計(jì)算系數(shù)矩陣的條件數(shù)，結(jié)果如下，cond(A)1

2557.500000,cond(A)2

1727.556025,cond(A)

2557.50000或A引起解的較大誤差；采用順序高斯消去法，計(jì)算結(jié)果為：......zz最終解為x=(1.0000,1.0000,1.0000, 1.0001,0.9998,1.0004,1.0012,0.9979,1.0028)T使用無窮數(shù)衡量誤差，得到XX* =2.0401e-14，可以發(fā)現(xiàn)，采用順序高斯消去法式并沒有對(duì)結(jié)果造成特別大的影響。若采用列主元高斯消元法，則結(jié)果為：最終解為x=(1.0000,1.0000,1.0000, 1.0000,1.0000,1.0000,1.0000,1.0000,1.0000)T同樣使用無窮數(shù)衡量誤差，有XX*=0；若使用完全主元高斯消元法，則結(jié)果為最終解x=(1.0000,1.0000,1.0000, 1.0000,1.0000,1.0000,1.0000,1.0000,1.0000)T同樣使用無窮數(shù)衡量誤差，有XX* =0；（2）若每步都選取模最小或盡可能小的元素為主元，則計(jì)算結(jié)果為最終解x=(1.0000 1.0000 1.0000 1.0001 0.9998 1.00040.9993 1.0012 0.9979 1.0028)T使用無窮數(shù)衡量誤差，有XX*為2.0401e-14；而完全主元消去法的誤差為XX*=0。從和的實(shí)驗(yàn)結(jié)果可以發(fā)現(xiàn)，列主元消去法和完全主元消去法都得4種法每步選取的主元數(shù)值，并列表進(jìn)行比較，結(jié)果如下：......zz第n次消元順序列主元完全主元模最小16.0000886.000024.6667884.666734.4286884.428644.3333884.333354.2258884.225864.6032884.603274.6063884.606384.4902884.490294.4853884.4853104.30990.54690.54694.30998，與4所以計(jì)算過程中的累積誤差也較小，最終4種法的輸出結(jié)果均較為精確。8，0認(rèn)的主元也就是該列最小的主元，因此兩種法所得到的計(jì)算結(jié)果是一致的。（元素相差不大并且維度不高，這個(gè)理論現(xiàn)象在這里并沒有充分體現(xiàn)出來。（3）選取模最小或盡可列主元高斯 完全主元高斯順序高斯消去法能小元素作為主元消去 消去消去X0.98860.98861 11.02270.9547選取模最小或盡可列主元高斯 完全主元高斯順序高斯消去法能小元素作為主元消去 消去消去X0.98860.98861 11.02270.95471.02270.95471 11.09020.82091.09020.82091 11.35240.31791.35240.31791 11.27320.81731.27320.81731 11.91021.91021 1111111111.00001.0000111.00001.0000111.00001.00011.00001.0001110.99981.0004110.99981.00040.99931.00140.99931.0014110.99721.00570.99721.005711111111112.9430e-11002.9430e-11n=10時(shí)候的誤差比相比，n=20時(shí)的誤差增長(zhǎng)了大約1000倍，這是由于計(jì)算過程中舍入誤差的（4）不同矩陣維度下的誤差如下，在這里，為便起見，選取2-條件數(shù)對(duì)不同維度的系數(shù)矩陣進(jìn)行比較。維度條件數(shù)順序消去列主元完全主元模盡量小n101.7e+32.84e-14002.84e-14n201.8e+62.91e-11002.91e-11n255.7e+79.31e-10009.31e-10n301.8e+92.98e-08002.98e-08n401.9e+123.05e-05003.05e-05n703.8e+163.28e+043.88e-123.88e-123.28e+04n1008.5e+163.52e+134.2e-34.2e-33.52e+13長(zhǎng)較快，這是由于舍入誤差逐步累計(jì)而造成的。不過，法二與法三在維度小于40結(jié)論本文矩陣中的元素差別不大，模最大和模最小的元素并沒有數(shù)量級(jí)上的差四種法都足夠精確。結(jié)果誤差均較小。算。附錄：程序代碼clearclc;formatlong;%法選擇n=input('矩陣A階數(shù):n=');disp('選取求解式');disp('1順序Gauss消元法，2列主元Gauss消元法，3完全選主元Gauss消元法，4模最小或近可能小的元素作為主元');a=input('求解式序號(hào)：');%賦值A(chǔ)和bA=zeros(n,n);b=zeros(n,1);fori=1:nA(i,i)=6;ifi>1A(i,i-1)=8;endifi<nA(i,i+1)=1;endendfori=1:nforj=1:nb(i)=b(i)+A(i,j);endenddisp(Adisp(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消元法ifa==1A1=A;b1=b;fork=1:nifA1(k,k)==0disp('主元為零，順序Gauss消元法無法進(jìn)行');breakendfprintf('第%d次消元所選取的主元：%g\n',k,A1(k,k))%disp('此次消元后系數(shù)矩陣為：');%A1forp=k+1:nl=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);endendx1(n)=b1(n)/A1(n,n);fork=n-1:-1:1forw=k+1:nb1(k)=b1(k)-A1(k,w)*x1(w);endx1(k)=b1(k)/A1(k,k);enddisp('順序Gauss消元法解為:');disp(x1);disp('所求解與精確解之差的無窮-數(shù)為');norm(x1-X,inf)end%列主元Gauss消元法ifa==2A2=A;b2=b;fork=1:n[max_i,max_j]=find(A2(:,k)==max(abs(A2(k:n,k))));ifmax_i~=kA2_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;endifA2(k,k)==0disp('主元為零，列主元Gauss消元法無法進(jìn)行');breakendfprintf('第%d次消元所選取的主元：%g\n',k,A2(k,k))%disp('此次消元后系數(shù)矩陣為：');%A2forp=k+1:nl=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);endendx2(n)=b2(n)/A2(n,n);fork=n-1:-1:1forb2(k)=b2(k)-A2(k,w)*x2(w);endx2(k)=b2(k)/A2(k,k);enddisp('列主元Gauss消元法解為:');disp(x2);disp('所求解與精確解之差的無窮-數(shù)為');norm(x2-X,inf)end%完全選主元Gauss消元法ifa==3A3=A;b3=b;fork=1:nVV=eye(n);[max_i,max_j]=find(A3(k:n,k:n)==max(max(abs(A3(k:n,k:n)))));ifnumel(max_i)==0[max_i,max_j]=find(A3(k:n,k:n)==-max(max(abs(A3(k:n,k:n)))));endW=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;ifA3(k,k)==0disp('主元為零，完全選主元Gauss消元法無法進(jìn)行');breakendfprintf('第%d次消元所選取的主元：%g\n',k,A3(k,k))%disp('此次消元后系數(shù)矩陣為：');%A3forp=k+1:nl=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);endendx3(n)=b3(n)/A3(n,n);fork=n-1:-1:1forw=k+1:nb3(k)=b3(k)-A3(k,w)*x3(w);endx3(k)=b3(k)/A3(k,k);enddisp('完全選主元Gauss消元法解為:');disp(x3);disp('所求解與精確解之差的無窮-數(shù)為');norm(x3-X,inf)end%模最小或近可能小的元素作為主元ifa==4A4=A;b4=b;fork=1:nAA=A4;AA(AA==0)=NaN;[min_i,j]=find(AA(k:n,k)==min(abs(AA(k:n,k))));ifnumel(min_i)==0[min_i,j]=find(AA(k:n,k)==-min(abs(AA(k:n,k:n))));endW=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;ifA4(k,k)==0disp('主元為零，模最小Gauss消元法無法進(jìn)行');breakendfprintf('第%d次消元所選取的主元：%g\n',k,A4(k,k))%A4forp=k+1:nl=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);endendx4(n)=b4(n)/A4(n,n);fork=n-1:-1:1forw=k+1:nb4(k)=b4(k)-A4(k,w)*x4(w);endx4(k)=b4(k)/A4(k,k);enddisp('模最小Gauss消元法解為:');disp(x4);disp('所求解與精確解之差的無窮-數(shù)為');norm(x4-X,inf)end二、實(shí)驗(yàn)3.3題目：考慮程組Hxb的解，其中系數(shù)矩陣H為Hilbert矩陣：Hh

,h

1 ,i,j1,2, ,ni,j

nn

i,j

ij1這是一個(gè)著名的病態(tài)問題。通過首先給定解（例如取為各個(gè)分量均為1）再計(jì)算出右端的辦法給出確定的問題。選擇問題的維數(shù)為6，分別用Gauss消去法（即LU分解、J迭代法、GS迭代法和SOR結(jié)論如。果說明的什么？討論病態(tài)問題求解的算法。算法設(shè)計(jì)對(duì)任意線性程組Axb，分析各種法的計(jì)算公式如下，（1）Gauss消去法：首先對(duì)系數(shù)矩陣進(jìn)行LUALULUxbUxy則原程可以分為兩步回代求解：yLyb具體法這里不再贅述。（2）J迭代法：首先分解ADLU ，再構(gòu)造迭代矩陣Xk1BXkf ，其中BD1LU,f，進(jìn)行迭代計(jì)算，直到誤差滿足要求。GS迭代法：首先分解ADLU，再構(gòu)造迭代矩陣Xk1GXkf ，其中GDL1U,fDL1b，進(jìn)行迭代計(jì)算，直到誤差滿足要求。SOR迭代法：首先分解ADLU，再構(gòu)造迭代矩陣XL

Xkf，其中L DL11DU,fDL1b，進(jìn)行迭代計(jì)算，直到誤差滿足要求。實(shí)驗(yàn)過程一、根據(jù)維度n確定矩陣H的各個(gè)元素和b的各個(gè)分量值；......zz二、選擇計(jì)算法（GaussJGSSOR迭代法法設(shè)定初值，此外SOR法還需要設(shè)定松弛因子；三、進(jìn)行計(jì)算，直至滿足誤差要求（窮數(shù)小于0.0001；對(duì)SOR100次之后的結(jié)果及誤差值，輸出實(shí)驗(yàn)結(jié)果。計(jì)算結(jié)果及分析n6時(shí)，問題可以具體定義為計(jì)算結(jié)果如下，Gauss消去法第1次消元所選取的主元是：1第2次消元所選取的主元是：0.0833333第3次消元所選取的主元是：0.00555556第4次消元所選取的主元是第5次消元所選取的主元是：2.26757e-05第6次消元所選取的主元是：1.43155e-06......zz解得X=(0.9228 1.1937 0.1792 1.5369 0.4584 1.7680)T使用無窮數(shù)衡量誤差，可得XX*=4.9847e-10；J迭代法J法的迭代矩陣B的譜半徑為以J法不收斂；GS迭代法設(shè)定迭代初值為零，計(jì)算得到GS法的迭代矩陣G的譜半徑為：0.999998＜1，故GS法收斂，經(jīng)過541次迭代計(jì)算后，結(jié)果為X=(1.27060.18600.49021.91621.02810.4608)T使用無窮數(shù)衡量誤差，有XX*=0.9162；SOR迭代法設(shè)定迭代初值為零向量，并設(shè)定0.5，計(jì)算得到SOR法迭代矩陣譜半徑為0.5223，經(jīng)過100次迭代后的計(jì)算結(jié)果為X=(1.50780.84231.45591.98811.71640.3527)T；使用無窮數(shù)衡量誤差，有XX*=0.6473；對(duì)SOR法，可變，改變值，計(jì)算結(jié)果可以列表如下0.50.81.31.8迭代次數(shù)100100100100迭代矩陣的譜半0.52230.31550.50130.2386徑1.46041.51471.39660.40960.98430.91561.46940.93531.24401.62651.40721.25101.18271.26810.62681.81481.78910.81191.53310.26021.82700.92550.44670.45880.37320.07450.40720.2510XX綜上，四種算法的結(jié)果列表如下：算法算法Gauss消去法Jacobi法GS法SOR法（取0.5）迭代--不收斂541100次數(shù)迭代矩陣--4.308530.9999980.5223的譜半徑0.92281.19371.27061.50780.84231.45591.988X0.1792 1.5369--0.18600.49021.91621.02811.71640.35270.45841.768010.46084.9847e-10--0.91620.6473H的條件數(shù)為1.50107>>1LU分解法的誤差存在主要是由于HilbertLU所以最后得到的結(jié)果不是程的精確解，但結(jié)果顯示該法的誤差非常?。籎acobi迭代矩陣的譜半徑為4.30853，故此迭代法不收斂；GS迭代法在迭代次數(shù)為5410.05較大。GS迭代矩陣的譜半徑為0.999998，很接近1，所以GS迭代法收斂速度較慢；SOR迭代法在迭代次數(shù)為100次時(shí)誤差約為0.08，誤差較大。SOR迭代矩陣的譜半徑為，所以0.5時(shí)SOR迭代法收斂速度不是很快，但是相比于GS，SOR法的迭代速度會(huì)相應(yīng)有變化，如果選用最佳松弛因子，可以實(shí)現(xiàn)更快的收斂；（2）算法Gauss消去J法計(jì)算結(jié)果0.62690.26761.9060算法Gauss消去J法計(jì)算結(jié)果0.62690.26761.90601.8103--GS法0.53491.68390.60151.5036SOR（w=0.5）1.56340.51710.30811.5154迭代次數(shù)譜半徑4.87107--1.031011.201010.41161.21641.25761.39740.71251.47471.66121.55140.53640.84690.55360.7655----356100--6.0421310.8776n11時(shí)，算法算法Gauss消去法Jacobi法GS法SOR（w=0.5）0.11971.63540.37000.11971.63540.37001.8631計(jì)算結(jié)果0.43291.78050.93701.15560.58121.90200.67511.65721.84090.44830.67390.79100.32370.27270.13531.3537--1.65911.64681.66081.50691.85331.11921.94220.24240.62500.41420.06970.1922----1019100--8.6496410.99669.73103--5.281021.471017.1018 -7.90817.0484-1.5142 1.1601 0.3362譜半徑7.1018 -7.90817.0484-1.5142 1.1601 0.3362譜半徑0.5643譜半徑算法Gauss消去法Jacobi法GS法SOR法（w=0.5）1.15901.75540.77290.37991.48960.37991.48960.97690.77061.92970.5194計(jì)算結(jié)果2.5890-2.21280.90530.84571.79031.8764--不收斂1.39441.97681.60561.51760.20410.80030.40650.3220----262100--6.04213＞11.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法迭代不收斂，無論n如取值，其譜半徑始終大于1，因而J法不收斂，所以J迭代法不能用于Hilbert矩陣的求解；對(duì)于GS迭代法和SOR迭代法兩種法均收斂迭代法是SOR迭代法松弛因子取值為1的特例法受到 取值的影響會(huì)有不同的收斂情況可以得出GS迭代矩陣的譜半徑小于1但是很接近1，收斂速度很慢。雖然隨著維數(shù)的增大，所需迭代的次數(shù)逐漸減少，但是當(dāng)維數(shù)達(dá)到15的時(shí)候法已經(jīng)不再收斂。因此可以得出結(jié)論迭代法在Hilbert矩陣維數(shù)較低時(shí)，能夠在一定程度上滿足迭代求解的需求，不過迭代的速度很慢。另外，隨著矩陣維數(shù)的增加SOR法的誤差水平基本穩(wěn)定，而且誤差在可以接受的圍之。經(jīng)過比較可以得出結(jié)論，如果求解較低維度的Hibert矩陣問題，Gauss消去法、GS迭代法和SOR迭代法均可使用，且Gauss消去法的結(jié)果精確度較高；如果需要求解較高維度的Hibert矩陣問題，只有采用SOR迭代法。（3）總體來看，對(duì)于一般病態(tài)程組的求解，可以采用以下式：低維度下采用GaussJacobi迭代法不適宜于求解病態(tài)問題；GS迭代法可以解決維數(shù)較低的病態(tài)問題，但其譜半徑非常趨近于1，導(dǎo)致迭代算法收斂速度很慢，維數(shù)較大的時(shí)候，GS法也不再收斂；SOR法較適合于求解病態(tài)問題，特別是矩陣維數(shù)較高的時(shí)候，其優(yōu)勢(shì)更為明顯。度；可以對(duì)原程組作某些預(yù)處理，從而有效降低系數(shù)矩陣的條件數(shù)。實(shí)驗(yàn)結(jié)論對(duì)HibertGS迭代法和SOR且可以優(yōu)先使用Gauss消去法；如果需要求解較高維度的Hibert有SOR迭代法能夠求解。SOR為明顯。從本次實(shí)驗(yàn)可以看出，隨著矩陣維數(shù)的增大，SOR法所需的迭代次數(shù)減少，而且誤差基本穩(wěn)定，是解決病態(tài)問題的適宜法。附錄：程序代碼clearallclc;formatlong;%矩陣賦值n=input('矩陣H的階數(shù):n=');fori=1:nforj=1:nH(i,j)=1/(i+j-1);endendb=H*ones(n,1);disp('H矩陣為：');Hdisp('向量b%法選擇disp('選取求解式');disp('1Gauss消去法，2J迭代法，3GS迭代法，4SOR迭代法');a=input('求解式序號(hào)：');%Gauss消去法ifa==1;H1=H;b1=b;fork=1:nifH1(k,k)==0disp('主元為零，Gauss消去法無法進(jìn)行');breakendfprintf('第%d次消元所選取的主元是：%g\n',k,H1(k,k))forp=k+1:nm5=-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);endendx1(n)=b1(n)/H1(n,n);fork=n-1:-1:1forv=k+1:nb1(k)=b1(k)-H1(k,v)*x1(v);endx1(k)=b1(k)/H1(k,k);enddisp('Gauss消去法解為:');disp(x1);disp('解與精確解之差的無窮數(shù)');norm((x1-a),inf)endD=diag(diag(H));L=-tril(H,-1);U=-triu(H,1);%J迭代法ifa==2;%給定初始x0ini=input('初始值設(shè)定:x0=');x0(:,1)=ini*diag(ones(n));disp('初始解向量為:');x0xj(:,1)=x0(:,1);B=(D^(-1))*(L+U);f=(D^(-1))*b;fprintf('(J法B矩陣譜半徑為：ifvrho(B)<1;form2=1:5000xj(:,m2+1)=B*xj(:,m2)+fj;ifnorm((xj(:,m2+1)-xj(:,m2)),inf)<0.0001breakendenddisp('J法計(jì)算結(jié)果為:');xj(:,m2+1)disp('解與精確解之差的無窮數(shù)');norm((xj(:,m2+1)-diag(ones(n))),inf)disp('J迭代法迭代次數(shù):');m2elsedisp('由于B矩陣譜半徑大于1，因而J法不收斂');endend%GS迭代法ifa==3;%給定初始x0ini=input('初始值設(shè)定:x0=');x0(:,1)=ini*diag(ones(n));disp('初始解向量為:');x0xG(:,1)=x0(:,1);G=inv(D-L)*U;fG=inv(D-L)*b;fprintf('GS法G矩陣譜半徑為：ifvrho(G)<1form3=1:5000xG(:,m3+1)=G*xG(:,m3)+fG;ifnorm((xG(:,m3+1)-xG(:,m3)),inf)<0.0001break;endenddisp('GS迭代法計(jì)算結(jié)果:');xG(:,m3+1)disp('解與精確解之差的無窮數(shù)');norm((xG(:,m3+1)-diag(ones(n))),inf)disp('GS迭代法迭代次數(shù):');m3elsedisp('由于G矩陣譜半徑大于1，因而GS法不收斂');endend%SOR迭代法ifa==4;%給定初始x0ini=input('初始值設(shè)定:x0=');x0(:,1)=ini*diag(ones(n));disp('初始解向量為:');x0A=H;fori=1:nb(i)=sum(A(i,:));endx_star=ones(n,1);formatlongw=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);%迭代初值fori=1:time_maxx=Lw*x+f;enddisp('SOR迭代法得到的解為');xdisp('解與精確解之差的無窮數(shù)');norm((x_star-x),inf)endpause三、實(shí)驗(yàn)4.1題目：對(duì)牛頓法和擬牛頓法。進(jìn)行非線性程組的數(shù)值求解較其迭代次數(shù)、CPU時(shí)間等。12xx24x

702 x x 11022 x

x01,1,T1 2 3x310x802 3 3xcosx

0.50 1 232x281x0.2sin

1.060,x00,0,0T1 2

31 exx

20x

30 12 3 3x，結(jié)果又如？反復(fù)選取不同的初值，比較其結(jié)果?？偨Y(jié)歸納你的實(shí)驗(yàn)結(jié)果，試說明各種法適用的問題。算法設(shè)計(jì)對(duì)需要求解的非線性程組Fx0而言，牛頓法和擬牛頓法的迭代公式如下，牛頓法：

Fxk

xk

F xk xkxkxk牛頓法為單步迭代法，需要取一個(gè)初值。（Broyden秩1法）

Xk1X

Xk

A1F xkk

TskA skkk

ykAk

sk

T sk sk

其中sk

xk1

x

;y

F xk

F xk ，擬牛頓法不需要求解F(x)的導(dǎo)數(shù)，因此節(jié)省了大量的運(yùn)算時(shí)間，但需要給定矩陣Ak

的初值，取為Fx0 。實(shí)驗(yàn)過程一、輸入初值；三、輸出數(shù)據(jù)；計(jì)算結(jié)果及分析首先求解程組（1，在這里，設(shè)定精度要求為x xkk2

104，法法牛頓法擬牛頓法初始值x(0)(1,1,1Tx(0)(1,1,1Tx10.59140.7151計(jì)算結(jié)x21.00311.4940果Xx30.83060.5304迭代次數(shù)313CPU計(jì)算時(shí)間/s3.7778152.739349迭代都需要求解矩陣的逆，所以牛頓法每次迭代的CPU計(jì)算時(shí)間更長(zhǎng)。之后求解程組（2，同樣設(shè)定精度要求為x x法初始值x(0)牛頓法(0,0,0)法初始值x(0)牛頓法(0,0,0)T擬牛頓法x(0)(0,0,0)T

104計(jì)算結(jié)x10.96990.3600果Xx20.34280.1856x3x3-0.0483-0.8871迭代次數(shù)412CPU計(jì)算時(shí)間/s2.7224373.920195同樣地，可以看出，在初始值相同情況下，牛頓法和擬牛頓法在達(dá)到同樣計(jì)的，由于每次迭代中有求解矩陣的逆的運(yùn)算，牛頓法每次迭代的CPU計(jì)算時(shí)間較長(zhǎng)。對(duì)程組（1，取其他初值，計(jì)算結(jié)果列表如下，同樣設(shè)定精度要求為x xkk2

104初始值 法計(jì)算結(jié)果x(0)(0,0,0)T迭代次數(shù)CPU計(jì)算時(shí)間/s計(jì)算結(jié)果x(0)(0.5,0.5,0.5)T迭代次數(shù)CPU計(jì)算時(shí)間/s計(jì)算結(jié)果x(0)(1,1,1T迭代次數(shù)CPU計(jì)算時(shí)間/s計(jì)算結(jié)果x(0)(2,2,2)T迭代次數(shù)CPU計(jì)算時(shí)間/sx(0)(10,10,10)T 計(jì)算結(jié)果

0.59141.00310.830543.9071640.59141.00310.830548.1272860.59141.00310.830633.7778150.59141.00310.830643.8356979.3722-5.4773

擬牛頓法9.7894-5.534618.1205584.8180199.4591-5.354918.280727355.6260230.71511.49400.5304132.7393490.57731.24580.01362.879070Matlab警告矩陣接近奇異值，程序進(jìn)入長(zhǎng)期x(0)(50,50,50)T

迭代次數(shù)CPU計(jì)算時(shí)間/s計(jì)算結(jié)果迭代次數(shù)CPU計(jì)算時(shí)間/s

18.8605194.0338680.73351.15360.49221312.243263

循環(huán)計(jì)算中----Matlab警告矩陣接近奇異值，程序進(jìn)入長(zhǎng)期循環(huán)計(jì)算中----從上表可以發(fā)現(xiàn)，程組（1）存在另一個(gè)在(9.2,-5.6,18.1)T附近的不動(dòng)點(diǎn)，初值的選取會(huì)直接影響到牛頓法和擬牛頓法最后的收斂點(diǎn)?，F(xiàn)迭代次數(shù)增加直至無法收斂的情況；斂情況不如牛頓法好（初值不夠接近時(shí)，甚至?xí)霈F(xiàn)奇異矩陣的情況，但由于牛頓法的求解比較復(fù)雜，計(jì)算時(shí)間較長(zhǎng)；同樣的，對(duì)程組2，取其他初值，計(jì)算結(jié)果列表如下，同樣設(shè)定精度要求為x xkk2

104初始值x(0)(0,0,0)Tx(0)(0.25,0.25,0.25)

法計(jì)算結(jié)果迭代次數(shù)CPU計(jì)算時(shí)間/s計(jì)算結(jié)果迭代次數(shù)CPU計(jì)算時(shí)間/s

牛頓法0.96990.3428-0.048342.7224370.1085

0.36000.1856-0.8871123.9201950.07530.5427-0.2266-0.650776.71295575.0471115.6197520.09160.5427-0.2266-0.650776.71295575.0471115.6197520.09161.0e+02*0.0410-0.4775x(0)(1,1,1Tx(0)(2,2,2)Tx(0)(10,10,10)Tx(0)(20,20,

迭代次數(shù)CPU計(jì)算時(shí)間/s計(jì)算結(jié)果迭代次數(shù)CPU計(jì)算時(shí)間/s計(jì)算結(jié)果迭代次數(shù)CPU計(jì)算時(shí)間/s計(jì)算結(jié)果迭代次數(shù)CPU計(jì)算時(shí)間/s計(jì)算結(jié)果迭代次數(shù)CPU計(jì)算時(shí)間/s

-0.567263.5406680.01520.6711-0.786272.2005710.00050.0503-0.828681.7190720.20220.1686-0.25001492.797116矩陣為奇異值，無法輸出準(zhǔn)確結(jié)果----

-0.28861.2843623.3878291.0e+04*0.0770-0.12951.8650552.640901出準(zhǔn)確結(jié)果----矩陣為奇異值，無法輸出準(zhǔn)確結(jié)果----矩陣為奇異值，無法輸出準(zhǔn)確結(jié)果----在這里，與前文類似的發(fā)現(xiàn)不再贅述。圍上選取初值并最終收斂到精確解附近；在初始值較接近于不動(dòng)點(diǎn)時(shí)，牛頓法和擬牛頓法計(jì)算所得到的結(jié)果是基本相同的，雖然迭代次數(shù)有所差別，但計(jì)算總的所需時(shí)間相近。（3）合應(yīng)用牛頓法，其對(duì)初始值敏感程度較低，算法具有很好的收斂性。另外，需要說明的是，每次計(jì)算給出的CPU時(shí)間與計(jì)算機(jī)當(dāng)時(shí)的運(yùn)行狀態(tài)參考價(jià)值。實(shí)驗(yàn)結(jié)論現(xiàn)象；1的2CPU耗時(shí)反而更少。另外，就程組壓縮映射區(qū)間來說，一般而言，對(duì)于在區(qū)間函數(shù)呈現(xiàn)單調(diào)或者所以對(duì)初始值的敏感程度較小。附錄：程序代碼%程1，牛頓法tic;formatlong;%%初值disp('第1個(gè)分量為：');b=input('第2個(gè)分量為：');c=input('第3個(gè)分量為：');disp(