黄河小浪底调水调沙工程数学实验实验报告(共10页).doc
精选优质文档-倾情为你奉上数学实验实验报告题目：黄河小浪底调水调沙工程姓名： 胡 迪 学号： 专业：信息与计算科学黄河小浪底调水调沙问题2004年6月至7月黄河进行了第三次调水调沙试验，特别是首次由小浪底、三门峡和万家寨三大水库联合调度，采用接力式防洪预泄放水，形成人造洪峰进行调沙试验获得成功。整个试验期为20多天，小浪底从6月19日开始预泄放水，至到7月13日恢复正常供水结束。小浪底水利工程按设计拦沙量为75.5亿m3,在这之前，小浪底共积泥沙达14.15亿t。这次调水调沙试验一个重要的目的就是由小浪底上游的三门峡和万家寨水库泄洪，在小浪底形成人造洪峰，冲刷小浪底库区沉积的泥沙，在小浪底水库开闸泄洪以后，从6月27日开始三门峡水库和万家寨水库陆续开闸放水，人造洪峰于29日先后到达小浪底，7月3日达到最大流量2700 ，使小浪底水库的排沙量也不断地增加。表1是由小浪底观测站从6月29日到7月10日检测到的试验数据。 表1 试验观测数据 （ 单位：水流为，含沙量为）日期6.296.307.17.27.37.4时间8：0020:008：0020:008：0020:008：0020:008：0020:008：0020:00水流量180019002100220023002400250026002650270027202650含沙量326075859098100102108112115116日期7.57.67.77.87.97.10时间8：0020:008：0020:008：0020:008：0020:008：0020:008：0020:00水流量26002500230022002000185018201800175015001000900含沙量11812011810580605030262085注：以上数据主要是根据媒体公开报道的结果整理而成。现在，根据试验数据建立数学模型研究下面的问题：（1）给出估算任意时刻的排沙量及总排沙量的方法；（2）确定排沙量与水流量的变化关系。关键词：拟合，SAS，Matlab，线性回归，调水调沙实验问题分析：1、对于问题一，所给数据中水流量x和含沙量h的乘积即为该时刻的排沙量y即：y=hx。2、对于问题二，研究排沙量与排水量的关系，从实验数据中可以看出，开始排沙量随水量增加而增加，而后随水流量的增加而减少，显然变化关系并非线性的关系，为此，把问题分为两部分，从水流量增加到最大值为第一阶段，从水流量最大值到结束为第二阶段，分别来研究水流量与排沙量之间的函数关系。模型假设：1、水流量和排沙量都是连续的，不考虑上游泄洪所带来的含沙量和外界带来的含沙量。2、时间是连续变化的，所取时间点依次为1,2,3，,24,单位时间为12h。模型的建立与求解：<一>对于问题一，因为排沙量与时间的散点图基本符合正态曲线，如图二所示。所以，排沙量的对数与时间的函数关系就应该符合二次函数关系，因而排沙量取对数后，再与时间t进行二次回归，排沙量取自然后的数据见表2.假设排沙量与时间函数关系的数学模型是两边取对数得Lny=at2+bt+c先由表二做出排沙量的自然对数lny与时间t的散点图见图一，并利用SAS软件进行拟合，得到排沙量的自然对数与时间的回归方程为：Lny=-0.0209t2+0.4298t+10.6321由回归拟合参数表可知回归方程是显著的，因为相关系数人R2=0.9629,误差均方S2=0.0543，说明回归曲线拟合效果很好。所以排沙量与时间之间的函数关系式为图二：排沙量对时间的曲线图时间点11800325760010.96128219006011.64395321007511.96718422008512.13886523009012.24047624009812.368197250010012.429228260010212.488249265010812.5644510270011212.6195111272011512.6533212265011612.6359113260011812.6339514250012012.6115415230011812.5113516220010512.350171720008011.982931818506011.61729191820509100011.41864201800305400010.89674211750264550010.72547221500203000010.30895231000880008.24900545008.最后对所求出的函数关系在区间0,24之间进行积分结果为总排沙量1.93962亿吨，此与媒体报道的排沙量几乎一样。<二>对于第二个问题，两个阶段的数据如表三、表四所示表三：第一阶段试验数据序号12345678水流量x18001900210022002300240025002600含沙量h326075859098100102表四：第二阶段的试验观测数据序号12345678水流量x26502600250023002200200018501820含沙量h116118120118105806050对于第一阶段，有表四用MATLAB作图（如图三）可以看出其变化趋势，我们用多项式做最小二乘拟合。设三次拟合函数关系h=a0+a1x+a2x2+a3x3其中a0,a1,a2,a3,为待定系数。四次拟合函数关系h= a0+a1x+a2x2+a3x3+a4x4其中a0,a1,a2,a3,a4为待定系数。图三：第一阶段水流量与排沙量之间的关系图三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=0.0032,a3=-2.4929.则拟合函数h=0.0032x2-2.4929x3,拟合效果如图四所示图四：三次多项式拟合效果，红线为拟合曲线类似的四次多项式拟合由MATLAB拟合函数求解出a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合函数h=0.0121x3-7.4347x4,拟合效果如图五所示图五：四次多项式拟合效果，蓝线线为拟合曲线对于第二阶段，有表五用MATLAB作图可以看出其变化趋势，我们用多项式做最小二乘拟合。设三次拟合函数关系h=a0+a1x+a2x2+a3x3其中a0,a1,a2,a3,为待定系数。四次拟合函数关系h= a0+a1x+a2x2+a3x3+a4x4其中a0,a1,a2,a3,a4为待定系数。三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-0.9475,a3= 464.9601.则拟合函数h=-0.9475x2+464.9601x3,拟合效果如图图六所示图六：三次拟合函数拟合效果类似的四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-354.5952则拟合函数h=-0.0013x2+1.1219x3-354.5952x4,拟合效果如图七所示图七：四次拟合函数拟合效果结论以及分析检验：<一>用SAS软件做线性回归得到排沙量与时间的函数关系式为:再利用所求函数在区间0,24上进行积分得到总排沙量1.93962亿吨，这与现实情况基本相符。<二>对于第一阶段三次多项式拟合由MATLAB拟合函数求解a0=a1=0,a2=0.0032,a3=-2.4929则拟合函数h=0.0032x2-2.4929x3对于第一阶段四次多项式拟合由MATLAB拟合函数求解出a0=a1=a2=0,a3=0.0121,a4= -7.4347则拟合函数h=0.0121x3-7.4347x4对于第二阶段三次多项式拟合由MATLAB拟合函数求解出a0=a1=0,a2=-0.9475,a3= 464.9601则拟合函数h=-0.9475x2+464.9601x3对于第二阶段四次多项式拟合由MATLAB拟合函数求解出a0=a1=0，a2=-0.0013，a3= 1.1219 a4=-354.5952则拟合函数h=-0.0013x2+1.1219x3-354.5952x4讨论与推广：1、对于第一个问题排沙量与时间不是严格的正态函数关系可能与实际有些偏差，此外还可以用SAS软件进行高次的多向式回归2、对于第二个问题，由于MATLAB软件的计算可能有些偏差导致拟合的函数关系可能与实际有稍微偏差，此外，还可以进行高次的拟合。附录：1、排沙量与时间的关系图像的MATLAB程序：>>t=1:1:24;>>y=57600,2400,91000,54000,45500,30000,8000,4500;>> plot(t,y,'r')2、对排沙量求自然对数的MATLAB程序与结果：y3=log(y)y3 = Columns 1 through 17 10.9613 11.6440 11.9672 12.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、第一阶段的排沙量与水流量之间的关系MATLAB程序：>>x=1800,1900,2100,2200,2300,2400,2500,2600,2650,2700,2720;>> h=32,60,75,85,90,98,100,102,108,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、第一阶段三次多项式拟合函数以及拟合效果程序与结果：>> 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、第一阶段四次多项式拟合函数以及拟合效果程序与结果：>>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、第二阶段三次多项式拟合函数以及拟合效果程序与结果：>>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、第二阶段四次多项式拟合函数以及拟合效果程序与结果：>>A4=polyfit(x1,h1,4)>> In polyfit at 80A4 = -0.0000 0.0000 -0.0013 1.1219 -354.5952>> z4=polyval(A4,x1);>> plot(x1,h1,'k*',x1,z4,'r:')参考文献：【1】姜启源，数学模型（第三版），高等教育出版社【2】楼顺天,MATLAB程序设计语言（第二版），西安电子科技大学出版社【3】刘娜，在SAS中拟合ARCH/GARCH模型专心-专注-专业
