matlab实现水塔水流量问题

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

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

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

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

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

 表2 各时刻水体积表

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

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

 

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

差分公式为：

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

 

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

 

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

 

 表3 各时刻水流量

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

方法一：直接积分法

 方法二：分段计算

 

clear,clc
c=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:n1
    if 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)));
    end
end
%计算第二段
n2=21;
for i=n1+1:n2
    if i<=n1+2
        f(i)=abs(-3*v(i)+4*v(i+1)-v(i+2))/(2*(t(i+1)-t(i)));
    elseif i<=n2-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>=n2-1
        f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));
    end
end
%计算第三段
n3=25;
for i=n2+1:n3
    if 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-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>=n3-1
        f(i)=abs(3*v(i)-4*v(i-1)+v(i-2))/(2*(t(i)-t(i-1)));
    end
end
plot(t,f,'r*');%原始图
tmin=min(t);tmax=max(t);
tt=tmin:0.1:tmax;%获取离散的时间点，以便做样条曲线
ff=spline(t,f,tt);%计算三次样条插值
hold on
plot(tt,ff,'b');%画样条曲线
xlabel('时间(小时)');
ylabel('水塔流量图');
hold off
dt=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);