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

數(shù)值分析上機(jī)實(shí)驗(yàn)報(bào)告學(xué)院：專業(yè)：班級(jí)：學(xué)號(hào)：學(xué)生姓名：指導(dǎo)教師：實(shí)驗(yàn)五解線性方程組的直接方法實(shí)驗(yàn)5.1〔主元的選取與算法的穩(wěn)定性〕問(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ù)值分析中十分典型的問(wèn)題。實(shí)驗(yàn)內(nèi)容：考慮線性方程組編制一個(gè)能自動(dòng)選取主元，又能手動(dòng)選取主元的求解線性方程組的Gauss消去過(guò)程。實(shí)驗(yàn)要求：〔1〕取矩陣，那么方程有解。取n=10計(jì)算矩陣的條件數(shù)。讓程序自動(dòng)選取主元，結(jié)果如何？〔2〕現(xiàn)選擇程序中手動(dòng)選取主元的功能。每步消去過(guò)程總選取按模最小或按模盡可能小的元素作為主元，觀察并記錄計(jì)算結(jié)果。假設(shè)每步消去過(guò)程總選取按模最大的元素作為主元，結(jié)果又如何？分析實(shí)驗(yàn)的結(jié)果?！?〕取矩陣階數(shù)n=20或者更大，重復(fù)上述實(shí)驗(yàn)過(guò)程，觀察記錄并分析不同的問(wèn)題及消去過(guò)程中選擇不同的主元時(shí)計(jì)算結(jié)果的差異，說(shuō)明主元素的選取在消去過(guò)程中的作用?！?〕選取其他你感興趣的問(wèn)題或者隨機(jī)生成矩陣，計(jì)算其條件數(shù)。重復(fù)上述實(shí)驗(yàn)，觀察記錄并分析實(shí)驗(yàn)結(jié)果。實(shí)驗(yàn)過(guò)程：程序：建立M文件：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)p=input('條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=')pp=cond(A,p)pause[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ù)值實(shí)驗(yàn)結(jié)果及分析：⑴取矩陣A的階數(shù):n=10，自動(dòng)選取主元：>>formatlong>>gauss請(qǐng)輸入矩陣A的階數(shù):n=10n=10條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=1p=1pp=2.557500000000000e+003請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=0r=0⑵取矩陣A的階數(shù):n=10，手動(dòng)選取主元：①選取絕對(duì)值最大的元素為主元：>>gauss請(qǐng)輸入矩陣A的階數(shù):n=10n=10條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=2p=2pp=1.727556024913903e+003請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=1r=1ans=1111111111②選取絕對(duì)值最小的元素為主元：>>gauss請(qǐng)輸入矩陣A的階數(shù):n=10n=10條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=2p=2pp=1.727556024913903e+003請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=1r=1ans=1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000000.999999999999991.000000000000010.999999999999981.00000000000003⑶取矩陣A的階數(shù):n=20，手動(dòng)選取主元：選取絕對(duì)值最大的元素為主元：>>gauss請(qǐng)輸入矩陣A的階數(shù):n=20條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=1p=1ans=11111111111111111111選取絕對(duì)值最小的元素為主元：>>gauss請(qǐng)輸入矩陣A的階數(shù):n=20.n=20條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=2p=2pp=1.789670565881683e+006請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=1r=1ans=1.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000001.000000000000010.999999999999971.000000000000060.999999999999891.000000000000230.999999999999551.000000000000900.999999999998211.000000000003520.999999999993181.000000000012730.999999999978171.00000000002910⑷將M文件中的第三行：A=diag(6*ones(1,n))+diag(ones(1,n-1),1)+diag(8*ones(1,n-1),-1)改為：A=hilb(n)①>>gauss請(qǐng)輸入矩陣A的階數(shù):n=7n=7條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=1p=1請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=1r=1ans=1.000000000000510.999999999972511.000000000313540.999999998641331.000000002688050.999999997541811.00000000084337②>>gauss請(qǐng)輸入矩陣A的階數(shù):n=7n=7條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=2p=2pp=4.753673569067072e+008請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=1r=1ans=0.999999999998691.000000000043370.999999999642991.000000001211430.999999998030381.000000001528250.99999999954491該問(wèn)題在主元選取與算出結(jié)果有著很大的關(guān)系，取絕對(duì)值大的元素作為主元比取絕對(duì)值小的元素作為主元時(shí)產(chǎn)生的結(jié)果比擬準(zhǔn)確，即選取絕對(duì)值小的主元時(shí)結(jié)果產(chǎn)生了較大的誤差，條件數(shù)越大產(chǎn)生的誤差就越大。討論：在gauss消去法解線性方程組時(shí)，主元的選擇與算法的穩(wěn)定性有密切的聯(lián)系，選取絕對(duì)值大的元素作為主元比絕對(duì)值小的元素作為主元時(shí)對(duì)結(jié)果產(chǎn)生的誤差較小。條件數(shù)越大對(duì)用gauss消去法解線性方程組時(shí)，對(duì)結(jié)果產(chǎn)生的誤差就越大。實(shí)驗(yàn)總結(jié)：對(duì)用gauss消去法解線性方程組時(shí)，主元的選取與算法的穩(wěn)定性有密切的聯(lián)系，選取適當(dāng)?shù)闹髟欣诘贸龇€(wěn)定的算法，在算法的過(guò)程中，選取絕對(duì)值較大的主元比選取絕對(duì)值較小的主元更有利于算法的穩(wěn)定，選取絕對(duì)值最大的元素作為主元時(shí)，得出的結(jié)果相對(duì)較準(zhǔn)確較穩(wěn)定。條件數(shù)越小，對(duì)用這種方法得出的結(jié)果更準(zhǔn)確。在算除法的過(guò)程中要盡量防止使用較小的數(shù)做為除數(shù)，以免發(fā)生結(jié)果數(shù)量級(jí)加大，使大數(shù)吃掉小數(shù)，產(chǎn)生舍入誤差。實(shí)驗(yàn)5.2〔線性代數(shù)方程組的性態(tài)與條件數(shù)的估計(jì)〕問(wèn)題提出：理論上，線性代數(shù)方程組的攝動(dòng)滿足矩陣的條件數(shù)確實(shí)是對(duì)矩陣病態(tài)性的刻畫(huà)，但在實(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é)果?！?〕選擇一系列維數(shù)遞增的矩陣〔可以是隨機(jī)生成的〕，比擬函數(shù)“condest”所需機(jī)器時(shí)間的差異.考慮假設(shè)干逆是的矩陣，借助函數(shù)“eig”很容易給出cond2(A)的數(shù)值。將它與函數(shù)“cond(A,2)”所得到的結(jié)果進(jìn)行比擬?！?〕利用“condest”給出矩陣A條件數(shù)的估計(jì)，針對(duì)〔1〕中的結(jié)果給出的理論估計(jì)，并將它與〔1〕給出的計(jì)算結(jié)果進(jìn)行比擬，分析所得結(jié)果。注意，如果給出了cond(A)和的估計(jì)，馬上就可以給出的估計(jì)?！?〕估計(jì)著名的Hilbert矩陣的條件數(shù)。實(shí)驗(yàn)過(guò)程：程序：⑴n=input('pleaseinputn:n=')%輸入矩陣的階數(shù)a=fix(100*rand(n))+1%隨機(jī)生成一個(gè)矩陣ax=ones(n,1)%假設(shè)知道方程組的解全為1b=a*x%用矩陣a和以知解得出矩陣bdata=rand(n)*0.00001%隨即生成擾動(dòng)矩陣datadatb=rand(n,1)*0.00001%隨即生成擾動(dòng)矩陣datbA=a+dataB=b+datbxx=geshow(A,B)%解擾動(dòng)后的解x0=norm(xx-x,1)/norm(x,1)%得出的理論結(jié)果保存為：fanshu.mfunctionx=geshow(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);end保存為：geshow.m⑵functioncond2(A)%自定義求二階條件數(shù)B=A'*A;[V1,D1]=eig(B);[V2,D2]=eig(B^(-1));cond2A=sqrt(max(max(D1)))*sqrt(max(max(D2)))end保存為：cond2.mformatlongforn=10:10:100n=n%n為矩陣的階A=fix(100*randn(n));%隨機(jī)生成矩陣AcondestA=condest(A)%用condest求條件數(shù)cond2(A)%用自定義的求條件數(shù)condA2=cond(A,2)%用cond求條件數(shù)pause%運(yùn)行一次暫停end保存為：shiyan52.m⑶n=input('pleaseinputn:n=')%輸入矩陣的階數(shù)a=fix(100*rand(n))+1;%隨機(jī)生成一個(gè)矩陣ax=ones(n,1);%假設(shè)知道方程組的解全為1b=a*x;%用矩陣a和以知解得出矩陣bdata=rand(n)*0.00001;%隨即生成擾動(dòng)矩陣datadatb=rand(n,1)*0.00001;%隨即生成擾動(dòng)矩陣datbA=a+data;B=b+datb;xx=geshow(A,B);%利用第一小問(wèn)的geshow.m求出解陣x0=norm(xx-x,1)/norm(x,1)%得出的理論結(jié)果x00=cond(A)/(1-norm(inv(A))*norm(xx-x))*(norm((xx-x))/(norm(A))+norm(datb)/norm(B))%得出的估計(jì)值datx=abs(x0-x00)%求兩者之間的誤差保存為：sy5_2.m⑷f(wàn)ormatlongforn=4:11n=n%n為矩陣的階數(shù)Hi=hilb(n);%生成Hilbert矩陣cond1Hi=cond(Hi,1)%求Hilbert矩陣得三種條件數(shù)cond2Hi=cond(Hi,2)condinfHi=cond(Hi,inf)pauseend數(shù)值實(shí)驗(yàn)結(jié)果及分析：⑴>>fanshupleaseinputn:n=6n=6a=142516881989329385489260144088501316235219292324010100737241437227701x=111111b=251410221157218187data=1.0e-005*datb=1.0e-005*0.335632943352170.275099821466210.04452752039203A=1.0e+002*0.320000064986810.930000023756510.850000054592420.480000097047240.920000035801710.60000002215793B=1.0e+002*xx=0.999998307797201.000000225695551.000000193415550.999999093880730.999999968940211.00000066032794x0=的計(jì)算結(jié)果為：〔2〕NcondestAcond2AcondA210e+00232.8905456307542132.89054563075420203.470959631940668e+002306.050503865112835e+002e+002e+002403.549487892582470e+00261.3753756968344861.3753756968336550e+002601.082004656409367e+0041.704830815154781e+0031.704830815108527e+003703.234679145192132e+0033.878481155980936e+0023.878481155978439e+00280e+002902.063634143407935e+003e+002e+0021001.536592818758897e+003e+002e+002⑶>>sy5_2pleaseinputn:n=8n=8x0=1.095033343195828e-006x00=1.705456352162135e-005datx=1.595953017842553e-005給出對(duì)的估計(jì)是：1.705456352162135e-005的理論結(jié)果是：1.095033343195828e-006結(jié)果相差：1.595953017842553e-005(4)ncond1Hicond2HicondinfHi42.837499999999738e+004e+0042.837499999999739e+00459.436559999999364e+0054.766072502414135e+0059.436559999999336e+00562.907027900294878e+007e+0072.907027900294064e+0077e+0084.753673565864586e+008e+00883.387279082022742e+0101.525757545841988e+0103.387279081949470e+01091.099650993366047e+012e+0111.099650991701052e+012103.535372424347474e+0131.602528637652488e+0133.535372455375642e+013111.230369955362001e+0155.223946340715823e+0141.230369938308720e+015討論：線性代數(shù)方程組的性態(tài)與條件數(shù)有著很重要的關(guān)系，既矩陣的條件數(shù)是刻畫(huà)矩陣性質(zhì)的一個(gè)重要的依據(jù)，條件數(shù)越大，矩陣“病態(tài)”性越嚴(yán)重，在解線性代數(shù)方程組的過(guò)程中較容易產(chǎn)生比擬大的誤差，那么在實(shí)際問(wèn)題的操作過(guò)程中，我們必須要減少對(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í)，在求解的過(guò)程中要盡量防止hilbert矩陣思考題一：〔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)系?！?〕你能由此解釋為什么不用插值函數(shù)存在定理直接求插值函數(shù)而要用拉格朗日或牛頓插值法的原因嗎？實(shí)驗(yàn)六解線性方程組的迭代法實(shí)驗(yàn)6.1〔病態(tài)的線性方程組的求解〕問(wèn)題提出：理論的分析說(shuō)明，求解病態(tài)的線性方程組是困難的。實(shí)際情況是否如此，會(huì)出現(xiàn)怎樣的現(xiàn)象呢？實(shí)驗(yàn)內(nèi)容：考慮方程組Hx=b的求解，其中系數(shù)矩陣H為Hilbert矩陣，這是一個(gè)著名的病態(tài)問(wèn)題。通過(guò)首先給定解〔例如取為各個(gè)分量均為1〕再計(jì)算出右端b的方法給出確定的問(wèn)題。實(shí)驗(yàn)要求：〔1〕選擇問(wèn)題的維數(shù)為6，分別用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程組，其各自的結(jié)果如何？將計(jì)算結(jié)果與問(wèn)題的解比擬，結(jié)論如何？〔2〕逐步增大問(wèn)題的維數(shù)，仍然用上述的方法來(lái)解它們，計(jì)算的結(jié)果如何？計(jì)算的結(jié)果說(shuō)明了什么？〔3〕討論病態(tài)問(wèn)題求解的算法第〔1〕問(wèn)：選擇問(wèn)題的維數(shù)為6，分別用Gauss消去法、J迭代法、GS迭代法和SOR迭代法求解方程組，其各自的結(jié)果如何？將計(jì)算結(jié)果與問(wèn)題的解比擬，結(jié)論如何？源程序：Gauss消去法functionx=gauss(n)formatshortdisp('請(qǐng)輸入當(dāng)前你要計(jì)算的Hilbert矩陣的階數(shù):')n=input('');disp('構(gòu)造出的Hilbert矩陣為:')A=hilb(n)d=diag((diag(A))')*eye(n,n);b=A*ones(n,1);[m,n]=size(A);n1=n+1;disp('該線性方程組的增廣矩陣為:')Ab=[Ab]fori=1:n-1[zhuyuan,p]=max(abs(Ab(i:n,i)));ip=p+i-1;ifip~=iAb([iip],:)=Ab([ipi],:);endzhuyuan=Ab(i,i);fork=i+1:nAb(k,i:n1)=Ab(k,i:n1)-(Ab(k,i)/zhuyuan)*Ab(i,i:n1);enddisp(Ab);pauseendx=zeros(n,1);x(n)=Ab(n,n1)/Ab(n,n);fori=n-1:-1:1x(i)=(Ab(i,n1)-Ab(i,i+1:n)*x(i+1:n))/Ab(i,i);endJ迭代法formatshortdisp('請(qǐng)輸入當(dāng)前你要計(jì)算的Hilbert矩陣的階數(shù):')n=input('');disp('構(gòu)造出的Hilbert矩陣為:')A=hilb(n)d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);x=[0;0;0;0;0;0];fori=1:1000x1=inv(d)*(l+u)*x+inv(d)*b;x=x1;enddisp('采用J迭代法結(jié)果為:')x1GS迭代法formatshortdisp('請(qǐng)輸入當(dāng)前你要計(jì)算的Hilbert矩陣的階數(shù):')n=input('');disp('構(gòu)造出的Hilbert矩陣為:')A=hilb(n)d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);x=[0;0;0;0;0;0];fori=1:1000x1=inv(d-l)*u*x+inv(d-l)*b;x=x1;enddisp('采用GS迭代法結(jié)果為:')x1SOR迭代法formatshortdisp('請(qǐng)輸入你當(dāng)前要計(jì)算的Hilbert矩陣的階數(shù):')n=input('');disp('構(gòu)造出的Hilbert矩陣為:')A=hilb(n)d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);disp(請(qǐng)輸入松弛因子w:')w=input('');x=[0;0;0;0;0;0];fori=1:1000x1=inv(d-w*l)*((1-w)*d+w*u)*x+inv(d-w*l)*w*b;x=x1;enddisp('采用SOR迭代法結(jié)果為:')x1實(shí)驗(yàn)結(jié)果：Gauss消去法請(qǐng)輸入當(dāng)前你要計(jì)算的Hilbert矩陣的階數(shù):6構(gòu)造出的Hilbert矩陣為:A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909該線性方程組的增廣矩陣為:Ab=1.00000.50000.33330.25000.20000.16672.45000.50000.33330.25000.20000.16670.14291.59290.33330.25000.20000.16670.14290.12501.21790.25000.20000.16670.14290.12500.11110.99560.20000.16670.14290.12500.11110.10000.84560.16670.14290.12500.11110.10000.09090.73651.00000.50000.33330.25000.20000.16672.450000.08330.08330.07500.06670.05950.367900.08330.08890.08330.07620.06940.401200.07500.08330.08040.07500.06940.383100.06670.07620.07500.07110.06670.355600.05950.06940.06940.06670.06310.32821.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.401200-0.0056-0.0083-0.0095-0.0099-0.0333000.00330.00540.00640.00690.0221000.00510.00830.01020.01110.0347000.00600.00990.01220.01350.04161.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.4012000.00600.00990.01220.01350.0416000-0.0002-0.0004-0.0006-0.0013000-0.0001-0.0003-0.0004-0.00090000.00090.00190.00270.00551.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.4012000.00600.00990.01220.01350.04160000.00090.00190.00270.00550000-0.0000-0.0000-0.00010000-0.0000-0.0001-0.00011.00000.50000.33330.25000.20000.16672.450000.08330.08890.08330.07620.06940.4012000.00600.00990.01220.01350.04160000.00090.00190.00270.00550000-0.0000-0.0001-0.000100000-0.0000-0.0000ans=1.00001.00001.00001.00001.00001.0000J迭代法請(qǐng)輸入當(dāng)前你要計(jì)算的Hilbert矩陣的階數(shù):6構(gòu)造出的Hilbert矩陣為:A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909采用J迭代法結(jié)果為:x1=NaNNaNNaNNaNNaNNaNGS迭代法disp('采用GS迭代法結(jié)果為:')x1請(qǐng)輸入當(dāng)前你要計(jì)算的Hilbert矩陣的階數(shù):6構(gòu)造出的Hilbert矩陣為:A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909采用GS迭代法結(jié)果為:x1=0.99931.01310.95371.03741.02960.9662SOR迭代法請(qǐng)輸入你當(dāng)前要計(jì)算的Hilbert矩陣的階數(shù):6構(gòu)造出的Hilbert矩陣為:A=1.00000.50000.33330.25000.20000.16670.50000.33330.25000.20000.16670.14290.33330.25000.20000.16670.14290.12500.25000.20000.16670.14290.12500.11110.20000.16670.14290.12500.11110.10000.16670.14290.12500.11110.10000.0909請(qǐng)輸入松弛因子w:1.5采用SOR迭代法結(jié)果為:x1=0.99851.02510.90761.10480.99180.9713結(jié)果分析整理以上采用四種方法求出的希爾伯特矩陣的解，如下表所示：GUESS迭代法J迭代法GS迭代法SOR迭代法1NAN0.99930.99851NAN1.01311.02511NAN0.95370.90761NAN1.03741.10481NAN1.02960.99181NAN0.96620.9713從上表可以看到，在希爾伯特矩陣為六階的時(shí)候。采用GUESS迭代法可以得到精確解。在J迭代法時(shí)候迭代法發(fā)散，不能得出解，改良GUESS的GS迭代法得到了近似解，在采用松弛因子為1.5時(shí)，用SOR迭代法也得到了近似解。但是我們必不能就此得出GUESS迭代法最好等結(jié)論。為了得到更加合理的結(jié)果，接下來(lái)繼續(xù)對(duì)不同階矩陣進(jìn)行討論。第〔2〕問(wèn)：逐步增大問(wèn)題的維數(shù)，仍然用上述的方法來(lái)解它們，計(jì)算的結(jié)果如何？計(jì)算的結(jié)果說(shuō)明了什么？源程序：GUESS迭代法依次改變高斯迭代法中矩陣的階數(shù)，分別對(duì)N等于7到12運(yùn)行程序，得到解。J迭代法formatshortforn=7:12A=hilb(n);d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);w=1.2;x=zeros(n,1);fori=1:1000x1=inv(d)*(l+u)*x+inv(d)*b;x=x1;endndisp('采用J迭代法結(jié)果為:')x1endGS迭代法formatshortforn=7:12A=hilb(n);d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);w=1.2;x=zeros(n,1);fori=1:1000x1=inv(d-l)*u*x+inv(d-l)*b;x=x1;endndisp('采用GS迭代法結(jié)果為:')x1endSOR迭代法formatshortforn=7:12A=hilb(n);d=diag((diag(A))')*eye(n,n);l=-(tril(A)-d);u=-(triu(A)-d);b=A*ones(n,1);w=1.2;x=zeros(n,1);fori=1:1000x1=inv(d-w*l)*((1-w)*d+w*u)*x+inv(d-w*l)*w*b;x=x1;endndisp('采用SOR迭代法結(jié)果為:')x1end實(shí)驗(yàn)結(jié)果GUESS迭代法采用GUESS迭代法，分別計(jì)算階數(shù)為7到12階時(shí)，對(duì)應(yīng)的解如下：7891011121.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00001.00000.99991.00001.00001.00001.00020.99981.00051.00001.00001.00000.99961.00070.99811.00001.00001.00001.00070.99841.00431.00001.00000.99931.00220.99401.00001.00040.99821.00520.99991.00080.99740.99981.00061.0000J迭代法采用J迭代法，分別計(jì)算階數(shù)為7到12階時(shí)，對(duì)應(yīng)的解如下：789101112NANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANNANGS迭代法采用GS迭代法，分別計(jì)算階數(shù)為7到12階時(shí)，對(duì)應(yīng)的解如下：7891011120.99810.99720.99670.99680.99730.99821.02521.03291.03441.02971.02021.00780.93210.92720.93820.96120.99031.02051.02381.00240.98130.96530.95620.95361.04851.04581.02991.00850.98710.96881.01651.04111.04581.03741.02211.00450.95451.00351.03081.04091.03931.03080.94790.99451.02261.03621.03970.94580.9891.01661.03230.94590.98531.01190.94650.98230.9467SOR迭代法采用SOR迭代法，分別計(jì)算階數(shù)為7到12階時(shí)，對(duì)應(yīng)的解如下：7891011120.99770.99680.99650.99680.99760.99881.03171.03871.03781.02951.01621.00060.9110.90920.92860.96251.00281.04281.04461.01830.98860.96270.94380.93251.04761.04631.03131.01030.98910.97061.01441.04011.04541.03741.02271.00560.95161.00161.02951.04041.03981.03240.94690.99381.0221.03621.04050.9460.98891.01651.03250.94660.98551.01190.94710.98230.9467結(jié)果分析在以上四種運(yùn)算中，我們可以看到最不穩(wěn)定的是采用J迭代法，在第一問(wèn)中我們看到，當(dāng)希爾伯特矩陣為六階的時(shí)候沒(méi)有得到近似解，隨著矩陣階數(shù)的繼續(xù)增加，依然沒(méi)有得到解。故此在解高度病態(tài)矩陣的時(shí)候我們應(yīng)該盡量防止使用J迭代法。在采用GUESS迭代法的時(shí)候我們看到，針對(duì)特定的希爾伯特矩陣求解，一直到九階以下都能夠得到精確解，但是隨著矩陣階數(shù)的提高，可以看到解出現(xiàn)了不穩(wěn)定的情況，在1的附近出現(xiàn)左右震蕩的情況。且從上表可以看出，隨著矩陣階數(shù)的增加，有震蕩越來(lái)越劇烈的傾向。故此在求解階數(shù)較小的矩陣的時(shí)候，可以優(yōu)先選擇高斯迭代法，但是在階數(shù)較大的矩陣的求解過(guò)程，應(yīng)該盡量選擇更有的數(shù)值解法。GS迭代法和SOR迭代法在階數(shù)逐漸增加的時(shí)候，解依然在一附近波動(dòng)，并且沒(méi)有出現(xiàn)誤差增大的情況。我們可以認(rèn)為這是兩種比擬優(yōu)秀的解高階病態(tài)矩陣的方法。特別值得提出的就是SOR迭代法的收斂性還和迭代因子的選擇有密切的關(guān)系。在實(shí)際計(jì)算過(guò)程中，需要通過(guò)簡(jiǎn)要計(jì)算選擇較好的迭代因子，或者是在MATLAB中編程，在1~2之間選擇微小的步長(zhǎng)，通過(guò)編程比擬求得其最正確的步長(zhǎng)因子，使得其迭代收斂并且收斂速度較快?？傊谇蠼獠煌母呓饩€性方程的時(shí)候，我們應(yīng)該根據(jù)具體情況選擇適宜的迭代法求解矩陣。而迭代矩陣并不僅限于上面提到的四種，我們可以自己構(gòu)造更加恰當(dāng)?shù)木仃嚕沟竭_(dá)事半功倍效果。第〔3〕問(wèn)：討論病態(tài)問(wèn)題求解的算法關(guān)于病態(tài)問(wèn)題的求解，主要的方法是對(duì)原方程作寫(xiě)預(yù)處理，以降低系數(shù)矩陣的條件數(shù)。例如選擇非奇異矩陣，使得方程化為等價(jià)方程組的，原方程的解。原那么上應(yīng)使PAQ的條件數(shù)比A有所改善。一般P和Q可選擇為三角形矩陣或者對(duì)角矩陣。理論上最好選擇對(duì)角陣P和Q，并且應(yīng)當(dāng)滿足：源代碼：cleark=500；i=i:k；forn=1:k；H=hilb(n)；K(n)=cond(H)；D1=diag(sqrt(diag(H)))；D=inv(D1)；Hy=D*H*D；Ky1(n)=cond(Hy)；Ky(n)=log10(Ky1(n)/K(n))；endplot(i,Ky)xlabel(‘n’,’fontsize’,20)ylabel(‘lg(cond(hn))’,’fontsize’,20)實(shí)驗(yàn)結(jié)果：實(shí)驗(yàn)分析:上圖中的〔Dn為Hn的對(duì)角線元素開(kāi)放構(gòu)成的對(duì)角矩陣〕，N為希爾伯特矩陣的階數(shù)。可以看出，隨著矩陣階數(shù)的增大，函數(shù)在之間波動(dòng)。還可以看出，對(duì)于上圖，大多數(shù)的函數(shù)值都是小于或者等于零的，這說(shuō)明經(jīng)過(guò)處理以后的矩陣的條件數(shù)較小，在一定程度上改變了矩陣的性質(zhì)。實(shí)驗(yàn)總結(jié)：對(duì)于希爾伯特矩陣病態(tài)的性質(zhì)，對(duì)于原來(lái)希爾伯特矩陣進(jìn)行據(jù)處理后的新的希爾伯特矩陣在一定區(qū)間里面呈現(xiàn)波動(dòng)的變化規(guī)律。從整體上觀察，對(duì)于大多說(shuō)的N，進(jìn)行預(yù)處理后的希爾伯特矩陣的條件數(shù)有了明顯的降低，這就說(shuō)明預(yù)處理改善了大多數(shù)階數(shù)的希爾伯特矩陣的性質(zhì)。實(shí)驗(yàn)七非線性方程求根實(shí)驗(yàn)7.1〔迭代法、初始值與收斂性〕實(shí)驗(yàn)?zāi)康模撼醪秸J(rèn)識(shí)非線性問(wèn)題的迭代法與線性問(wèn)題迭代法的差異，探討迭代法及初始值與迭代收斂性的關(guān)系。問(wè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)，記錄各算法的迭代過(guò)程。實(shí)驗(yàn)要求：〔1〕取定某個(gè)初始值，分別計(jì)算〔7.1〕-〔7.3〕迭代結(jié)果，它們的收斂性如何？重復(fù)選取不同的初始值，反復(fù)實(shí)驗(yàn)。請(qǐng)自選設(shè)計(jì)一種比擬形象的記錄方式〔如利用MATLAB的圖形功能〕，分析三種迭代法的收斂性與初值選取的關(guān)系?！?〕對(duì)三個(gè)迭代法中的某個(gè)，取不同的初始值進(jìn)行迭代，結(jié)果如何？試分析迭代法對(duì)不同的初值是否有差異？〔3〕線性方程組迭代法的收斂性是不依賴初始值選取的。比擬線性與非線性問(wèn)題迭代的差異，有何結(jié)論和問(wèn)題。實(shí)驗(yàn)過(guò)程：程序：clearclcs=input('請(qǐng)輸入要運(yùn)行的方程，運(yùn)行第幾個(gè)輸入幾s=');clfifs==1%決定坐標(biāo)軸的范圍和初始值a=-1.5;b=2.5;y00