ODE初值问题实验报告
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
然后根据上两个阶段分别写函数rocket2.m和rocket3.m。具体实现和测试程序可参见代 码。不再赘述。
3.3实验结果
2 / 25
3 / 25
上面三种颜色代表三种不同的状态。 所有的数据都储存在T1.mat，T2.mat，T3.mat三个matlab数据文件中。第一列为时间， 第二列为高度，第三列为速度，第四列为加速度。得到如下近似结果： 关闭引擎时的速度为267.257m/s，加速度a=0.9170 m / s2 ，高度为12190m。 火箭达到最高点即第二阶段结束时，高度为13115m，加速度为-9.8 m / s2 。
6.2实验过程 照着前面列出的常微分方程组写出函数测试即可。需要注意的是由于这个过程计算时 间过程，故使用解刚性方程的方法来解。得到结果为：
5 / 25
蓝色为横向方向，红色为竖直方向。
6 / 25
由数据可以得到当t=66.66s时，x=0.00697m,y=-1.9451e-6m。这时小船已经到达了目 的地B。 解析解得话,由前面的解析公示带入数据
dh v dt ,初始条件为上一阶段结束状态。 2 dv 9.8 0.4v 1400 1080 dt
最后一个阶段为火箭下降阶段。这一阶段火箭受力为向上的空气阻力和重力。常微分方程 组如下：
dh v dt ,初值条件为上阶段结束状态。 0.4v 2 dv 9.8 1400 1080 dt
列出常微分方程组：
1 / 25
dh v(0) 0 v dt ，初值为 h(0) 0 2 dv Fp 0.4v 9.8 t 60 s dt 1400 18t
3.2 实验过程 列出常微分方程组，就可以使用龙格-库塔方法实现算法。编写函数。具体为:
（2）1、r1 = 2;r2 = 1;n1 = 100;n2 = 100;s1 = 0.5;s2 = 2;
可见甲的数量上升更快，乙灭绝得也更快。
14 / 25
2、r1 = 1;r2 = 2;n1 = 100;n2 = 100;s1 = 0.5;s2 = 2;
乙刚开始上升阶段的最大种群数量有所上升,但不影响最后的结果.
V1=2m/s 时
可以看出小船永远也到不了对岸 B 点。解析解得:
11 / 25
dx v2 x x v1 v2 (1 ) 0 。小船最后静止在(50,0)处。 2 2 2 dt x y x y2
题目 9： 两种群相互竞争模型如下：
x y (t ) r1 x(1 s1 ), x n1 n2 x y y (t ) r2 y (1 s2 ), n1 n2
7、s1=1.5,s2=0.7
21 / 25
情况几乎翻转,乙趋向饱和,甲趋向灭绝. 这说明 r1 和 r2 是甲乙种群的固有增长率， 他们决定了种群在开始一段时间内数量 增长的快慢，固有增长率大的种群增长更快。同时决定了到达稳态的时间快慢，固有增 长率大的种群到达稳态时间短。n1 和 n2 是甲乙两种群的最大容量，决定了某一种群最 终的稳态值。x0 和 y0 是甲乙两种群的初始值，决定了物种曲线的起点，但不影响最终 的稳态。s1 是对于供养甲的资源而言，单位数量乙的消耗为单位数量甲消耗的 s1 倍。 s2 的意义与此相反。若 s1<1<s2，则甲物种会抢夺乙物种的资源，乙物种并不抢夺甲 物种的资源，最终乙没有充足的食物，从而灭绝，甲乙的地位是对称的，反过来灭绝的 是甲，而乙会达到饱和。这些都非常符合生态学常识。说明了实验的正确性。
1 1 1 1 p (cy ) k (cy ) k ,代入前式整理得 x c k y1k c k y1 k ,c是待定系数。代入初值 2 2 2 2 1 1 1 k 1k 1 1 k 1 k y （- ）y 。 ( x, y) (0, d ) ，得 c , x （- ） d 2 d 2 d
附代码列表： T3 Rocket1.m Rocket2.m Rocket3.m
24 / 25
Rocket_test.m T1.mat T2.mat T3.mat T6 Boat.m Boat_inl.m Boat_test.m Boat_test_3.m T9 Comp.m Comp_test.m compA.mat
题目6 一只小船渡过宽为 d 的河流，目标是起点 A 正对着的另一岸 B 点。已知河水流速 v1 与船在静水中的速度 v2 之比为 k。 （1） 建立描述小船航线的数学模型，求其解析解； （2） 设 d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需的时间、任意时刻小船 的位置及航行曲线，作图，并与解析解比较。
17 / 25
将乙的种群数量上限改为 200 后，只是在刚开始很短的时间内数量超过了甲，但最后仍然 是灭绝，甲仍然是饱和。
5、更改初值,让甲改为原来的 10 倍 100
18 / 25
乙都没有上升一直减少，很快灭绝。
6、将乙初值改为 100
19 / 25
乙一直下降最后灭绝。甲一直上升最后饱和。
20 / 25
（4）s1=0.8(<1)，s2=0.7(<1)
22 / 25
S1<1,s2<1。甲乙不会互相抢夺对方资源，从而基本稳定，但明显乙对资源更有竞 争力。所以乙最后上升到更大的稳定值。
s1=1.5(>1)，s2=1.7(>1)
23 / 25
当 s1>1 且 s2>1 时，两物种对彼此的资源都有抢夺性。因为 s1 略小于 s2，所以 甲物种对乙物种资源的掠夺性更大， 导致在一段时间后甲物种数量最终饱和， 乙物种几 乎灭绝。 从上面的实验可以看出,对于两种处于竞争关系的物种来说,s1 和 s2 起着非常重要 的作用。
( x, y) (0, 100) ，得 c=-0.01,k=0.5 故
1 1 x (0.01) 0.5 y 0.5 (0.01) 0.5 y 0.5 2 2
写出函数 boat_inl.m 计算得结果：
通过比较可以看出我们前面用的数值方法精确性还是非常不错的。
(3)改变 v1=0 结果如下：
替换变量
x 2 y 2 x v1 x x ( )2 1 y y v2 y y
dp x v dx dp p 0, 1 k则 y p, 得 y p k p 2 1 p 化简为 y v2 dy dy dy
p 2 1 (cy ) 2 k 2 p (cy ) k p 2 ,这个方程可以解得：
9.1 题目分析: 这道题目已经给出常微分方程组直接写出程序即可。代码在附件中可以使用。得到的
12 / 25
结果储存在 compA.mat 中。 （1）
由得到的结果可知： 在[0,2]中，乙种群数量一度上升，但之后衰减，t=10 之后，两种种群数量都趋于稳定，
13 / 25
甲达到最大容量，乙趋于 0，属于灭绝。
15 / 25
3、r1 = 1;r2 = 2;n1 = 200;n2 = 100;s1 = 0.5;s2 = 2;
将甲的种群数量上限改为 200 后，依然可以在很短的时间内达到饱和。结果乙仍然是灭 绝。
16 / 25
4、r1 = 1;r2 = 2;n1 = 100;n2 = 200;s1 = 0.5;s2 = 2;
4 / 25
（3） 若流速 v1=0，0.5，1.5，2（m/s） ，结果将如何？
6.1模型分析 对速度进行两个方向的分解。得ห้องสมุดไป่ตู้微分方程组：
v2 x dx dt v1 2 x y2 dy v2 y dt x2 y 2
两个式子相除得到：
dx v 1 dy v2
数学实验
ODE 初值问题
计 03 林智慧 2010011292
Lab2_ODE 初值问题实验报告
实验目的： 1.用matlab软件掌握求微分方程数值解的方法，并对结果作初步分 析； 2.通过实例学习用微分方程模型解决简化的实际问题。
实验内容 题目 3 小型火箭初始重量为 1400kg， 其中包括 1080kg 燃料。 火箭竖直向上发射时燃料燃烧 率为 18kg/s，由此产生 32000N 的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气 阻力正比于速度的平方，比例系数为 0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速 度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
7 / 25
小船 50s 沿直线过河。小船位置解析解为 x
8 / 25
1 1 y y 0 。与实验结果一致。 2 2
V1=0.5m/s 时
9 / 25
与前面求解过程一致，经过 53.33s 小船过河。
V1=1.5m/s 时
可以看出小船需要更多地时间过河。经计算为 114.28s。
10 / 25
3.1 问题分析 火箭质量 m,高度为 h，速度为 v，加速度 a，火箭推力为 Fp,重力加速度 g，空气阻力 f。 首先从题目条件中易得火箭燃料燃烧 60s。阻力 f 0.4v 。火箭质量 m(t ) m0 18t 。
2
dh d dv 32000 0.4 v 2 dt 。故根据牛顿第二定律可以得到 a 加速度 a 9.8 dt dt 1400 18t
其中 x(t)，y(t)分别为甲乙两种群的数量；r1，r2 为它们的固有增长率；n1，n2 为 它们的最大容量。s1 的含义是，对于供养甲的资源而言，单位数量乙（相对 n2）的消 耗为单位数量甲（相对 n1）消耗的 s1 倍，对 s2 可作相应的解释。 该模型无解析解，试用数值解法研究以下问题： （1）设 r1=r2=1，n1=n2=100，s1=0.5，s2=2，初值 x0=y0=10，计算 x(t)， y(t)，画出它们的图形及相图（x，y） ，说明时间 t 充分大以后 x(t)，y(t)的 变化趋势（人们今天看到的已经是自然界长期演变的结局） 。 （2）改变 r1,r2,n1,n2,x0,y0，但 s1,s2 不变（或保持 s1<1,s2>1） ，计算分析所 得结果；若 s1=1.5(>1)，s2=0.7(<1)，再分析结果。由此你能得到什么结 论，请用各参数生态学上的含义做出解释。 （3）试 验 当 s1=0.8(<1) ， s2=0.7(<1) 时 会 有 什 么 结 果 ； 当 s1=1.5(>1) ， s2=1.7(>1)时又会有什么结果。能解释这些结果吗？
function dx = rocket1(t,x) %初始条件,记x(1)=h,x(2)=v Fp = 32000;g = 9.8;k = 0.4;m0 = 1400;u = 18; dx = [x(2);(Fp-0.4*x(2)^2)/(m0-u*t)-g]; end
测试实例可以参见代码。主要就是画图。但这只是分段函数的第一部分。第二部分为燃料 用尽，但仍然上升的阶段.这一阶段,火箭受重力和空气阻力。常微分方程组如下：
25 / 25
3.3实验结果
2 / 25
3 / 25
上面三种颜色代表三种不同的状态。 所有的数据都储存在T1.mat，T2.mat，T3.mat三个matlab数据文件中。第一列为时间， 第二列为高度，第三列为速度，第四列为加速度。得到如下近似结果： 关闭引擎时的速度为267.257m/s，加速度a=0.9170 m / s2 ，高度为12190m。 火箭达到最高点即第二阶段结束时，高度为13115m，加速度为-9.8 m / s2 。
6.2实验过程 照着前面列出的常微分方程组写出函数测试即可。需要注意的是由于这个过程计算时 间过程，故使用解刚性方程的方法来解。得到结果为：
5 / 25
蓝色为横向方向，红色为竖直方向。
6 / 25
由数据可以得到当t=66.66s时，x=0.00697m,y=-1.9451e-6m。这时小船已经到达了目 的地B。 解析解得话,由前面的解析公示带入数据
dh v dt ,初始条件为上一阶段结束状态。 2 dv 9.8 0.4v 1400 1080 dt
最后一个阶段为火箭下降阶段。这一阶段火箭受力为向上的空气阻力和重力。常微分方程 组如下：
dh v dt ,初值条件为上阶段结束状态。 0.4v 2 dv 9.8 1400 1080 dt
列出常微分方程组：
1 / 25
dh v(0) 0 v dt ，初值为 h(0) 0 2 dv Fp 0.4v 9.8 t 60 s dt 1400 18t
3.2 实验过程 列出常微分方程组，就可以使用龙格-库塔方法实现算法。编写函数。具体为:
（2）1、r1 = 2;r2 = 1;n1 = 100;n2 = 100;s1 = 0.5;s2 = 2;
可见甲的数量上升更快，乙灭绝得也更快。
14 / 25
2、r1 = 1;r2 = 2;n1 = 100;n2 = 100;s1 = 0.5;s2 = 2;
乙刚开始上升阶段的最大种群数量有所上升,但不影响最后的结果.
V1=2m/s 时
可以看出小船永远也到不了对岸 B 点。解析解得:
11 / 25
dx v2 x x v1 v2 (1 ) 0 。小船最后静止在(50,0)处。 2 2 2 dt x y x y2
题目 9： 两种群相互竞争模型如下：
x y (t ) r1 x(1 s1 ), x n1 n2 x y y (t ) r2 y (1 s2 ), n1 n2
7、s1=1.5,s2=0.7
21 / 25
情况几乎翻转,乙趋向饱和,甲趋向灭绝. 这说明 r1 和 r2 是甲乙种群的固有增长率， 他们决定了种群在开始一段时间内数量 增长的快慢，固有增长率大的种群增长更快。同时决定了到达稳态的时间快慢，固有增 长率大的种群到达稳态时间短。n1 和 n2 是甲乙两种群的最大容量，决定了某一种群最 终的稳态值。x0 和 y0 是甲乙两种群的初始值，决定了物种曲线的起点，但不影响最终 的稳态。s1 是对于供养甲的资源而言，单位数量乙的消耗为单位数量甲消耗的 s1 倍。 s2 的意义与此相反。若 s1<1<s2，则甲物种会抢夺乙物种的资源，乙物种并不抢夺甲 物种的资源，最终乙没有充足的食物，从而灭绝，甲乙的地位是对称的，反过来灭绝的 是甲，而乙会达到饱和。这些都非常符合生态学常识。说明了实验的正确性。
1 1 1 1 p (cy ) k (cy ) k ,代入前式整理得 x c k y1k c k y1 k ,c是待定系数。代入初值 2 2 2 2 1 1 1 k 1k 1 1 k 1 k y （- ）y 。 ( x, y) (0, d ) ，得 c , x （- ） d 2 d 2 d
附代码列表： T3 Rocket1.m Rocket2.m Rocket3.m
24 / 25
Rocket_test.m T1.mat T2.mat T3.mat T6 Boat.m Boat_inl.m Boat_test.m Boat_test_3.m T9 Comp.m Comp_test.m compA.mat
题目6 一只小船渡过宽为 d 的河流，目标是起点 A 正对着的另一岸 B 点。已知河水流速 v1 与船在静水中的速度 v2 之比为 k。 （1） 建立描述小船航线的数学模型，求其解析解； （2） 设 d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需的时间、任意时刻小船 的位置及航行曲线，作图，并与解析解比较。
17 / 25
将乙的种群数量上限改为 200 后，只是在刚开始很短的时间内数量超过了甲，但最后仍然 是灭绝，甲仍然是饱和。
5、更改初值,让甲改为原来的 10 倍 100
18 / 25
乙都没有上升一直减少，很快灭绝。
6、将乙初值改为 100
19 / 25
乙一直下降最后灭绝。甲一直上升最后饱和。
20 / 25
（4）s1=0.8(<1)，s2=0.7(<1)
22 / 25
S1<1,s2<1。甲乙不会互相抢夺对方资源，从而基本稳定，但明显乙对资源更有竞 争力。所以乙最后上升到更大的稳定值。
s1=1.5(>1)，s2=1.7(>1)
23 / 25
当 s1>1 且 s2>1 时，两物种对彼此的资源都有抢夺性。因为 s1 略小于 s2，所以 甲物种对乙物种资源的掠夺性更大， 导致在一段时间后甲物种数量最终饱和， 乙物种几 乎灭绝。 从上面的实验可以看出,对于两种处于竞争关系的物种来说,s1 和 s2 起着非常重要 的作用。
( x, y) (0, 100) ，得 c=-0.01,k=0.5 故
1 1 x (0.01) 0.5 y 0.5 (0.01) 0.5 y 0.5 2 2
写出函数 boat_inl.m 计算得结果：
通过比较可以看出我们前面用的数值方法精确性还是非常不错的。
(3)改变 v1=0 结果如下：
替换变量
x 2 y 2 x v1 x x ( )2 1 y y v2 y y
dp x v dx dp p 0, 1 k则 y p, 得 y p k p 2 1 p 化简为 y v2 dy dy dy
p 2 1 (cy ) 2 k 2 p (cy ) k p 2 ,这个方程可以解得：
9.1 题目分析: 这道题目已经给出常微分方程组直接写出程序即可。代码在附件中可以使用。得到的
12 / 25
结果储存在 compA.mat 中。 （1）
由得到的结果可知： 在[0,2]中，乙种群数量一度上升，但之后衰减，t=10 之后，两种种群数量都趋于稳定，
13 / 25
甲达到最大容量，乙趋于 0，属于灭绝。
15 / 25
3、r1 = 1;r2 = 2;n1 = 200;n2 = 100;s1 = 0.5;s2 = 2;
将甲的种群数量上限改为 200 后，依然可以在很短的时间内达到饱和。结果乙仍然是灭 绝。
16 / 25
4、r1 = 1;r2 = 2;n1 = 100;n2 = 200;s1 = 0.5;s2 = 2;
4 / 25
（3） 若流速 v1=0，0.5，1.5，2（m/s） ，结果将如何？
6.1模型分析 对速度进行两个方向的分解。得ห้องสมุดไป่ตู้微分方程组：
v2 x dx dt v1 2 x y2 dy v2 y dt x2 y 2
两个式子相除得到：
dx v 1 dy v2
数学实验
ODE 初值问题
计 03 林智慧 2010011292
Lab2_ODE 初值问题实验报告
实验目的： 1.用matlab软件掌握求微分方程数值解的方法，并对结果作初步分 析； 2.通过实例学习用微分方程模型解决简化的实际问题。
实验内容 题目 3 小型火箭初始重量为 1400kg， 其中包括 1080kg 燃料。 火箭竖直向上发射时燃料燃烧 率为 18kg/s，由此产生 32000N 的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气 阻力正比于速度的平方，比例系数为 0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速 度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
7 / 25
小船 50s 沿直线过河。小船位置解析解为 x
8 / 25
1 1 y y 0 。与实验结果一致。 2 2
V1=0.5m/s 时
9 / 25
与前面求解过程一致，经过 53.33s 小船过河。
V1=1.5m/s 时
可以看出小船需要更多地时间过河。经计算为 114.28s。
10 / 25
3.1 问题分析 火箭质量 m,高度为 h，速度为 v，加速度 a，火箭推力为 Fp,重力加速度 g，空气阻力 f。 首先从题目条件中易得火箭燃料燃烧 60s。阻力 f 0.4v 。火箭质量 m(t ) m0 18t 。
2
dh d dv 32000 0.4 v 2 dt 。故根据牛顿第二定律可以得到 a 加速度 a 9.8 dt dt 1400 18t
其中 x(t)，y(t)分别为甲乙两种群的数量；r1，r2 为它们的固有增长率；n1，n2 为 它们的最大容量。s1 的含义是，对于供养甲的资源而言，单位数量乙（相对 n2）的消 耗为单位数量甲（相对 n1）消耗的 s1 倍，对 s2 可作相应的解释。 该模型无解析解，试用数值解法研究以下问题： （1）设 r1=r2=1，n1=n2=100，s1=0.5，s2=2，初值 x0=y0=10，计算 x(t)， y(t)，画出它们的图形及相图（x，y） ，说明时间 t 充分大以后 x(t)，y(t)的 变化趋势（人们今天看到的已经是自然界长期演变的结局） 。 （2）改变 r1,r2,n1,n2,x0,y0，但 s1,s2 不变（或保持 s1<1,s2>1） ，计算分析所 得结果；若 s1=1.5(>1)，s2=0.7(<1)，再分析结果。由此你能得到什么结 论，请用各参数生态学上的含义做出解释。 （3）试 验 当 s1=0.8(<1) ， s2=0.7(<1) 时 会 有 什 么 结 果 ； 当 s1=1.5(>1) ， s2=1.7(>1)时又会有什么结果。能解释这些结果吗？
function dx = rocket1(t,x) %初始条件,记x(1)=h,x(2)=v Fp = 32000;g = 9.8;k = 0.4;m0 = 1400;u = 18; dx = [x(2);(Fp-0.4*x(2)^2)/(m0-u*t)-g]; end
测试实例可以参见代码。主要就是画图。但这只是分段函数的第一部分。第二部分为燃料 用尽，但仍然上升的阶段.这一阶段,火箭受重力和空气阻力。常微分方程组如下：
25 / 25