实验4-常微分方程数值解(2012080076)2.0分

得到表：
时间 高度 /s h/m 0 0 0.1 0.06528 0.2 0.26133 0.3 0.58841 0.4 1.04681 0.5 1.63678 0.6 2.35857 0.7 3.21244 0.8 4.1986 0.9 5.3173 时间 /s 10 10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9 高度 h/m 659.3598 672.2493 685.2485 698.3566 711.5732 724.8975 738.329 751.867 765.5109 779.26 时间 /s 30 30.1 30.2 30.3 30.4 30.5 30.6 30.7 30.8 30.9 时间 /s 4609.958 50 4633.462 50.1 4656.982 50.2 4680.519 50.3 4704.073 … 50.4 4727.642 50.5 4751.228 50.6 4774.83 50.7 4798.447 50.8 4822.08 50.9 高度 h/m 时间 /s 9560.388 59.1 9586.187 59.2 9611.995 59.3 9637.812 59.4 9663.639 … 59.5 9689.476 59.6 9715.322 59.7 9741.178 59.8 9767.043 59.9 9792.918 60 高度 h/m 高度 h/m 11946.71 11973.36 12000.02 12026.69 12053.36 12080.05 12106.74 12133.45 12160.16 12186.88
得到表：
时间 /s 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
加速度 a/ （m/s^2)
时间 /s 10 10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9
加速度 a/ （m/s^2)
时间 /s 30 30.1 30.2 30.3 30.4 30.5 30.6 30.7 30.8 30.9
t=0:0.1:60; %设定起终点及时间间隔 v0=[0,0]; %输入 v 与 h 的初始值 [t,v]=ode45(@huojian,t,v0);% 调用 ode45 计算 a=(32000-0.4*v(:,1).^2)./(1400-18*t)-9.81; %解出各个时刻的加速度 plot(t,v(:,2));xlabel('t/s'),ylabel('h/m');title('引擎关闭前 h-t 关系图'); title('引擎关闭前 h-t 关系图');%画出高度时间关系图 [t,v(:,2)]
引 擎 关 闭 前 h-t 关 系 图 14000
ห้องสมุดไป่ตู้
12000
10000
8000
h/m
6000 4000 2000
0
0
10
20
30 t/s
40
50
60
引 擎 关 闭 前 v-t 关 系 图 300
250
200
v/(m/s)
150
100
50
0
0
10
20
引 擎 关 闭 前 a-t 关 系 图 14
30 t/s
…
11.0177 10.95547 10.89279 10.82969 10.76619 10.70228 10.638 10.57335 10.50835 10.44301
…
1.745256 1.731881 1.71871 1.705734 1.692945 … 1.680334 1.667892 1.65561 1.643478 1.631486
0.95384 0.952645 0.951442 0.950235 0.949034 … 0.947844 0.946671 0.945522 0.944402 0.943316
0.917757 0.918349 0.918864 0.919278 0.919563 0.919692 0.919632 0.919352 0.918817 0.91799
得到图：
引 擎 关 闭 前 a-t 关 系 图 14
12
10
a/(m/s 2)
8
6
4
2
0 0
10
20
30 t/s
40
50
60
数学实验 2012080076 朱燚豪
从上面的分析知道：在引擎关闭的瞬间（前） ，火箭的高度、速度及加速度为：
h=12186.88m, v=267.2407 m/s ,a=0.91799m/s^2 三者随时间变化图形：
问题分析及其求解：
整个过程可以分为两个过程： 引擎关闭前以及引擎关闭后， 问题中所求的引擎关闭瞬间应该 是第一个阶段恰好结束时，按照物理中只是，火箭达到最高点时应该为速度为0时。
（1） 引擎关闭前
假设火箭在上升过程中，重力加速度 g 不随高度而变化，即固定 g = 9.8m/s^2。 (1)从火箭开始上升到引擎关闭：设火箭质量为 m，高度为 h，速度为 v，加 速度为 a，阻力为 f： dh m = m(t) = 1400 − 18t， = v，f = 0.4v 2 dt 由牛顿第二定律可得： dv F总 F − mg − f 32000 − 9.8(1400 − 18t) − 0.4v 2 a= = = = dt m m 1400 − 18t 综上可得： dh = v； dt dv 32000 − 0.4v 2 = − 9.8 dt 1400 − 18t 初值条件为：v = 0 m/s，h = 0 m；定义域为：0 ≤ t ≤
数学实验 2012080076 朱燚豪
实验 4 常微分方程数值解
化学工程系 化22班 朱燚豪 2012080076
实验目的：
1. 掌握用 MATLAB 软件求微分方程初值问题数值解的方法； 2. 通过实例学习用微分方程模型解决简化的实际问题； 3. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
function d=huojian2(t,v) %建立命名为huojian2的函数M文件 g=9.81;m0=1400;k1=18;k2=0.4; d=[(-k2*v(1)^2)/(m0-k1*t)-g;v(1)]; %以向量形式表示微分方程 end
v0=[v(601,1),v(601,2)]; %设定起终点及时间间隔 t=60:0.01:80; %设定起终点及时间间隔 [t,v]=ode45(@huojian2,t,v0); % 调用 ode45 计算 [t,v(:,1)]
得到表：
时间 速度 v/ /s （m/s) 60 267.241 60.01 266.253 60.02 265.271 60.03 264.296 60.04 263.326 60.05 262.362 60.06 261.404 60.07 260.452 60.08 259.505 60.09 258.564 时间 /s 64 64.01 64.02 64.03 64.04 64.05 64.06 64.07 64.08 64.09 速度 v/ （m/s) 86.0164 85.79919 85.5825 85.36633 85.15067 84.93552 84.72088 84.50675 84.29312 84.08 时间 /s 70.39 70.4 70.41 70.42 70.43 70.44 70.45 70.46 70.47 70.48 速度 v/ 时间 速度 v/ （m/s) /s （m/s) 0.45481 76 -472.786 0.356728 76.01 -502.476 0.258649 76.02 -536.444 0.160571 76.03 -575.744 0.062494 … 76.04 -621.629 -0.03558 76.05 -675.714 -0.13366 76.06 -740.462 -0.23174 76.07 -818.999 -0.32982 76.08 -916.976 -0.4279 76.09 -1042.78
加速度 a/ （m/s^2)
时间 /s 50 50.1 50.2 50.3 50.4 50.5 50.6 50.7 50.8 50.9
加速度 a/ （m/s^2)
时间 /s 59.1 59.2 59.3 59.4 59.5 59.6 59.7 59.8 59.9 60
加速度 a/ （m/s^2)
13.0471 13.0761 13.1041 13.1312 13.1574 13.1827 13.207 13.2303 13.2527 13.2741
…
…
得到图：
引 擎 关 闭 前 v-t 关 系 图 300
250
200
v/(m/s)
150
100
50
0
0
10
20
30 t/s
40
50
60
数学实验 2012080076 朱燚豪
再输入：
plot(t,a);xlabel('t/s'),ylabel('a/(m/s^2)'); title('引擎关闭前 a-t 关系图'); %画出加速度时间关系图 [t,a]
40
50
60
12
10
a/(m/s 2)
8
6
4
2
0 0
10
20
30 t/s
40
50
60
数学实验 2012080076 朱燚豪
（2） 引擎关闭后（到达最高点前）:
0.4v 2 a = −9.8 − ； 320 可得到新的微分方程组： dv 0.4v 2 = −9.8 − ； dt 320 dh = v； dt 相应的初值为：h = 12190m，v = 267.26m/s，经过估计，火箭速度为0 的时刻不会超过80s，故定义域为60s ≤ t ≤ 80s。 编写 MATLAB 程序：
1080 18
= 60 s。
编写 MATLAB 程序：
数学实验 2012080076 朱燚豪 function d=huojian(t,v) %建立命名为huojian的函数M文件 g=9.81;F=32000;m0=1400;k1=18;k2=0.4; d=[(F-k2*v(1)^2)/(m0-k1*t)-g;v(1)]; %以向量形式表示微分方程 end
实验内容：
问题 1： 问题陈述：小型火箭初始重量为 1400kg，其中包括 1080kg 燃料。火箭竖直向
上发射时燃料燃烧率为 18kg/s，由此产生 32000N 的推力，火箭引擎在燃烧用尽 时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为 0.4kg/m，求引 擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度， 并画出高度、速度、加速度随时间变化的图形。
…
…
得到图：
引 擎 关 闭 前 h-t 关 系 图 14000
12000
10000
8000
h/m
6000 4000 2000
0
0
10
20
30 t/s
40
50
60
数学实验 2012080076 朱燚豪
再输入：
plot(t,v(:,1));xlabel('t/s'),ylabel('v/(m/s)'); title('引擎关闭前 v-t 关系图'); %画出速度时间关系图 [t,v(:,1)]
得到表：
时间 速度 v/ /s （m/s) 0 0 0.1 1.30617 0.2 2.61519 0.3 3.92696 0.4 5.2414 0.5 6.55841 0.6 7.8779 0.7 9.19977 0.8 10.5239 0.9 11.8503 时间 /s 10 10.1 10.2 10.3 10.4 10.5 10.6 10.7 10.8 10.9 速度 v/ （m/s) 128.357 129.4557 130.5481 131.6342 132.714 133.7874 134.8544 135.9149 136.969 138.0165 时间 /s 30 30.1 30.2 30.3 30.4 30.5 30.6 30.7 30.8 30.9 速度 v/ 时间 （m/s) /s 234.8536 50 235.0253 50.1 235.1957 50.2 235.3649 50.3 235.5328 … 50.4 235.6996 50.5 235.8653 50.6 236.0299 50.7 236.1935 50.8 236.3561 50.9 速度 v/ 时间 （m/s) /s 257.9636 59.1 258.0603 59.2 258.157 59.3 258.2537 59.4 258.3503 … 59.5 258.4468 59.6 258.5432 59.7 258.6395 59.8 258.7357 59.9 258.8317 60 速度 v/ （m/s) 266.427 266.5166 266.6064 266.6962 266.7863 266.8766 266.9671 267.058 267.1491 267.2407