數(shù)值分析實(shí)驗(yàn)報(bào)告-2

實(shí)用文檔實(shí)驗(yàn)五解線性方程組的直接方法實(shí)驗(yàn)5.1（主元的選取與算法的穩(wěn)定性）問題提出：Gauss消去法是我們?cè)诰€性代數(shù)中已經(jīng)熟悉的。但由于計(jì)算機(jī)的數(shù)值運(yùn)算是在一個(gè)有限的浮點(diǎn)數(shù)集合上進(jìn)行的，如何才能確保Gauss消去法作為數(shù)值算法的穩(wěn)定性呢？Gauss消去法從理論算法到數(shù)值算法，其關(guān)鍵是主元的選擇。主元的選擇從數(shù)學(xué)理論上看起來(lái)平凡，它卻是數(shù)值分析中十分典型的問題。實(shí)驗(yàn)內(nèi)容：考慮線性方程組編制一個(gè)能自動(dòng)選取主元，又能手動(dòng)選取主元的求解線性方程組的Gauss消去過程。實(shí)驗(yàn)要求：（1）取矩陣，則方程有解。取n=10計(jì)算矩陣的條件數(shù)。讓程序自動(dòng)選取主元，結(jié)果如何？（2）現(xiàn)選擇程序中手動(dòng)選取主元的功能。每步消去過程總選取按模最小或按模盡可能小的元素作為主元，觀察并記錄計(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é)果。思考題一：（Vadermonde矩陣）設(shè)，其中，，（1）對(duì)n=2,5,8，計(jì)算A的條件數(shù)；隨n增大，矩陣性態(tài)如何變化？（2）對(duì)n=5，解方程組Ax=b；設(shè)A的最后一個(gè)元素有擾動(dòng)10-4，再求解Ax=b（3）計(jì)算（2）擾動(dòng)相對(duì)誤差與解的相對(duì)偏差，分析它們與條件數(shù)的關(guān)系。（4）你能由此解釋為什么不用插值函數(shù)存在定理直接求插值函數(shù)而要用拉格朗日或牛頓插值法的原因嗎？相關(guān)MATLAB函數(shù)提示：zeros(m,n)生成m行，n列的零矩陣ones(m,n)生成m行，n列的元素全為1的矩陣eye(n)生成n階單位矩陣rand(m,n)生成m行,n列(0,1)上均勻分布的隨機(jī)矩陣diag(x)返回由向量x的元素構(gòu)成的對(duì)角矩陣tril(A)提取矩陣A的下三角部分生成下三角矩陣triu(A)提取矩陣A的上三角部分生成上三角矩陣rank(A)返回矩陣A的秩det(A)返回方陣A的行列式inv(A)返回可逆方陣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實(shí)驗(yàn)過程：5.1.1程序：functionx=gauss(n,r)n=input('請(qǐng)輸入矩陣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)fori=1:4p=input('條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=')pp=cond(A,p)endpause[m,n]=size(A);nb=n+1;Ab=[Ab]r=input('請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=')fori=1:n-1ifr==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);disp(Ab);pauseendendifr==1i=iip=input('輸入i列所選元素所處的行數(shù)：ip=');Ab([iip],:)=Ab([ipi],:);disp(Ab);pauseendpivot=Ab(i,i);fork=i+1:nAb(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end5.1.2實(shí)驗(yàn)結(jié)果如下：1.按照實(shí)驗(yàn)要求一：取矩陣A的階數(shù):n=10且自動(dòng)選取主元，程序結(jié)果運(yùn)行如下：現(xiàn)選擇程序中手動(dòng)選取主元的功能，觀察并記錄計(jì)算結(jié)果。①選取絕對(duì)值最大的元素為主元：程序運(yùn)行開始如第一問的截圖也是求范數(shù)故這里不在給出。Theansweris：1111111111②選取絕對(duì)值最小的元素為主元：Theansweris：1.0e+003*（INF0.0070.0057-0.0410-0.03030.34300.2577-2.7290-2.04632.7308）⑶取矩陣A的階數(shù):n=20，手動(dòng)選取主元：選取絕對(duì)值最大的元素為主元：Theansweris：11111111111111111111選取絕對(duì)值最小的元素為主元：Theansweris：1.0e+007*（-Inf0.00000.0000-0.0000-0.00000.00000.0000-0.0003-0.00020.00220.0016-0.0175-0.01310.13980.1049-1.1185-0.83898.94786.7109-8.9478）⑷修改程序如下：functionx=gaussong(n,r)n=input('請(qǐng)輸入矩陣A的階數(shù):n=')A=hilb(n)b=A*ones(n,1)fori=1:4p=input('條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=')pp=cond(A,p)endpause[m,n]=size(A);nb=n+1;Ab=[Ab]r=input('請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=')fori=1:n-1ifr==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);disp(Ab);pauseendendifr==1i=iip=input('輸入i列所選元素所處的行數(shù)：ip=');Ab([iip],:)=Ab([ipi],:);disp(Ab);pauseendpivot=Ab(i,i);fork=i+1:nAb(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end①所求范數(shù)為：自動(dòng)輸入結(jié)果為：ans=1.00001.00001.00001.00001.00020.99961.00070.99931.00040.9999②選取絕對(duì)值最大的元素為主元結(jié)果為：Theansweris：NaNNaNNaNNaNNaNInf-Inf-Inf281.3945-283.3708③選取絕對(duì)值最小的元素為主元結(jié)果為：Theansweris：NaNNaNNaN-Inf-5.8976-1.9243-2.0291-4.997223.4548-11.10125.1.3對(duì)實(shí)驗(yàn)結(jié)果進(jìn)行分析：5.1.3.1對(duì)實(shí)驗(yàn)要求一的結(jié)果進(jìn)行分析：對(duì)于Gauss消去法就是用行的初等變換將原線性方程組系數(shù)矩陣轉(zhuǎn)化為簡(jiǎn)單形式，從而進(jìn)行求解，缺點(diǎn)是迭代次數(shù)可能較多，效率不高，且在消去過程中不可以將主元素很小的做除數(shù)，否則將導(dǎo)致其他元素?cái)?shù)量級(jí)的嚴(yán)重增長(zhǎng)和舍入誤差的擴(kuò)散，使得計(jì)算解不可靠。5.1.3.2對(duì)實(shí)驗(yàn)要求二的結(jié)果進(jìn)行分析：通過每次選取最大或最小的主元可以發(fā)現(xiàn)取絕對(duì)值大的元素作為主元比取絕對(duì)值小的元素作為主元時(shí)產(chǎn)生的結(jié)果比較準(zhǔn)確，即選取絕對(duì)值小的主元時(shí)結(jié)果產(chǎn)生了較大的誤差，條件數(shù)越大產(chǎn)生的誤差就越大。所以應(yīng)盡量避免很小的數(shù)作為除數(shù)。5.1.3.3對(duì)實(shí)驗(yàn)要求三的結(jié)果進(jìn)行分析：此要求是對(duì)要求一和要求二的一個(gè)延續(xù)，通過實(shí)驗(yàn)結(jié)果可以看出若采用很小的數(shù)作為主元迭代次數(shù)越多導(dǎo)致的結(jié)果越不可靠，甚至出現(xiàn)錯(cuò)誤。5.1.3.4對(duì)實(shí)驗(yàn)要求四的結(jié)果進(jìn)行分析：對(duì)新矩陣進(jìn)行實(shí)驗(yàn)發(fā)現(xiàn)依然符合上述規(guī)律，可以知道，在進(jìn)行迭代時(shí)主元的選擇與算法的穩(wěn)定性有密切的聯(lián)系選取絕對(duì)值大的元素作為主元比絕對(duì)值小的元素作為主元時(shí)對(duì)結(jié)果產(chǎn)生的誤差較小。條件數(shù)越大對(duì)用gauss消去法解線性方程組時(shí)，對(duì)結(jié)果產(chǎn)生的誤差就越大。5.1.4實(shí)驗(yàn)總結(jié)：在用gauss消去法解線性方程組時(shí)，主元的選取與算法的穩(wěn)定性有密切的聯(lián)系，選取適當(dāng)?shù)闹髟欣诘贸龇€(wěn)定的算法，在算法的過程中，選取絕對(duì)值較大的主元比選取絕對(duì)值較小的主元更有利于算法的穩(wěn)定，選取絕對(duì)值最大的元素作為主元時(shí)，得出的結(jié)果相對(duì)較準(zhǔn)確較穩(wěn)定。條件數(shù)越小，對(duì)用這種方法得出的結(jié)果更準(zhǔn)確。在算除法的過程中要盡量避免使用較小的數(shù)做為除數(shù)，以免發(fā)生結(jié)果數(shù)量級(jí)加大，使大數(shù)吃掉小數(shù)，產(chǎn)生舍入誤差。實(shí)驗(yàn)5.2（線性代數(shù)方程組的性態(tài)與條件數(shù)的估計(jì)）問題提出：理論上，線性代數(shù)方程組的攝動(dòng)滿足矩陣的條件數(shù)確實(shí)是對(duì)矩陣病態(tài)性的刻畫，但在實(shí)際應(yīng)用中直接計(jì)算它顯然不現(xiàn)實(shí)，因?yàn)橛?jì)算通常要比求解方程還困難。實(shí)驗(yàn)內(nèi)容：Matlab中提供有函數(shù)“condest”可以用來(lái)估計(jì)矩陣的條件數(shù)，它給出的是按1-范數(shù)的條件數(shù)。首先構(gòu)造非奇異矩陣A和右端，使得方程是可以精確求解的。再人為地引進(jìn)系數(shù)矩陣和右端的攝動(dòng)，使得充分小。實(shí)驗(yàn)要求：（1）假設(shè)方程Ax=b的解為x，求解方程，以1-范數(shù)，給出的計(jì)算結(jié)果。（2）選擇一系列維數(shù)遞增的矩陣（可以是隨機(jī)生成的），比較函數(shù)“condest”所需機(jī)器時(shí)間的差別.考慮若干逆是已知的矩陣，借助函數(shù)“eig”很容易給出cond2(A)的數(shù)值。將它與函數(shù)“cond(A,2)”所得到的結(jié)果進(jìn)行比較。（3）利用“condest”給出矩陣A條件數(shù)的估計(jì)，針對(duì)（1）中的結(jié)果給出的理論估計(jì)，并將它與（1）給出的計(jì)算結(jié)果進(jìn)行比較，分析所得結(jié)果。注意，如果給出了cond(A)和的估計(jì)，馬上就可以給出的估計(jì)。（4）估計(jì)著名的Hilbert矩陣的條件數(shù)。5.2實(shí)驗(yàn)過程如下：5.2.1.1實(shí)驗(yàn)要求一的程序如下：functionn=jisuan(n)a=fix(100*rand(n))+1x=ones(n,1)b=a*xdata=rand(n)*0.00001datb=rand(n,1)*0.00001A=a+dataB=b+datbx0=get(A,B)x1=norm(x0-x,1)/norm(x,1)functionx=get(A,B)[m,n]=size(A);nb=n+1;AB=[AB];fori=1:n-1pivot=AB(i,i);fork=i+1:nAB(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);fori=n-1:-1:1x(i)=(AB(i,nb)-AB(i,i+1:n)*x(i+1:n))/AB(i,i);End5.2.1.2實(shí)驗(yàn)要求一程序運(yùn)行結(jié)果如下：系數(shù)矩陣a為：7018142952536398247346580289074421962620992338538830595879897467437769加擾動(dòng)后的系數(shù)矩陣A為：70.000004618.000004214.000009929.000003252.000004453.000001363.000005798.00000302.000007947.000009634.000009365.000002180.000007928.000008790.00000447.000007344.000006821.000006196.000000626.000000220.000005099.000004123.000002138.000006353.000006088.000007730.000002159.000007458.000008479.000003789.000000574.000009767.000006443.000002777.000006369.0000058b值為：236309270302367419加擾動(dòng)后的b值為：236.00000451309.00000044270.00000027302.00000313367.00000013419.00000384data的值為：4.610951267E-064.153748604E-069.900825926E-063.200355775E-064.399243096E-061.337727485E-065.678287124E-063.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.445657831E-068.392382403E-063.704768261E-065.02688037E-079.708449393E-066.434922879E-062.679472507E-066.287846E-065.75147779E-06datb的值為：4.514248268E-064.38953253E-072.7185123E-073.126850481E-061.28625747E-073.839672885E-06xx的值為：0.999999684491.00000038780.999999143481.00000065171.00000221560.9999975453x0的值為：1.146958549E-06x1的值為：6.8990e-0075.2.1.3實(shí)驗(yàn)結(jié)果為：的計(jì)算結(jié)果為：6.8990e-0075.2.2.1實(shí)驗(yàn)要求二的程序如下：functioncond2(A)B=A'*A;[V1,D1]=eig(B);[V2,D2]=eig(B^(-1));cond2A=sqrt(max(max(D1)))*sqrt(max(max(D2)))endforn=10:10:100n=nA=fix(100*randn(n));condestA=condest(A)cond2(A)condA2=cond(A,2)pauseend5.2.2.2實(shí)驗(yàn)要求二的程序運(yùn)行結(jié)果如下：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+025.2.3.1實(shí)驗(yàn)要求三的程序如下：functionbijiao(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實(shí)驗(yàn)要求三的程序運(yùn)行結(jié)果如下：給出對(duì)的估計(jì)是：7.310559817408125e-007的理論結(jié)果是：3.828481757617297e-007結(jié)果相差：3.482078059790828e-0075.2.4.1實(shí)驗(yàn)要求四的程序如下：forn=4:11n=nHi=hilb(n);cond1Hi=cond(Hi,1)cond2Hi=cond(Hi,2)condinfHi=cond(Hi,inf)pauseend5.2.4.2實(shí)驗(yàn)要求四的程序運(yùn)行結(jié)果如下：ncond1Hicond2HicondinfHi42.84E+041.55E+042.84E+0459.44E+054.77E+059.44E+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ù)有著很重要的關(guān)系，既矩陣的條件數(shù)是刻畫矩陣性質(zhì)的一個(gè)重要的依據(jù)，條件數(shù)越大，矩陣“病態(tài)”性越嚴(yán)重，在解線性代數(shù)方程組的過程中較容易產(chǎn)生比較大的誤差，則在實(shí)際問題的操作過程中，我們必須要減少對(duì)條件數(shù)來(lái)求解，把條件數(shù)較大的矩陣化成條件數(shù)較小的矩陣來(lái)進(jìn)行求解。實(shí)驗(yàn)總結(jié)：在本次實(shí)驗(yàn)中，使我們知道了矩陣條件數(shù)對(duì)線性代數(shù)方程組求解的影響，條件數(shù)越大，對(duì)最后解的影響的越大，hilbert矩陣是一個(gè)很”病態(tài)”的矩陣，他的條件數(shù)隨著階數(shù)的增加而增大，每增加一階，條件數(shù)就增大一個(gè)數(shù)量級(jí)，在求解的過程中要盡量避免hilbert矩陣實(shí)驗(yàn)六解線性方程組的迭代法實(shí)驗(yàn)6.1（病態(tài)的線性方程組的求解）問題提出：理論的分析表明，求解病態(tài)的線性方程組是困難的。實(shí)際情況是否如此，會(huì)出現(xiàn)怎樣的現(xiàn)象呢？實(shí)驗(yàn)內(nèi)容：考慮方程組Hx=b的求解，其中系數(shù)矩陣H為Hilbert矩陣，這是一個(gè)著名的病態(tài)問題。通過首先給定解（例如取為各個(gè)分量均為1）再計(jì)算出右端b的辦法給出確定的問題。實(shí)驗(yàn)要求：（1）選擇問題的維數(shù)為6，分別用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程組，其各自的結(jié)果如何？將計(jì)算結(jié)果與問題的解比較，結(jié)論如何？（2）逐步增大問題的維數(shù)，仍然用上述的方法來(lái)解它們，計(jì)算的結(jié)果如何？計(jì)算的結(jié)果說明了什么？（3）討論病態(tài)問題求解的算法由于本實(shí)驗(yàn)用到迭代法，故先給出迭代法相關(guān)資料：迭代法也稱輾轉(zhuǎn)法，是一種不斷用變量的舊值遞推新值的過程,跟迭代法相對(duì)應(yīng)的是直接法（或者稱為一次解法），即一次性解決問題。迭代法又分為精確迭代和近似迭代。“二分法”和“牛頓迭代法”屬于近似迭代法。迭代算法是用計(jì)算機(jī)解決問題的一種基本方法。它利用計(jì)算機(jī)運(yùn)算速度快、適合做重復(fù)性操作的特點(diǎn)，讓計(jì)算機(jī)對(duì)一組指令（或一定步驟）進(jìn)行重復(fù)執(zhí)行，在每次執(zhí)行這組指令（或這些步驟）時(shí)，都從變量的原值推出它的一個(gè)新值。迭代是數(shù)值分析中通過從一個(gè)初始估計(jì)出發(fā)尋找一系列近似解來(lái)解決問題（一般是解方程或者方程組）的過程，為實(shí)現(xiàn)這一過程所使用的方法統(tǒng)稱為迭代法（IterativeMethod）。一般可以做如下定義：對(duì)于給定的線性方程組（這里的同為矩陣，任意線性方程組都可以變換成此形式），用公式(括號(hào)中為上標(biāo)，代表迭代k次得到的x，初始時(shí)k=0)逐步帶入求近似解的方法稱為迭代法（或稱一階定常迭代法）。如果k趨向無(wú)窮大時(shí)存在，記為，稱此迭代法收斂。顯然就是此方程組的解，否則稱為迭代法發(fā)散。跟迭代法相對(duì)應(yīng)的是直接法（或者稱為一次解法），即一次性的快速解決問題，例如通過開方解決方程。一般如果可能，直接解法總是優(yōu)先考慮的。但當(dāng)遇到復(fù)雜問題時(shí)，特別是在未知量很多，方程為非線性時(shí)，我們無(wú)法找到直接解法（例如五次以及更高次的代數(shù)方程沒有解析解，參見阿貝爾定理），這時(shí)候或許可以通過迭代法尋求方程（組）的近似解。最常見的迭代法是牛頓法。其他還包括最速下降法、共軛迭代法、變尺度迭代法、最小二乘法、線性規(guī)劃、非線性規(guī)劃、單純型法、懲罰函數(shù)法、斜率投影法、遺傳算法、模擬退火法等等。利用迭代算法解決問題，需要做好以下三個(gè)方面的工作：（1）確定迭代變量在可以用迭代算法解決的問題中，至少存在一個(gè)直接或間接地不斷由舊值遞推出新值的變量，這個(gè)變量就是迭代變量。（2）建立迭代關(guān)系式所謂迭代關(guān)系式，指如何從變量的前一個(gè)值推出其下一個(gè)值的公式（或關(guān)系）。迭代關(guān)系式的建立是解決迭代問題的關(guān)鍵，通?？梢皂樛苹虻雇频姆椒▉?lái)完成。（3）對(duì)迭代過程進(jìn)行控制在什么時(shí)候結(jié)束迭代過程？這是編寫迭代程序必須考慮的問題。不能讓迭代過程無(wú)休止地重復(fù)執(zhí)行下去。迭代過程的控制通常可分為兩種情況：一種是所需的迭代次數(shù)是個(gè)確定的值，可以計(jì)算出來(lái)；另一種是所需的迭代次數(shù)無(wú)法確定。對(duì)于前一種情況，可以構(gòu)建一個(gè)固定次數(shù)的循環(huán)來(lái)實(shí)現(xiàn)對(duì)迭代過程的控制；對(duì)于后一種情況，需要進(jìn)一步分析出用來(lái)結(jié)束迭代過程的條件。實(shí)驗(yàn)過程：6.1.1.1實(shí)驗(yàn)要求一的gauss實(shí)現(xiàn)程序如下：Gauss消去法程序：functionx=gaussong(n,r)A=hilb(n)b=A*ones(n,1)fori=1:4p=input('條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=')pp=cond(A,p)endpause[m,n]=size(A);nb=n+1;Ab=[Ab]r=input('請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=')fori=1:n-1ifr==0[pivot,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);disp(Ab);pauseendendifr==1i=iip=input('輸入i列所選元素所處的行數(shù)：ip=');Ab([iip],:)=Ab([ipi],:);disp(Ab);pauseendpivot=Ab(i,i);fork=i+1:nAb(k,i:nb)=Ab(k,i:nb)-(Ab(k,i)/pivot)*Ab(i,i:nb);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,nb)/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,nb)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);end6.1.1.2實(shí)驗(yàn)要求一的gauss實(shí)現(xiàn)程序運(yùn)行結(jié)果如下：6.1.2.1實(shí)驗(yàn)要求一的J程序如下：functionn=JD(n)A=hilb(n);fori=1:n;a0(i)=1;x(i)=0;endb=A*a0';fori=1:100;y=x;forj=1:n;x(j)=b(j)/A(j,j);fork=1:j-1;x(j)=x(j)-A(j,k)*y(k)/A(j,j);endfork=j+1:n;x(j)=x(j)-A(j,k)*y(k)/A(j,j);endendendx6.1.2.2實(shí)驗(yàn)要求一的J程序運(yùn)行結(jié)果如下：6.1.3.1實(shí)驗(yàn)要求一的GS程序如下：functionn=GS(n)A2=hilb(n);fori=1:n;a02(i)=1;x2(i)=0;endb2=A2*a02';fori=1:100000;forj=1:n;x2(j)=b2(j)/A2(j,j);fork=1:j-1;x2(j)=x2(j)-A2(j,k)*x2(k)/A2(j,j);endfork=j+1:n;x2(j)=x2(j)-A2(j,k)*x2(k)/A2(j,j);endendendx26.1.3.2實(shí)驗(yàn)要求一的GS程序運(yùn)行結(jié)果如下：6.1.4.1實(shí)驗(yàn)要求一的SOR程序如下：function[nss]=SOR(n,ss)A3=hilb(n);fori=1:n;a03(i)=1;x3(i)=0;endb3=A3*a03';fori=1:100000;forj=1:n;rc=x3(j);x3(j)=b3(j)/A3(j,j);fork=1:j-1;x3(j)=x3(j)-A3(j,k)*x3(k)/A3(j,j);endfork=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實(shí)驗(yàn)要求一的SOR程序運(yùn)行結(jié)果如下：（注意由于運(yùn)行結(jié)果迭代過程很長(zhǎng)故不體現(xiàn)在實(shí)驗(yàn)報(bào)告中）6.2.1實(shí)驗(yàn)要求二的gauss程序運(yùn)行結(jié)果如下：x=00000000000000000000.99756.2.2實(shí)驗(yàn)要求二的J程序運(yùn)行結(jié)果如下：6.2.3實(shí)驗(yàn)要求二的GS程序運(yùn)行結(jié)果如下：6.2.4實(shí)驗(yàn)要求二的SOR程序運(yùn)行結(jié)果如下：實(shí)驗(yàn)要求二的結(jié)果分析：選擇問題的維數(shù)為20時(shí)：1.用Gauss消去法求得的解與精確解相差很大，說明能否得到優(yōu)秀的解取決于算法的穩(wěn)定性，如果算法不夠穩(wěn)定產(chǎn)生的結(jié)果將變的無(wú)法理喻。2.取初始向量為0，用J迭代方法迭代發(fā)散，無(wú)法求解；3.取初始向量為0，用GS迭代方法迭代不發(fā)散，能求得解，但收斂非常緩慢，而且迭代次數(shù)越多，與準(zhǔn)確解的偏差就越大，說明GS迭代適合迭代次數(shù)少的，但是通常我們無(wú)法得知需要迭代的次數(shù)。4.取初始向量為0，用SOR迭代方法迭代不發(fā)散，能求得解，但同樣收斂非常緩慢?？傊瑥纳厦娴慕Y(jié)果可以看出當(dāng)病態(tài)問題的階數(shù)升高時(shí)作為直接法的Gauss消去法又能求解變成不能求解。而GS和SOR迭代法在階數(shù)升高時(shí)仍能求解。但在階數(shù)較低時(shí)直接法能求得精確解而迭代發(fā)卻總存在一定的誤差?？梢娭苯臃ㄅc迭代法各有各的優(yōu)勢(shì)與不足。關(guān)于病態(tài)問題的求解，主要的方法是對(duì)原方程作某些預(yù)處理，以降低系數(shù)矩陣的條件數(shù)?？梢圆扇⑾禂?shù)矩陣A的每一行本別乘上適當(dāng)常數(shù)的方法。即找到可逆的對(duì)角陣和使方程組化為理論上最好選擇對(duì)角陣滿足：。補(bǔ)充：本實(shí)驗(yàn)用到了Gauss消去法、J迭代法、GS迭代法、SOR迭代，故對(duì)其重新理解學(xué)習(xí)了一下：Gauss消去法：1將其增廣矩陣化為行階梯形2若最后有形如[00...0a]（a<>0)的行則無(wú)解3若含有自由變量則有無(wú)窮組解4原方程有唯一解。采用回代求解。至于有無(wú)窮組解的方程組的求解，需將其化為行最簡(jiǎn)形矩陣，其方法稱為高斯-若爾當(dāng)消元法。J迭代法：一、算法理論 迭代格式的引出是依據(jù)迭代法的基本思想：構(gòu)造一個(gè)向量系列，使其收斂至某個(gè)極限，則就是要求的方程組的準(zhǔn)確解。迭代將方程組：在假設(shè)，改寫成如果引用系數(shù)矩陣，及向量方程組（1）和（2）分別可寫為：及，這樣就得到了迭代格式用迭代解方程組時(shí)，就可任意取初值帶入迭代可知式，然后求。但是，比較大的時(shí)候，寫方程組（1）和（2）是很麻煩的，如果直接由，能直接得到，就是矩陣與向量的運(yùn)算了，那么如何得到，呢 實(shí)際上，如果引進(jìn)非奇異對(duì)角矩陣將分解成：要求的解，實(shí)質(zhì)上就有，而是非奇異的，所以存在，，從而有我們?cè)谶@里不妨令，就得到迭代格式：現(xiàn)在考慮迭代法的計(jì)算程序floata[3][3]={{10,-2,-1},{-2,10,-1},{-1,-2,-5}};floatb[3]={3,15,10};分別代表的系數(shù)和等號(hào)右邊的常數(shù)項(xiàng),即先輸入方程，運(yùn)行main函數(shù)，如果first不為null，則執(zhí)行if括號(hào)里的，否則執(zhí)行else里面的，最后會(huì)調(diào)用方法sum（）。在sum（）中sum=a[m][n]*x[n]+sum;y=(b[m]-sum)/a[m][m]，之后運(yùn)行for循環(huán)，最后輸出結(jié)果，算法結(jié)束。二、算法框圖調(diào)用方法開始讀入調(diào)用方法開始讀入輸出輸出結(jié)束結(jié)束GS迭代法：1.高斯-塞德爾迭代法公式的矩陣形式首先將高斯-塞德爾迭代法的公式表示為矩陣形式，為此設(shè)這里是系數(shù)矩陣的對(duì)角部分,是嚴(yán)格下三角部分，是嚴(yán)格上三角部分，則高斯-塞德爾迭代法的公式可表示為用矩陣乘等式兩邊得再用矩陣乘等式兩邊得其中矩陣稱為高斯—塞德爾迭代矩陣。由此可見，高斯-塞德爾迭代法是一般迭代法中迭代矩陣為的特殊情形。需要指出的是，由于矩陣難于計(jì)算，所以式（2）多用在理論分析中。高斯—塞德爾迭代法計(jì)算框圖（見圖）高斯—塞德爾迭代法計(jì)算框圖松弛法：其中相當(dāng)于在的基礎(chǔ)上加個(gè)余項(xiàng)生成，下面令，希望通過選取來(lái)加速收斂，這就是松弛法。實(shí)驗(yàn)七非線性方程求根實(shí)驗(yàn)7.1（迭代法、初始值與收斂性）實(shí)驗(yàn)?zāi)康模撼醪秸J(rèn)識(shí)非線性問題的迭代法與線性問題迭代法的差別，探討迭代法及初始值與迭代收斂性的關(guān)系。問題提出：迭代法是求解非線性方程的基本思想方法，與線性方程的情況一樣，其構(gòu)造方法可以有多種多樣，但關(guān)鍵是怎樣才能使迭代收斂且有較快的收斂速度。實(shí)驗(yàn)內(nèi)容：考慮一個(gè)簡(jiǎn)單的代數(shù)方程針對(duì)上述方程，可以構(gòu)造多種迭代法，如在實(shí)軸上取初始值x0，請(qǐng)分別用迭代（7.1）-（7.3）作實(shí)驗(yàn)，記錄各算法的迭代過程。實(shí)驗(yàn)要求：（1）取定某個(gè)初始值，分別計(jì)算（7.1）-（7.3）迭代結(jié)果，它們的收斂性如何？重復(fù)選取不同的初始值，反復(fù)實(shí)驗(yàn)。請(qǐng)自選設(shè)計(jì)一種比較形象的記錄方式（如利用Matlab的圖形功能），分析三種迭代法的收斂性與初值選取的關(guān)系。（2）對(duì)三個(gè)迭代法中的某個(gè)，取不同的初始值進(jìn)行迭代，結(jié)果如何？試分析迭代法對(duì)不同的初值是否有差異？（3）線性方程組迭代法的收斂性是不依賴初始值選取的。比較線性與非線性問題迭代的差異，有何結(jié)論和問題。實(shí)驗(yàn)過程：7.1程序：7.1.1對(duì)于第一個(gè)迭代方程的程序：保存為：diedai71clearclca=-1.5;b=2.5;y00=0;x00=input('請(qǐng)輸入第一個(gè)函數(shù)的初值：x00=');x=linspace(a,b,80);y0=x;%計(jì)算直線y=xy1=diedai7f1(x);%計(jì)算迭代函數(shù)y=f(x)cleary;y=[y0;y1];plot(x,y,'linewidth',1)legend('y=x','y=f1')title('x(n+1)=[x(n)]^2-1')%輸出標(biāo)題holdonplot([ab],[0,0],'k-',[00],[ab],'k-')axis([a,b,a,b])%畫坐標(biāo)軸z=[];fori=1:15%畫蛛網(wǎng)圖，迭代過程為n=15次xt(1)=x00;yt(1)=y00;%決定始點(diǎn)坐標(biāo)xt(2)=diedai7f1(xt(1));%決定終點(diǎn)坐標(biāo)yt(2)=diedai7f1(xt(1));diedaiplot71(xt,yt,0.6);%畫蛛網(wǎng)圖ifi<=5pause%按任意鍵逐次觀察前5次迭代的蛛網(wǎng)圖endx00=xt(2);y00=yt(2);%將本次迭代的終點(diǎn)作為下次的始點(diǎn)z=[z,xt(1)];%保存迭代點(diǎn)end保存為：diedai7f1functiony=diedai7f1(x)y=(x.*x-1);保存為：diedaiplot72functionout=diedaiplot72(x,y,p)%畫一次迭代的蛛網(wǎng)圖，改變p調(diào)節(jié)箭頭的大小u(1)=0;v(1)=(y(2)-y(1));%畫出始點(diǎn)(x(1),y(1))終點(diǎn)(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');holdonu(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對(duì)于第一個(gè)迭代方程的程序運(yùn)行結(jié)果：（注：由于可以估計(jì)出方程的根，故選取根附近的值開始迭代）當(dāng)x=1.5時(shí)程序運(yùn)行如下圖：當(dāng)x=-1時(shí)程序運(yùn)行如下圖：當(dāng)x=0.5時(shí)程序運(yùn)行如下圖：當(dāng)x=-0.5時(shí)程序運(yùn)行如下圖：當(dāng)x=0時(shí)程序運(yùn)行如下圖：7.2.1對(duì)于第二個(gè)迭代方程的程序：保存為diedai72：clearclca=0.1;b=6.5;y00=0;x00=input('請(qǐng)輸入第二個(gè)函數(shù)的初值：x00=');x=linspace(a,b,80);y0=x;%計(jì)算直線y=xy1=diedai7f2(x);%計(jì)算迭代函數(shù)y=f(x)cleary;y=[y0;y1];plot(x,y,'linewidth',1)legend('y=x','y=f1')title('x(n+1)=[x(n)]^2-1')%輸出標(biāo)題holdonplot([ab],[0,0],'k-',[00],[ab],'k-')axis([a,b,a,b])%畫坐標(biāo)軸z=[];fori=1:15%畫蛛網(wǎng)圖，迭代過程為n=15次xt(1)=x00;yt(1)=y00;%決定始點(diǎn)坐標(biāo)xt(2)=diedai7f2(xt(1));%決定終點(diǎn)坐標(biāo)yt(2)=diedai7f2(xt(1));diedaiplot72(xt,yt,0.6);%畫蛛網(wǎng)圖ifi<=5pause%按任意鍵逐次觀察前5次迭代的蛛網(wǎng)圖endx00=xt(2);y00=yt(2);%將本次迭代的終點(diǎn)作為下次的始點(diǎn)z=[z,xt(1)];%保存迭代點(diǎn)end保存為迭代7f2：functiony=diedai7f2(x)y=(1+1./x);保存為：diedaiplot72functionout=diedaiplot72(x,y,p)%畫一次迭代的蛛網(wǎng)圖，改變p調(diào)節(jié)箭頭的大小u(1)=0;v(1)=(y(2)-y(1));%畫出始點(diǎn)(x(1),y(1))終點(diǎn)(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');holdonu(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.2.2對(duì)于第二個(gè)迭代方程的程序運(yùn)行結(jié)果：（注：由于可以估計(jì)出方程的根，故選取根附近的值開始迭代）當(dāng)x=1.5時(shí)程序運(yùn)行如下圖：當(dāng)x=5時(shí)程序運(yùn)行如下圖：當(dāng)x=2.5時(shí)程序運(yùn)行如下圖：當(dāng)x=0.5時(shí)程序運(yùn)行如下圖：7.3.1對(duì)于第三個(gè)迭代方程的程序：保存為diedai73：clearclca=0;b=2;y00=0;x00=input('請(qǐng)輸入第三個(gè)函數(shù)的初值：x00=');x=linspace(a,b,80);y0=x;%計(jì)算直線y=xy1=diedai7f3(x);%計(jì)算迭代函數(shù)y=f(x)cleary;y=[y0;y1];plot(x,y,'linewidth',1)legend('y=x','y=f1')title('x(n+1)=[x(n)]^2-1')%輸出標(biāo)題holdonplot([ab],[0,0],'k-',[00],[ab],'k-')a