數(shù)值分析實驗二

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

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

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

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

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

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

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

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

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

10、t(x,y,2)擬合圖像圖13 polyfit(x,y,8)擬合圖像我們發(fā)現(xiàn)直接用polyfit命令非常方便，只需要改變維數(shù)就行了，結果和多項式擬合的圖像一樣，因此我們可以用這個命令來判斷上述的結論，當n=100時擬合圖像如下：圖14 polyfit(x,y,100)擬合圖像結論：并非多項式擬合的次數(shù)越高越好，比如當n=13時圖像就如下圖：圖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題的實驗過程一、 實驗分析（1）、要求構造龍格現(xiàn)象，則采用多項式插值的方法呈現(xiàn)。（2）、針對于（1）中的龍格現(xiàn)象，采用分段線性插值和分段三次插值的兩種算法進行改進。二、實驗求解過程（代碼見附錄）（1）、利用多項式插值的方法對函數(shù)進行插值，分別取5、10個點進行插值進而構造出Runge現(xiàn)象，并呈現(xiàn)圖像。（2）、分別采用分段線性插值和分段三次插值的兩種算法對（1）中的Runge現(xiàn)象進行改進計算。并呈現(xiàn)圖像作比較。三、實驗結果（1）MATLAB編程實現(xiàn)Runge現(xiàn)象的呈現(xiàn)如下圖：（2）分別采用分段線性插值和分段三次插值的兩種算法對（1）中的Ru

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

13、度m，n判斷m與n是否相等，如果相等用for循環(huán)求插值基函數(shù)的值在命令窗口中調用M文件，求出插值點對應的函數(shù)值結 束三、實驗求解結果時的函數(shù)近似值分別為：0.8265 ,0.8268 ,0.8231。附 錄A題實驗程序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'); %散點數(shù)據(jù)的插值函數(shù)griddata()figure(1)plot(x0,y0,'r*')figure(2)plot3(x0,y0,z0,'r*'); %已知節(jié)點的三維散點圖hold onsurf(x,y,z) %畫出經插值的圖形figure(3)c,h=contour(x,y,z); %畫等高線clabel(c,h) %標明等高線的高程，高程為5的區(qū)域避免進入B題實驗程序clcclear% 數(shù)據(jù)文件（dat

15、a）的讀取v d = xlsread('data.xls');x = v(2,:);y = v(3,:);n = length(x); % 讀取文件中的點的個數(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')% 誤差的計算y1 = a+b.*x;format longw1 = sum(abs(y1-y)/sum(y)% 指數(shù)模型的擬合% y=A*exp(B/x)% matlab內部函數(shù)的命令a

16、1 = polyfit(x,y,1);% 誤差的計算y2 = a1(1).*x+a1(2);w2 = sum(abs(y2-y)/sum(y)C題實驗程序function = Runge(a1,a2,n)% (1)構造龍格現(xiàn)象% n:構造n次的多項式插值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% 插值多項式x = -5:0.01:5;yy = eval(y);plot(x,yy)% 龍格現(xiàn)象圖 Endclcclear% c題問題主程序% 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% 畫圖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('分段線性圖 紅色：插值 藍色：原值')figure % 分段三次插值for i = 1:a(2) x = x1(i):(a2-a1)/(10*a(2):x1