实验4常微分方程数值解.docx
《实验4常微分方程数值解.docx》由会员分享，可在线阅读，更多相关《实验4常微分方程数值解.docx（33页珍藏版）》请在冰豆网上搜索。
实验4常微分方程数值解
实验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文件
functiondx=Rocket（t,x）
dx=[x
（2）;（32000-0.4*x
（2）^2）/（1400-18*t）-9.8];%以向量形式表示微分方程
end
ts=0:
60%终点时间为60s，步长定义为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]
由MATLAB计算得到从开始上升到关闭引擎瞬间的情况如下表：
t
h
v
a
0
0
0
13.057
1
6.5737
13.189
13.305
2
26.444
26.577
13.453
3
59.762
40.062
13.497
4
106.57
53.535
13.433
5
166.79
66.89
13.261
6
240.27
80.021
12.985
7
326.72
92.829
12.612
8
425.79
105.22
12.152
9
536.99
117.11
11.617
10
659.8
128.43
11.021
11
793.63
139.14
10.38
12
937.85
149.18
9.7083
13
1091.8
158.55
9.0209
14
1254.7
167.23
8.3309
15
1425.9
175.22
7.6502
16
1604.8
182.55
6.9901
17
1790.8
189.22
6.3593
18
1983.1
195.27
5.7646
19
2181.2
200.75
5.2095
20
2384.5
205.7
4.6946
21
2592.4
210.18
4.222
22
2804.5
214.19
3.7943
23
3020.6
217.79
3.412
24
3240.1
221.01
3.073
25
3462.7
223.92
2.7726
26
3687.9
226.56
2.5044
27
3915.6
228.97
2.2677
28
4145.6
231.14
2.0633
29
4377.8
233.11
1.8898
30
4611.9
234.91
1.7433
31
4847.7
236.57
1.6178
32
5085
238.14
1.5062
33
5323.8
239.61
1.4095
34
5564.1
240.99
1.3293
35
5805.8
242.28
1.265
36
6048.7
243.5
1.2139
37
6292.9
244.68
1.1708
38
6538.1
245.83
1.1303
39
6784.5
246.96
1.0947
40
7032
248.05
1.0663
41
7280.5
249.1
1.0456
42
7530.2
250.12
1.0308
43
7780.9
251.14
1.0178
44
8032.5
252.15
1.0024
45
8285.1
253.16
0.98757
46
8538.8
254.15
0.97632
47
8793.4
255.12
0.96964
48
9049
256.07
0.96629
49
9305.6
257.03
0.96239
50
9563.1
257.99
0.95272
51
9821.5
258.95
0.94118
52
10081
259.9
0.93373
53
10341
260.83
0.93278
54
10603
261.75
0.93633
55
10865
262.67
0.93695
56
11128
263.61
0.92585
57
11392
264.54
0.91379
58
11657
265.46
0.91063
59
11923
266.35
0.91612
60
12190
267.26
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文件
functiondy=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]
forn=1:
2000%此处利用一个循环语句找出时间终点
T=80-0.01*n;%估计出时间不会晚于80
tss=60:
0.02:
T;
y0=[x（61,1）,x（61,2）];
option=odeset（'reltol',1e-3,'abstol',1e-6）;
[t2,y]=ode45（@Rocket2,tss,y0,option）;
[t2,y];
ify（:
2）>=0
break
end
end
plot（t,x（:
1）,'b',t2,y（:
1）,'r'）,grid,
title（'图1.高度-时间'）
xlabel（'t/s'）
ylabel（'h/m'）
pause
plot（t,x（:
2）,'b',t2,y（:
2）,'r'）,grid,
title（'图2.速度-时间'）
xlabel（'t/s'）
ylabel（'v/（m/s）'）
pause
a=（32000-0.4*x（:
2）.^2）./（1400-18*t）-9.8;
a2=-9.8-0.4*y（:
2）.^2/320;
plot（t,a,'b',t2,a2,'r'）,grid,
title（'图3.加速度-时间'）
xlabel（'t/s'）
ylabel（'a/（m/s^2）'）
由MATLAB数据计算可知，当t=71.31s时，火箭上升到最大高度h=13115m，此时火箭的速度v=0.019874，几乎可认为已经停止，加速度a=-9.8m/s2。
【小结】
由以上三图的曲线可以看出，起初火箭迅速加速，而且速度值不大时，燃料动力的贡献要超过空气阻力，而质量的减少又会使加速度变大，故而两个因素共同作用使加速度较为稳定。
此后速度增加，阻力变大，加速度渐渐变为0.当燃料耗尽时，推动力消失，加速度突变负，而此时速度由增大转为减少。
此时速度迅速下降，但是由于它还是正值，故高度上升。
速度为0时高度最大，无阻力作用，加速度等于重力加速度。
这就是全过程中高度、速度、加速度的定性分析过程。
题目6
一只小船渡过宽为d的河流，目标是起点A正对着的另一岸B点。
已知河水流速v1与船在静水中的速度v2之比为k。
（1）建立描述小船航线的数学模型，求其解析解；
（2）设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需的时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。
（3）若流速v1=0，0.5，1.5，2（m/s），结果将如何？
6.1建立数学模型及分析
此种模型的前提是船并不事先知道水速，否则只要调整一个合适的角度，直接沿直线通过即可。
现小船不知道水速，则它的策略应为始终使船头瞄准B点。
对速度进行xy两个方向的分解，可列出常微分方程如下：
其初值条件为x（0）=0,y（0）=-d.
6.2常微分方程解析解
两式相除，得
变量代换，令
则有dxd
代入有
分离变量积分可得
代入初值，可解出
整理得解析解最终表达式为
6.3MATLAB数值解
先建立M文件函数
functiondx=Guohe（t,x,v1,v2）
s=（x
（1）^2+x
（2）^2）^0.5;
dx=[v1-x
（1）/s*v2;-x
（2）/s*v2];
end
ts=0:
0.01:
80;
x0=[0,-100];
option=odeset（'reltol',1e-6,'abstol',1e-9）;
[t,x]=ode23s（@Guohe,ts,x0,option,1,2）;
plot（t,x）,grid
title（'图1.分位移-时间'）
xlabel（'t/s'）
ylabel（'x,y/m'）
gtext（'x（t）'）
gtext（'y（t）'）
pause
plot（x（:
1）,x（:
2））,grid,
title（'图2.过河轨迹'）
xlabel（'x/m'）
ylabel（'y/m'）
[t,x（:
1）,x（:
2）]
节选部分数据如下:
t
x
y
0
0
0
0.01
0.009999
0.009999
0.02
0.019996
0.019996
0.03
0.029991
0.029991
0.04
0.039984
0.039984
0.05
0.049975
0.049975
0.06
0.059964
0.059964
0.07
0.069951
0.069951
0.08
0.079936
0.079936
0.09
0.089919
0.089919
0.1
0.0999
0.0999
0.11
0.10988
0.10988
0.12
0.11986
0.11986
0.13
0.12983
0.12983
0.14
0.1398
0.1398
0.15
0.14977
0.14977
0.16
0.15974
0.15974
0.17
0.16971
0.16971
0.18
0.17968
0.17968
0.19
0.18964
0.18964
0.2
0.1996
0.1996
……
……
……
66.5
0.16685
-0.0011133
66.51
0.15685
-0.00098384
66.52
0.14685
-0.00086239
66.53
0.13685
-0.000748
66.54
0.12685
-0.00064349
66.55
0.11685
-0.00054603
66.56
0.10685
-0.00045657
66.57
0.096853
-0.00037511
66.58
0.086853
-0.00030164
66.59
0.076853
-0.00023618
66.6
0.066853
-0.00017871
66.61
0.056853
-0.00012924
66.62
0.046853
-8.7773e-005
66.63
0.036853
-5.4301e-005
66.64
0.026853
-2.8827e-005
66.65
0.016853
-1.1351e-005
66.66
0.0068532
-1.8739e-006
66.67
6.7389e-011
0
【小结】
小船起初速度较快，相当于在顺水被冲下去。
后来快到达终点时较慢，相当于逆水而行。
6.4解析解的图像
先建立函数
functionx=Guohe2（y,k）
x=0.5*（-0.01）^（-k）*y.^（1-k）-0.5*（-0.01）^k*y.^（1+k）;
end
y=0:
-0.01:
-100;
x=Guohe2（y,0.5）;
plot（x,y），grid;
title（'图4.过河轨迹，解析解'）
xlabel（'x/m'）
ylabel（'y/m'）
可以看出，由数值解法得到的数据绘成的图与解析解的结果几乎一样。
6.5改变流速v1=0
ts=0:
0.01:
80;
x0=[0,-100];
option=odeset（'reltol',1e-6,'abstol',1e-9）;
[t,x]=ode23s（@Guohe,ts,x0,option,0,2）;
plot（t,x）,grid
title（'图1.分位移-时间'）
xlabel（'t/s'）
ylabel（'x,y/m'）
gtext（'x（t）'）
gtext（'y（t）'）
pause
plot（x（:
1）,x（:
2））,grid,
title（'图2.过河轨迹'）
xlabel（'x/m'）
ylabel（'y/m'）
6.6改变流速v1=0.5m/s
程序基本相同，在此略去，只给出图像。
6.7改变流速v1=1.5m/s
6.8改变流速v1=2m/s
题目9
两种群相互竞争模型如下：
其中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）时又会有什么结果。
能解释这些结果吗？
9.1模型分析
此题的模型很类似教材中的“弱肉强食”模型，只是多考虑了竞争资源的问题。
微分方程描述了竞争关系和捕食关系。
通过对微分方程组进行数值求解，即得出竞争关系的直观表示。
9.2程序代码
设x
（1）=x,x
（2）=y
编写M文件如下:
functiondx=jing（t,x,r1,r2,n1,n2,s1,s2）
dx=[r1*x
（1）*（1-x
（1）/n1-s1*x
（2）/n2）;r2*x
（2）*（1-s2*x
（1）/n1-x
（2）/n2）];
end
ts=0:
0.1:
15;
x0=[10,10];
option=odeset（'reltol',1e-6,'abstol',1e-9）;
[t,x]=ode45（@jing,ts,x0,option,1,1,100,100,0.5,2）;
[t,x]
plot（t,x）,grid,
title（'图1.时间-种群数量'）
xlabel（'t'）
ylabel（'x,y'）
gtext（'x（t）'）
gtext（'y（t）'）
pause
plot（x（:
1）,x（:
2））,grid
title（'图2.相图'）
xlabel（'x'）
ylabel（'y'）
【小结】
如图，可以清晰地看出，t充分大之后，y起初小有增长，随后衰减，直至灭亡，而x则一直增长。
二者数量在大约10年时已经趋于稳定，此时x的数量约稳定在100，而y已经彻底灭亡。
9.3改变参数
1.改变r1与r2,即改变固有增长率
（1）.r1=0.4,r2=0.4
jing函数中设计了各参数的接口，只要改变输入即可，程序代码省略。
数据分析：
由图中数据可以看出，改变二者增长率，最终结果并没有变化。
仍然是甲达到100，而乙彻底灭绝。
与之前变化不同的是，达到稳定需要更长的时间，约25年。
可以看出，增长率减少，种群数量变化的速度会减慢。
（2）.r1=5,r2=0.4
数据分析：
最终结果仍然没有改变！
但不同的是，由于此时甲物种的固有增长率远高于乙，即竞争优势十分巨大，因此乙物种在初期甚至几乎没有增长就迅速衰减到0.二者约10年时间就达到了平衡。
（2）.r1=0.4,r2=5
数据分析：
至此我们基本可以下结论，改变增长率的值不会改变它们最终的演化结果。
图中数据表明，甲物种增长率小，增长速度很慢，相当于减缓了乙物种的灭绝进程。
乙初期数量迅速增长，但是后来随着甲的增多渐渐灭绝。
不同的是，乙大约在11年时灭绝，但是甲的数量持续增长到约22年才稳定。
2.改变n1与n2,即改变最大容量
（1）.n1=1000,n2=100
数据分析：
初期甲物种的数量占最大容量的比例较小（x/n1）即威胁小，所以乙物种增长迅速但最终仍然灭绝。
物种容量的改变并不能影响最终谁会灭绝。
（2）.n1=100,n2=1000
数据分析：
在此数据下，乙物种完全没有机会增长到容量值就早已灭亡.
3.改变x0,y0,即改变初值
（1）.x0=10,y0=100.
数据分析：
乙物种的初始数量大使其灭绝时间稍稍延后，但它灭绝的趋势不变。
（2）.x0=100,y010
数据分析：
乙就这么毫无悬念地死掉了……才用了5年时间啊
【小结】由以上分析知，不论怎么改变以上6个参数，都不会改变甲乙最终的宿命。
甲总会优胜，达到环境容量值，而乙永远逃脱不了灭绝的命运。
9.4改变s1,s2
1.s1=1.5,s2=0.7
数据分析：
情况正好相反，最后是甲物种灭绝，而乙物种达到环境容量值。
【小结】
由s1,s2的生态学意义，当某个s1和s2一个大于1，另一个小于1时，其中一物种将过量消耗与其竞争的物种的生存资源，从而导致另一物种灭绝。
2.s1=0.8,s2=0.7
数据分析：
数值计算结果最后稳定在x=45.47，y=68.168，说明此时两种物种竞争会达到一个相对平衡的值，并以此值稳定共存。
3.s1=1.5,s2=1.7
数据分析：
此时s1、s2都大于1，都在过量地消耗着对方的资源，但是更大的一方会被消耗更多的资源，最终导致灭亡。
【本题总结】
s1,s2小于1时，互相消耗程度较轻，因此可以达到平衡共存的状态。
但两者都无法达到容量值，因为互相在制约着对方的生长。
当其中之一大于1时，它就会因为被消耗资源而灭亡；当s1,s2都大于1时，两物种竞争激烈，最后s1,s2中更小者在竞争中获胜，另一物种灭绝。
正所谓物竞天择，适者生存。
自然中能供养物种的资源是十分有限的，谁能更好地适应环境，谁就能在竞争中取得最终的胜利。
这对我们人类也是一个警示。
【作业总结】
如果说上次作业让我简单地了解了数学建模，那么这次作业是让我更加深刻地理解了数学建模的重要意义。
当我编完程序，图像显示了种群数量的变化时，我震惊了！
一个简单的计算就可以反映出生活
升级会员 