数学实验课程设计【常微分方程数值解】【数学建模倒葫芦流水问题】

4 4
h 2 dh d02 gh

建立模型：
（1） 在试验中我们不考虑圆锥的缺省对流水的影响， 以及其他外界因素和玻璃 的毛细作用，试验中水可以顺利流完。实验中重力加速度 g=9.8 m / s ；倒 圆锥的液面最初高度为 H=1.2m， 液面直径 D=1.2m=0.03， 小孔的直径为 d0 =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 =
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))
/ 4)d0 2 gh ； 则

4 2 间内液面下降高度为 h ， 可得到关系式：0.6 gh (

由此可知水下降 h 时需要的时间： t
0.6

4 4
h 2 dh d 0 2 gh

h h 0.6d 0 2 g
；
根据此关系式知道。 （2） 在第二问中，考虑倒葫芦形容器时因为他的高度 h 不同容器直径 D 变化 没有规律可循，同第一题相比我们只知道他的一些数值，这就需要我们建 立高度 h 和容器直径 D 之间的关系矩阵，然后再欧拉方程和龙格—库塔 方法找出时间 t 和液面高度之间的分量关系。 由 （ 1 ）可 同 理推 知 ：假 设 在时 间 t 时 ，液 面 高度 为 h， 此 时流 量 为
2
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 程序编写及计算：
表： x/m 0 D/m
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1.0
1.1
1.2
0.03 0.05 0.08 0.14 0.19 0.33 0.45 0.68 0.98 0.10 1.20 1.13 1.00
问题分析：
（1） 倒圆锥形容器流水问题中随时间 t 液面高度 h 也在变化，同时水的流速也 在变化，再写变化难以用普通的方程进行模拟求解，考虑建立常微分方程 竟而代入数值求解。水面的直径等于液面的高度。可以建立容器中水流失 的液面高度对时间 t 的变化率。 假设 t 时，液面的高度 h，此时水的流速流量 Q 为： 0.6(
安徽工业大学 数学实验课程设计
班级： 任课老师： 姓名： 学号：
数学实验 课程设计
问题提出：
某容器盛满水后，低端直径为 d 0 的小孔开启（图） 。根据水力学知识，当水面 高度 h 时，水冲小孔中流出的速度 v 0.6 gh (g 为重力加速度，0.6 为孔口的收 缩系数)。 ⑴若容器为倒圆锥形（如图 1） ，现测得容器高和上底面直径均为 1.2m，小孔直 径为 3cm，问水从小孔中流完需要多长时间；2min 水面高度是多少。 ⑵若容器为倒葫芦形（如图 2） ，现测得容器高为 1.2m，小孔直径为 3cm，有低 端（记作 x=0）向上每隔 0.1m 测出容器的直径 D（ m）如表所示，问水从小孔中 流完需要多少时间；2min 时水面的高度是多少。 图1 ： 图 2：
运行结果：
>> a3 ans = 1.0278 Warning: Imaginary parts of complex X and/or Y arguments ignored > In a3 at 14 T= 396.5254
1.2
1
0.8
0.6
0.4
0.2
0
0
50
100
150
200
250
300
函数没有数值解怎么回事。 是思路错了还是关系建立错了此题是不是存在一些问 题？
0.6( / 4)d0 2
gh
；经过 t 时，液面下降 h ，若我们取的 t 是在 t(n)和 t(n+1)
之 间 的某 一时 刻， 于是 就可 在误 差范 围内 得到
t (n 1) t (n) t
h h ； 0.6d 0 2 g
； 可 以得 到

dt t (n 1) t (n) 0.6
1 H 1 H h( t) 2 3 ) H ( H h( t)) 3 ； V= V H ( ) 2 ( H h(t )) ( 3 2 3 2 12 12
dV/dt=

4
( H h(t )) 2 h(t )’= 0.6 gh(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)
350
400
运行结果得知： 倒葫芦形容器的水大约396.53s流完， 2mine后液面高度为1.03m。
实验总结：
在研究实际问题时， 建立模型找出微分方程是分析客观对象变化规律的有力工具， 但不是所有问题多可以建立解析表达式（如（2）问） ，即使获得了解析表达式也 十分复杂， 因此就必须通过数值求解的方法算出微分方程的某些离散点处得近似 解析，进而分析微分方程反映客观问题。 指点迷津：在上述问题的（1）中仍有疑问，之前思考时建立模型是倒圆锥内剩 余的水的体积变化量dV对时间dt的关系，即dV/dt=Q=vt( v 0.6 gh ); 假设t时间内水流了V