安徽工业大学
大学数学实验课程设计
姓名：
班级：
任课老师：
数学实验 课程设计
问题提出：
某容器盛满水后，低端直径为0d 的小孔开启（图）。
根据水力学知识，当水面
高度h 时，水冲小孔中流出的速度v =(g 为重力加速度，0.6为孔口的收缩系数)。
⑴若容器为倒圆锥形（如图1），现测得容器高和上底面直径均为1.2m ，小孔直径为3cm ，问水从小孔中流完需要多长时间；2min 水面高度是多少。
⑵若容器为倒葫芦形（如图2），现测得容器高为1.2m ，小孔直径为3cm ，有低端（记作x=0）向上每隔0.1m 测出容器的直径D （m ）如表所示，问水从小孔中流完需要多少时间；2min 时水面的高度是多少。
图1 ： 图2：
问题分析：
（1） 倒圆锥形容器流水问题中随时间t 液面高度h 也在变化，同时水的流速也
在变化，再写变化难以用普通的方程进行模拟求解，考虑建立常微分方程竟而代入数值求解。
水面的直径等于液面的高度。
可以建立容器中水流失的液面高度对时间t 的变化率。
假设t 时，液面的高度h ，此时水的流速流量Q 为：00.6(/4)d π
； 则
在t ∆时间内液面下降高度为h ∆，可得到关系式：220(
)2
4
d dt h dh π
=
；
由此可知水下降h ∆
时需要的时间：20
40.6
4
h dh
t d π
π
∆=
=
根据此关系式知道。
（2） 在第二问中，考虑倒葫芦形容器时因为他的高度h 不同容器直径D 变化
没有规律可循，同第一题相比我们只知道他的一些数值，这就需要我们建立高度h 和容器直径D 之间的关系矩阵，然后再欧拉方程和龙格—库塔方法找出时间t 和液面高度之间的分量关系。
由（1）可同理推知：假设在时间t 时，液面高度为h ，此时流量
为
2
00.6(/4)d π；经过t ∆时，液面下降h ∆，若我们取的t 是在t(n)和t(n+1)
之间的某一时刻，于是就可在误差范围内得到
(1)()t n t n t +=+∆；可以得
到
204
(1)()0.64
h d h
dt t n t n d π
π
=+-=-
= ；
建立模型：
（1） 在试验中我们不考虑圆锥的缺省对流水的影响，以及其他外界因素和玻璃
的毛细作用，试验中水可以顺利流完。
实验中重力加速度g=9.82
/m s ；倒圆锥的液面最初高度为H=1.2m ，液面直径D=1.2m=0.03，小孔的直径为
0d =0.03m ；
接上文中分析结论代入数据：即在T 时间内将1.2m 的液面高度放完， （matlab 不支持一些运算符号，故用matlab 运算格式）
dt=-((pi/4)h^2*dh)/(0.6*(pi/4)*d^2*sqrt(gh))=-(h^1.5*dh)/(0.6*d^2*sqrt(g)) h 是由0→1.2m 对t 积分 用matlab 计算上式 编辑文件：a1.m ， d0=0.03; g=9.8; syms h
t=(h^1.5)/(0.6*d0^2*sqrt(g)); T=int(t,0,1.2); eval(T)
运行结果： >> a1
ans =
373.2556
结果：水从倒圆锥中流完需要373.26s；
2mine之后液面的高度为h1;
373.26-120= h1^2.5/(1.5*d^2*(g)^0.5)=153.26
可知h1=((1.5*d^2*sqrt(g))*153.26)^(2/5);
Malab计算
>> g=9.8
g =
9.8000
>> d0=0.03
d0 =
0.0300
>> h1=((1.5*d0^2*sqrt(g))*153.26)^(2/5)
h1 =
0.8405
即2mine之后液面的高度为0.84m；
上述运行结果可知：谁需要373.26s流完，2mine之后液面高度为0.84m；
（2）在与（1）同样的条件下，倒葫芦形容器的液面最初高度H=1.2m，小孔的直径为d0=0.03m，液面直径和液面高度关系如表。
在分析中已经讨论出Δt和Δh的关系。
dt=t(k+1)-t(k)=-((pi/4)h^2*dh)/( 0.6(pi/4)d^2*sqrt(gh))=-(h^1.5*dh)/(0.6d^2*sqrt(g)) k1=0.15*sqrt(g*(x(n)))*d^2/(-43.6359*x(n)^8+213.0457*x(n)^7-414.873*x(n)^6+41 0.2075*x(n)^5-218.8936*x(n)^4+62.553*x(n)^3-8.3215*x(n)^2+0.49619*x(n)+.0148 92)^2;
k2=0.15*sqrt(g*(x(n))-h*k1)*d^2/(-43.6359*(x(n)-h*k1)^8+213.0457*(x(n)-h*k1)^7 -414.873*(x(n)-h*k1)^6+410.2075*(x(n)-h*k1)^5-218.8936*(x(n)-h*k1)^4+62.553*( x(n)-h*k1)^3-8.3215*(x(n)-h*k1)^2+0.49619*(x(n)-h*k1)+.014892)^2;
x(n+1)=x(n)-h*(k1+k2)/2;
matlab程序编写及计算：
g=9.8;
d=0.03;
k1=0;
k2=0;
h=0.4;
x(1)=1.2;
for n=1:1000
k1=0.15*sqrt(g*(x(n)))*d^2/(-43.6359*x(n)^8+213.0457*x(n)^7-414.873*x (n)^6+410.2075*x(n)^5-218.8936*x(n)^4+62.553*x(n)^3-8.3215*x(n)^2+0.4 9619*x(n)+.014892)^2;
k2=0.15*sqrt(g*(x(n))-h*k1)*d^2/(-43.6359*(x(n)-h*k1)^8+213.0457*(x(n )-h*k1)^7-414.873*(x(n)-h*k1)^6+410.2075*(x(n)-h*k1)^5-218.8936*(x(n) -h*k1)^4+62.553*(x(n)-h*k1)^3-8.3215*(x(n)-h*k1)^2+0.49619*(x(n)-h*k1 )+.014892)^2;
x(n+1)=x(n)-h*(k1+k2)/2;
end
x(300)
t=0:h:1000*h;
plot(t,x);
axis([0,400,0,1.21]);
p=polyfit(x,t,5);
T=real(polyval(p,0))
运行结果：
>> a3
ans =
1.0278
Warning: Imaginary parts of complex X and/or Y arguments ignored
> In a3 at 14
T =
396.5254
运行结果得知：倒葫芦形容器的水大约396.53s 流完，2mine 后液面高度为1.03m 。
实验总结：
在研究实际问题时，建立模型找出微分方程是分析客观对象变化规律的有力工具，但不是所有问题多可以建立解析表达式（如（2）问），即使获得了解析表达式也十分复杂，因此就必须通过数值求解的方法算出微分方程的某些离散点处得近似解析，进而分析微分方程反映客观问题。
指点迷津：在上述问题的（1）中仍有疑问，之前思考时建立模型是倒圆锥内剩余的水的体积变化量dV 对时间dt 的关系，即dV/dt=Q=vt(v =
假设t 时间内水流了V
V=22331H 1()()(())(
)(())32321212H h t V H H h t H H h t ππ
ππ-=--=--； dV/dt=2(())()4
H h t h t π
-’=求次微分方程：
>> h=dsolve('Dh=(0.6*sqrt(9.8*h))/((pi/4)*(1.2-h)^2)','h(0)=1.2','t') h =
RootOf(X224^5 - 8*X224^4 + (152*X224^3)/5 - (288*X224^2)/5 + (1296*X224)/25 - (55296*pi^2 + 1102500*t^2 - 201600*pi*6^(1/2)*t)/(3125*pi^2), X224)
RootOf(X233^5 - 8*X233^4 + (152*X233^3)/5 - (288*X233^2)/5 + (1296*X233)/25 - (55296*pi^2 + 1102500*t^2 + 201600*pi*6^(1/2)*t)/(3125*pi^2), X233)
函数没有数值解怎么回事。
是思路错了还是关系济啉错了此题是不是存在一些问题？。