水塔流量问题资料.doc

水塔流量问题

精品资料

本科生课程设计报告

实习课程
数值分析
学院名称
管理科学学院
专业名称

学生姓名

学生学号

指导教师

实验地点

实验成绩

二〇 一 六 年 六 月 二〇 一 六 年 六 月　　　　　

估计水塔的水流量
摘要
水塔流量的估计是一个较为经典的数学建模问题，本问题最大的困难在于不知泵启动时水位的变化和向外水流的速度.解决该问题，先确定近似流速，利用中点数值求导公式计算出每个时间点出的流速，再利用插值与拟合计算出流速与时间的函数，对0到24小时积分可得总用水量，这是第一种方法.第二种方法，水泵没有开动时利用高度差计算用水量，水泵开动时利用积分，这样计算出的结果较为准确，2种方法比较，可得出误差.
关键词：中点数值求导；插值与拟合；积分

目录
第1章 前言 1
1.1 内容及要求 1
1.2 研究思路及结构安排 2
第2章 模型建立与求解 3
2.1模型假设 3
2.2确定近似流速 3
2.3 确定水泵启动时的流量及总流量曲线 4
2.4确定总用水量 4
第3章 算法步骤 6
3.1 中点数值求导函数步骤及流程图 6
3.2 三次样条插值函数步骤及流程图 7
第4章 算法实现 7
4.1 程序总体结构 7
4.2 源程序清单 8
4.3 程序运行 12
第5章 误差分析 14
第6章 模型的评价和改进 15
6.1 优点 15
6.2 缺点 15
6.3 模型的改进方向 15
参考文献 15

仅供学习与交流，如有侵权请联系网站删除 谢谢15

第1章 前言
1.1 内容及要求
某地的用水管理机构要求各社区提供各个时刻的用水率以及每天所用的总用水量。但许多社区并没有测量流入或流出当地水塔的水量的设备，他们只能代之以每小时测量水塔中的水位，其误差不超过0.5%。更为重要的是，无论什么时候，只要水塔中的水位下降到最低水位L时，水泵就自动启动向水塔重新充水直到最高水位H时水泵自动停止，但也无法得到水泵的供水量的测量数据。因此，在水泵正在工作时，人们不容易建立水塔中水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两次，每次约二小时。下表为某地一天中的真实的数据。
表1某天水塔水位测量记录
时刻t（秒）
0
3316
6635
10619
13937
17921
21240
水位（0.01英尺）
3175
3110
3054
2994
2947
2892
2850
时刻t（秒）
25223
28543
32284
35932
39332
39435
43318
水位（0.01英尺）
2795
2752
2697
水泵启动
水泵启动
3550
3445
时刻t（秒）
46636
49953
53936
57254
60574
64554
68535
水位（0.01英尺）
3350
3260
3167
3087
3012
2927
2842
时刻t（秒）
71854
75021
79254
82649
85968
89953
93270
水位（0.01英尺）
2767
2697
水泵启动
水泵启动
3475
3397
3340
水塔是一个高40英尺、直径57英尺的圆柱。按照设计，水塔水位降至约L=27英尺时，水泵自动启动加水；当水位升高到约H=35.5英尺米时，水泵自动停止工作。
试估计在任何时刻（包括水泵正在供水时）水从水塔流出的流量f(t)，并估计一天的总用水量。
1.2 研究思路及结构安排
本问题最大的困难在于不知道水泵启动时水位的变化和向外水流的速度。用水量等于向外水流速度乘以时间。因此，如何确定流速是解决该问题的关键。

第2章 模型建立与求解
2.1模型假设
（1）除了问题中特别说明的数据以外，其他给定的数据其测量误差不超过0.5%；
（2）一天之中，任意从一个开始时刻，如从t=0起到t=24小时结束，一天开始时刻的不同不影响一天总水量；
（3）管理部门不考虑水流速度的瞬间值，感兴趣的是整个一天中的用水总趋势；
（4）水泵抽水的速度是均匀的；
（5）假设水流的速度是连续变化的，流速可以用一条光滑的曲线近似表示，即在给定的数据点上有连续的二阶导数。
2.2确定近似流速
水泵不供水时,水流速度不难用数值方法描述，但如何描述清楚水泵向水塔供水期间的流速?
我们首先估计水泵不启动时,每个已知数据点上的流速。流速等于水量的变化除以时间的变化。假设为时刻水塔中水的容量，表示时刻水塔中水的容量为。
使用中点数值求导公式，求处的导数的近似值。水泵不启动时段共有3段，对于每段的前两个点

对于每段的后两个点

对于每段的其余点

如此，可得各时刻点上的流速。 (负号保证流速为正数)
2.3 确定水泵启动时的流量及总流量曲线
在2.3.1中已经得到了已知数据点(水泵不启动时)上的近似流速，根据假设（5），水流的速度可以表示为一条光滑的曲线。因此考虑使用已经求出的水泵未启动时的流速来近似模拟水泵启动时的流速，这种模拟有两种方式:插值和拟合。
由于假设流速具有连续的二阶导数，因此采用三次样条插值多项式来模拟任意时刻的流速。用表示时刻 的流速。
2.4确定总用水量
前面已经求出水流的流速曲线,如何求总用水量呢? 如果画出流速曲线，就可以从流速曲线图看出,流速在0时刻和24时的流速近似相等,因此流速是近似地以一天为周期的。
可以考虑以下两种方式：
1）一种方式是使用流速在时刻0至24小时上积分

2)直接根据高度的变化求出水泵不启动时的用水量
用积分求出水泵启动时的用水量近似值,然后相加。选择024小时内的范围，08.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<stdio.h>
#include<math.h>
#define N 24
double Q1=0,Q2=0,q1=0,q2=0,n[N]={0},m[N]={0};//依次为法一积分、法二积分、两次启动水泵期间用水量、每次积分
void zhuigan2(double a[],double b[],double c[],double d[],double m[]) //第二类边界追赶法求解方程组的解
{
double e[N],f[N];
e[0]=c[1]/b[1];
f[0]=d[1]/b[1];
for(int i=1;i<N-3;i++)
e[i]=c[i+1]/(b[i+1]-a[i+1]*e[i-1]);
for(i=1;i<N-2;i++)
f[i]=(d[i+1]-a[i+1]*f[i-1])/(b[i+1]-a[i+1]*e[i-1]);
m[N-2]=f[N-3];
for(i=N-3;i>0;i--)
m[i]=f[i-1]-e[i-1]*m[i+1];
}

void chazhi2(double x[],double y[],double s[]) //第二类边界条件及自然插值
{
double h[N],a[N],b[N],c[N],d[N],m[N];
for(int i=0;i<N;i++)
b[i]=2;
for(i=0;i<N-1;i++) //计算h[N]
h[i]=x[i+1]-x[i];
for(i=1;i<N-1;i++) //计算a[N],c[N]
{
a[i]=h[i-1]/(h[i]+h[i-1]);
c[i]=1-a[i];
}
for(i=1;i<N-1;i++) //计算d[N]
d[i]=(6/(h[i]+h[i-1]))*(((y[i+1]-y[i])/h[i])-((y[i]-y[i-1])/h[i-1]));
d[1]=d[1]-a[1]*s[0];
d[N-2]=d[N-2]-c[N-2]*s[1];
m[0]=s[0];
m[N-1]=s[1];
zhuigan2(a,b,c,d,m); //追赶法求m[N]
printf("\n追赶法求出m:\n");
for(i=0;i<N;i++) //输出m[N]
printf("m[%d]=%.2f\n",i,m[i]);
//printfcz(h,x,y,m);

}

void median(double *t,double *f,double *v)
{
v[0]=1/(t[2]-t[0])*(-3*f[0]+4*f[1]-f[2]);//第一段数据中值求导
v[1]=1/(t[2]-t[0])*(f[2]-f[0]);
v[8]=1/(t[9]-t[7])*(f[9]-f[7]);
v[9]=1/(t[9]-t[7])*(3*f[9]-4*f[8]+f[7]);

for(int i=2;i<8;i++)
v[i]=1/(3*(t[i+2]-t[i-2]))*(f[i-2]-8*f[i-1]+8*f[i+1]-f[i+2]);

v[10]=1/(t[12]-t[10])*(-3*f[10]+4*f[11]-f[12]);//第二段数据中值求导
v[11]=1/(t[12]-t[10])*(f[12]-f[10]);
v[19]=1/(t[20]-t[18])*(f[20]-f[18]);
v[20]=1/(t[20]-t[18])*(3*f[20]-4*f[19]+f[18]);

for(i=12;i<19;i++)
v[i]=1/(3*(t[i+2]-t[i-2]))*(f[i-2]-8*f[i-1]+8*f[i+1]-f[i+2]);

v[21]=1/(t[23]-t[21])*(-3*f[21]+4*f[22]-f[23]);//第三段数据中值求导
v[22]=1/(t[23]-t[21])*(f[23]-f[21]);
v[23]=1/(t[23]-t[21])*(3*f[23]-4*f[22]+f[21]);
}

void Integral(double h[],double x[],double y[]) //积分法
{
int i;
double k[3],c,d; //k数组为三次样条插值函数各项系数
printf("*******************************************************************************\n");
printf("计算一天总水量:\n");
printf("\t\t方法1：0~24小时函数直接积分计算总用水量\n");
printf("积分结果为:\n");
for(i=0;i<N-1;i++) //计算每段函数的积分
{
k[0]=(m[i+1]-m[i])/(6*h[i]);
k[1]=m[i]/2;
k[2]=(y[i+1]-y[i])/h[i]-h[i]*(2*m[i]+m[i+1])/6;
if(i<21)
{
c=pow(x[i+1]-x[i],4)*k[0]/4+pow(x[i+1]-x[i],3)*k[1]/3+pow(x[i+1]-x[i],2)*k[2]/2+y[i]*x[i+1]; //代上限
d=y[i]*x[i]; //代下限
n[i]=c-d;
}
else
{
c=pow(24-x[i],4)*k[0]/4+pow(24-x[i],3)*k[1]/3+pow(24-x[i],2)*k[2]/2+y[i]*24; //代上限
d=y[i]*x[i]; //代下限
n[i]=c-d;
}
}
for(i=0;i<N-1;i++)
{
Q1+=n[i]; //每段积分求和,即0到24小时用水量
}
printf("%f\n",Q1);
printf("*******************************************************************************\n");
}

void way2(double t[],double f[]) //高度变化结合积分求水量
{
double PI=3.14159265,DIA=57;
q1=fabs(PI*(DIA/2.0)*(DIA/2.0)*(f[9]-f[0]))*0.01;
q2=fabs(PI*(DIA/2.0)*(DIA/2.0)*(f[19]-f[10]))*0.01;
Q2+=q1+q2+n[9]+n[20]+n[21];
printf("\t\t方法2：积分与高度变化结合求用水量");
printf("\n第1次水泵启动时用水量为：%f\n",n[9]);
printf("第2次水泵启动时用水量为：%f\n",n[20]);
printf("\n%f~%f用水量为：%f\n",t[0],t[9],q1);
printf("%f~%f用水量为：%f\n",t[10],t[20],q2);
printf("总用水量为：%f\n",Q2);
printf("*******************************************************************************\n");
}
void CountError() //计算误差
{
double w;
w=((Q1-n[9]-n[20]-n[21])-(Q2-n[9]-n[20]-n[21]))/(Q2-n[9]-n[20]-n[21]);
printf("用水量误差为：%f%%\n",w*100);
}

void main()
{
double t[N]={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 f[N]={3175,3110,3054,2994,2947,2892,2850,2795,2752,2697,3550,3445,3350,3260,3167,3087,3012,2927,2842,2767,2697,3475,3397,3340};
double v[N]={0};
double S=3.1415*28.5*28.5;
for(int i=0;i<N;i++)
f[i]=f[i]*S*0.01;

printf("数据已录入，现在计算\n");

median(t,f,v);//计算导数
for(i=0;i<N;i++)
v[i]=-v[i];
printf("中值求导计算完毕，各时间点流量导数为:\n");
for(i=0;i<N;i++)
printf("%10.2f",v[i]);

printf("现在开始三次样条插值近似出导数函数\n");
double s[2]={0};

chazhi2(t,v,s);//计算拟合函数

for(i=0;i<N;i++)
f[i]=f[i]/S/0.01;

Integral(f,t,v);
way2(t,f); //方法2计算总水量
CountError(); //计算误差
}

4.3 程序运行

图4 水塔流量程序运行结果

图5 水塔流量程序运行结果

运行结果说明：通过积分的方法算出的用水量为37879，通过综合计算的方法计算出的用水量为36363，两次结果都与实验说明相差8000，原因应该是计算过程中，曲线的拟合处出现了较大误差；
而两次计算的相对误差为4.7%，结果相对比较接近实验说明的要求。所以，此次编程虽然有些许误差，但是实验结果仍然表明整个方法没有根本性错误。

第5章 误差分析
考虑前面两种求总用水量的结果，由数值积分求出，结合高度和数值积分求出。根据高度求出的水泵未启动时的用水量可以认为是精确的。由数值积分求出的水泵未启动时的用水量。根据高度求出的水泵未启动时的用水量，可以认为是精确值。我们使用水泵未启动时的用水量的误差近似替代总的误差。

第6章 模型的评价和改进
6.1 优点
1）模型的灵活性较好，可以用到那些拥有垂直圆柱型水塔的社区或单位，模型中的数据只要是离散的时间和水位的高度；
2）)模型所用的数学概念简单，并且也容易理解，一阶导数、三次样条插值和积分都是被广泛使用的概念；
3）只要提供了数据和计算机，模型很容易实现；
4）模型为一天里的流速和总用水量提供了较为精确的估计。
6.2 缺点
1)模型受到水塔形状限制；
2)光滑曲线不能模拟真实流速观测中的微小变化。
6.3 模型的改进方向
1) 可以增加一些参数,以使模型适用于不同形状的水塔；
2) 也可以加入随机因素，并进行方差分析，使模型更适应于真实的情况。

参考文献

学生学习心得
通过对三次样条插值法的应用，更加深入地了解了其在实际生产生活中的应用；通过对中值求导法的应用，能够更加熟练地对实际数据进行处理，以得到我们想要的数据。
在对水塔问题的处理中，充分体会到了一些简单的数值处理方法的有机组合能够达到的强大功能，能够处理的问题的复杂性。

学生（签名）：
2016 年 6月 16日
诚信承诺
本人郑重声明所呈交的课程报告是本人在指导教师指导下进行的研究工作及取得的研究成果。据我所知，除了文中特别加以标注的地方外，论文中不包含其他人已经发表或撰写过的研究成果。与我一同工作的同学对本文研究所做的贡献均已在报告中作了明确的说明并表示谢意。
学生（签名）：
任课
教师
评语

成绩评定：
任课教师（签名）：
年 月 日