數(shù)值分析實(shí)驗(yàn)二

1、內(nèi)江師范學(xué)院數(shù) 值 分 析實(shí) 驗(yàn) 報(bào) 告 冊(cè)編制 張莉 審定 牟廉明專(zhuān)業(yè)： 班級(jí)： 級(jí) 班學(xué)號(hào)： 姓名： 數(shù)學(xué)與信息科學(xué)學(xué)院2013年9月說(shuō) 明一、學(xué)生在做實(shí)驗(yàn)之前必須要準(zhǔn)備實(shí)驗(yàn)，主要包括預(yù)習(xí)與本次實(shí)驗(yàn)相關(guān)的理論知識(shí)，熟練與本次實(shí)驗(yàn)相關(guān)的軟件操作，收集整理相關(guān)的實(shí)驗(yàn)參考資料，要求學(xué)生在做實(shí)驗(yàn)時(shí)能帶上充足的參考資料；若準(zhǔn)備不充分，則學(xué)生不得參加本次實(shí)驗(yàn)，不得書(shū)寫(xiě)實(shí)驗(yàn)報(bào)告；二、要求學(xué)生要認(rèn)真做實(shí)驗(yàn)，主要是指不得遲到、早退和曠課，在做實(shí)驗(yàn)過(guò)程中要嚴(yán)格遵守實(shí)驗(yàn)室規(guī)章制度，認(rèn)真完成實(shí)驗(yàn)內(nèi)容，極積主動(dòng)地向?qū)嶒?yàn)教師提問(wèn)等； 三、各個(gè)實(shí)驗(yàn)按照學(xué)生水平分別設(shè)置了A、B、C、D四個(gè)等級(jí)，其中對(duì)應(yīng)的難度系數(shù)為1、0

2、.8、0.7、0.6，也可根據(jù)實(shí)際完成情況制定相應(yīng)地的難度系數(shù)，但總體保證難度排序?yàn)锳級(jí)難度最大，B級(jí)次之，C級(jí)較易，D級(jí)最簡(jiǎn)單。四、學(xué)生可以根據(jù)自己對(duì)各個(gè)實(shí)驗(yàn)涉及到的知識(shí)點(diǎn)掌握的程度自由選取A、B、C、D等級(jí)的實(shí)驗(yàn)題目。五、學(xué)生要認(rèn)真工整地書(shū)寫(xiě)實(shí)驗(yàn)報(bào)告，實(shí)驗(yàn)報(bào)告的內(nèi)容要緊扣實(shí)驗(yàn)的要求和目的，不得抄襲他人的實(shí)驗(yàn)報(bào)告；四、根據(jù)實(shí)驗(yàn)準(zhǔn)備、實(shí)驗(yàn)態(tài)度、實(shí)驗(yàn)報(bào)告的書(shū)寫(xiě)、實(shí)驗(yàn)報(bào)告的內(nèi)容進(jìn)行綜合評(píng)定，并給出實(shí)驗(yàn)成績(jī)?cè)u(píng)定分。實(shí)驗(yàn)名稱(chēng)： 實(shí)驗(yàn)二 插值與擬合 指導(dǎo)教師： 吳開(kāi)騰 張莉 實(shí)驗(yàn)時(shí)數(shù)： 4 實(shí)驗(yàn)設(shè)備：安裝了Matlab、C+、VF軟件的計(jì)算機(jī) 實(shí)驗(yàn)日期：2013年 10 月 23、30 日 實(shí)驗(yàn)地點(diǎn)：

3、第五教學(xué)樓北902 實(shí)驗(yàn)?zāi)康模?. 掌握插值方法的基本思想和基本步驟，能夠根據(jù)實(shí)際問(wèn)題選用適當(dāng)?shù)夭逯捣椒ㄟM(jìn)行數(shù)值實(shí)驗(yàn)，并從實(shí)驗(yàn)過(guò)程中理解各類(lèi)插值方法之間的聯(lián)系與區(qū)別。2. 理解各類(lèi)插值方法優(yōu)缺點(diǎn)，并能自行編程求解。3. 理解插值方法與數(shù)據(jù)擬合的區(qū)別，掌握數(shù)據(jù)擬合方法解決實(shí)際問(wèn)題的基本步驟和求解理論，并能通過(guò)數(shù)值實(shí)驗(yàn)進(jìn)行驗(yàn)證。實(shí)驗(yàn)準(zhǔn)備：1. 在開(kāi)始本實(shí)驗(yàn)之前，請(qǐng)回顧教科書(shū)的相關(guān)內(nèi)容；2. 需要一臺(tái)準(zhǔn)備安裝Windows XP Professional操作系統(tǒng)和裝有數(shù)學(xué)軟件的計(jì)算機(jī)。實(shí)驗(yàn)內(nèi)容及要求A題 在某海域測(cè)得一些點(diǎn)處的水深由表1給出，在矩形區(qū)域內(nèi)畫(huà)出海底曲面的圖形；若船的吃水深度為5英尺，請(qǐng)

4、問(wèn)在矩形區(qū)域里的哪些地方船要避免進(jìn)入。表1 海域數(shù)據(jù)表129140103.588185.51951057.5141.52314722.5137.585.54868688157.5107.57781162162117.5-6.5-81356.5-66.584-33.59988949B題 設(shè)從某一實(shí)驗(yàn)中測(cè)得2個(gè)變量和的一組數(shù)據(jù)表2所示。表2 變量關(guān)系數(shù)據(jù)表12345678913456789101054211234求一代數(shù)多項(xiàng)式曲線，使其最好地?cái)M合這組給定數(shù)據(jù)。要求：請(qǐng)用多項(xiàng)式和指數(shù)模型進(jìn)行曲線擬合，觀察其結(jié)果的精度，再直接調(diào)用polyfit命令進(jìn)行比較。C題 考察函數(shù)（1）設(shè)計(jì)插值方案，構(gòu)造龍格（

5、Runge）現(xiàn)象。（2）設(shè)計(jì)插值方法消除龍格（Runge）現(xiàn)象，要求至少用分段線性插值和分段三次插值兩種算法進(jìn)行設(shè)計(jì)，并比較兩種算法的精度。D題 已知數(shù)據(jù)（表3）表3 數(shù)據(jù)表0.561600.562800.564010.565210.827410.826590.825770.81495要求：試用Lagrange插值多項(xiàng)式求時(shí)的函數(shù)近似值說(shuō)明：實(shí)驗(yàn)過(guò)程應(yīng)包括對(duì)問(wèn)題的簡(jiǎn)要分析、求解方法、求解步驟、程序及其必要的圖表等內(nèi)容。實(shí)驗(yàn)過(guò)程：A題的實(shí)驗(yàn)過(guò)程1、實(shí)驗(yàn)中問(wèn)題的重述在某海域測(cè)得一些點(diǎn)處的水深由表1給出，在矩形區(qū)域內(nèi)畫(huà)出海底曲面的圖形；若船的吃水深度為5英尺，請(qǐng)問(wèn)在矩形區(qū)域里的哪些地方船要避免進(jìn)入。

6、2、對(duì)實(shí)驗(yàn)的分析2.1 對(duì)軟件的選擇對(duì)于此題來(lái)說(shuō)，利用Matlab10.0軟件編程計(jì)算。2.2 對(duì)“如何畫(huà)海底曲面圖形”的分析要畫(huà)出海底曲面圖形，首先就應(yīng)該有該海域海底是平滑的假設(shè)；然后，由于測(cè)量點(diǎn)是散亂分布的，所以在平面上先作出測(cè)量點(diǎn)的分布圖；其次，利用二維插值方法griddadt函數(shù)補(bǔ)充一些點(diǎn)的水深；最后，再利用surf函數(shù)作出海底曲面圖。2.3 對(duì)“哪些地方需避免進(jìn)入”的分析由于船只的吃水深度為5，所以在水深為5以下的區(qū)域都是危險(xiǎn)區(qū)域，船只都需避免。于是，先作出等高線圖，再通過(guò)等高線圖得到水深小于5的區(qū)域。2.4 對(duì)結(jié)果呈現(xiàn)的分析由于題中明確要求畫(huà)出海底曲面，所以利用插值，以曲面圖顯示結(jié)

7、果。題目中還要求給出哪些地方需要避免進(jìn)入，所以可以利用圖形顯示出那一部分區(qū)域，也可以給出該危險(xiǎn)區(qū)域中點(diǎn)的坐標(biāo)范圍。3、實(shí)驗(yàn)求解過(guò)程（程序見(jiàn)附錄）首先，利用plot函數(shù)，畫(huà)出所給14個(gè)點(diǎn)的位置分布圖，分布圖見(jiàn)下圖1；然后，利用meshgrid函數(shù)，生成一些采樣點(diǎn)，便于接下來(lái)根據(jù)這些采樣點(diǎn)繪制出整個(gè)海底曲面圖形；其次，利用griddadt函數(shù)進(jìn)行插值；再次，利用surf函數(shù)，畫(huà)出經(jīng)過(guò)插值后的海底曲面圖形，見(jiàn)下圖2。圖1 測(cè)量點(diǎn)的分布圖圖2 海底曲面圖最后，利用contour函數(shù)，畫(huà)出海水的等高線圖，見(jiàn)下圖3。并且，利用clabel函數(shù)標(biāo)記出各個(gè)等高線的水深高度。圖3 等高線圖于是得到該危險(xiǎn)區(qū)域?yàn)椋?/p>

8、圖3中水深為-5的等高線區(qū)域內(nèi)部所有的位置。B題的實(shí)驗(yàn)過(guò)程一、用多項(xiàng)式擬合曲線觀察這一實(shí)驗(yàn)中測(cè)得的變量x、y的數(shù)據(jù)并進(jìn)行分析，當(dāng)i=1,5時(shí)y隨x的增大而減小，當(dāng)i=6,9時(shí)，y隨x的增大呈線性增大。那么我們就可以采取n>=2次多項(xiàng)式進(jìn)行曲線擬合,然后觀察圖像的變化及其精度，并分析。用matlab多項(xiàng)式擬合曲線圖像如下：圖4 二次多項(xiàng)式擬合圖像圖5 三次多項(xiàng)式擬合圖像圖6 四次多項(xiàng)式擬合圖像圖7 五次多項(xiàng)式擬合圖像圖8 六次多項(xiàng)式擬合圖像圖9 七次多項(xiàng)式擬合圖像通過(guò)上面六個(gè)擬合圖像我們可以觀察到當(dāng)擬合多項(xiàng)式的次數(shù)越高，我們就猜想：擬合多項(xiàng)式的次數(shù)越高那么擬合曲線的精度就越高，這一結(jié)論是否

9、正確呢？如果一直這樣下去會(huì)產(chǎn)生什么樣的情況？圖10八次多項(xiàng)式擬合圖像從上一圖像我們可以看出，在x=1,3這個(gè)區(qū)間時(shí)，圖像幅度變化的非常大。但是我們可以觀察到下圖是當(dāng)n=12時(shí)的圖像，可以看出當(dāng)x=9 10范圍內(nèi)時(shí)變化幅度大。圖11次多項(xiàng)式擬合圖像盡管曲線會(huì)因n的取值的增大，會(huì)不同程度的發(fā)生變化。但是通過(guò)上面的圖像我們可以總結(jié)出當(dāng)多項(xiàng)式的次數(shù)越高時(shí)，曲線穿過(guò)點(diǎn)的個(gè)數(shù)是越多的，當(dāng)n=8時(shí)，基本上曲線能穿過(guò)每一個(gè)點(diǎn)，當(dāng)n=12時(shí)也一樣。二、用polyfit命令擬合曲線我們直接調(diào)用命令進(jìn)行作圖： p=polyfit(x,y,2)和p=polyfit(x,y,8)，分別得到以下結(jié)果：圖12 polyfi

10、t(x,y,2)擬合圖像圖13 polyfit(x,y,8)擬合圖像我們發(fā)現(xiàn)直接用polyfit命令非常方便，只需要改變維數(shù)就行了，結(jié)果和多項(xiàng)式擬合的圖像一樣，因此我們可以用這個(gè)命令來(lái)判斷上述的結(jié)論，當(dāng)n=100時(shí)擬合圖像如下：圖14 polyfit(x,y,100)擬合圖像結(jié)論：并非多項(xiàng)式擬合的次數(shù)越高越好，比如當(dāng)n=13時(shí)圖像就如下圖：圖15 polyfit(x,y,13)擬合圖像三、用指數(shù)模型擬合曲線同樣的道理，指數(shù)模型擬合也是一樣，得到以下擬合圖像：Matlab編程：M文件：function yy=linefit2(x,y,A)n=length(x);y=reshape(y,n,1);

11、A=A'yy=Ay;yy=yy'程序：C題的實(shí)驗(yàn)過(guò)程一、 實(shí)驗(yàn)分析（1）、要求構(gòu)造龍格現(xiàn)象，則采用多項(xiàng)式插值的方法呈現(xiàn)。（2）、針對(duì)于（1）中的龍格現(xiàn)象，采用分段線性插值和分段三次插值的兩種算法進(jìn)行改進(jìn)。二、實(shí)驗(yàn)求解過(guò)程（代碼見(jiàn)附錄）（1）、利用多項(xiàng)式插值的方法對(duì)函數(shù)進(jìn)行插值，分別取5、10個(gè)點(diǎn)進(jìn)行插值進(jìn)而構(gòu)造出Runge現(xiàn)象，并呈現(xiàn)圖像。（2）、分別采用分段線性插值和分段三次插值的兩種算法對(duì)（1）中的Runge現(xiàn)象進(jìn)行改進(jìn)計(jì)算。并呈現(xiàn)圖像作比較。三、實(shí)驗(yàn)結(jié)果（1）MATLAB編程實(shí)現(xiàn)Runge現(xiàn)象的呈現(xiàn)如下圖：（2）分別采用分段線性插值和分段三次插值的兩種算法對(duì)（1）中的Ru

12、nge現(xiàn)象進(jìn)行改進(jìn)計(jì)算。并呈現(xiàn)圖像如下圖：由圖像中的關(guān)系可以看出，分段三次插值方法的誤差相對(duì)低些。D題的實(shí)驗(yàn)過(guò)程一、實(shí)驗(yàn)分析針對(duì)問(wèn)題試用Lagrange插值多項(xiàng)式求時(shí)的函數(shù)近似值。首先求出Lagrange插值基函數(shù)的值，然后輸入表格的數(shù)據(jù)及插值點(diǎn)，最后計(jì)算插值點(diǎn)所對(duì)應(yīng)的函數(shù)值。二、實(shí)驗(yàn)求解過(guò)程用Matlab編程可以求出插值點(diǎn)的函數(shù)值。Step 1：在M文件中編寫(xiě)程序，求出和的字符長(zhǎng)度；Step 2：判斷m與n是否相等，如果不等，程序結(jié)束，反之繼續(xù)編寫(xiě)；Step 3：用for循環(huán)求插值基函數(shù)的值；Step 4：在命令窗口中調(diào)用函數(shù)求出插值點(diǎn)對(duì)應(yīng)的函數(shù)值；具體求解的流程圖如下：開(kāi) 始求出x，y的長(zhǎng)

13、度m，n判斷m與n是否相等，如果相等用for循環(huán)求插值基函數(shù)的值在命令窗口中調(diào)用M文件，求出插值點(diǎn)對(duì)應(yīng)的函數(shù)值結(jié) 束三、實(shí)驗(yàn)求解結(jié)果時(shí)的函數(shù)近似值分別為：0.8265 ,0.8268 ,0.8231。附 錄A題實(shí)驗(yàn)程序clcclearx0 = 129 140 103.5 88 185.5 195 105 157.5 107.5 77 81 162 162 117.5; y0 = 7.5 141.5 23 147 22.5 137.5 85.5 -6.5 -81 3 56.5 -66.5 84 -33.5; z0 = -4 8 6 8 6 8 8 9 9 8 8 9 4 9; %原始數(shù)據(jù)x0、y

14、0、z0x,y=meshgrid(75:5:200,-50:5:150);z = griddata(x0,y0,z0,x,y,'v4'); %散點(diǎn)數(shù)據(jù)的插值函數(shù)griddata()figure(1)plot(x0,y0,'r*')figure(2)plot3(x0,y0,z0,'r*'); %已知節(jié)點(diǎn)的三維散點(diǎn)圖hold onsurf(x,y,z) %畫(huà)出經(jīng)插值的圖形figure(3)c,h=contour(x,y,z); %畫(huà)等高線clabel(c,h) %標(biāo)明等高線的高程，高程為5的區(qū)域避免進(jìn)入B題實(shí)驗(yàn)程序clcclear% 數(shù)據(jù)文件（dat

15、a）的讀取v d = xlsread('data.xls');x = v(2,:);y = v(3,:);n = length(x); % 讀取文件中的點(diǎn)的個(gè)數(shù)% 最小二乘法syms a b y1 = a*n+b*sum(x)-sum(y);y2 = a*sum(x)+b*sum(x.2)-sum(x.*y);a b = solve(y1,y2,'a','b')% 誤差的計(jì)算y1 = a+b.*x;format longw1 = sum(abs(y1-y)/sum(y)% 指數(shù)模型的擬合% y=A*exp(B/x)% matlab內(nèi)部函數(shù)的命令a

16、1 = polyfit(x,y,1);% 誤差的計(jì)算y2 = a1(1).*x+a1(2);w2 = sum(abs(y2-y)/sum(y)C題實(shí)驗(yàn)程序function = Runge(a1,a2,n)% (1)構(gòu)造龍格現(xiàn)象% n:構(gòu)造n次的多項(xiàng)式插值x1 = a1:(a2-a1)/n:a2;y1 = 1./(1+x1.2);syms x y = 0;for i = 1:n+1 ss = 1; for j = 1:n+1 if j=i c = (x-x1(j)/(x1(i)-x1(j); ss = ss.*c; else continue end end y = y+ss.*y1(i);en

17、d% 插值多項(xiàng)式x = -5:0.01:5;yy = eval(y);plot(x,yy)% 龍格現(xiàn)象圖 Endclcclear% c題問(wèn)題主程序% Runge現(xiàn)象的呈現(xiàn)的程序a = 5 10;a1=-5;a2=5;for i = 1:length(a) Runge(a1,a2,a(i); hold onendxlabel('Runge現(xiàn)象圖')% 分段線性插值x1 = a1:(a2-a1)/a(2):a2;y1 = 1./(1+x1.2); syms x s1 = ;for i = 1:a(2) b0 = (1-(x-x1(i)/(x1(i+1)-x1(i)*y1(i); b

18、1 = (x-x1(i)/(x1(i+1)-x1(i)*y1(i+1); s11,i = b0+b1;end% 分段三次插值x1 = a1:(a2-a1)/a(2):a2;y1 = 1./(1+x1.2);y2 = -2.*(x1)./(1+x1.2).2);s3 = ;for i = 1:a(2) h = x1(i+1)-x1(i); c0 = (x-x1(i)/(x1(i+1)-x1(i)-1)2)*(2*(x-x1(i)/(x1(i+1)-x1(i)+1)*y1(i); c1 = (x-x1(i)/(x1(i+1)-x1(i)2*(-2*(x-x1(i)/(x1(i+1)-x1(i)+3

19、)*y1(i+1); c2 = (x-x1(i)/(x1(i+1)-x1(i)*(x-x1(i)/(x1(i+1)-x1(i)-1)2*(y2(i); c3 = (x-x1(i)/(x1(i+1)-x1(i)2*(x-x1(i)/(x1(i+1)-x1(i)-1)*y2(i+1); s31,i = c0+c1+c2+c3;end% 畫(huà)圖figure x11 = a1:(a2-a1)/(a(2)*100):a2;y11 = 1./(1+x11.2);% 分段線性插值圖for i = 1:a(2) x = x1(i):(a2-a1)/(10*a(2):x1(i+1); yy = eval(s11,i); plot(x,yy,'r'); hold on endplot(x11,y11)% axis(a1,a2,-0.5,2);xlabel('分段線性圖 紅色：插值 藍(lán)色：原值')figure % 分段三次插值for i = 1:a(2) x = x1(i):(a2-a1)/(10*a(2):x1