數(shù)值計算方法I實(shí)驗(yàn)報告作業(yè)

1、16 重 慶 交 通 大 學(xué)學(xué) 生 實(shí) 驗(yàn) 報 告實(shí)驗(yàn)課程名稱 數(shù)值計算方法I 開課實(shí)驗(yàn)室 數(shù)學(xué)實(shí)驗(yàn)室 學(xué) 院 理學(xué)院 年級12專業(yè)班 信息與計算科學(xué)1班學(xué) 生 姓 名 浦中正 學(xué) 號 631222020101 開 課 時 間 2012 至 2013 學(xué)年第 2 學(xué)期評分細(xì)則評分報告表述的清晰程度和完整性（20分）程序設(shè)計的正確性（40分）實(shí)驗(yàn)結(jié)果的分析（30分）實(shí)驗(yàn)方法的創(chuàng)新性（10分）總成績教師簽名 實(shí)驗(yàn)一 誤差分析試驗(yàn)1.1(病態(tài)問題）問題提出：考慮一個高次的代數(shù)多項(xiàng)式顯然該多項(xiàng)式的全部根為1,2,20共計20個，且每個根都是單重的?，F(xiàn)考慮該多項(xiàng)式的一個擾動其中是一個非常小的數(shù)。這相當(dāng)于

2、是對（1.1）中的系數(shù)作一個小的擾動。我們希望比較（1.1）和（1.2）根的差別，從而分析方程（1.1）的解對擾動的敏感性。實(shí)驗(yàn)內(nèi)容：為了實(shí)現(xiàn)方便，我們先介紹兩個MATLAB函數(shù)：“roots”和“poly”。其中若變量a存儲n+1維的向量，則該函數(shù)的輸出u為一個n維的向量。設(shè)a的元素依次為，則輸出u的各分量是多項(xiàng)式方程的全部根；而函數(shù) 的輸出b是一個n+1維向量，它是以n維向量v的各分量為根的多項(xiàng)式的系數(shù)。可見“roots”和“poly”是兩個互逆的運(yùn)算函數(shù)。上述簡單的MATLAB程序便得到（1.2）的全部根，程序中的“ess”即是（1.2）中的。實(shí)驗(yàn)要求：（1）選擇充分小的ess，反復(fù)進(jìn)行

3、上述實(shí)驗(yàn)，記錄結(jié)果的變化并分析它們。如果擾動項(xiàng)的系數(shù)很小，我們自然感覺（1.1）和（1.2）的解應(yīng)當(dāng)相差很小。計算中你有什么出乎意料的發(fā)現(xiàn)？表明有些解關(guān)于如此的擾動敏感性如何？ （2）將方程（1.2）中的擾動項(xiàng)改成或其它形式，實(shí)驗(yàn)中又有怎樣的現(xiàn)象？ （3）（選作部分）請從理論上分析產(chǎn)生這一問題的根源。注意我們可以將方程（1.2）寫成展開的形式， 同時將方程的解x看成是系數(shù)的函數(shù)，考察方程的某個解關(guān)于的擾動是否敏感，與研究它關(guān)于的導(dǎo)數(shù)的大小有何關(guān)系？為什么？你發(fā)現(xiàn)了什么現(xiàn)象，哪些根關(guān)于的變化更敏感？實(shí)驗(yàn)過程：對ess取不同的值，帶入函數(shù)得出相同i的情況下，結(jié)果的差異；對i取不同的值，帶入函數(shù)得出

4、相同ess的情況下，結(jié)果的不同。程序：function effect(ess)ve=zeros(1,21);ve(2)=ess;format shortroots(poly(1:20)+ve)實(shí)驗(yàn)結(jié)果:essansi0.0000010.00000010.000000010.0000000010.0000000001 2 21.3025 + 1.5672i 21.3025 - 1.5672i 18.5028 + 3.6004i 18.5028 - 3.6004i 15.1651 + 3.7612i 15.1651 - 3.7612i 12.4866 + 2.8828i 12.4866 - 2.8

5、828i 10.5225 + 1.7196i 10.5225 - 1.7196i 9.0448 + 0.5947i 9.0448 - 0.5947i 7.9489 7.0025 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.4220 + 0.9992i 20.4220 - 0.9992i 18.1572 + 2.4702i 18.1572 - 2.4702i 15.3147 + 2.6987i 15.3147 - 2.6987i 12.8460 + 2.0622i 12.8460 - 2.0622i 10.9206 + 1.1013i 10.9206

6、 - 1.1013i 9.5767 9.1074 7.9948 7.0001 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 19.8692 + 0.4830i 19.8692 - 0.4830i 17.8767 + 1.5257i 17.8767 - 1.5257i 15.4031 + 1.7507i 15.4031 - 1.7507i 13.1302 + 1.2796i 13.1302 - 1.2796i 11.2511 + 0.4806i 11.2511 - 0.4806i 9.9302 9.0101 7.9990 7.0001 6.0000 5.00

7、00 4.0000 3.0000 2.0000 1.0000 19.9513 19.2322 17.6604+ 0.7003i 17.6604 - 0.7003i 15.4537 + 0.8891i 15.4537 - 0.8891i 13.3455 + 0.4972i 13.3455 - 0.4972i 11.8733 11.0239 10.0014 8.9984 8.0004 6.9999 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 19.9952 19.0330 17.8638 17.2299 15.5044 + 0.1556i 15.5044 -

8、 0.1556i 13.6984 13.2332 11.9006 11.0500 9.9829 9.0050 7.9989 7.0002 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 3 20.2296 + 0.8314i 20.2296 - 0.8314i 18.1285 + 2.1925i 18.1285 - 2.1925i 15.4044 + 2.5140i 15.4044 - 2.5140i 12.9354 + 1.9888i 12.9354 - 1.9888i 10.9628 + 1.0902i 10.9628 - 1.0902i 9.5623

9、9.1226 7.9934 7.0003 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 19.7653 + 0.2995i 19.7653 - 0.2995i 17.8425 + 1.3013i 17.8425 - 1.3013i 15.4531 + 1.5826i 15.4531 - 1.5826i 13.1873 + 1.1911i 13.1873 - 1.1911i 11.2774 + 0.4266i 11.2774 - 0.4266i 9.9429 9.0061 8.0000 6.9999 6.0000 5.0000 4.0000 3.0000 2

10、.0000 1.0000 19.9765 19.1384 17.6369 + 0.5206i 17.6369 - 0.5206i 15.4896 + 0.7750i 15.4896 - 0.7750i 13.3819 + 0.5003i 13.3819 - 0.5003i 11.8068 11.0737 9.9859 9.0022 7.9998 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 19.9976 19.0177 17.9330 17.1227 15.8172 15.1710 13.9089 13.0372 11.9976 10.99

11、44 10.0038 8.9987 8.0003 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 19.9994 19.0055 17.9744 17.0653 15.8551 15.1864 13.7981 13.1723 11.9091 11.0465 9.9846 9.0041 7.9992 7.0001 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 4 19.8184 19.5494 17.8047 + 1.0917i 17.8047 - 1.0917i 15.4956 + 1.4207i 15.4

12、956 - 1.4207i 13.2439 + 1.1128i 13.2439 - 1.1128i 11.3003 + 0.4128i 11.3003 - 0.4128i 9.9340 9.0099 7.9992 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 19.9881 19.0816 17.6096 + 0.3292i 17.6096 - 0.3292i 15.5064 + 0.6411i 15.5064 - 0.6411i 13.4056 + 0.4138i 13.4056 - 0.4138i 11.8365 11.0604 9.98

13、88 9.0014 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 19.9984 19.0139 17.9374 17.1353 15.6615 15.3816 13.6778 13.2631 11.8919 11.0513 9.9847 9.0038 7.9993 7.0001 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9979 18.0104 16.9680 16.0585 14.9185 14.0786 12.9473 12.0286 10.9894 10.0

14、031 8.9993 8.0001 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15.9995 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 519.9948 19.0405 17.7775 17.3662 15.5112 + 0.4154i 15.5112 - 0.4154i 13.4374 + 0.17

15、97i 13.4374 - 0.1797i 11.8905 11.0414 9.9903 9.0020 7.9997 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 19.9998 19.0013 17.9951 17.0102 15.9866 15.0100 13.9989 12.9933 12.0089 10.9933 10.0034 8.9988 8.0003 6.9999 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0003 18.9976 18.0083 16.9825 16.0228 1

16、4.9808 14.0104 12.9970 12.0005 10.9994 10.0007 8.9995 8.0002 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15.9995 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15.9995

17、15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 6 20.0007 18.9932 18.0283 16.9158 16.1375 14.7973 14.1994 12.8695 12.0816 10.9679 10.0111 8.9973 8.0005 6.9999 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15.99

18、95 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.000020.0002 18.9987 18.0038 16.9952 15.9995 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15.999

19、5 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15.9995 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 721 20.0002 18.9987 18.0038 16.9952 15

20、.9995 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15.9995 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15

21、.9995 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15.9995 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000 20.0002 18.9987 18.0038 16.9952 15

22、.9995 15.0117 13.9791 13.0207 11.9868 11.0058 9.9985 9.0002 8.0000 7.0000 6.0000 5.0000 4.0000 3.0000 2.0000 1.0000實(shí)驗(yàn)分析：（1）對于同一項(xiàng)即(i=21,20,19,18,)的擾動，隨著擾動系數(shù)的減小，擾動也越來越小（2）當(dāng)相同的擾動系數(shù)作用于不同的項(xiàng)式，隨著冪數(shù)的降低，擾動越來越小，擾動的項(xiàng)的冪數(shù)越高，擾動越敏感，誤差越大。實(shí)驗(yàn)總結(jié)： 通過本次試驗(yàn)，了解了病態(tài)問題的產(chǎn)生原因和產(chǎn)生的影響。在以后的學(xué)習(xí)中，對于病態(tài)問題的處理應(yīng)該給予相當(dāng)大的重視。實(shí)驗(yàn)二 插值法實(shí)驗(yàn)2.1（多項(xiàng)式插值

23、的振蕩現(xiàn)象）問題提出：考慮一個固定的區(qū)間上用插值逼近一個函數(shù)。顯然拉格朗日插值中使用的節(jié)點(diǎn)越多，插值多項(xiàng)式的次數(shù)就越高。我們自然關(guān)心插值多項(xiàng)式的次數(shù)增加時，是否也更加靠近被逼近的函數(shù)。龍格（Runge）給出一個例子是極著名并富有啟發(fā)性的。設(shè)區(qū)間-1,1上函數(shù) 實(shí)驗(yàn)內(nèi)容：考慮區(qū)間-1,1的一個等距劃分，分點(diǎn)為 則拉格朗日插值多項(xiàng)式為 其中的是n次拉格朗日插值基函數(shù)。實(shí)驗(yàn)要求：（1） 選擇不斷增大的分點(diǎn)數(shù)目n=2,3.，畫出原函數(shù)f(x)及插值多項(xiàng)式函數(shù)在-1，1上的圖像，比較并分析實(shí)驗(yàn)結(jié)果。（2）選擇其他的函數(shù)，例如定義在區(qū)間-5，5上的函數(shù)重復(fù)上述的實(shí)驗(yàn)看其結(jié)果如何。（3）區(qū)間a,b上切比雪夫

24、點(diǎn)的定義為 以為插值節(jié)點(diǎn)構(gòu)造上述各函數(shù)的拉格朗日插值多項(xiàng)式，比較其結(jié)果。實(shí)驗(yàn)2.2（樣條插值的收斂性）問題提出：多項(xiàng)式插值是不收斂的，即插值的節(jié)點(diǎn)多，效果不一定就好。對樣條函數(shù)插值又如何呢？理論上證明樣條插值的收斂性是比較困難的，但通過本實(shí)驗(yàn)可以驗(yàn)證這一理論結(jié)果。實(shí)驗(yàn)內(nèi)容：請按一定的規(guī)則分別選擇等距或者非等距的插值節(jié)點(diǎn)，并不斷增加插值節(jié)點(diǎn)的個數(shù)?？紤]實(shí)驗(yàn)2.1中的函數(shù)或選擇其他你有興趣的函數(shù)，可以用MATLAB的函數(shù)“spline”作此函數(shù)的三次樣條插值。實(shí)驗(yàn)要求：（1）隨節(jié)點(diǎn)個數(shù)增加，比較被逼近函數(shù)和樣條插值函數(shù)誤差的變化情況。分析所得結(jié)果并與拉格朗日多項(xiàng)式插值比較。（2）樣條插值的思想是早

25、產(chǎn)生于工業(yè)部門。作為工業(yè)應(yīng)用的例子考慮如下問題：某汽車制造商用三次樣條插值設(shè)計車門的曲線，其中一段的數(shù)據(jù)如下：xk012345678910yk0.00.791.532.192.713.033.272.893.063.193.29yk0.80.2 思考題一：（二維插值）在一丘陵地帶測量高程，x和y方向每隔100米測一個點(diǎn)，得高程數(shù)據(jù)如下。試用MATLAB的二維插值函數(shù)“interp2”進(jìn)行插值，并由此找出最高點(diǎn)和該點(diǎn)的高程。yx100200300400100636697624478200698712630478300680674598412400662626552334相關(guān)MATLAB函數(shù)提示：

26、plot(x,y) 作出以數(shù)據(jù)(x(i),y(i)為節(jié)點(diǎn)的折線圖，其中x,y為同長度的向量subplot(m,n,k) 將圖形窗口分為m*n個子圖，并指向第k幅圖yi=interp1(x,y,xi) 根據(jù)數(shù)據(jù)(x,y)給出在xi的分段線性插值結(jié)果yipp=spline(x,y) 返回樣條插值的分段多項(xiàng)式(pp)形式結(jié)構(gòu)pp=csape(x,y,邊界類型,邊界值) 生成各種邊界條件的三次樣條插值yi=ppval(pp,xi) pp樣條在xi的函數(shù)值ZI=interp2(x,y,z,xi,yi) x,xi為行向量，y,yi為列向量，z為矩陣的雙線性二維插值ZI=interp2(,'spli

27、ne') 使用二元三次樣條插值ZI=griddata(x,y,z,xi,yi) x,y,z均為向量（不必單調(diào)）表示數(shù)據(jù)，xi,yi為網(wǎng)格向量的三角形線性插值（不規(guī)則數(shù)據(jù)的二維插值）實(shí)驗(yàn)過程：牛頓插值源程序:function y=j(xx,k)n=length(xx);s=1;for m=1:n if m=k s=s*(xx(k)-xx(m); end endy=s; 均差函數(shù)function jj=juncha(x0,y0)jj(1)=y0(1);n=length(x0);for m=2:n s=0; for mm=1:m s=s+y0(mm)/j(x0(1:m),mm);end jj

28、(m)=s;end基函數(shù)function b=basic(x0,x)b(1)=1;n=length(x0);for m=2:n s=1; for mm=1:(m-1) s=s*(x-x0(mm); end b(m)=s;end主程序n=10;%n次插值多項(xiàng)式（n+1）個等距插值點(diǎn)m=200;%要看m+1個等距點(diǎn)上的插值函數(shù)值a=-5;b=5;x0=a:(b-a)/n):b;y0=atan(x0);x=a:(b-a)/m):b;%插值函數(shù)在這些點(diǎn)上的值y=Lagrange(x0,y0,x);x0=a:(b-a)/m):b;y0=atan(x0);plot(x0,y0,'bo-',

29、x,y,'rs-');實(shí)驗(yàn)結(jié)果:原函數(shù)圖象和牛頓插值多項(xiàng)式圖象figure: f(x)=x/(1+x4) 插值點(diǎn)個數(shù)n=5 n=10 n=15 n=20figure: f(x)=1/(1+25*x2) n=5 n=10 n=15 n=20figure: f(x)=arctan(x) n=5 n=10 n=15 n=20實(shí)驗(yàn)分析：從圖形可以看的出，Ln(x)在兩端點(diǎn)處發(fā)生劇烈的抖動，造成了較大的誤差。即插值點(diǎn)越多，龍格現(xiàn)象越明顯，兩端誤差越大。實(shí)驗(yàn)總結(jié)： 通過本次實(shí)驗(yàn)我對龍格現(xiàn)象有了更加深刻的認(rèn)識。龍格現(xiàn)象表明高階插值多項(xiàng)式Ln（x）近似被逼近函數(shù)效果并不好，需要用其他的方法來避

30、免這樣的錯誤。挺有意義的。實(shí)驗(yàn)2.2（1）（樣條插值的收斂性）程序n=10;%n次插值多項(xiàng)式（n+1）個等距插值點(diǎn)m=200;%要看m+1個等距點(diǎn)上的插值函數(shù)值a=-5;b=5;x0=a:(b-a)/n):b;y0=x0./(1+x0.*x0.*x0.*x0);x=a:(b-a)/m):b;%插值函數(shù)在這些點(diǎn)上的值y=Lagrange(x0,y0,x);x0=a:(b-a)/m):b;y0=x0./(1+x0.*x0.*x0.*x0);plot(x0,y0,'-',x,y);實(shí)驗(yàn)結(jié)果圖：虛線為原函數(shù)圖象，實(shí)線為三次樣條插值函數(shù)的圖象figure: f(x)=x/(1+x4)結(jié)果

31、分析： 通過對牛頓插值函數(shù)的圖象與樣條插值函數(shù)的圖象的比較，隨著n的增大，使用Newton插值多項(xiàng)式會出現(xiàn)龍格現(xiàn)象。隨著n的增大，三次樣條插值多項(xiàng)式將越來越接近被插值的函數(shù)（2）樣條插值的思想是早產(chǎn)生于工業(yè)部門。作為工業(yè)應(yīng)用的例子考慮如下問題：某汽車制造商用三次樣條插值設(shè)計車門的曲線，其中一段的數(shù)據(jù)如下：xk012345678910yk0.00.791.532.192.713.033.272.893.063.193.29yk0.80.2程序：x=0 1 2 3 4 5 6 7 8 9 10;y=0.0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29;pp=csape(x,y,'second',0.8,0.2);pp.coefs;xi=0:0.1:10;yi=ppval(pp,xi);plot(x,y,'o',xi,yi);實(shí)驗(yàn)結(jié)果圖：實(shí)驗(yàn)總結(jié)： 對于函數(shù)逼近由于高次插值存在龍格現(xiàn)象，所以要用分段低次插值（特別是三次樣條差值），在理論上和實(shí)際應(yīng)用中更有意義。思考題一：（二維插值）一實(shí)驗(yàn)內(nèi)容：在一丘陵地帶測量高程，x和y方向每隔100米測一個點(diǎn)，得高程數(shù)據(jù)如下。試用MATLAB的二維插值函數(shù)“