1、数学建模实验报告-流水问题一 问题描述三个横截面积为常数 A，高分别为 H1，H2，H3 的水池内都盛满了水，都由池底一横截面积为 B的小孔放水。设水从小孔流出的速度为 v（i）=sqrt（2*g*h（i） ） ，求水流空所需的时间。二 前提假设1. 假设在一段极微小的时间间隔 dt内，三个浴缸的高度变化速率以及三个排水口的排水速率是一个不变化的定值。 2. 排水速率仅与水池高度有关。3. 排水口的高度为水池最低处，即不会出现因水位低于排水口而无法排完水的现象。三 问题分析将此问题抽象成数学问题：求初值分别为 H1，H2，H3 的函数 h1，h2，h3 随时间变化的函数，以及他们变为 0所需的
2、时间（设水池 1，2，3 的流出速度为 v1，v2，v3） 。根据任何一个水池内水量在一个时间微元内的减少量等于流出量减去流入量可以得到如下关系：浴缸 1： 。(/2)*hHBAgt水池 2： 。两边取极限后得-dh2*A= *1shds2*B- dh1*A。注意到|dh1*A|= |ds1*B|,除以 dt可得(-dh2/dt)*A= (v2-v1)*B，化简并带入 v1的函数表达式可得，再代入 h1的表达式(2/)*(21)*dhtAghB可以得到如下的常微分方程/)HBt。这是一个非(/)(/2)*t HAgtB线性常微分方程，难以得到解析解（并非不可求，可用待定系数法等求解，但是在此处
3、求出解析解并不是建模的重点，因为即使 h2存在解析解，h3 也不一定存在，而对于求出排水时间图，求出近似解更为重要） ，在这里我们采用计算方法中的一些数值计算方法求出几组 h2和 t2，v2 的近似解。浴缸 3： ，由此递推公式可得如下微分3*2*hASBhA方程： ，可以用 Euler法求出近(/)/(3)dtg似解 。四 问题求解对于 h1，书上已给出解法，即用一个微分方程求解，在此不做累述。在代码中，我们用一个数组 arrayt记录时间，一个数组 h1记录水位高度，然后对应画图。本部分Matlab代码如下： A=2;B=1;H=10;g=10;h=H;temp=0;t=0;tf=0i=1
4、;h1=;h2=;h3=;temp=;tempt=;arrayt=;arrayt2=;arrayt3=;sqrtH=sqrt(H);h1=h1,H;arrayt=arrayt,t;while (h0)t=t+0.01;arrayt=arrayt,t;h=sqrt(H)-(B/2/A)*(sqrt(2*g )*t);h1=h1,h*h;end对于 h2，由于我们得到了一个递推公式，并得到了一个常微分方程，本来可以直接求解，但由于此方程属于不可解的一类，故仅求出近似解。由于第二步的数值会对后续的浴缸产生影响，我们选用了精度较高的 Runge-Kutta法（龙格-库塔法，以下简称 RK法，具体方法参
5、见数值计算方法第六章），其原理是改进了的欧拉法，将 x的区间分成很多小间隔来分段应用欧拉法，matlab 中给出了实现变步长 Runge-Kutta法的函数：t,y=ODE45(odefun,tspan,y0)其中 t，y 用于输出广义时间与相空间向量，odefun 为要求的微分方程，tspan为广义时间区间，y0 为初值向量。由于 Runge-Kutta法是一种近似解法，仅在局部区间收敛，所以会有虚解及负解出现，而求解 h3时会需要用到 h2的解，所以我们要对得到的解进行筛选，只留下非负实数解。本部分 Matlab代码如下：funt.M：function y=funt(t,x)A=2;B=1
6、;g=10;H=10;if(x0)arrayt2=arrayt2,tempt(i);h2=h2,temp(i);i=i+1;end对于 h3，浴缸 3： -Dh3*A= Ds3*B- Dh2*A ，由此递推公式可得如下微分方程：由于只有三步，且得到的 h2的点是离散的，Runge-Kutta 法使用不便。此处我们采用简单一些的欧拉法(Euler 法，具体方法见数值计算方法第六章)进行数值求解。采用此方法可以得到 h3的数值解，即 h3关于时间的离散函数。 h3 的解不需要筛选。本部分 Matlab代码如下：Euler法的实现： j=1;h3(1)=H/2;arrayt3(1)=0;temp3=
7、0;while(ji-1)temp3=h3(j)+(h2(j)-h2(j+1)*(-1*B/A*(sqrt(2*g*h3(j)-sqrt(2*g*h2(j);h3=h3,temp3;arrayt3=arrayt3,arrayt2(j+1);j=j+1;end至于求解水流完的时间，可采用遍历 h1，h2，h3 数组的方法，直到找到 0点所在，然后取出对应的 t即可。在计算 h3时要注意对模型进行修正，由于 h3所有的分量均由 h2参与迭代运算产生，故当 h2为 0之后，h3 也停止了迭代。故在我们绘制的图中 h3没有到 0就停止了。可以改进为，当 t时刻 h2为 0之后，h3 采用 h3(t)为
8、初值，带入 h1的表达式进行计算来计算零点。本部分 Matlab代码如下：t1=0;t2=0;t3=0;i=length(arrayt);t1=arrayt(i);i=length(arrayt2);t2=arrayt2(i);i=length(arrayt3);t3=arrayt3(i)h3(i)t_extra=A/B*sqrt(2*h3(i)/g);t3=t3+t_extra;t1t2t3五 数值模拟结果取 H0=10，讨论以下初始情况下三个浴缸水位的关系：分别是三个浴缸初始量相等的情况，一多两少的情况和一少二多的情况.H1,H2,H3,t1,t2,t3 分别表示三浴缸的初始水量以及全部流出时间排水时间统计如下： H1=10,H2=10,H3=10H1=10,H2=10,H3=5H1=10,H2=5,H3=10H1=5,H2=10,H3=10H1=5,H2=5,H3=10H1=5,H2=10,H3=5H1=10,H2=5,H3=5
