数学建模流水问题的实验报告

数学建模实验报告-流水问题

一． 问题描述

三个横截面积为常数A ，高分别为H1，H2，H3的水池内都盛满了水，都由池底一横截面积为B 的小孔放水。设水从小孔流出的速度为v （i ）=sqrt （2*g*h （i ）），求水流空所需的时间。

二． 前提假设

1. 假设在一段极微小的时间间隔dt 内，三个浴缸的高度变化

速率以及三个排水口的排水速率是一个不变化的定值。

2. 排水速率仅与水池高度有关。

3. 排水口的高度为水池最低处，即不会出现因水位低于排水口

而无法排完水的现象。

三． 问题分析

将此问题抽象成数学问题：求初值分别为H1，H2，H3的函数h1，h2，h3随时间变化的函数，以及他们变为0所需的时间（设水池1，2，3的流出速度为v1，v2，v3）。根据任何一个水池内水量在一个时间微元内的减少量等于流出量减去流入量可以得到如下关系：

浴缸 1(/2B A =。

水池2：2*2*1*h A s B h A -∆=∆-∆。两边取极限后得-dh2*A=

ds2*B- dh1*A 。注意到|dh1*A|= |ds1*B|,除以dt 可得(-dh2/dt)*A= (v2-v1)*B ，化简并带入v1的函数表达式可得

(2/)*dh dt A B -=，再代入h1的表达式

(/2B A -可以得到如下的常微分方程

(2/)*(/2))*dh dt A B A B -=。这是一个非线性常微分方程，难以得到解析解（并非不可求，可用待定系数法等求解，但是在此处求出解析解并不是建模的重点，因为即使h2存在解析解，h3也不一定存在，而对于求出排水时间图，求出近似解更为重要），在这里我们采用计算方法中的一些数值计算方法求出几组h2和t2，v2的近似解。

浴缸 3：3*3*2*h A S B h A -∆=∆-∆，由此递推公式可得如下微分

方程：(3/)/dh dt B -=，可以用Euler 法求出近似 解 。

四． 问题求解

对于h1，书上已给出解法，即用一个微分方程求解，在

此不做累述。在代码中，我们用一个数组arrayt 记录时间，一个数组h1记录水位高度，然后对应画图。

本部分Matlab 代码如下：

A=2;

B=1;

H=10;

g=10;

h=H;

temp=0;

t=0;

tf=0

i=1;

h1=[];

h2=[];

h3=[];

temp=[];

tempt=[];

arrayt=[];

arrayt2=[];

arrayt3=[];

sqrtH=sqrt(H);

h1=[h1,H];

arrayt=[arrayt,t];

while (h>0)

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法，具体方法参见数值计算方法第六章），其原理是改进了的欧拉法，将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;g=10;H=10;

if(x<0)

x=0

end

y=-1*B/A*(sqrt(2*g*x)-sqrt(2*g*H)+B/A*g*t)

end

变步长Runge-Kutta法的实现：

tf=2*t;

[tempt,temp]=ode45('funt',[0,tf],H);

解的筛选：

i=1;

while (temp(i)>0)

arrayt2=[arrayt2,tempt(i)];

h2=[h2,temp(i)];

i=i+1;

end

对于h3，浴缸3： -D h3*A= D s3*B- D h2*A，由此递推公式可得如下微分方程：

由于只有三步，且得到的h2的点是离散的，Runge-Kutta法使

用不便。此处我们采用简单一些的欧拉法(Euler法，具体方法

见数值计算方法第六章)进行数值求解。采用此方法可以得到

h3的数值解，即h3关于时间的离散函数。 h3的解不需要筛选。

本部分Matlab代码如下：

Euler法的实现：

j=1;

h3(1)=H/2;

arrayt3(1)=0;

temp3=0;

while(j

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)为初值，带入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);