.

实验4 常微分方程数值解

化工系毕啸天 2010011811

【实验目的】

1. 练习数值微分的计算；

2. 掌握用MATLAB 软件求微分方程初值问题数值解的方法； 3. 通过实例学习用微分方程模型解决简化的实际问题；

4. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。

【实验内容】 题目3

小型火箭初始重量为1400kg，其中包括1080kg 燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N 的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。

3.1 燃料燃烧过程物理模型分析

设火箭质量为m，高度为h，速度为v，加速度为a，火箭推力为F，重力加速度为g，阻力为f。

1.由火箭上总共携带燃料1080kg，燃料燃烧率为18kg/s，可知火箭上升时间t=60s时，燃料全部烧尽。

2.由阻力正比于速度的平方，比例系数0.4kg/m，可知阻力表达式为f=0.4v2。 3.由于燃料燃烧，火箭的质量是时间的函数，易知m(t)=m0-18t

4.

5.根据牛顿第二运动定律，有 。代入数据有

解出

由以上5条分析，我们得到了一个常微分方程组：

初值条件为：v0=0,h(0)=0,t 60s.

3.2 程序代码

根据常微分方程组的初值问题，在MATLAB中计算数值解。 记x(1) = h，x(2)= v，x =(x(1), x(2))T 首先编写M文件

function dx = Rocket(t,x)

dx=[x(2);(32000-0.4*x(2)^2)/(1400-18*t)-9.8];%以向量形式表示微分方程 end

ts=0:60 %终点时间为60s，步长定义为1即可

1页

.

x0=[0,0];

[t,x]=ode45(@Rocket,ts,x0); [t,x]

plot(t,x(:,1)),grid, title('图1.高度-时间') xlabel('t/s') ylabel('h/m') pause

plot(t,x(:,2)),grid,

title('图2.速度-时间') xlabel('t/s') ylabel('v/(m/s)') pause

a=(32000-0.4*x(:,2).^2)./(1400-18*t)-9.8 plot(t,a),grid,

title('图3.加速度-时间') xlabel('t/s')

ylabel('a/(m/s2)') [t,x(:,1),x(:,2),a]

2页

.

由MATLAB计算得到从开始上升到关闭引擎瞬间的情况如下表：

3页

.

t 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 h 0 6.5737 26.444 59.762 106.57 166.79 240.27 326.72 425.79 536.99 659.8 793.63 937.85 1091.8 1254.7 1425.9 1604.8 1790.8 1983.1 2181.2 2384.5 2592.4 2804.5 3020.6 3240.1 3462.7 3687.9 3915.6 4145.6 4377.8 4611.9 4847.7 5085 5323.8 5564.1 5805.8 6048.7 6292.9 6538.1 6784.5 7032 7280.5 v 0 13.189 26.577 40.062 53.535 66.89 80.021 92.829 105.22 117.11 128.43 139.14 149.18 158.55 167.23 175.22 182.55 189.22 195.27 200.75 205.7 210.18 214.19 217.79 221.01 223.92 226.56 228.97 231.14 233.11 234.91 236.57 238.14 239.61 240.99 242.28 243.5 244.68 245.83 246.96 248.05 249.1 4页

a 13.057 13.305 13.453 13.497 13.433 13.261 12.985 12.612 12.152 11.617 11.021 10.38 9.7083 9.0209 8.3309 7.6502 6.9901 6.3593 5.7646 5.2095 4.6946 4.222 3.7943 3.412 3.073 2.7726 2.5044 2.2677 2.0633 1.8898 1.7433 1.6178 1.5062 1.4095 1.3293 1.265 1.2139 1.1708 1.1303 1.0947 1.0663 1.0456 .

42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 7530.2 7780.9 8032.5 8285.1 8538.8 8793.4 9049 9305.6 9563.1 9821.5 10081 10341 10603 10865 11128 11392 11657 11923 12190 250.12 251.14 252.15 253.16 254.15 255.12 256.07 257.03 257.99 258.95 259.9 260.83 261.75 262.67 263.61 264.54 265.46 266.35 267.26 1.0308 1.0178 1.0024 0.98757 0.97632 0.96964 0.96629 0.96239 0.95272 0.94118 0.93373 0.93278 0.93633 0.93695 0.92585 0.91379 0.91063 0.91612 0.91701 由数据可以看出，引擎关闭瞬间，火箭的高度为h=12190m，速度v=267.26m/s，加速度 a=0.91701m/s2。

3.3燃烧耗尽后物理模型分析

当引擎关闭后，火箭由重力作用上升，燃料耗尽后火箭质量为320kg。由牛顿第二定律 即 可再次列出微分方程如下:

初值条件即为(1)中的末态条件,t>60s.

3.4程序代码

记y(1) = h，y(2)= v，y=(y(1), y(2))T

首先编写M文件

function dy = Rocket2(t,y)

dy=[y(2);-9.8-0.4*y(2).^2/320];%以向量形式表示微分方程 end

ts=0:60 %终点时间为60s，步长定义为1即可 x0=[0,0];

[t,x]=ode45(@Rocket,ts,x0); [t,x]

5页