數(shù)值計(jì)算方法實(shí)驗(yàn)報(bào)告

1、重 慶 交 通 大 學(xué)學(xué) 生 實(shí) 驗(yàn) 報(bào) 告實(shí)驗(yàn)課程名稱 數(shù)值計(jì)算方法i 開(kāi)課實(shí)驗(yàn)室 數(shù)學(xué)實(shí)驗(yàn)室 學(xué) 院 理學(xué)院 年級(jí)11專業(yè)班 信息與計(jì)算科學(xué) 學(xué) 生 姓 名 李偉凱 學(xué) 號(hào) 631122020203 開(kāi) 課 時(shí) 間 2013 至 2014 學(xué)年第 1 學(xué)期評(píng)分細(xì)則評(píng)分報(bào)告表述的清晰程度和完整性（20分）程序設(shè)計(jì)的正確性（40分）實(shí)驗(yàn)結(jié)果的分析（30分）實(shí)驗(yàn)方法的創(chuàng)新性（10分）總成績(jī)教師簽名鄒昌文實(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)行的，如何才能確保g

2、auss消去法作為數(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é)果。若每步消去過(guò)程總選取按模最大的元素作為主元，結(jié)果又如何？分析實(shí)驗(yàn)的結(jié)果。（3）取矩陣階數(shù)n=20或者更大，重復(fù)上述實(shí)驗(yàn)過(guò)程，觀

3、察記錄并分析不同的問(wèn)題及消去過(guò)程中選擇不同的主元時(shí)計(jì)算結(jié)果的差異，說(shuō)明主元素的選取在消去過(guò)程中的作用。（4）選取其他你感興趣的問(wèn)題或者隨機(jī)生成矩陣，計(jì)算其條件數(shù)。重復(fù)上述實(shí)驗(yàn)，觀察記錄并分析實(shí)驗(yàn)結(jié)果。實(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)，使得充

4、分小。實(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ù)。思考題一：（vadermonde矩陣）設(shè)

5、，其中，（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)

6、提取矩陣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矩陣實(shí)驗(yàn)程序：m文件程序?yàn)椋篺unction x=gauss(n,r)n=input('請(qǐng)輸入矩陣a的階數(shù):n=')a=diag(6*ones(1,n)+diag(

7、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)pausem,n=size(a);nb=n+1;ab=a br=input('請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入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(&

8、#39;輸入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);end（1）取矩陣a的階數(shù):n=10，自動(dòng)選取主元：>> f

9、ormat long>> 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= 1 1 1 1 1 1 1 1 1 1選取絕對(duì)值最小的元

10、素為主元：>> 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.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 0.99999999999999 1.00000000000001 0.99999999999998 1.00000000000003（2）取矩陣a的

11、階數(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= 1 1 1 1 1 1 1 1 1 1選取絕對(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.00

12、000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 0.99999999999999 1.00000000000001 0.99999999999998 1.00000000000003（3）取矩陣a的階數(shù):n=20，手動(dòng)選取主元： 選取絕對(duì)值最大的元素為主元：>> gauss請(qǐng)輸入矩陣a的階數(shù):n=20條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=1p = 1pp = 2.621437500000000e+006ans = 1 1 1 1 1 1 1

13、1 1 1 1 1 1 1 1 1 1 1 1 1 選取絕對(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.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000000 1.00000000000001 0.99999999999997 1.000

14、00000000006 0.99999999999989 1.00000000000023 0.99999999999955 1.00000000000090 0.99999999999821 1.00000000000352 0.99999999999318 1.00000000001273 0.99999999997817 1.00000000002910（4）該題目的m程序如下所示function x=gauss(n,r)n=input('請(qǐng)輸入矩陣a的階數(shù):n=')a= hilb(n)b=a*ones(n,1)p=input('條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=&

15、#39;)pp=cond(a,p)pausem,n=size(a);nb=n+1;ab=a br=input('請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入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(

16、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);end >>gauss請(qǐng)輸入矩陣a的階數(shù)：n=7n = 7a = 6 1 0 0 0 0 0 8 6 1 0 0 0 0 0 8 6 1 0 0 0 0 0 8 6 1 0 0 0 0 0 8 6 1 0 0 0 0 0

17、 8 6 1 0 0 0 0 0 8 6b = 7 15 15 15 15 15 14pp = 3.174999999999999e+002ab = 6 1 0 0 0 0 0 7 8 6 1 0 0 0 0 15 0 8 6 1 0 0 0 15 0 0 8 6 1 0 0 15 0 0 0 8 6 1 0 15 0 0 0 0 8 6 1 15 0 0 0 0 0 8 6 14請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入0：r=1r = 1條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=1p = 1i = 1ans = 0.99999999999869 1.00000000004337 0.9999999996

18、4299 1.00000000121143 0.99999999803038 1.00000000152825 0.99999999954491顯然的是，該問(wèn)題在主元選取與算出結(jié)果有著很大的關(guān)系，取絕對(duì)值大的元素作為主元比取絕對(duì)值小的元素作為主元時(shí)產(chǎn)生的結(jié)果比較準(zhǔn)確，即選取絕對(duì)值小的主元時(shí)結(jié)果產(chǎn)生了較大的誤差，條件數(shù)越大產(chǎn)生的誤差就越大實(shí)驗(yàn)體會(huì)：運(yùn)用高斯消去法求解線性方程組問(wèn)題的時(shí)候，主元的選取和相應(yīng)的消去法的選取決定了該算法的穩(wěn)定性，選取絕對(duì)值大的元素比選取絕對(duì)值比較小的元素作為主元時(shí)對(duì)結(jié)果產(chǎn)生的誤差影響比較小。并且增加條件數(shù)反而對(duì)結(jié)果的誤差產(chǎn)生更大的影響。并且在運(yùn)算中要盡量避免出現(xiàn)運(yùn)用小數(shù)

19、作為除數(shù)，使數(shù)量級(jí)加大，令大數(shù)吃掉小數(shù)的情況發(fā)生。實(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é)果。（2）選擇一系列維數(shù)遞增的矩陣（可以是隨機(jī)生成的），比較函數(shù)“c

20、ondest”所需機(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ù)。程序代碼：1保存m文件名為：fanshu.mn=input('please input n:n=') a=fix(100*rand(n)+1 x=ones(n,1) b=a*x

21、 data=rand(n)*0.00001 datb=rand(n,1)*0.00001 a=a+datab=b+datbxx=geshow(a,b) x0=norm(xx-x,1)/norm(x,1) 保存m文件名為：geshow.m function x=geshow(a,b) m,n=size(a);nb=n+1;ab=a b;for i=1:n-1 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); endendx=zeros(n,1);x(n)=ab(n,nb)/ab(n,n);fo

22、r i=n-1:-1:1 x(i)=(ab(i,nb)-ab(i,i+1:n)*x(i+1:n)/ab(i,i);end2保存m文件名為：cond2.mfunction cond2(a) b=a'*a;v1,d1=eig(b);v2,d2=eig(b(-1);cond2a=sqrt(max(max(d1)*sqrt(max(max(d2)end保存m文件為：shiyan52.mformat longfor n=10:10:100n=n a=fix(100*randn(n); condesta=condest(a) cond2(a) conda2=cond(a,2) pause end3

23、保存m文件為：sy5_2.mn=input('please input n:n=') %輸入矩陣的階數(shù)a=fix(100*rand(n)+1; %隨機(jī)生成一個(gè)矩陣ax=ones(n,1); %假設(shè)知道方程組的解全為1b=a*x; data=rand(n)*0.00001; datb=rand(n,1)*0.00001; a=a+data;b=b+datb;xx=geshow(a,b); x0=norm(xx-x,1)/norm(x,1) x00=cond(a)/(1-norm(inv(a)*norm(xx-x)*(norm(xx-x)/(norm(a)+norm(datb)/n

24、orm(b) datx=abs(x0-x00) （4）format longfor n=4:11 n=n hi=hilb(n); cond1hi=cond(hi,1) cond2hi=cond(hi,2) condinfhi=cond(hi,inf) pauseend實(shí)驗(yàn)結(jié)果及其分析：(1)>>fanshuplease input n:n=6n=6 a = 82 28 96 80 68 71 91 55 49 96 76 4 13 96 81 66 75 28 92 97 15 4 40 5 64 16 43 85 66 10 10 98 92 94 18 83x = 1 1 1

25、1 1 1b = 425 371 359 253 284 395data = 1.0e-005 * columns 1 through 4 0.694828622975817 0.765516788149002 0.709364830858073 0.118997681558377 0.317099480060861 0.795199901137063 0.754686681982361 0.498364051982143 0.950222048838355 0.186872604554379 0.276025076998578 0.959743958516081 0.034446080502

26、909 0.489764395788231 0.679702676853675 0.340385726666133 0.438744359656398 0.445586200710900 0.655098003973841 0.585267750979777 0.381558457093009 0.646313010111265 0.162611735194631 0.223811939491137 columns 5 through 6 0.751267059305653 0.547215529963803 0.255095115459269 0.138624442828679 0.5059

27、57051665142 0.149294005559058 0.699076722656686 0.257508254123737 0.890903252535799 0.840717255983663 0.959291425205445 0.254282178971531datb = 1.0e-005 * 0.814284826068817 0.243524968724989 0.929263623187228 0.349983765984809 0.196595250431208 0.251083857976031a = columns 1 through 4 82.00000694828

28、6227 28.000007655167881 96.000007093648307 80.000001189976814 91.000003170994802 55.000007951999009 49.000007546866819 96.000004983640522 13.000009502220488 96.000001868726045 81.000002760250766 66.000009597439586 92.000000344460801 97.000004897643961 15.000006797026769 4.000003403857266 64.00000438

29、7443596 16.000004455862008 43.000006550980039 85.000005852677504 10.000003815584570 98.000006463130106 92.000001626117353 94.000002238119393 columns 5 through 6 68.000007512670592 71.000005472155294 76.000002550951152 4.000001386244429 75.000005059570512 28.000001492940054 40.000006990767226 5.00000

30、2575082541 66.000008909032530 10.000008407172560 18.000009592914253 83.000002542821790b = 1.0e+002 * 4.250000081428483 3.710000024352497 3.590000092926362 2.530000034998376 2.840000019659525 3.950000025108386xx = 1.000000428381561 0.999999831084732 1.000002591353244 0.999999600566705 0.9999982265399

31、37 0.999997826103164x0 =1.255906711054392e-006的計(jì)算結(jié)果為：1.255906711054392e-006(2)ncondestacond2aconda2101.152530883943102e+00232.8905456307542132.89054563075420203.470959631940668e+00265.5412238417896665.54122384178720306.050503865112835e+0021.126539755706398e+0021.126539755706322e+002403.5494878925824

32、70e+00261.3753756968344861.37537569683365506.855018184779408e+00281.1213899375359481.12138993753482601.082004656409367e+0041.704830815154781e+0031.704830815108527e+003703.234679145192132e+0033.878481155980936e+0023.878481155978439e+002808.318226153918658e+00286.2381429985251386.23814299853018902.063

33、634143407935e+0032.120696380331705e+0022.120696380331079e+0021001.536592818758897e+0031.559132035738491e+0021.559132035738373e+002(3)sy5_2please input n:n=8n = 8x0 = 3.064664488849900e-007x00 = 5.064016823313296e-007datx = 1.999352334463396e-007給出對(duì)的估計(jì)是5.064016823313296e-007的理論結(jié)果是：3.064664488849900e-

34、007結(jié)果相差：1.999352334463396e-007(4)ncond1hicond2hicondinfhi42.837499999999738e+0041.551373873892786e+0042.837499999999739e+00459.436559999999364e+0054.766072502414135e+0059.436559999999336e+00562.907027900294878e+0071.495105864009243e+0072.907027900294064e+00779.851948897194700e+0084.753673565864586e+

35、0089.851948897198483e+00883.387279082022742e+0101.525757545841988e+0103.387279081949470e+01091.099650993366047e+0124.931544439891016e+0111.099650991701052e+012103.535372424347474e+0131.602528637652488e+0133.535372455375642e+013111.230369955362001e+0155.223946340715823e+0141.230369938308720e+015討論：線性

36、代數(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)體會(huì)：在本次實(shí)驗(yàn)中，使我們知道了矩陣條件數(shù)對(duì)線性代數(shù)方程組求解的影響，條件數(shù)越大，對(duì)最后解的影響的越大，hilbert矩陣是一個(gè)很”病態(tài)”的矩陣，他的條件數(shù)隨著階數(shù)的增加而增大，每增加一階，條件數(shù)就增大一個(gè)數(shù)量級(jí)，在求解的過(guò)程中要盡量避免hilbert矩陣實(shí)驗(yàn)六 解線性方程組的迭代法實(shí)驗(yàn)6.1（病態(tài)的線性方程組的求解

37、）問(wèn)題提出：理論的分析表明，求解病態(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)題求解的算法程序代碼：gauss消去法程序：function x

38、=gauss(n,r)n=input('請(qǐng)輸入矩陣的階數(shù):n=')a=hilb(n);%構(gòu)造hilbert矩陣p=input('條件數(shù)對(duì)應(yīng)的范數(shù)是p-范數(shù)：p=')pp=cond(a,p)pausem,n=size(a);nb=n+1;ab=a br=input('請(qǐng)輸入是否為手動(dòng)，手動(dòng)輸入1，自動(dòng)輸入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

39、=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)=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);endj迭代法程序：n=input('

40、系數(shù)矩陣的階數(shù):');a=hilb(n);%構(gòu)造hilbert矩陣for i=1:n; a0(i)=1; %給定解 x(i)=0;end;b=a*a0' %由給定的解算出相應(yīng)的b%進(jìn)行迭代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);end;for k=j+1:n; x(j)=x(j)-a(j,k)*y(k)/a(j,j);end;end;end;xgs迭代程序：n=input('系數(shù)矩陣的階數(shù):');%對(duì)題中給定的矩陣進(jìn)行消元a2

41、=hilb(n);for i=1:n; a02(i)=1; x2(i)=0;end;b2=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);end;for k=j+1:n; x2(j)=x2(j)-a2(j,k)*x2(k)/a2(j,j);end;end;end;x2sor迭代程序：n=input('系數(shù)矩陣的階數(shù):');ss=input('松弛因子：');%對(duì)題中給定的矩陣進(jìn)行消元a3=hilb(n

42、);for i=1:n; a03(i)=1; x3(i)=0;end;b3=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);end;for k=j+1:n; x3(j)=x3(j)-a3(j,k)*x3(k)/a3(j,j);end;x3(j)=(1-ss)*rc+ss*x3(j);end;end;x3實(shí)驗(yàn)結(jié)果及其分析：給定各分量為1的解，計(jì)算出右端作為問(wèn)題。1 選擇問(wèn)題的維數(shù)為6時(shí)：用gauss消去法求得的解與

43、精確解一致；取初始向量為0，用j迭代方法迭代出現(xiàn)發(fā)散的不穩(wěn)定現(xiàn)象，無(wú)法求解；取初始向量為0，用gs迭代方法迭代不發(fā)散，能求得解，但收斂非常緩慢，當(dāng)?shù)螖?shù)取得相當(dāng)大（100000次）時(shí)解仍在精確解附近波動(dòng)；取初始向量為0，用sor迭代方法迭代不發(fā)散，能求得解，當(dāng)松弛因子去1.25左右時(shí)收斂較gs迭代快一些，但仍非常緩慢。2 選擇問(wèn)題的維數(shù)為20時(shí)：用gauss消去法求得的解與精確解相差很大，相差10的量級(jí)。取初始向量為0，用j迭代方法迭代發(fā)散，無(wú)法求解；取初始向量為0，用gs迭代方法迭代不發(fā)散，能求得解，但收斂非常緩慢，迭代100000次后，算得的值與精確值1相差0.001的量級(jí)。取初始向量為

44、0，用sor迭代方法迭代不發(fā)散，能求得解，但收斂非常緩慢。從上面的結(jié)果可以看出當(dāng)病態(tài)問(wèn)題的階數(shù)升高時(shí)作為直接法的gauss消去法又能求解變成不能求解。而gs和sor迭代法在階數(shù)升高時(shí)仍能求解。但在階數(shù)較低時(shí)直接法能求得精確解而迭代發(fā)卻總存在一定的誤差。可見(jiàn)直接法與迭代法各有各的優(yōu)勢(shì)與不足。關(guān)于病態(tài)問(wèn)題的求解，主要的方法是對(duì)原方程作某些預(yù)處理，以降低系數(shù)矩陣的條件數(shù)。可以采取將系數(shù)矩陣a的每一行本別乘上適當(dāng)常數(shù)的方法。即找到可逆的對(duì)角陣和使方程組化為 理論上最好選擇對(duì)角陣滿足：實(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)鍵是怎樣才能使迭代收斂且有較快的