數(shù)學(xué)模型與數(shù)學(xué)建模實驗五

1、實驗報告五學(xué)院名稱：理學(xué)院 專業(yè)年級： 姓 名： 學(xué) 號：課 程：數(shù)學(xué)模型與數(shù)學(xué)建模 報告日期：2015年12月8日一、實驗題目例2.2.1 水庫庫容量與高程 設(shè)一水庫將河道分為上、下游兩個河段，降雨的開始時刻為8時，這是水位的高程為168m，水庫容量為，預(yù)測上游的流量，d取值如表2.2.1所示。 表2.2.1 上有流量的預(yù)測812162430444856360054007800920010100350025001600已知水庫中水的容量與水位高程H（m）的數(shù)值關(guān)系為表2.2.2 表2.2.2 水庫庫容量與水位高程的關(guān)系23.9324.0624.1224.3324.4724.624.75168

2、.75168.8168.85168.9168.95169169.05 如果當日從8時開始，水一直保持的泄流量，根據(jù)所給數(shù)據(jù)，預(yù)報從降雨時刻到56h以內(nèi)每小時整點時刻水庫中水的庫容量與水位高程。例2.2.2 地下含沙量 某地區(qū)有優(yōu)質(zhì)細沙埋在地下，某公司擬在此處采沙，已得到該地區(qū)鉆探資料圖的一角如下表，在每個格點上有三個數(shù)字列，都是相對于選定基點的高度（m），最上面的數(shù)字是覆蓋表面的標高，中間的數(shù)字是沙層頂部的標高最下面的數(shù)字是沙層底部的標高，每個格子都是正方形，邊長50m。畫星號處，即沼澤表層地帶，沒有鉆探數(shù)據(jù)。試估計整個矩形區(qū)域內(nèi)的含沙量。ABCDEFGH022.420.05.8*22.518

3、.40.523.017.80.413.218.00.4123.019.96.023.120.03.223.220.01.623.419.81.023.519.91.124.020.01.024.019.80.824.019.40.9223.119.82.223.219.71.423.419.40.623.420.00.523.520.10.324.220.3-0.224.120.3-0.124.120.50.0二、實驗?zāi)康?插值模型是數(shù)據(jù)挖掘的另一類模型，插值（Interpolation）的目的是根據(jù)能夠獲得的觀測數(shù)據(jù)推測缺損的數(shù)據(jù)，此時觀測數(shù)據(jù)被視為精確的基準數(shù)據(jù)，尋找一個至少滿足條件的函數(shù)

4、，使得，在本節(jié)我們強調(diào)的是插值模型的應(yīng)用，而不是插值方法的構(gòu)造。三、問題陳述2.2.1 一維插值 例2.2.1 水庫庫容量與高程2.2.2 二維插值 例2.2.2 地下含沙量2.2.3 泛克里金插值四、模型及求解結(jié)果 2.2.1 一維插值 一元函數(shù)差值公式為其中是滿足條件的函數(shù)，依據(jù)插值的公式，如最近鄰差值，線性插值、分段三次Hermite插值等，分別取階梯函數(shù)、線性函數(shù)、三次多項式函數(shù)等，相應(yīng)的數(shù)學(xué)表達式可以查閱本科生數(shù)值計算教材。下面先通過簡單的Matlab一維插值令interpret1了解相應(yīng)的計算結(jié)果。例：2.2.1 水庫庫容量與高程 為了給出每小時的報告，需要補充每小時整點時刻上有流

5、量的數(shù)據(jù)，以及相應(yīng)不同庫容量的水位高程。 假設(shè)（a）已知數(shù)據(jù)準確（b）相鄰兩個時刻之間的流量變化是現(xiàn)行的（c）相鄰兩個水位高程之間的高程對水的庫容量的變化也是線性的 首先，利用Matlab線性插值令，確定每小時的上游流量q（t）.由程序在Matlab中運行的結(jié)果為：q = 1.0e+004 * Columns 1 through 9 0.3600 0.4050 0.4500 0.4950 0.5400 0.6000 0.6600 0.7200 0.7800 Columns 10 through 18 0.7975 0.8150 0.8325 0.8500 0.8675 0.8850 0.902

6、5 0.9200 0.9350 Columns 19 through 27 0.9500 0.9650 0.9800 0.9950 1.0100 0.9629 0.9157 0.8686 0.8214 Columns 28 through 36 0.7743 0.7271 0.6800 0.6329 0.5857 0.5386 0.4914 0.4443 0.3971 Columns 37 through 45 0.3500 0.3000 0.2500 0.2410 0.2320 0.2230 0.2140 0.2050 0.1960 Columns 46 through 49 0.1870

7、0.1780 0.1690 0.1600 然后確定每個時刻t的水庫容量，因為，水庫容量=原庫存量+流入量-泄流量，即：這里我們遇到數(shù)值積分，被積函數(shù)沒有解析表達式，只有一個數(shù)列表示，表示在i整點時刻的流量，利用Matlab逐點積分指令，可以得到水庫容量在每一刻的值。 最后確定每時刻t水庫的水位高程h（t），因為最大水庫容量已經(jīng)遠遠超出了已知數(shù)據(jù)范圍，需要利用外插方法補充數(shù)據(jù)，確定水庫高程對水庫容量的依賴關(guān)系h=H（v）。最后利用函數(shù)復(fù)合得到水位高程確定每小時的上游流量q（t）利用函數(shù)復(fù)合得到水位高程2.2.2 二維插值二維插值大致可以分為兩種，規(guī)則點插值和散亂點插值，前者即利用在方形網(wǎng)格點給定

8、的數(shù)據(jù)，在加密網(wǎng)格點上補充相應(yīng)函數(shù)值，插值公式為其中是滿足的二元函數(shù)，如雙線性函數(shù)、雙三次多項式函數(shù)、樣條函數(shù)等，一句不同的插值方式選取，相應(yīng)的表達式可以從本科生數(shù)值計算教材中查閱。這里先通過Matlab二維插值指令interp2了解相應(yīng)的計算結(jié)果散亂點插值是要利用在不規(guī)則排列的觀測點上的數(shù)據(jù)確定其他點上的函數(shù)值，插值公式形式為常用到的插值方法，如距離加權(quán)反比插值，Kriging插值，需要查閱專業(yè)書籍，下面先通過Matlab二維插值指令griddata了解相應(yīng)的計算結(jié)果例2.2.2 地下含沙量 假設(shè)地下沙層連續(xù)變化，于是可以利用積分運算矩形區(qū)域內(nèi)的含沙量。首先通過插值補充缺失的數(shù)據(jù)，采用一維樣

9、條插值。為了提高梯形積分精度，利用二維樣條插值加密數(shù)據(jù)，得到邊長為0.5m的長方形網(wǎng)格上的數(shù)據(jù)。最后，將二重積分化為累次積分，利用梯形公式得到積分的近似值。v = 6.4290e+005如果在一維插值補齊缺失數(shù)據(jù)中用線性插值代替上面的樣條插值，將得到沙層體積634525，和以上結(jié)果有些差別。 對于加密二維網(wǎng)格，應(yīng)該加到多密？注意，從理論上講，網(wǎng)格越密計算精度越高，但是從計算角度看，網(wǎng)格越密計算誤差累計也越大，所以需要通過逐步加密網(wǎng)格，確定精度最優(yōu)的積分值。 如果直接利用二維插值，例如采用Matlab的散亂點插值指令，仍用二次梯形公式，會得到相近的結(jié)果。v = 6.3572e+0052.2.3

10、泛克里金插值考慮到沙子是一種特殊的地址，下面試用地質(zhì)學(xué)常用的泛克里金插值方法計算。先介紹克里金插值方法，這本身就是運用統(tǒng)計學(xué)方法建模的一個有趣的例子。1951年南非地質(zhì)學(xué)家Krige講處礦藏的貯藏量看成是隨機函數(shù)的一個實現(xiàn)，提出了依據(jù)觀測值，尋找基函數(shù)，以獲得的具有形式 的最小方差無偏估計，取的條件期望 給出未知點估值的插值方法，即考慮形如 的隨機函數(shù)，其中在地質(zhì)學(xué)中稱為飄移，在一些隨機分析的場合也成為趨勢，是均值為的未知的隨機變量，是多項式基函數(shù)，例如這里取，是滿足，的隨機函數(shù)是給定的且滿足當當?shù)暮撕瘮?shù)。這里按通常取法，取高斯核函數(shù)根據(jù)無偏估計要求， 即： 從而得到無偏條件，對任意的有： 在

11、這個約束條件下極小化方差 =：其中，。引入拉格朗日乘子，由拉格朗日函數(shù)的極小值點必為其駐點的結(jié)論，得到關(guān)系式 其中，于是，在未知點上的隨機函數(shù)值可以用它的最小方差無偏估計的條件期望近似，從而得到泛克里金插值函數(shù)： = =其中練習(xí)：證明上面構(gòu)造的泛克里金函數(shù)是連續(xù)的，且通過已知點，即插值函數(shù)滿足： 用泛克里金插值得到的邊長為0.5m的加密網(wǎng)格上的沙層厚度，通過二次梯形公式，得到沙層體積近似值為637320五、程序代碼2.2.1 一維插值>> x=0:4:20;>> y=37 51 45 74 83 88;>> xx=0:1:20;>> y1=int

12、erp1(x,y,xx,'nearest'); %最近鄰插值，間斷函數(shù)>> y2=interp1(x,y,xx); %線性插值，連續(xù)函數(shù)>> y3=interp1(x,y,xx,'cubic'); %分段三次Hermite插值，一階連續(xù)函數(shù)>> y4=interp1(x,y,xx,'splinet'); %樣條插值，二階倒數(shù)連續(xù)的函數(shù)>> subplot(2,2,1),plot(x,y,'kd',xx,y1),title('nearest')>>

13、 subplot(2,2,2),plot(x,y,'kd',xx,y2),title('linear')>> subplot(2,2,3),plot(x,y,'kd',xx,y3),title('cubic')>> subplot(2,2,4),plot(x,y,'kd',xx,y4),title('spline')例2.2.1 水庫庫容量與高程利用Matlab線性插值令，確定每小時的上游流量q（t）>> T=8,12,16,24,3

14、0,44,46,56;>> Q=3600,5400,7800,9200,10100,3500,2500,1600;>> t=8:56;>> q=interp1(T,Q,t,'linear')%得到每小時上游流量>> plot(T,Q,'kd',t,q),title('time-flow')水庫容量在每一刻的值>> v=21.9+36*10(-6)*cumtrapz(t,q-103*ones(size(t); %每小時水庫容量>> vmax=max(v); %得到最大的水庫容量

15、為30.7524（108立方米）利用函數(shù)復(fù)合得到水位高程>> V=21.9 23.93 24.06 24.12 24.33 24.47 24.3 24.75;>> H=168 168.75 168.8 168.85 168.9 168.95 169 169.05;>> h=interp1(V,H,v,'linear','extrap'); %得到每小時水位高程>> hmax=max(h);%最大水位高程171.0508米>> plot(V,H,'kd',v,h),title('v

16、olume-atitude')2.2.2 二維插值二維插值指令>> x=0:4:16;y=0:4:16;>> z=620 730 800 850 870;760 880 970 720 1050880 1080 630 1250 1280;980 1180 1320 1450 14201060 1230 1390 1500 1500;>> X,Y=meshgrid(0:16,0:16);>> Z1=interp2(x,y,z,X,Y,'nearest');Z2=interp2(x,y,z,X,Y);>> Z3=

17、interp2(x,y,z,X,Y,'cubic');Z4=interp2(x,y,z,X,Y,'spline');>> subplot(2,2,1),surfc(X,Y,Z1),title('nearest')>> subplot(2,2,2),surfc(X,Y,Z2),title('linear')>> subplot(2,2,3),surfc(X,Y,Z3),title('cubic')>> subplot(2,2,4),surfc(X,Y,Z4),title

18、('spline')散亂點插值>> x=1 2 3 4 5;>> y=4 1 3 2 5;>> z=4 1 6 2 5;>> X,Y=meshgrid(0:0.2:5);>> Z1=griddata(x,y,z,X,Y,'nearest');>> Z2=griddata(x,y,z,X,Y);>> Z3=griddata(x,y,z,X,Y,'cubic');>> Z4=griddata(x,y,z,X,Y,'v4');>>

19、 subplot(2,2,1),surfc(X,Y,Z1),title('nearest')>> subplot(2,2,2),surfc(X,Y,Z2),title('linear')>> subplot(2,2,3),surfc(X,Y,Z3),title('cubic')>> subplot(2,2,4),surfc(X,Y,Z4),title('v4')例2.2.2 地下含沙量通過一維插值補充缺失的數(shù)據(jù)>> D=23.0,23.1,23.2,23.4,23.5,24.0,24

20、.0,24.0;23.1,23.3,23.4,23.4,23.5,24.2,24.1,24.1;>> E=19.9,20.0,20.0,19.8,19.9,20.0,19.8,19.6;19.8,19.7,19.4,20.0,20.1,20.3,20.3,20.5;>> F=6.0,3.2,1.6,1.0,1.1,1.0,0.8,0.9;2.2,1.4,0.6,0.5,0.3,-0.2,-0.1,0.0;>> P=1,6,7,8;K=1:8;>> A=22.4,22.5,23.0,23.2;>> A1=interp1(P,A,K,&#

21、39;spline');a=A1;D>> B=20.0,18.4,17.8,18.0;>> B1=interp1(P,B,K,'spline');b=B1;E>> C=5.8,0.5,0.4,0.4;>> C1=interp1(P,C,K,'spline');c=C1;F二維樣條插值加密數(shù)據(jù)，得到邊長為0.5m的長方形網(wǎng)格上的數(shù)據(jù)>> M,N=meshgrid(1:8,0:2);>> delta=0.01;X,Y=meshgrid(1:delta:8,0:delta:2);>&

22、gt; a1=interp2(M,N,a,X,Y,'spline');>> b1=interp2(M,N,b,X,Y,'spline');>> c1=interp2(M,N,c,X,Y,'spline');>> mesh(X,Y,a1),hold on,>> mesh(X,Y,b1),hold on,mesh(X,Y,c1)最后，將二重積分化為累次積分，利用梯形公式得到積分的近似值。>> f=b1-c1;>> v=trapz(trapz(f,1),2)*delta2*2500

23、v = 6.4290e+005（v為沙層體積） 如果直接利用二維插值，例如采用Matlab的散亂點插值指令，仍用二次梯形公式，會得到相近的結(jié)果。>> x0=1 6 7 8 1:8 1:8;>> y0=zeros(1,4) ones(1,8) 2*ones(1,8);>> f0=B E'-C F'>> f=griddata(x0,y0,f0,X,Y,'cubic');>>>> v=trapz(trapz(f,1),2)*delta2*2500v = 6.3572e+005（v為沙層體積）2.2

24、.3 泛克里金插值泛克里金插值的MATLAB指令為：function f=krige2(x0,y0,f0,X,Y)%x0,y0,f0·Ö±ðÊÇÒÑÖªÊý¾ÝµÄºá¡¢×Ý×ø±êºÍ¹Û²âÖµ£¬±í´

25、ï³ÉÐÐÏòÁ¿%X,Y²åÖµµãµÄºá¡¢×Ý×ø±ê£¬¿ÉÒÔÊÇÐÐÏòÁ¿£¬Ò²¿ÉÒÔ&#

26、202;ǾØÕó±í´ïÍø¸ñµã×ø±ê¡£%Êä³öfÊǶÔÓ¦X,YµÄº¯ÊýÖµn=length(x0);syms x y;%¶¨Òå·ûºÅ±äÁ¿h0=(x-x0').2+(y-y0').2;%¾àÀ뺯