水塔流量问题

1、本科生课程设计报告实习课程数值分析学院名称管理科学学院专业名称学生姓名学生学号指导教师实验地点实验成绩 二 一 六 年 六 月 二 一 六 年 六 月估计水塔的水流量摘要水塔流量的估计是一个较为经典的数学建模问题，本问题最大的困难在于不知泵启动时水位的变化和向外水流的速度.解决该问题，先确定近似流速，利用中点数值求导公式计算出每个时间点出的流速，再利用插值与拟合计算出流速与时间的函数，对0到24小时积分可得总用水量，这是第一种方法.第二种方法，水泵没有开动时利用高度差计算用水量，水泵开动时利用积分，这样计算出的结果较为准确，2种方法比较，可得出误差.关键词：中点数值求导；插值与拟合；积分目录第

2、1章 前言11.1 内容及要求11.2 研究思路及结构安排2第2章 模型建立与求解32.1模型假设32.2确定近似流速32.3 确定水泵启动时的流量及总流量曲线42.4确定总用水量4第3章 算法步骤63.1 中点数值求导函数步骤及流程图63.2 三次样条插值函数步骤及流程图7第4章 算法实现74.1 程序总体结构74.2 源程序清单84.3 程序运行12第5章 误差分析14第6章 模型的评价和改进156.1 优点156.2 缺点156.3 模型的改进方向15参考文献15第1章 前言1.1 内容及要求某地的用水管理机构要求各社区提供各个时刻的用水率以及每天所用的总用水量。但许多社区并没有测量流入

3、或流出当地水塔的水量的设备，他们只能代之以每小时测量水塔中的水位，其误差不超过0.5%。更为重要的是，无论什么时候，只要水塔中的水位下降到最低水位L时，水泵就自动启动向水塔重新充水直到最高水位H时水泵自动停止，但也无法得到水泵的供水量的测量数据。因此，在水泵正在工作时，人们不容易建立水塔中水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两次，每次约二小时。下表为某地一天中的真实的数据。表1某天水塔水位测量记录时刻t（秒）03316663510619139371792121240水位（0.01英尺）3175311030542994294728922850时刻t（秒）252232854

4、33228435932393323943543318水位（0.01英尺）279527522697水泵启动水泵启动35503445时刻t（秒）46636499535393657254605746455468535水位（0.01英尺）3350326031673087301229272842时刻t（秒）71854750217925482649859688995393270水位（0.01英尺）27672697水泵启动水泵启动347533973340水塔是一个高40英尺、直径57英尺的圆柱。按照设计，水塔水位降至约L=27英尺时，水泵自动启动加水；当水位升高到约H=35.5英尺米时，水泵自动停止工作。试

5、估计在任何时刻（包括水泵正在供水时）水从水塔流出的流量f(t)，并估计一天的总用水量。1.2 研究思路及结构安排本问题最大的困难在于不知道水泵启动时水位的变化和向外水流的速度。用水量等于向外水流速度乘以时间。因此，如何确定流速是解决该问题的关键。第2章 模型建立与求解2.1模型假设（1）除了问题中特别说明的数据以外，其他给定的数据其测量误差不超过0.5%；（2）一天之中，任意从一个开始时刻，如从t=0起到t=24小时结束，一天开始时刻的不同不影响一天总水量；（3）管理部门不考虑水流速度的瞬间值，感兴趣的是整个一天中的用水总趋势；（4）水泵抽水的速度是均匀的；（5）假设水流的速度是连续变化的，流

6、速可以用一条光滑的曲线近似表示，即在给定的数据点上有连续的二阶导数。2.2确定近似流速水泵不供水时,水流速度不难用数值方法描述，但如何描述清楚水泵向水塔供水期间的流速?我们首先估计水泵不启动时,每个已知数据点上的流速。流速等于水量的变化除以时间的变化。假设为时刻水塔中水的容量，表示时刻水塔中水的容量为。使用中点数值求导公式，求处的导数的近似值。水泵不启动时段共有3段，对于每段的前两个点对于每段的后两个点 对于每段的其余点如此，可得各时刻点上的流速。 (负号保证流速为正数)2.3 确定水泵启动时的流量及总流量曲线在2.3.1中已经得到了已知数据点(水泵不启动时)上的近似流速，根据假设（5），水流

7、的速度可以表示为一条光滑的曲线。因此考虑使用已经求出的水泵未启动时的流速来近似模拟水泵启动时的流速，这种模拟有两种方式:插值和拟合。由于假设流速具有连续的二阶导数，因此采用三次样条插值多项式来模拟任意时刻的流速。用表示时刻 的流速。2.4确定总用水量前面已经求出水流的流速曲线,如何求总用水量呢? 如果画出流速曲线，就可以从流速曲线图看出,流速在0时刻和24时的流速近似相等,因此流速是近似地以一天为周期的。可以考虑以下两种方式：1）一种方式是使用流速在时刻0至24小时上积分 2)直接根据高度的变化求出水泵不启动时的用水量用积分求出水泵启动时的用水量近似值,然后相加。选择024小时内的范围，08.

8、9678小时内,高度从31.75下降至26.97英尺，用水量为在10.954220.8392小时内,高度从35.50降至2697英尺，用水量为 然后用数值积分求出两次水泵启动期间的用水量 总用水量为第3章 算法步骤3.1 中点数值求导函数步骤及流程图图 1 中点数值求导函数流程图3.2 三次样条插值函数步骤及流程图图2 三次样条插值函数流程图第4章 算法实现4.1 程序总体结构图 3 程序函数调用关系4.2 源程序清单#include#include#define N 24double Q1=0,Q2=0,q1=0,q2=0,nN=0,mN=0;/依次为法一积分、法二积分、两次启动水泵期间用水

9、量、每次积分void zhuigan2(double a,double b,double c,double d,double m) /第二类边界追赶法求解方程组的解double eN,fN;e0=c1/b1;f0=d1/b1;for(int i=1;iN-3;i+)ei=ci+1/(bi+1-ai+1*ei-1);for(i=1;i0;i-)mi=fi-1-ei-1*mi+1;void chazhi2(double x,double y,double s) /第二类边界条件及自然插值double hN,aN,bN,cN,dN,mN; for(int i=0;iN;i+)bi=2;for(i=0

10、;iN-1;i+) /计算hNhi=xi+1-xi;for(i=1;iN-1;i+) /计算aN,cNai=hi-1/(hi+hi-1);ci=1-ai;for(i=1;iN-1;i+) /计算dNdi=(6/(hi+hi-1)*(yi+1-yi)/hi)-(yi-yi-1)/hi-1);d1=d1-a1*s0;dN-2=dN-2-cN-2*s1;m0=s0;mN-1=s1;zhuigan2(a,b,c,d,m); /追赶法求mNprintf(n追赶法求出m:n);for(i=0;iN;i+) /输出mNprintf(m%d=%.2fn,i,mi);/printfcz(h,x,y,m);voi

11、d median(double *t,double *f,double *v)v0=1/(t2-t0)*(-3*f0+4*f1-f2);/第一段数据中值求导v1=1/(t2-t0)*(f2-f0);v8=1/(t9-t7)*(f9-f7);v9=1/(t9-t7)*(3*f9-4*f8+f7);for(int i=2;i8;i+)vi=1/(3*(ti+2-ti-2)*(fi-2-8*fi-1+8*fi+1-fi+2);v10=1/(t12-t10)*(-3*f10+4*f11-f12);/第二段数据中值求导v11=1/(t12-t10)*(f12-f10);v19=1/(t20-t18)*(

12、f20-f18);v20=1/(t20-t18)*(3*f20-4*f19+f18);for(i=12;i19;i+)vi=1/(3*(ti+2-ti-2)*(fi-2-8*fi-1+8*fi+1-fi+2);v21=1/(t23-t21)*(-3*f21+4*f22-f23);/第三段数据中值求导v22=1/(t23-t21)*(f23-f21);v23=1/(t23-t21)*(3*f23-4*f22+f21);void Integral(double h,double x,double y) /积分法int i;double k3,c,d; /k数组为三次样条插值函数各项系数printf

13、(*n);printf(计算一天总水量:n);printf(tt方法1：024小时函数直接积分计算总用水量n);printf(积分结果为:n);for(i=0;iN-1;i+) /计算每段函数的积分 k0=(mi+1-mi)/(6*hi);k1=mi/2;k2=(yi+1-yi)/hi-hi*(2*mi+mi+1)/6;if(i21)c=pow(xi+1-xi,4)*k0/4+pow(xi+1-xi,3)*k1/3+pow(xi+1-xi,2)*k2/2+yi*xi+1; /代上限d=yi*xi; /代下限ni=c-d; elsec=pow(24-xi,4)*k0/4+pow(24-xi,3)

14、*k1/3+pow(24-xi,2)*k2/2+yi*24; /代上限d=yi*xi; /代下限ni=c-d; for(i=0;iN-1;i+)Q1+=ni; /每段积分求和,即0到24小时用水量printf(%fn,Q1);printf(*n);void way2(double t,double f) /高度变化结合积分求水量double PI=3.14159265,DIA=57;q1=fabs(PI*(DIA/2.0)*(DIA/2.0)*(f9-f0)*0.01;q2=fabs(PI*(DIA/2.0)*(DIA/2.0)*(f19-f10)*0.01;Q2+=q1+q2+n9+n20+

15、n21;printf(tt方法2：积分与高度变化结合求用水量);printf(n第1次水泵启动时用水量为：%fn,n9);printf(第2次水泵启动时用水量为：%fn,n20);printf(n%f%f用水量为：%fn,t0,t9,q1);printf(%f%f用水量为：%fn,t10,t20,q2);printf(总用水量为：%fn,Q2);printf(*n);void CountError() /计算误差double w;w=(Q1-n9-n20-n21)-(Q2-n9-n20-n21)/(Q2-n9-n20-n21);printf(用水量误差为：%f%n,w*100);void ma

16、in()double tN=0,0.9211,1.8431,2.9497,3.8714,4.9781,5.9000,7.0064,7.9286,8.9678,10.9256,12.0328,12.9544,13.8758,14.9822,15.9039,16.8261,17.9317,19.0375,19.9594,22.9581,23.8800,24.9869,25.9083;double fN=3175,3110,3054,2994,2947,2892,2850,2795,2752,2697,3550,3445,3350,3260,3167,3087,3012,2927,2842,2767

17、,2697,3475,3397,3340;double vN=0;double S=3.1415*28.5*28.5;for(int i=0;iN;i+)fi=fi*S*0.01;printf(数据已录入，现在计算n);median(t,f,v);/计算导数for(i=0;iN;i+)vi=-vi;printf(中值求导计算完毕，各时间点流量导数为:n);for(i=0;iN;i+)printf(%10.2f,vi);printf(现在开始三次样条插值近似出导数函数n);double s2=0;chazhi2(t,v,s);/计算拟合函数for(i=0;iN;i+)fi=fi/S/0.01;I

18、ntegral(f,t,v);way2(t,f); /方法2计算总水量CountError(); /计算误差4.3 程序运行图4 水塔流量程序运行结果图5 水塔流量程序运行结果运行结果说明：通过积分的方法算出的用水量为37879，通过综合计算的方法计算出的用水量为36363，两次结果都与实验说明相差8000，原因应该是计算过程中，曲线的拟合处出现了较大误差；而两次计算的相对误差为4.7%，结果相对比较接近实验说明的要求。所以，此次编程虽然有些许误差，但是实验结果仍然表明整个方法没有根本性错误。 第5章 误差分析考虑前面两种求总用水量的结果，由数值积分求出，结合高度和数值积分求出。根据高度求出的水泵未启动时的用水量可以认为是精确的。由数值积分求出的水泵未启动时的用水量。根据高度求出的水泵未启动时的用水量，可以认为是精确值。我们使用水泵未启动时的用水量的误差近似替代总的误差。第6章 模型的评价和改进6.1 优点1）模型的灵活性较好，可以用到那些拥有垂直圆柱型水塔的社区或单位，模