數(shù)值分析實驗報告

1、實驗五 解線性方程組的直接方法實驗5.1 （主元的選取與算法的穩(wěn)定性）問題提出：gauss消去法是我們在線性代數(shù)中已經(jīng)熟悉的。但由于計算機的數(shù)值運算是在一個有限的浮點數(shù)集合上進行的，如何才能確保gauss消去法作為數(shù)值算法的穩(wěn)定性呢？gauss消去法從理論算法到數(shù)值算法，其關鍵是主元的選擇。主元的選擇從數(shù)學理論上看起來平凡，它卻是數(shù)值分析中十分典型的問題。實驗內容：考慮線性方程組 編制一個能自動選取主元，又能手動選取主元的求解線性方程組的gauss消去過程。實驗要求：（1）取矩陣，則方程有解。取n=10計算矩陣的條件數(shù)。讓程序自動選取主元，結果如何？（2）現(xiàn)選擇程序中手動選取主元的功能。每步消

2、去過程總選取按模最小或按模盡可能小的元素作為主元，觀察并記錄計算結果。若每步消去過程總選取按模最大的元素作為主元，結果又如何？分析實驗的結果。（3）取矩陣階數(shù)n=20或者更大，重復上述實驗過程，觀察記錄并分析不同的問題及消去過程中選擇不同的主元時計算結果的差異，說明主元素的選取在消去過程中的作用。（4）選取其他你感興趣的問題或者隨機生成矩陣，計算其條件數(shù)。重復上述實驗，觀察記錄并分析實驗結果。思考題一：（vadermonde矩陣）設 ，其中，（1）對n=2,5,8，計算a的條件數(shù)；隨n增大，矩陣性態(tài)如何變化？（2）對n=5，解方程組ax=b；設a的最后一個元素有擾動10-4，再求解ax=b（3

3、）計算（2）擾動相對誤差與解的相對偏差，分析它們與條件數(shù)的關系。（4）你能由此解釋為什么不用插值函數(shù)存在定理直接求插值函數(shù)而要用拉格朗日或牛頓插值法的原因嗎？相關matlab函數(shù)提示：zeros(m,n) 生成m行，n列的零矩陣ones(m,n) 生成m行，n列的元素全為1的矩陣eye(n) 生成n階單位矩陣rand(m,n) 生成m行,n列(0,1)上均勻分布的隨機矩陣diag(x) 返回由向量x的元素構成的對角矩陣tril(a) 提取矩陣a的下三角部分生成下三角矩陣triu(a) 提取矩陣a的上三角部分生成上三角矩陣rank(a) 返回矩陣a的秩det(a) 返回方陣a的行列式inv(a)

4、 返回可逆方陣a的逆矩陣v，d=eig(a) 返回方陣a的特征值和特征向量norm(a,p) 矩陣或向量的p范數(shù)cond(a,p) 矩陣的條件數(shù)l，u，p=lu(a) 選列主元lu分解r=chol(x) 平方根分解hi=hilb(n) 生成n階hilbert矩陣5.1實驗過程：5.1.1程序：function x=gauss(n,r)n=input(請輸入矩陣a的階數(shù):n=)a=diag(6*ones(1,n)+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)b=a*ones(n,1)for i=1:4p=input(條件數(shù)對應的范數(shù)是p-范數(shù)：p=)pp=

5、cond(a,p)endpausem,n=size(a);nb=n+1;ab=a br=input(請輸入是否為手動，手動輸入1，自動輸入0：r=)for i=1:n-1 if r=0 pivot,p=max(abs(ab(i:n,i); ip=p+i-1; if ip=i ab(i ip,:)=ab(ip i,:);disp(ab); pause end end if r=1 i=i ip=input(輸入i列所選元素所處的行數(shù)：ip=); ab(i ip,:)=ab(ip i,:);disp(ab); pause end pivot=ab(i,i); for k=i+1:n ab(k,i:

6、nb)=ab(k,i:nb)-(ab(k,i)/pivot)*ab(i,i:nb); end disp(ab); pauseendx=zeros(n,1);x(n)=ab(n,nb)/ab(n,n);for i=n-1:-1:1 x(i)=(ab(i,nb)-ab(i,i+1:n)*x(i+1:n)/ab(i,i);end5.1.2實驗結果如下： 1.按照實驗要求一：取矩陣a的階數(shù):n=10且自動選取主元，程序結果運行如下：（2） 現(xiàn)選擇程序中手動選取主元的功能，觀察并記錄計算結果。 選取絕對值最大的元素為主元：程序運行開始如第一問的截圖也是求范數(shù)故這里不在給出。the answer is ：

7、1 1 1 1 1 1 1 1 1 1 選取絕對值最小的元素為主元：the answer is： 1.0e+003*（inf 0.007 0.0057 -0.0410 -0.0303 0.3430 0.2577 -2.7290 -2.0463 2.7308）取矩陣a的階數(shù):n=20，手動選取主元： 選取絕對值最大的元素為主元：the answer is ：1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 選取絕對值最小的元素為主元：the answer is： 1.0e+007*（-inf 0.0000 0.0000 -0.0000 -0.0000 0.0000

8、0.0000 -0.0003 -0.0002 0.0022 0.0016 -0.0175 -0.0131 0.1398 0.1049 -1.1185 -0.8389 8.9478 6.7109 -8.9478）修改程序如下：function x=gaussong(n,r)n=input(請輸入矩陣a的階數(shù):n=)a=hilb(n)b=a*ones(n,1)for i=1:4p=input(條件數(shù)對應的范數(shù)是p-范數(shù)：p=)pp=cond(a,p)endpausem,n=size(a);nb=n+1;ab=a br=input(請輸入是否為手動，手動輸入1，自動輸入0：r=)for i=1:n-

9、1 if r=0 pivot,p=max(abs(ab(i:n,i); ip=p+i-1; if ip=i ab(i ip,:)=ab(ip i,:);disp(ab); pause end end if r=1 i=i ip=input(輸入i列所選元素所處的行數(shù)：ip=); ab(i ip,:)=ab(ip i,:);disp(ab); pause end pivot=ab(i,i); for k=i+1:n ab(k,i:nb)=ab(k,i:nb)-(ab(k,i)/pivot)*ab(i,i:nb); end disp(ab); pauseendx=zeros(n,1);x(n)=a

10、b(n,nb)/ab(n,n);for i=n-1:-1:1 x(i)=(ab(i,nb)-ab(i,i+1:n)*x(i+1:n)/ab(i,i);end所求范數(shù)為：自動輸入結果為：ans = 1.0000 1.00001 .0000 1.0000 1.0002 0.9996 1.0007 0.9993 1.0004 0.9999選取絕對值最大的元素為主元結果為：the answer is ：nan nan nan nan nan inf -inf -inf 281.3945 -283.3708選取絕對值最小的元素為主元結果為：the answer is ：nan nan nan -inf

11、-5.8976 -1.9243 -2.0291 -4.9972 23.4548 -11.10125.1.3 對實驗結果進行分析：5.1.3.1 對實驗要求一的結果進行分析：對于gauss消去法就是用行的初等變換將原線性方程組系數(shù)矩陣轉化為簡單形式，從而進行求解，缺點是迭代次數(shù)可能較多，效率不高，且在消去過程中不可以將主元素很小的做除數(shù)，否則將導致其他元素數(shù)量級的嚴重增長和舍入誤差的擴散，使得計算解不可靠。5.1.3.2 對實驗要求二的結果進行分析：通過每次選取最大或最小的主元可以發(fā)現(xiàn)取絕對值大的元素作為主元比取絕對值小的元素作為主元時產(chǎn)生的結果比較準確，即選取絕對值小的主元時結果產(chǎn)生了較大的誤

12、差，條件數(shù)越大產(chǎn)生的誤差就越大。所以應盡量避免很小的數(shù)作為除數(shù)。5.1.3.3 對實驗要求三的結果進行分析：此要求是對要求一和要求二的一個延續(xù)，通過實驗結果可以看出若采用很小的數(shù)作為主元迭代次數(shù)越多導致的結果越不可靠，甚至出現(xiàn)錯誤。5.1.3.4 對實驗要求四的結果進行分析：對新矩陣進行實驗發(fā)現(xiàn)依然符合上述規(guī)律，可以知道，在進行迭代時主元的選擇與算法的穩(wěn)定性有密切的聯(lián)系選取絕對值大的元素作為主元比絕對值小的元素作為主元時對結果產(chǎn)生的誤差較小。條件數(shù)越大對用gauss消去法解線性方程組時，對結果產(chǎn)生的誤差就越大。5.1.4實驗總結：1. 在用gauss消去法解線性方程組時，主元的選取與算法的穩(wěn)定

13、性有密切的聯(lián)系，選取適當?shù)闹髟欣诘贸龇€(wěn)定的算法，2. 在算法的過程中，選取絕對值較大的主元比選取絕對值較小的主元更有利于算法的穩(wěn)定，選取絕對值最大的元素作為主元時，得出的結果相對較準確較穩(wěn)定。3. 條件數(shù)越小，對用這種方法得出的結果更準確。4. 在算除法的過程中要盡量避免使用較小的數(shù)做為除數(shù)，以免發(fā)生結果數(shù)量級加大，使大數(shù)吃掉小數(shù)，產(chǎn)生舍入誤差。實驗5.2（線性代數(shù)方程組的性態(tài)與條件數(shù)的估計）問題提出：理論上，線性代數(shù)方程組的攝動滿足 矩陣的條件數(shù)確實是對矩陣病態(tài)性的刻畫，但在實際應用中直接計算它顯然不現(xiàn)實，因為計算通常要比求解方程還困難。實驗內容：matlab中提供有函數(shù)“condest

14、”可以用來估計矩陣的條件數(shù)，它給出的是按1-范數(shù)的條件數(shù)。首先構造非奇異矩陣a和右端，使得方程是可以精確求解的。再人為地引進系數(shù)矩陣和右端的攝動，使得充分小。實驗要求：（1）假設方程ax=b的解為x，求解方程，以1-范數(shù)，給出的計算結果。（2）選擇一系列維數(shù)遞增的矩陣（可以是隨機生成的），比較函數(shù)“condest”所需機器時間的差別.考慮若干逆是已知的矩陣，借助函數(shù)“eig”很容易給出cond2(a)的數(shù)值。將它與函數(shù)“cond(a,2)”所得到的結果進行比較。（3）利用“condest”給出矩陣a條件數(shù)的估計，針對（1）中的結果給出的理論估計，并將它與（1）給出的計算結果進行比較，分析所得結

15、果。注意，如果給出了cond(a)和的估計，馬上就可以給出的估計。（4）估計著名的hilbert矩陣的條件數(shù)。5.2 實驗過程如下： 5.2.1.1 實驗要求一的程序如下：function n=jisuan(n)a=fix(100*rand(n)+1 x=ones(n,1) b=a*x data=rand(n)*0.00001 datb=rand(n,1)*0.00001 a=a+datab=b+datbx0=get(a,b) x1=norm(x0-x,1)/norm(x,1) function x=get(a,b)m,n=size(a);nb=n+1;ab=a b;for i=1:n-1 p

16、ivot=ab(i,i); for k=i+1:n ab(k,i:nb)=ab(k,i:nb)-(ab(k,i)/pivot)*ab(i,i:nb); endendx=zeros(n,1);x(n)=ab(n,nb)/ab(n,n);for i=n-1:-1:1 x(i)=(ab(i,nb)-ab(i,i+1:n)*x(i+1:n)/ab(i,i);end5.2.1.2 實驗要求一程序運行結果如下：系數(shù)矩陣a為：7018142952536398247346580289074421962620992338538830595879897467437769加擾動后的系數(shù)矩陣a為：70.0000046

17、 18.0000042 14.0000099 29.0000032 52.0000044 53.0000013 63.0000057 98.0000030 2.0000079 47.0000096 34.0000093 65.0000021 80.0000079 28.0000087 90.0000044 7.0000073 44.0000068 21.0000061 96.0000006 26.0000002 20.0000050 99.0000041 23.0000021 38.0000063 53.0000060 88.0000077 30.0000021 59.0000074 58.0

18、000084 79.0000037 89.0000005 74.0000097 67.0000064 43.0000027 77.0000063 69.0000058 b值為：236309270302367419加擾動后的b值為：236.00000451309.00000044270.00000027302.00000313367.00000013419.00000384 data的值為：4.610951267e-064.153748604e-069.900825926e-063.200355775e-064.399243096e-061.337727485e-065.678287124e-0

19、63.04998677e-067.888616922e-069.600986004e-069.333801082e-062.071327296e-067.942106514e-068.743671716e-064.386585338e-067.266317666e-066.833323243e-066.071989445e-065.91825935e-071.50094987e-074.983113034e-064.119532082e-062.125598643e-066.298878488e-066.028690857e-067.6795039e-062.139633316e-067.44

20、5657831e-068.392382403e-063.704768261e-065.02688037e-079.708449393e-066.434922879e-062.679472507e-066.287846e-065.75147779e-06 datb的值為：4.514248268e-064.38953253e-072.7185123e-073.126850481e-061.28625747e-073.839672885e-06 xx的值為：0.999999684491.00000038780.999999143481.00000065171.00000221560.99999754

21、53 x0的值為：1.146958549e-06 x1的值為：6.8990e-0075.2.1.3實驗結果為：的計算結果為：6.8990e-0075.2.2.1 實驗要求二的程序如下：function cond2(a) b=a*a;v1,d1=eig(b);v2,d2=eig(b(-1);cond2a=sqrt(max(max(d1)*sqrt(max(max(d2)endfor n=10:10:100n=n a=fix(100*randn(n); condesta=condest(a) cond2(a) conda2=cond(a,2) pause end5.2.2.2 實驗要求二的程序運行

22、結果如下：ncondestacond2aconda2101.46e+0242.50783044142.507830441204.59e+021.23e+021.23e+02304.05e+0279.26495488579.264954885402.21e+034.26e+024.26e+02503.02e+034.08e+024.08e+02604.75e+037.78e+027.78e+02704.69e+035.14e+025.14e+02805.47e+034.89e+024.89e+02905.66e+035.50e+025.50e+021004.47e+035.06e+025.06e

23、+025.2.3.1 實驗要求三的程序如下：function bijiao(n) a=fix(100*rand(n)+1; x=ones(n,1); b=a*x; data=rand(n)*0.00001; datb=rand(n,1)*0.00001; a=a+data;b=b+datb;xx=geshow(a,b); x1=norm(xx-x,1)/norm(x,1) x2=cond(a)/(1-norm(inv(a)*norm(xx-x)*(norm(xx-x)/(norm(a)+norm(datb)/norm(b) datx=abs(x1-x2) 5.2.3.2 實驗要求三的程序運行結

24、果如下：給出對的估計是：7.310559817408125e-007的理論結果是： 3.828481757617297e-007結果相差： 3.482078059790828e-0075.2.4.1 實驗要求四的程序如下：for n=4:11 n=n hi=hilb(n); cond1hi=cond(hi,1) cond2hi=cond(hi,2) condinfhi=cond(hi,inf) pauseend5.2.4.2 實驗要求四的程序運行結果如下：ncond1hicond2hicondinfhi42.84e+041.55e+042.84e+0459.44e+054.77e+059.44

25、e+0562.91e+071.50e+072.91e+0779.85e+084.75e+089.85e+0883.39e+101.53e+103.39e+1091.10e+124.93e+111.10e+12103.54e+131.60e+133.54e+13111.23e+155.22e+141.23e+15討論：線性代數(shù)方程組的性態(tài)與條件數(shù)有著很重要的關系，既矩陣的條件數(shù)是刻畫矩陣性質的一個重要的依據(jù)，條件數(shù)越大，矩陣“病態(tài)”性越嚴重，在解線性代數(shù)方程組的過程中較容易產(chǎn)生比較大的誤差，則在實際問題的操作過程中，我們必須要減少對條件數(shù)來求解，把條件數(shù)較大的矩陣化成條件數(shù)較小的矩陣來進行求解。

26、實驗總結：在本次實驗中，使我們知道了矩陣條件數(shù)對線性代數(shù)方程組求解的影響，條件數(shù)越大，對最后解的影響的越大，hilbert矩陣是一個很”病態(tài)”的矩陣，他的條件數(shù)隨著階數(shù)的增加而增大，每增加一階，條件數(shù)就增大一個數(shù)量級，在求解的過程中要盡量避免hilbert矩陣實驗六 解線性方程組的迭代法實驗6.1（病態(tài)的線性方程組的求解）問題提出：理論的分析表明，求解病態(tài)的線性方程組是困難的。實際情況是否如此，會出現(xiàn)怎樣的現(xiàn)象呢？實驗內容：考慮方程組hx=b的求解，其中系數(shù)矩陣h為hilbert矩陣， 這是一個著名的病態(tài)問題。通過首先給定解（例如取為各個分量均為1）再計算出右端b的辦法給出確定的問題。實驗要求

27、：（1）選擇問題的維數(shù)為6，分別用gauss消去法、j迭代法、gs迭代法和sor迭代法求解方程組，其各自的結果如何？將計算結果與問題的解比較，結論如何？（2）逐步增大問題的維數(shù)，仍然用上述的方法來解它們，計算的結果如何？計算的結果說明了什么？（3）討論病態(tài)問題求解的算法由于本實驗用到迭代法，故先給出迭代法相關資料： 迭代法也稱輾轉法，是一種不斷用變量的舊值遞推新值的過程,跟迭代法相對應的是直接法（或者稱為一次解法），即一次性解決問題。迭代法又分為精確迭代和近似迭代?！岸址ā焙汀芭ｎD迭代法”屬于近似迭代法。迭代算法是用計算機解決問題的一種基本方法。它利用計算機運算速度快、適合做重復性操作的特點

28、，讓計算機對一組指令（或一定步驟）進行重復執(zhí)行，在每次執(zhí)行這組指令（或這些步驟）時，都從變量的原值推出它的一個新值。 迭代是數(shù)值分析中通過從一個初始估計出發(fā)尋找一系列近似解來解決問題（一般是解方程或者方程組）的過程，為實現(xiàn)這一過程所使用的方法統(tǒng)稱為迭代法（iterative method）。一般可以做如下定義：對于給定的線性方程組（這里的同為矩陣，任意線性方程組都可以變換成此形式），用公式(括號中為上標，代表迭代k次得到的x，初始時k=0)逐步帶入求近似解的方法稱為迭代法（或稱一階定常迭代法）。如果k趨向無窮大時存在，記為，稱此迭代法收斂。顯然就是此方程組的解，否則稱為迭代法發(fā)散。跟迭代法相對

29、應的是直接法（或者稱為一次解法），即一次性的快速解決問題，例如通過開方解決方程。一般如果可能，直接解法總是優(yōu)先考慮的。但當遇到復雜問題時，特別是在未知量很多，方程為非線性時，我們無法找到直接解法（例如五次以及更高次的代數(shù)方程沒有解析解，參見阿貝爾定理），這時候或許可以通過迭代法尋求方程（組）的近似解。最常見的迭代法是牛頓法。其他還包括最速下降法、共軛迭代法、變尺度迭代法、最小二乘法、線性規(guī)劃、非線性規(guī)劃、單純型法、懲罰函數(shù)法、斜率投影法、遺傳算法、模擬退火法等等。利用迭代算法解決問題，需要做好以下三個方面的工作： （1）確定迭代變量在可以用迭代算法解決的問題中，至少存在一個直接或間接地不斷由舊

30、值遞推出新值的變量，這個變量就是迭代變量。 （2）建立迭代關系式所謂迭代關系式，指如何從變量的前一個值推出其下一個值的公式（或關系）。迭代關系式的建立是解決迭代問題的關鍵，通?？梢皂樛苹虻雇频姆椒▉硗瓿?。 （3）對迭代過程進行控制在什么時候結束迭代過程？這是編寫迭代程序必須考慮的問題。不能讓迭代過程無休止地重復執(zhí)行下去。迭代過程的控制通?？煞譃閮煞N情況：一種是所需的迭代次數(shù)是個確定的值，可以計算出來；另一種是所需的迭代次數(shù)無法確定。對于前一種情況，可以構建一個固定次數(shù)的循環(huán)來實現(xiàn)對迭代過程的控制；對于后一種情況，需要進一步分析出用來結束迭代過程的條件。實驗過程：6.1.1.1 實驗要求一的ga

31、uss實現(xiàn)程序如下：gauss消去法程序：function x=gaussong(n,r)a=hilb(n)b=a*ones(n,1)for i=1:4p=input(條件數(shù)對應的范數(shù)是p-范數(shù)：p=)pp=cond(a,p)endpausem,n=size(a);nb=n+1;ab=a br=input(請輸入是否為手動，手動輸入1，自動輸入0：r=)for i=1:n-1 if r=0 pivot,p=max(abs(ab(i:n,i); ip=p+i-1; if ip=i ab(i ip,:)=ab(ip i,:);disp(ab); pause end end if r=1 i=i i

32、p=input(輸入i列所選元素所處的行數(shù)：ip=); ab(i ip,:)=ab(ip i,:);disp(ab); pause end pivot=ab(i,i); for k=i+1:n ab(k,i:nb)=ab(k,i:nb)-(ab(k,i)/pivot)*ab(i,i:nb); end disp(ab); pauseendx=zeros(n,1);x(n)=ab(n,nb)/ab(n,n);for i=n-1:-1:1 x(i)=(ab(i,nb)-ab(i,i+1:n)*x(i+1:n)/ab(i,i);end6.1.1.2實驗要求一的gauss實現(xiàn)程序運行結果如下：6.1.2

33、.1實驗要求一的j程序如下：function n=jd(n)a=hilb(n);for i=1:n; a0(i)=1; x(i)=0;endb=a*a0; for i=1:100; y=x; for j=1:n; x(j)=b(j)/a(j,j); for k=1:j-1; x(j)=x(j)-a(j,k)*y(k)/a(j,j); endfor k=j+1:n; x(j)=x(j)-a(j,k)*y(k)/a(j,j);end endendx6.1.2.2實驗要求一的j程序運行結果如下：6.1.3.1實驗要求一的gs程序如下：function n=gs(n)a2=hilb(n);for i=

34、1:n; a02(i)=1; x2(i)=0;endb2=a2*a02;for i=1:100000; for j=1:n; x2(j)=b2(j)/a2(j,j); for k=1:j-1; x2(j)=x2(j)-a2(j,k)*x2(k)/a2(j,j); endfor k=j+1:n; x2(j)=x2(j)-a2(j,k)*x2(k)/a2(j,j);end endendx26.1.3.2實驗要求一的gs程序運行結果如下：6.1.4.1實驗要求一的sor程序如下：function n ss=sor(n,ss)a3=hilb(n);for i=1:n; a03(i)=1; x3(i)=

35、0;endb3=a3*a03;for i=1:100000; for j=1:n; rc=x3(j); x3(j)=b3(j)/a3(j,j); for k=1:j-1; x3(j)=x3(j)-a3(j,k)*x3(k)/a3(j,j); endfor k=j+1:n; x3(j)=x3(j)-a3(j,k)*x3(k)/a3(j,j);endx3(j)=(1-ss)*rc+ss*x3(j); endendx36.1.4.2實驗要求一的sor程序運行結果如下：（注意由于運行結果迭代過程很長故不體現(xiàn)在實驗報告中）6.2.1實驗要求二的gauss程序運行結果如下： x = 0 0 0 0 0 0

36、 0 0 0 0 0 0 0 0 0 0 0 0 0 0.99756.2.2實驗要求二的j程序運行結果如下：6.2.3實驗要求二的gs程序運行結果如下：6.2.4實驗要求二的sor程序運行結果如下：實驗要求二的結果分析：選擇問題的維數(shù)為20時：1. 用gauss消去法求得的解與精確解相差很大，說明能否得到優(yōu)秀的解取決于算法的穩(wěn)定性，如果算法不夠穩(wěn)定產(chǎn)生的結果將變的無法理喻。2. 取初始向量為0，用j迭代方法迭代發(fā)散，無法求解；3. 取初始向量為0，用gs迭代方法迭代不發(fā)散，能求得解，但收斂非常緩慢，而且迭代次數(shù)越多，與準確解的偏差就越大，說明gs迭代適合迭代次數(shù)少的，但是通常我們無法得知需要迭

37、代的次數(shù)。4. 取初始向量為0，用sor迭代方法迭代不發(fā)散，能求得解，但同樣收斂非常緩慢?？傊?，從上面的結果可以看出當病態(tài)問題的階數(shù)升高時作為直接法的gauss消去法又能求解變成不能求解。而gs和sor迭代法在階數(shù)升高時仍能求解。但在階數(shù)較低時直接法能求得精確解而迭代發(fā)卻總存在一定的誤差?？梢娭苯臃ㄅc迭代法各有各的優(yōu)勢與不足。關于病態(tài)問題的求解，主要的方法是對原方程作某些預處理，以降低系數(shù)矩陣的條件數(shù)。可以采取將系數(shù)矩陣a的每一行本別乘上適當常數(shù)的方法。即找到可逆的對角陣和使方程組化為 理論上最好選擇對角陣滿足：。補充：本實驗用到了gauss消去法、j迭代法、gs迭代法、sor迭代，故對其重新

38、理解學習了一下：gauss消去法：1 將其增廣矩陣化為行階梯形 2 若最后有形如0 0 . 0 a （a0)的行則無解3 若含有自由變量則有無窮組解 4 原方程有唯一解。采用回代求解。 至于有無窮組解的方程組的求解，需將其化為行最簡形矩陣，其方法稱為高斯-若爾當消元法。j迭代法：一、算法理論迭代格式的引出是依據(jù)迭代法的基本思想：構造一個向量系列，使其收斂至某個極限，則就是要求的方程組的準確解。迭代將方程組： 在假設，改寫成 如果引用系數(shù)矩陣，及向量方程組（1）和（2）分別可寫為：及，這樣就得到了迭代格式用迭代解方程組時，就可任意取初值帶入迭代可知式，然后求。但是，比較大的時候，寫方程組（1）和

39、（2）是很麻煩的，如果直接由，能直接得到，就是矩陣與向量的運算了，那么如何得到，呢實際上，如果引進非奇異對角矩陣將分解成：要求的解，實質上就有，而是非奇異的，所以存在，從而有我們在這里不妨令，就得到迭代格式：現(xiàn)在考慮迭代法的計算程序 float a33=10,-2,-1,-2,10,-1,-1,-2,-5;float b3=3,15,10;分別代表的系數(shù)和等號右邊的常數(shù)項,即 先輸入方程，運行main函數(shù)，如果first不為null，則執(zhí)行if括號里的，否則執(zhí)行else里面的，最后會調用方法sum（）。在sum（）中sum=amn*xn+sum;y=(bm-sum)/amm，之后運行for循環(huán)

40、，最后輸出結果，算法結束。二、算法框圖 調用方法 開始讀入輸出 結束gs迭代法：1.高斯 - 塞德爾迭代法 公式的矩陣形式 首先將高斯 - 塞德爾迭代法的公式表示為矩陣形式，為此設 這里 是系數(shù)矩陣的對角部分,是嚴格下三角部分，是嚴格上三角部分，則高斯 - 塞德爾迭代法的公式可表示為 用矩陣 乘等式兩邊得 再用矩陣 乘等式兩邊得 其中矩陣 稱為高斯塞德爾迭代矩陣。由此可見，高斯 - 塞德爾迭代法是一般迭代法中迭代矩陣為 的特殊情形。需要指出的是，由于矩陣 難于計算，所以式（2）多用在理論分析中。2. 高斯塞德爾迭代法計算框圖（見圖）高斯塞德爾迭代法計算框圖松弛法：其中相當于在的基礎上加個余項生

41、成，下面令，希望通過選取來加速收斂，這就是松弛法。實驗七 非線性方程求根實驗7.1（迭代法、初始值與收斂性）實驗目的：初步認識非線性問題的迭代法與線性問題迭代法的差別，探討迭代法及初始值與迭代收斂性的關系。問題提出：迭代法是求解非線性方程的基本思想方法，與線性方程的情況一樣，其構造方法可以有多種多樣，但關鍵是怎樣才能使迭代收斂且有較快的收斂速度。實驗內容：考慮一個簡單的代數(shù)方程針對上述方程，可以構造多種迭代法，如 在實軸上取初始值x0，請分別用迭代（7.1）-（7.3）作實驗，記錄各算法的迭代過程。實驗要求：（1）取定某個初始值，分別計算（7.1）-（7.3）迭代結果，它們的收斂性如何？重復選

42、取不同的初始值，反復實驗。請自選設計一種比較形象的記錄方式（如利用matlab的圖形功能），分析三種迭代法的收斂性與初值選取的關系。（2）對三個迭代法中的某個，取不同的初始值進行迭代，結果如何？試分析迭代法對不同的初值是否有差異？（3）線性方程組迭代法的收斂性是不依賴初始值選取的。比較線性與非線性問題迭代的差異，有何結論和問題。實驗過程：7.1程序：7.1.1 對于第一個迭代方程的程序：保存為：diedai71clearclca=-1.5;b=2.5; y00=0; x00=input(請輸入第一個函數(shù)的初值：x00=);x=linspace(a,b,80);y0=x; %計算直線y=xy1=

43、diedai7f1(x); %計算迭代函數(shù)y=f(x)clear y;y=y0;y1;plot(x,y,linewidth,1)legend(y=x,y=f1)title(x(n+1)=x(n)2-1) %輸出標題hold onplot(a b,0,0,k-,0 0,a b,k-)axis(a,b,a,b) %畫坐標軸z=;for i=1:15 %畫蛛網(wǎng)圖，迭代過程為n=15次 xt(1)=x00;yt(1)=y00; %決定始點坐標 xt(2)=diedai7f1(xt(1); %決定終點坐標 yt(2)=diedai7f1(xt(1); diedaiplot71(xt,yt,0.6); %

44、畫蛛網(wǎng)圖 if i=5 pause %按任意鍵逐次觀察前5次迭代的蛛網(wǎng)圖 end x00=xt(2);y00=yt(2); %將本次迭代的終點作為下次的始點 z=z,xt(1); %保存迭代點end保存為：diedai7f1function y=diedai7f1(x)y=(x.*x-1);保存為：diedaiplot72function out=diedaiplot72(x,y,p)%畫一次迭代的蛛網(wǎng)圖，改變p調節(jié)箭頭的大小u(1)=0;v(1)=(y(2)-y(1); %畫出始點(x(1),y(1)終點(x(2),y(2)的有向折線段u(2)=eps;v(2)=eps;x=x(1) x(1

45、);y=y(1) y(2);h=quiver(x,y,u,v,p);set(h,color,red);hold onu(1)=(x(2)-x(1);v(1)=0;u(2)=eps;v(2)=eps;h=quiver(x(1) x(2),y(2) y(2),u,v,p);set(h,color,red);plot(x(1) x(1) x(2),y(1) y(2) y(2),r.-)7.1.1 對于第一個迭代方程的程序運行結果：（注：由于可以估計出方程的根，故選取根附近的值開始迭代）當x=1.5時程序運行如下圖：當x=-1時程序運行如下圖：當x=0.5時程序運行如下圖：當x=-0.5時程序運行如下

46、圖：當x=0時程序運行如下圖：7.2.1 對于第二個迭代方程的程序：保存為diedai72：clearclca=0.1;b=6.5; y00=0; x00=input(請輸入第二個函數(shù)的初值：x00=);x=linspace(a,b,80);y0=x; %計算直線y=xy1=diedai7f2(x); %計算迭代函數(shù)y=f(x)clear y;y=y0;y1;plot(x,y,linewidth,1)legend(y=x,y=f1)title(x(n+1)=x(n)2-1) %輸出標題hold onplot(a b,0,0,k-,0 0,a b,k-)axis(a,b,a,b) %畫坐標軸z=

47、;for i=1:15 %畫蛛網(wǎng)圖，迭代過程為n=15次 xt(1)=x00;yt(1)=y00; %決定始點坐標 xt(2)=diedai7f2(xt(1); %決定終點坐標 yt(2)=diedai7f2(xt(1); diedaiplot72(xt,yt,0.6); %畫蛛網(wǎng)圖 if i=5 pause %按任意鍵逐次觀察前5次迭代的蛛網(wǎng)圖 end x00=xt(2);y00=yt(2); %將本次迭代的終點作為下次的始點 z=z,xt(1); %保存迭代點end保存為迭代7f2：function y=diedai7f2(x)y=(1+1./x);保存為：diedaiplot72function out=diedaiplot72(x,y,p)%畫一次迭代的蛛網(wǎng)圖，改變p調節(jié)箭頭的大小u(1)=0;v(1)=(y(2)-y(1); %畫出始點(x(1),y(1)終點(x(2),y(2)的有向折線段u(2)=eps;v(2)=eps;x=x(1) x(1);y=y(1) y(2);h=quiver(x,y,u,v,p);set(h,color,red);hold onu(1)