实验一-水塔流量问题

实验一水塔流量问题1.1问题的提出与要求某地的用水管理机构要求各社区提供各个时刻的用水率以及每天所用的总用水量。但许多社区并没有测量流入或流出当地水塔的水量的设备，他们只能代之以每小时测量水塔中的水位，其误差不超过0.5%。更为重要的是，无论什么时候，只要水塔中的水位下降到最低水位L时，水泵就自动启动向水塔重新充水直到最高水位H时水泵自动停止，但也无法得到水泵的供水量的测量数据。因此，在水泵正在工作时，人们不容易建立水塔中水位与水泵工作时的用水量之间的关系。水泵每天向水塔充水一次或两次，每次约二小时。表1-1为某地一天中的真实数据。表1-1某天水塔水位测量记录时刻t（秒）03316663510619139371792121240水位（0.01英尺）3175311030542994294728922850时刻t（秒）25223285433228435932393323943543318水位（0.01英尺）279527522697水泵启动水泵启动35503445时刻t（秒）46636499535393657254605746455468535水位（0.01英尺）3350326031673087301229272842时刻t（秒）71854750217925482649859688995393270水位（0.01英尺）27672697水泵启动水泵启动347533973340水塔是一个高40英尺、直径57英尺的圆柱。按照设计，水塔水位降至约L=27英尺时，水泵自动启动加水；当水位升高到约H=35.5英尺米时，水泵自动停止工作。试估计在任何时刻（包括水泵正在供水时）水从水塔流出的流量，并估计一天的总用水量。实验要求和提示:（1）本问题要用到插值拟合、数值微分和数值积分等知识。（2）在数值微分、插值拟合、数值积分阶段都有多种方法，多去尝试不同方法，合理选择算法。（3）对一天总用水量做误差估计。（4）根据流速曲线，可以向管理机构提交一份用水量变化的报告，以便管理机构掌握向居民供水的强度，如用水高峰时的安排等。（5）注意单位要统一。1.2模型的建立1.2.1问题分析为了分析的方便，先将数据中的时间单位化为小时，如表1-2表1-2时间h（小时）t1t2t3t4t5t6t700.921111.843052.949723.871384.978055.9000t8t9t10水泵启动水泵启动t11t127.006387.928618.967779.9811110.9255510.9541612.03277t13t14t15t16t17t18t1912.9544413.8758314.9822215.9038816.8261117.9316619.03750t20t21水泵启动水泵启动t22t23t2419.9594420.8391622.0150022.9580523.8799924.9869425.90833本问题最大的困难在于不知道水泵启动时水位的变化和向外水流的速度。用水量等于向外水流速度乘以时间，即流速的积分。因此如何确定流速是解决该问题的关键。1.2.2模型的假设在给出问题解决方法之前，需要做下面假设：（1）除了问题中特别说明的数据以外，其他给定的数据其测量误差不超过0.5%；（2）一天之中，任意从一个开始时刻，如从起到小时结束，一天开始时刻的不同不影响一天总水量；（3）管理部门不考虑水流速度的瞬间值，感兴趣的是整个一天中的用水总趋势；（4）水泵抽水的速度是均匀的；（5）假设水流的速度是连续变化的，流速可以用一条光滑的曲线近似表示；（6）水塔中水流量是时间的连续光滑函数，与水泵工作无关；（7）水泵工作与否完全取决于水塔内水位的高低，且每次加水的工作时间大约为2小时；（8）水塔为标准圆柱体,水塔截面积是常数3.14*28.5*28.5;利用水塔截面积是常数,得到不同时刻水塔中水的体积如表1-3表1-3不同时刻水塔中水的体积时刻t(小时)00.921111.843052.949723.871384.978055.9000体积v（立方英尺）8.097727.931947.789127.636097.516227.375947.26882时刻t(小时)7.006387.928618.967779.9811110.9255510.9541612.03277体积V（立方英尺）7.128547.018876.87860水泵启动水泵启动9.054158.78635时刻t(小时)12.9544413.8758314.9822215.9038816.8261117.9316619.03750体积V（立方英尺）8.544058.314518.077327.873287.682007.465217.24842时刻t(小时)19.9594420.8391622.0150022.9580523.8799924.9869425.90833体积V（立方英尺）7.057136.87860水泵启动水泵启动8.862868.663928.518551.2.3模型的建立1、确定各时刻点的近似流速水流速度应该是水塔中水的体积对时间的导数。由于没有水的体积关于时间的函数表达式，而只是一个离散的函数值表表1-3，因此考虑用差商代替导数,，这也是离散反映连续的常用思想。因为所有数据被水泵工作分割成三组数据,对每组数据的中间数据采用中心差商,前后两个数据采用向前或向后差商。（采用数值微分的五点公式）算法原理步骤：对于每组的前两个点（1-1）对于每组的后两个点（1-2）对于每组的其他点（1-3）计算出不同时刻水塔中水的流速,见表1-4表1-4水塔中水的流速时刻t(小时)00.921111.843052.949723.871384.978055.9000流速v1.923511.674421.537191.336961.282471.209211.21743时刻t(小时)7.006387.928618.967779.9811110.9255510.9541612.03277流速v1.255071.274321.58641水泵启动水泵启动2.805122.55011时刻t(小时)12.9544413.8758314.9822215.9038816.8261117.9316619.03750流速v2.428942.296841.901182.002072.180122.044322.42894时刻t(小时)19.9594420.8391622.0150022.9580523.8799924.9869425.90833流速v2.052641.91107水泵启动水泵启动2.225631.697512.150182、确定水泵启动时的流速及总流速曲线问题已经转变为根据流速的一个函数值表（表1-4），产生函数在整个区间(024小时)上的函数。插值和拟合是最常用的方法……三次样条插值函数算法原理步骤：先计算出，其中：（1-4）可得(1-5)根据边界条件，建立方程组，求解未知向量。然后带入下面的式子（1-6），就得到三次样条插值函数的分段表达式：（1-6）公式（1-7）式也可以写成(1-7）其中，，，，，(1-8)3、确定总用水量（1）用构造出来的三次样条插值函数表示流速的函数，对流速进行积分，区间为0到24小时，求出一天的总用水量。（2）根据表1-3可得三组用水量的变化量，由于第三组的时间超过了24小时，我们只代入第一、二组的用水变化量（两组变化量之和表示为c），水泵启动时用水量未知，我们对两段时间的流速进行积分得出估计用水量a和b，总用水量可表示为s1=a+b+c.1.3模型的求解1.3.1函数的定义各个函数的功能如表1-5所示：函数名输入参数返回值函数功能名称类型说明function[v,t,h]=ls()--------v,t,h计算对应时间点的流速functionS=csfit1()X,Y,dx0,dxnintX为时间t，Y为流速v，dx0=0,dxn=0S根据时间t与流速v得出三次样条插值函数v(t)function[f]=hs()S，t,t1doubleS为v(t)的系数矩阵，t为总的时间数组，t1为时刻f计算某时刻的流速functionf1=z()S,,t,t1,t2doubleS为v(t)的系数矩阵，t为总的时间数组，t1，t2为起始，终止时刻f1计算某个时间段的水量表1-5函数表流程图：1.函数function[v,t,h]=ls()，流程图如图1-1所示图1-1函数ls流程图2.函数functionS=csfit1(X,Y,dx0,dxn)，流程图如图1-2所示图1-2函数csfit1流程图3.函数function[f]=hs(S,t,t1)，流程图如图1-3所示图1-3函数hs流程图4.函数functionf1=z(S,t,t1,t2)，流程图如图1-4所示图1-4函数z流程图1.3.2程序总体结构函数间的调用关系如图1-5所示：图1-5函数间调用关系图主函数流程图如图1-6所示：图1-6主函数流程图1.3.3源程序function[v,t]=ls()%计算对应时间点的流速r=28.5;t11=[03316663510619139371792121240252232854332284];h1=[3175311030542994294728922850279527522697];t1=t11/3600;a=3.14*r*r*h1/100;%对应时刻水塔中的水量n1=10;v1(1)=abs(-3*a(1)+4*a(2)-a(3))/(t1(3)-t1(1));v1(2)=abs(a(3)-a(1))*1/(t1(3)-t1(1));fori=3:n1-2v1(i)=abs(a(i-2)-8*a(i-1)+8*a(i+1)-a(i+2))*1/(3*(t1(i+2)-t1(i-2)));endv1(n1-1)=abs(a(n1)-a(n1-2))*1/(t1(n1)-t1(n1-2));v1(n1)=abs(3*a(n1)-4*a(n1-1)+a(n1-2))*1/(t1(n1)-t1(n1-2));%第一组中对应时间的流速t22=[3943543318466364995353936572546057464554685357185475021];h2=[35503445335032603167308730122927284227672697];t2=t22/3600;b=3.14*r*r*h2/100;n2=11;v2(1)=abs(-3*b(1)+4*b(2)-b(3))*1/(t2(3)-t2(1));v2(2)=abs(b(3)-b(1))*1/(t2(3)-t2(1));fori=3:n2-2v2(i)=abs(b(i-2)-8*b(i-1)+8*b(i+1)-b(i+2))*1/(3*(t2(i+2)-t2(i-2)));endv2(n2-1)=abs(b(n2)-b(n2-2))*1/(t2(n2)-t2(n2-2));v2(n2)=abs(3*b(n2)-4*b(n2-1)+b(n2-2))*1/(t2(n2)-t2(n2-2));%第二组中对应时间的流速t33=[859688995393270];h3=[347533973340];t3=t33/3600;c=3.14*r*r*h3/100;v3(1)=abs(-3*c(1)+4*c(2)-c(3))*1/(t3(3)-t3(1));v3(2)=abs(c(3)-c(1))*1/(t3(3)-t3(1));v3(3)=abs(3*c(3)-4*c(2)+c(2))*1/(t3(3)-t3(1));%第三组中对应时间的流速v=[v1,v2,v3];t=[t1,t2,t3];endfunctionS=csfit1(X,Y,dx0,dxn)%Input-Xisthe1xnabscissavector%-Yisthe1xnordinatevector%-dxo=S''(x0)firstderivativeboundarycondition%-dxn=S''(xn)firstderivativeboundarycondition%Output-S:rowsofSarethecoefficientsforthecubicinterpolantsN=length(X);H=diff(X);D=diff(Y)./H;fork=1:N-2u(k)=H(k)/(H(k)+H(k+1));endfork=1:N-2L(k)=1-u(k);endfork=2:N-1g(k)=6*(D(k)-D(k-1))/(H(k-1)+H(k));endg(1)=6*(D(1)-dx0)/H(1);g(N)=6*(dxn-D(N-1))/H(N-1);g(2)=g(2)-u(1)*dx0;g(N-1)=g(N-1)-L(N-2)*dxn;fork=1:N-2R(k,1)=g(k+1);endA=2*eye(N-2);fork=2:N-2A(k,k-1)=u(k);A(k-1,k)=L(k-1);endm=inv(A)*R;M(1)=dx0;M(N)=dxn;fork=2:N-1M(k)=m(k-1);endfork=0:N-2S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));S(k+1,2)=M(k+1)/2;S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6;S(k+1,4)=Y(k+1);%计算三次样条插值函数的系数矩阵endfunction[f]=hs(S,t,t1)%计算某时刻的流速n=length(t);i=1;whilet1>=t(i)i=i+1;endf=S(i,4)+S(i,3)*(t1-t(i))+S(i,2)*((t1-t(i))^2)+S(i,1)*((t1-t(i))^3);endfunctionf1=z(S,t,t1,t2)%计算某段时间内的用水量n=length(t);k1=1;k2=1;m=1;whilet1>=t(k1)m=k1;k1=k1+1;endwhilet2>=t(k2)n=k2;k2=k2+1;endf1=0;ifm==n%该时间段在同一分段函数内k=m;f=S(k,4)*t2+(1/2)*S(k,3)*((t2-t1)^2)+(1/3)*S(k,2)*((t2-t1)^3)+(1/4)*S(k,1)*((t2-t1)^4)-S(k,4)*t1;f1=f;elsefork=m+1:n-1%该时间段在不同分段函数内f=S(k,4)*t(k+1)+(1/2)*S(k,3)*((t(k+1)-t(k))^2)+(1/3)*S(k,2)*((t(k+1)-t(k))^3)+(1/4)*S(k,1)*((t(k+1)-t(k))^4)-S(k,4)*t(k);f1=f1+f;endendk=m;%计算边界的值f=S(k,4)*t(k+1)+(1/2)*S(k,3)*((t(k+1)-t1)^2)+(1/3)*S(k,2)*((t(k+1)-t1)^3)+(1/4)*S(k,1)*((t(k+1)-t1)^4)-S(k,4)*t1;f1=f1+f;k=n;f=S(k,4)*t2+(1/2)*S(k,3)*((t2-t(k))^2)+(1/3)*S(k,2)*((t2-t(k))^3)+(1/4)*S(k,1)*((t2-t(k))^4)-S(k,4)*t(k);f1=f1+f;end%主函数>>[v,t]=ls();X=[t];Y=[v];dx0=0;dxn=0;S=csfit1(X,Y,dx0,dxn);a=z(S,t,32284/3600,39435/3600);b=z(S,t,75021/3600,24);c=3.1415*28.5*28.5*(3550-2697-2697+3175)/100;s1=a+b+c;s2=z(S,t,0,24);1.3.4程序运行调用函数func