黃河小浪底調(diào)水調(diào)沙工程數(shù)學(xué)實(shí)驗(yàn)實(shí)驗(yàn)報(bào)告

1、數(shù)學(xué)實(shí)驗(yàn)實(shí)驗(yàn)報(bào)告題目：黃河小浪底調(diào)水調(diào)沙工程姓名： 胡 迪 學(xué)號(hào)： 專業(yè)：信息與計(jì)算科學(xué)黃河小浪底調(diào)水調(diào)沙問(wèn)題2004年6月至7月黃河進(jìn)行了第三次調(diào)水調(diào)沙試驗(yàn)，特別是首次由小浪底、三門峽和萬(wàn)家寨三大水庫(kù)聯(lián)合調(diào)度，采用接力式防洪預(yù)泄放水，形成人造洪峰進(jìn)行調(diào)沙試驗(yàn)獲得成功。整個(gè)試驗(yàn)期為20多天，小浪底從6月19日開始預(yù)泄放水，至到7月13日恢復(fù)正常供水結(jié)束。小浪底水利工程按設(shè)計(jì)攔沙量為75.5億m3,在這之前，小浪底共積泥沙達(dá)14.15億t。這次調(diào)水調(diào)沙試驗(yàn)一個(gè)重要的目的就是由小浪底上游的三門峽和萬(wàn)家寨水庫(kù)泄洪，在小浪底形成人造洪峰，沖刷小浪底庫(kù)區(qū)沉積的泥沙，在小浪底水庫(kù)開閘泄洪以后，從6月27日

2、開始三門峽水庫(kù)和萬(wàn)家寨水庫(kù)陸續(xù)開閘放水，人造洪峰于29日先后到達(dá)小浪底，7月3日達(dá)到最大流量2700 ，使小浪底水庫(kù)的排沙量也不斷地增加。表1是由小浪底觀測(cè)站從6月29日到7月10日檢測(cè)到的試驗(yàn)數(shù)據(jù)。 表1 試驗(yàn)觀測(cè)數(shù)據(jù) （ 單位：水流為，含沙量為）日期6.296.307.17.27.37.4時(shí)間8：0020:008：0020:008：0020:008：0020:008：0020:008：0020:00水流量180019002100220023002400250026002650270027202650含沙量326075859098100102108112115116日期7.57.67.77.

3、87.97.10時(shí)間8：0020:008：0020:008：0020:008：0020:008：0020:008：0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085注：以上數(shù)據(jù)主要是根據(jù)媒體公開報(bào)道的結(jié)果整理而成?，F(xiàn)在，根據(jù)試驗(yàn)數(shù)據(jù)建立數(shù)學(xué)模型研究下面的問(wèn)題：（1）給出估算任意時(shí)刻的排沙量及總排沙量的方法；（2）確定排沙量與水流量的變化關(guān)系。關(guān)鍵詞：擬合，SAS，Matlab，線性回歸，調(diào)水調(diào)沙實(shí)驗(yàn)問(wèn)題分析：1、對(duì)于問(wèn)題一，所給數(shù)據(jù)中水流量x和含沙量h的乘積即為該時(shí)刻的排沙

4、量y即：y=hx。2、對(duì)于問(wèn)題二，研究排沙量與排水量的關(guān)系，從實(shí)驗(yàn)數(shù)據(jù)中可以看出，開始排沙量隨水量增加而增加，而后隨水流量的增加而減少，顯然變化關(guān)系并非線性的關(guān)系，為此，把問(wèn)題分為兩部分，從水流量增加到最大值為第一階段，從水流量最大值到結(jié)束為第二階段，分別來(lái)研究水流量與排沙量之間的函數(shù)關(guān)系。模型假設(shè)：1、水流量和排沙量都是連續(xù)的，不考慮上游泄洪所帶來(lái)的含沙量和外界帶來(lái)的含沙量。2、時(shí)間是連續(xù)變化的，所取時(shí)間點(diǎn)依次為1,2,3，,24,單位時(shí)間為12h。模型的建立與求解：對(duì)于問(wèn)題一，因?yàn)榕派沉颗c時(shí)間的散點(diǎn)圖基本符合正態(tài)曲線，如圖二所示。所以，排沙量的對(duì)數(shù)與時(shí)間的函數(shù)關(guān)系就應(yīng)該符合二次函數(shù)關(guān)系，因

5、而排沙量取對(duì)數(shù)后，再與時(shí)間t進(jìn)行二次回歸，排沙量取自然后的數(shù)據(jù)見表2.假設(shè)排沙量與時(shí)間函數(shù)關(guān)系的數(shù)學(xué)模型是兩邊取對(duì)數(shù)得Lny=at2+bt+c先由表二做出排沙量的自然對(duì)數(shù)lny與時(shí)間t的散點(diǎn)圖見圖一，并利用SAS軟件進(jìn)行擬合，得到排沙量的自然對(duì)數(shù)與時(shí)間的回歸方程為：Lny=-0.0209t2+0.4298t+10.6321由回歸擬合參數(shù)表可知回歸方程是顯著的，因?yàn)橄嚓P(guān)系數(shù)人R2=0.9629,誤差均方S2=0.0543，說(shuō)明回歸曲線擬合效果很好。所以排沙量與時(shí)間之間的函數(shù)關(guān)系式為圖二：排沙量對(duì)時(shí)間的曲線圖時(shí)間點(diǎn)11800325760010.96128219006011.643953210075

6、11.96718422008512.13886523009012.24047624009812.368197250010012.429228260010212.488249265010812.5644510270011212.6195111272011512.6533212265011612.6359113260011812.6339514250012012.6115415230011812.5113516220010512.350171720008011.982931818506011.61729191820509100011.41864201800305400010.89674211750

7、264550010.72547221500203000010.30895231000880008.24900545008.最后對(duì)所求出的函數(shù)關(guān)系在區(qū)間0,24之間進(jìn)行積分結(jié)果為總排沙量1.93962億噸，此與媒體報(bào)道的排沙量幾乎一樣。對(duì)于第二個(gè)問(wèn)題，兩個(gè)階段的數(shù)據(jù)如表三、表四所示表三：第一階段試驗(yàn)數(shù)據(jù)序號(hào)12345678水流量x18001900210022002300240025002600含沙量h326075859098100102表四：第二階段的試驗(yàn)觀測(cè)數(shù)據(jù)序號(hào)12345678水流量x26502600250023002200200018501820含沙量h11611812011810580

8、6050對(duì)于第一階段，有表四用MATLAB作圖（如圖三）可以看出其變化趨勢(shì)，我們用多項(xiàng)式做最小二乘擬合。設(shè)三次擬合函數(shù)關(guān)系h=a0+a1x+a2x2+a3x3其中a0,a1,a2,a3,為待定系數(shù)。四次擬合函數(shù)關(guān)系h= a0+a1x+a2x2+a3x3+a4x4其中a0,a1,a2,a3,a4為待定系數(shù)。圖三：第一階段水流量與排沙量之間的關(guān)系圖三次多項(xiàng)式擬合由MATLAB擬合函數(shù)求解出a0=a1=0,a2=0.0032,a3=-2.4929.則擬合函數(shù)h=0.0032x2-2.4929x3,擬合效果如圖四所示圖四：三次多項(xiàng)式擬合效果，紅線為擬合曲線類似的四次多項(xiàng)式擬合由MATLAB擬合函數(shù)求解

9、出a0=a1=a2=0,a3=0.0121,a4= -7.4347則擬合函數(shù)h=0.0121x3-7.4347x4,擬合效果如圖五所示圖五：四次多項(xiàng)式擬合效果，藍(lán)線線為擬合曲線對(duì)于第二階段，有表五用MATLAB作圖可以看出其變化趨勢(shì)，我們用多項(xiàng)式做最小二乘擬合。設(shè)三次擬合函數(shù)關(guān)系h=a0+a1x+a2x2+a3x3其中a0,a1,a2,a3,為待定系數(shù)。四次擬合函數(shù)關(guān)系h= a0+a1x+a2x2+a3x3+a4x4其中a0,a1,a2,a3,a4為待定系數(shù)。三次多項(xiàng)式擬合由MATLAB擬合函數(shù)求解出a0=a1=0,a2=-0.9475,a3= 464.9601.則擬合函數(shù)h=-0.9475x

10、2+464.9601x3,擬合效果如圖圖六所示圖六：三次擬合函數(shù)擬合效果類似的四次多項(xiàng)式擬合由MATLAB擬合函數(shù)求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-354.5952則擬合函數(shù)h=-0.0013x2+1.1219x3-354.5952x4,擬合效果如圖七所示圖七：四次擬合函數(shù)擬合效果結(jié)論以及分析檢驗(yàn)：用SAS軟件做線性回歸得到排沙量與時(shí)間的函數(shù)關(guān)系式為:再利用所求函數(shù)在區(qū)間0,24上進(jìn)行積分得到總排沙量1.93962億噸，這與現(xiàn)實(shí)情況基本相符。對(duì)于第一階段三次多項(xiàng)式擬合由MATLAB擬合函數(shù)求解a0=a1=0,a2=0.0032,a3=-2.4929則擬合

11、函數(shù)h=0.0032x2-2.4929x3對(duì)于第一階段四次多項(xiàng)式擬合由MATLAB擬合函數(shù)求解出a0=a1=a2=0,a3=0.0121,a4= -7.4347則擬合函數(shù)h=0.0121x3-7.4347x4對(duì)于第二階段三次多項(xiàng)式擬合由MATLAB擬合函數(shù)求解出a0=a1=0,a2=-0.9475,a3= 464.9601則擬合函數(shù)h=-0.9475x2+464.9601x3對(duì)于第二階段四次多項(xiàng)式擬合由MATLAB擬合函數(shù)求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-354.5952則擬合函數(shù)h=-0.0013x2+1.1219x3-354.5952x4討論與推廣：

12、1、對(duì)于第一個(gè)問(wèn)題排沙量與時(shí)間不是嚴(yán)格的正態(tài)函數(shù)關(guān)系可能與實(shí)際有些偏差，此外還可以用SAS軟件進(jìn)行高次的多向式回歸2、對(duì)于第二個(gè)問(wèn)題，由于MATLAB軟件的計(jì)算可能有些偏差導(dǎo)致擬合的函數(shù)關(guān)系可能與實(shí)際有稍微偏差，此外，還可以進(jìn)行高次的擬合。附錄：1、排沙量與時(shí)間的關(guān)系圖像的MATLAB程序：t=1:1:24;y=57600,2400,91000,54000,45500,30000,8000,4500; plot(t,y,r)2、對(duì)排沙量求自然對(duì)數(shù)的MATLAB程序與結(jié)果：y3=log(y)y3 = Columns 1 through 17 10.9613 11.6440 11.9672 12.

13、1389 12.2405 12.3682 12.4292 12.4882 12.5644 12.6195 12.6533 12.6359 12.6340 12.6115 12.5113 12.3502 11.9829 Columns 18 through 24 11.6173 11.4186 10.8967 10.7255 10.3090 8.9872 8.41183、第一階段的排沙量與水流量之間的關(guān)系MATLAB程序：x=1800,1900,2100,2200,2300,2400,2500,2600,2650,2700,2720; h=32,60,75,85,90,98,100,102,10

14、8,112,115; x1=2650,2600,2500,2300,2200,2000,1850,1820,1800,1750,1500,1000,900; h1=116,118,120,118,105,80,60,50,40,32,20,8,5; plot(x,h,r:)4、第一階段三次多項(xiàng)式擬合函數(shù)以及擬合效果程序與結(jié)果： A1=polyfit(x,h,3) in polyfit at 80A1 = 0.0000 -0.0000 0.0032 -2.4929 z1=polyval(A1,x);plot(x,h,k+,x,z1,r)5、第一階段四次多項(xiàng)式擬合函數(shù)以及擬合效果程序與結(jié)果：A2=polyfit(x,h,4) In polyfit at 80A2 = -0.0000 0.0000 -0.0000 0.0121 -7.4347 z2=polyval(A2,x); plot(x,h,*,x,z2,r)6、第二階段三次多項(xiàng)式擬合函數(shù)以及擬合效果程序與結(jié)果：A3=polyfit(x1,h1,3) In polyfit at 80A3 = -0.0000 0.0006 -0.9475 464.9601 z3=polyval(A3,x1); plot(x,h,*,x1,z3,b)7、第二階段