实验4-常微分方程数值解

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