數(shù)學(xué)模型與數(shù)學(xué)建模 實(shí)驗(yàn)五

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

2、.75168.8168.85168.9168.95169169.05 如果當(dāng)日從8時(shí)開始，水一直保持的泄流量，根據(jù)所給數(shù)據(jù)，預(yù)報(bào)從降雨時(shí)刻到56h以內(nèi)每小時(shí)整點(diǎn)時(shí)刻水庫(kù)中水的庫(kù)容量與水位高程。例2.2.2 地下含沙量 某地區(qū)有優(yōu)質(zhì)細(xì)沙埋在地下，某公司擬在此處采沙，已得到該地區(qū)鉆探資料圖的一角如下表，在每個(gè)格點(diǎn)上有三個(gè)數(shù)字列，都是相對(duì)于選定基點(diǎn)的高度（m），最上面的數(shù)字是覆蓋表面的標(biāo)高，中間的數(shù)字是沙層頂部的標(biāo)高最下面的數(shù)字是沙層底部的標(biāo)高，每個(gè)格子都是正方形，邊長(zhǎng)50m。畫星號(hào)處，即沼澤表層地帶，沒有鉆探數(shù)據(jù)。試估計(jì)整個(gè)矩形區(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二、實(shí)驗(yàn)?zāi)康?插值模型是數(shù)據(jù)挖掘的另一類模型，插值（Interpolation）的目的是根據(jù)能夠獲得的觀測(cè)數(shù)據(jù)推測(cè)缺損的數(shù)據(jù)，此時(shí)觀測(cè)數(shù)據(jù)被視為精確的基準(zhǔn)數(shù)據(jù)，尋找一個(gè)至少滿足條件的函數(shù)

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

5、量的數(shù)據(jù)，以及相應(yīng)不同庫(kù)容量的水位高程。 假設(shè)（a）已知數(shù)據(jù)準(zhǔn)確（b）相鄰兩個(gè)時(shí)刻之間的流量變化是現(xiàn)行的（c）相鄰兩個(gè)水位高程之間的高程對(duì)水的庫(kù)容量的變化也是線性的 首先，利用Matlab線性插值令，確定每小時(shí)的上游流量q（t）.由程序在Matlab中運(yùn)行的結(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 然后確定每個(gè)時(shí)刻t的水庫(kù)容量，因?yàn)?，水?kù)容量=原庫(kù)存量+流入量-泄流量，即：這里我們遇到數(shù)值積分，被積函數(shù)沒有解析表達(dá)式，只有一個(gè)數(shù)列表示，表示在i整點(diǎn)時(shí)刻的流量，利用Matlab逐點(diǎn)積分指令，可以得到水庫(kù)容量在每一刻的值。 最后確定每時(shí)刻t水庫(kù)的水位高程h（t），因?yàn)樽畲笏畮?kù)容量已經(jīng)遠(yuǎn)遠(yuǎn)超出了已知數(shù)據(jù)范圍，需要利用外插方法補(bǔ)充數(shù)據(jù)，確定水庫(kù)高程對(duì)水庫(kù)容量的依賴關(guān)系h=H（v）。最后利用函數(shù)復(fù)合得到水位高程確定每小時(shí)的上游流量q（t）利用函數(shù)復(fù)合得到水位高程2.2.2 二維插值二維插值大致可以分為兩種，規(guī)則點(diǎn)插值和散亂點(diǎn)插值，前者即利用在方形網(wǎng)格點(diǎn)給定

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

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

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

11、這個(gè)約束條件下極小化方差 =：其中，。引入拉格朗日乘子，由拉格朗日函數(shù)的極小值點(diǎn)必為其駐點(diǎn)的結(jié)論，得到關(guān)系式 其中，于是，在未知點(diǎn)上的隨機(jī)函數(shù)值可以用它的最小方差無(wú)偏估計(jì)的條件期望近似，從而得到泛克里金插值函數(shù)： = =其中練習(xí)：證明上面構(gòu)造的泛克里金函數(shù)是連續(xù)的，且通過已知點(diǎn)，即插值函數(shù)滿足： 用泛克里金插值得到的邊長(zhǎng)為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 水庫(kù)庫(kù)容量與高程利用Matlab線性插值令，確定每小時(shí)的上游流量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')%得到每小時(shí)上游流量>> plot(T,Q,'kd',t,q),title('time-flow')水庫(kù)容量在每一刻的值>> v=21.9+36*10(-6)*cumtrapz(t,q-103*ones(size(t); %每小時(shí)水庫(kù)容量>> vmax=max(v); %得到最大的水庫(kù)容量

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'); %得到每小時(shí)水位高程>> 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')散亂點(diǎn)插值>> 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 地下含沙量通過一維插值補(bǔ)充缺失的數(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ù)，得到邊長(zhǎng)為0.5m的長(zhǎng)方形網(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的散亂點(diǎn)插值指令，仍用二次梯形公式，會(huì)得到相近的結(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;%¾àÀ뺯