matlab实现水塔水流量问题

代码以及过程学习来源于：西北工业大学 肖华用老师课程 代码在文章最后

附上链接：9.5 插值与拟合模型(二)-----水塔流量问题_哔哩哔哩_bilibili

也可以去中国mooc观看肖老师视频

问题描述：美国某洲的各用水管理机构要求各社区提供以每小时多少加仑计的用水率以及每天总的用水量,但许多社区并没有测量水流入或流出水塔水量的设备,他们只能每小时测量水塔中的水位，精度在0.5%以内,更为重要的是,无论什么时候,只要水塔中的水位下降到某一最低水位L时，水泵就启动向水塔重新充水至某一最高水位H，但也无法得到水泵的供水量的测量数据。水泵每天向水塔充水一次或两次，每次约两小时。水塔是一个垂直圆形柱体，高为40英尺，直径57英尺。，当水位升至35.50英尺时停止充水。试估计在任何时刻，甚至包括水泵正在工作期间内，水从水塔流出的流量f(t)，并估计一天的总用水量。

表1小镇某天的水塔水位(0.01英尺)

表2各时刻水体积表

要从水的流速得到水的流量，可以使用直接积分法，对0到24小时的水的流速进行积分计算即可得到一天的总用水量。或者使用分部积分法，对第一次充水前，第一次充水到第二次期间，第二次充水后，这三段进行分部积分，然后将两种积分方法进行比较。

计算各时刻点的水流量(加仑/小时)：

25个时刻处的水流量采用差分的方法得到，共分三段分别处理

差分公式为：

对每段前两点采用向前三点差分公式：

对每段最后两点采用向后三点差分公式：

对每段中间点采用中心差分公式：

表3各时刻水流量

使用数值积分的方式，求解一天的总用水量

方法一：直接积分法

方法二：分段计算

clear,clcc=0.3048;%英尺和米的换算p=1.0/3.785;%升和加仑的换算d=57*c;%直径h=31.75*c;%v=pi*d*d*h/4*1000*p;data=[0,3175;3316,3110;6635,3054;10619,2994;13937,2947;17921,2892;21240,2850;25223,2797;28543,2752;32284,2697;39435,3550;43318,3445;46636,3350;49953,3260;53936,3167;57254,3087;60574,3012;64554,2927;68535,2842;71854,2767;75021,2697;82649,3550;85968,3475;89953,3397;93270,3340];%数据导入t=data(:,1)/3600;%计算时间v=pi*d*d*data(:,2)/100*c/4*1000*p;%计算体积%计算差分n=length(v);f=zeros(n,1);%储存差分值%计算第一段n1=10;for i=1:n1if i<=2 %前两点用向前三点差分f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));elseif i<=n1-2%中心差分公式f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));elseif i>=n1-1 %后两点用向后三点差分f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));endend%计算第二段n2=21;for i=n1+1:n2if i<=n1+2f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));elseif i<=n2-2f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));elseif i>=n2-1f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));endend%计算第三段n3=25;for i=n2+1:n3if i<=n2+2 f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));elseif i<=n3-2f(i)=abs(-v(i+2)+8*v(i+1)-8*v(i-1)+v(i-2))/(12*(t(i+1)-t(i)));elseif i>=n3-1f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));endendplot(t,f,'r*');%原始图tmin=min(t);tmax=max(t);tt=tmin:0.1:tmax;%获取离散的时间点，以便做样条曲线ff=spline(t,f,tt);%计算三次样条插值hold onplot(tt,ff,'b');%画样条曲线xlabel('时间(小时)');ylabel('水塔流量图');hold offdt=0.05;t2=0.5:dt:24.5;%获取离散时间点，以便积分nn=length(t2);f2=spline(t,f,t2);%计算24小时用水量，使用复化梯形公式s=(f2(1)+f2(nn)+2*sum(f2(2:nn-1)))*dt/2;fprintf('(全部积分法)1天总水量s=%8.2f\n',s);%第一次水塔增加的水 v10=v(11)-v(10);dt1=t(11)-t(10);%第一次充本其间流出的水tp=t(10):dt:t(11);nn=length(tp); yp=spline(t,f,tp); %计算三次样条插值v11=(yp(1)+yp(nn)+2*sum(yp(2:nn-1)))*dt/2;v1=v10+v11;%第一次充水总量 p1=v1/dt1; %第一次充水的平均水流量%第一次充水前总流量vv1=v(1)-v(10);%两次充水间总流量vv2=v(11)-v(21);% t=83649到85968其间流量vv3=v(22)-v(23);%第一次充水期间流量ta=t(10):dt:t(11);%获得离散的时间点，用于积分nn=length(ta);fa=spline(t,f,ta);s1=(fa(1)+fa(nn)+2*sum(fa(2:nn-1)))*dt/2;%第二次充水期间流量tb=t(21):dt:t(22); %获得离散的时间点用于积分nn=length(tb);fb=spline(t,f,tb);s2=(f(1)+fb(nn)+2*sum(fb(2:nn-1)))*dt/2;%t=85968到86400流量tc=t(23):dt:24; %获得离散的时间点，用于积分nn=length(tc);fc=spline(t,f,tc);s3=(fc(1)+fc(nn)+2*sum(fc(2:nn-1)))*dt/2;ss=vv1 +vv2+vv3+s1+s2+s3;fprintf('部分积分法)1天总水流量ss= %8.2f\n',ss);err=abs((s-ss)/s);fprintf('两种计算法总量相对误差%6.2f%%\n',err* 100);