常微分方程数值解Word文档下载推荐.docx
《常微分方程数值解Word文档下载推荐.docx》由会员分享，可在线阅读，更多相关《常微分方程数值解Word文档下载推荐.docx（35页珍藏版）》请在冰豆网上搜索。
时间t/s
高度h/m
速度v/（m/s）
加速度a/（m/s2）
0.000
13.05
1.000
6.569
13.18
13.29
2.000
26.42
26.56
13.44
3.000
59.72
40.03
13.49
4.000
106.5
53.50
13.42
5.000
166.7
66.84
13.25
6.000
240.1
79.97
12.98
7.000
326.5
92.77
12.61
8.000
425.5
105.2
12.15
9.000
536.6
117.0
11.61
10.00
659.3
128.4
11.02
11.00
792.9
139.1
10.38
12.00
937.0
149.1
9.704
13.00
1090.8
158.5
9.017
14.00
1253.4
167.2
8.325
15.00
1424.4
175.2
7.642
16.00
1603.3
182.5
6.982
17.00
1788.9
189.2
6.347
18.00
1980.7
195.4
5.742
19.00
2178.3
200.9
5.176
20.00
2381.4
205.9
4.659
21.00
2589.2
210.3
4.184
22.00
2800.7
214.5
3.735
23.00
3015.6
218.3
3.318
24.00
3234.2
221.7
2.945
25.00
3456.3
224.7
2.624
26.00
3681.9
227.2
2.367
27.00
3909.9
229.5
2.141
28.00
4139.8
231.7
1.928
29.00
4371.6
233.8
1.734
30.00
4605.3
235.7
1.567
31.00
4841.0
237.3
1.435
32.00
5078.9
238.8
1.346
33.00
5318.2
1.276
34.00
5558.7
241.4
1.207
35.00
5800.4
242.7
1.144
36.00
6043.3
243.9
1.090
37.00
6287.6
245.1
1.051
38.00
6533.1
246.2
1.031
39.00
6779.7
247.2
1.019
40.00
7027.4
248.2
1.005
41.00
7276.0
249.2
0.991
42.00
7525.7
250.3
0.978
43.00
7776.4
251.3
0.969
44.00
8028.1
252.2
0.965
45.00
8280.8
253.2
0.963
46.00
8534.5
254.2
0.958
47.00
8789.1
255.1
0.953
48.00
9044.7
256.1
0.948
49.00
9301.2
257.0
0.943
50.00
9558.8
258.0
0.941
51.00
9817.2
258.9
0.937
52.00
10076.6
259.9
0.930
53.00
10336.9
260.8
0.922
54.00
10598.2
261.8
0.914
55.00
10860.4
262.7
0.909
56.00
11123.5
263.6
57.00
11387.6
264.5
0.912
58.00
11652.6
265.4
59.00
11918.5
266.3
60.00
12185.3
267.3
0.911
从表中可以看出，在引擎关闭即燃料燃尽的瞬间，火箭的高度为h=12185.3m，速度v=267.3m/s，加速度a=0.911m/s2。
3.3燃料耗尽后物理模型分析
当引擎关闭后，火箭继续上升，但由于受重力加速度的影响，速度减小，直至速度减为0时达到最高点。
（1）燃料耗尽后火箭质量为320kg。
（2）火箭所受合力F=-f-mg;
（3）火箭加速度为a=F/m=-f/m-g=
;
于是可以列出以下微分方程：
初值条件为h=12185.3m，v=267.3m/s
3.4matlab程序及图像数据
记y
（1）=h,y
（2）=v,y=（y
（1）,y
（2））T。
%-----------------题目3燃料燃尽后函数rocketM文件----------------------
functiondx=rocket（t,y）;
dx=dy=[y
（2）;
-9.81-0.4*y
（2）^2/320]%以向量形式表示上述方程组
先运行以下程序，以求得v=0时的时间t：
%---------------题目3燃料燃尽后试验rocket1源程序-------------------
%先算出引擎关闭前的速度和高度
y0=[x（61,1）,x（61,2）];
%提取初值
fori=1:
1500;
T=75-0.01*i;
%可以先估算出火箭达到最高点的时间小于75s
tss=60:
0.01:
T;
[t2,y]=ode23（@rocket,tss,y0）;
ify（:
2）>
=0;
%速度小于0时循环终止
break
end
1）,'
b'
t2,y（:
r'
）,grid;
%绘制全部的高度-时间曲线
图4.高度-时间'
2）,'
%绘制全部的速度-时间曲线
图5.速度-时间'
%计算燃料燃尽后的加速度
a2=-9.81-0.4*y（:
2）.^2/320;
plot（t,a,'
t2,a2,'
%绘制全部的加速度-时间曲线
图6.加速度-时间'
v/（m/s2）'
[t2,y（:
1）,y（:
2）,a2];
以下为高度、速度、加速度随时间变化的图形：
由MATLAB数据计算可知，当t=71.30s时，火箭上升到最大高度h=13110m，此时火箭的速度v=0.00309m/s，几乎可认为已经停止，加速度a=-9.81m/s2。
3.5小结
全过程中高度、速度、加速度的定性分析：
由以上三图的曲线可以看出，起初火箭迅速加速，而且速度值不大时，燃料动力的贡献要超过空气阻力，而质量的减少又会使加速度变大，故而两个因素共同作用使加速度较为稳定。
此后速度增加，阻力变大，加速度渐渐变为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常微分方程解析解
两式相除，得
变量代换，令
则有
代入有
分离变量积分可得
代入初值，可解出
整理得解析解最终表达式为
6.3MATLAB数值解
记x
（1）=x,x
（2）=y,x=（x
（1）,x
（2））T。
%-------------------------题目6函数guohe.m源文件----------------------
functiondx=guohe（t,x）
v1=1;
v2=2;
s=（x
（1）^2+x
（2）^2）^0.5;
dx=[v1-x
（1）*v2/s;
-x
（2）*v2/s];
%----------------------题目6求数值解guohev.m源文件--------------------
x0=[0,-100];
0.1:
70;
%经过试验设定终值为70
[t,x]=ode15s（@guohe,ts,x0）;
%求解方程组
plot（t,x）,grid;
%画图同时做出x、y方向的位移
图1.分位移-时间'
）
x,y/m'
gtext（'
x（t）'
y（t）'
plot（x（:
%画出航行曲线
航行曲线'
x'
y'
2）]%输出数据表
得到如下图像：
由图7可以看出，在前30s中，y（t）有较好的的线性，因此，在这个时间段内，小船沿y轴方向的为匀速。
从图8可以看出，船合速度的方向为轨迹的切线方向。
在刚开始阶段，可以认为传顺水而行，小船在刚开始的速度方向的确变化不大，与上图相吻合。
小船在达到x最大值时开始转向，着这一时刻附近的合速度最小。
下面摘取部分数据：
t/s
x/m
y/m
-100.0000
0.1000
0.0999
-99.8000
0.2000
0.1996
-99.6000
0.3000
0.2991
-99.4000
0.4000
0.3984
-99.2000
0.5000
0.4975
-99.0000
0.6000
0.5964
-98.8000
0.7000
0.6951
-98.6000
0.8000
0.7936
-98.4000
0.9000
0.8919
-98.2000
1.0000
0.9900
-98.0000
1.1000
1.0879
-97.8000
1.2000
1.1856
-97.6000
1.3000
1.2831
-97.4000
1.4000
1.3804
-97.2000
1.5000
1.4775
-97.0001
1.6000
1.5744
-96.8001
1.7000
1.6710
-96.6001
1.8000
1.7675
-96.4001
……
65.1000
1.5882
-0.0987
65.2000
1.4886
-0.0866
65.3000
1.3889
-0.0754
65.4000
1.2892
-0.0649
65.5000
1.1894
-0.0552
65.6000
1.0896
-0.0463
65.7000
0.9898
-0.0382
65.8000
0.8899
-0.0308
65.9000
0.7900
-0.0243
66.0000
0.6901
-0.0185
66.1000
0.5902
-0.0135
66.2000
0.4902
-0.0093
66.3000
0.3902
-0.0059
66.4000
0.2903
-0.0033
66.5000
0.1903
-0.0014
66.6000
0.0903
-0.0003
66.7000
0.0000
66.8000
6.4解析解与数值解的比较
下面是matlab源文件程序：
%--------------------题目6guohebijiao.m源程序-------------------------
y0=0:
-1:
-100;
x0=50.*（（y0./（-100））.^0.5-（y0./（-100））.^1.5）;
x0,y0,'
legend（'
数值解'
'
解析解'
航行曲线比较'
由图9可以看出，数值解得到的曲线与解析解的到的曲线几乎完全重合，说明数值解得到的结果还是很准确的。
6.5改变流速
（1）v1=0m/s时,更改程序中的数值，可得如下分位移-时间图像及航行曲线：
由图像可知，v1=0m/s时，水的流速为0，即小船的合速度等于v2，因此小船沿y轴方向匀速前进。
此时小船的位置的解析解为x=0.5*y-0.5*y=0，所以小船过河轨迹为沿y轴方向的一条直线。
小船沿y方向的速度
，故小船匀速过河。
过河所用时间为50s。
（2）v1=0.5s时，
由数据表可知，t=53.3s时，y=-0.0131,可近似认为已到达岸边。
将V1=0.5m/s与v1=1m/s时相比，河水流速减小后，船在x方向上的分速度减小了，因此小船在x方向上行驶的最大值减少。
这与途中x最大值减小了一半左右相符。
与v1=0时比较，小船过河花的时间并没有增加很多，但路程长了很多，这说明小船的平均合速度大于v2。
这是因为实际速度为河水与船水速度的矢量叠加结果。
（3）v1=1.5m/s时，
由数据表可知，t=114.3s时，x=0.0007m,y=0.000m，可认为小船到达对岸。
从图中可以看出，小船在[30s，40s]时间范围内达到x方向的最大值，然后逆流驶向B点。
所以v1越小，船回到B的时间就越短，这一点与常识相符。
将解析式
求导并令其为零，可以求出x在
处，取得最大值，所以k在[0,1]区间内越大，x取得最大值的点就越靠近x轴。
（4）v1=2m/s时
利用程序可以求出，小船最终在下游50m处靠岸，不能到达正对岸。
也就是说，无论多长时间小船都无法到达正对岸。
由解析解也可以证明这一点。
，且不恒等于0，所以，对于任意的T＞0，都有x=
，所以小船不可能到达正对岸。
将k=1代入
得到y=0，所以小船将在y=0时到达x方向的最大值，及当小船到达对岸时，向下游漂的距离最大。
由解析式可得，
，为一条开口向左的抛物线，其轨迹与由matlab作出的轨迹相符。
所以小船无法到达B，最终静止在（50,0）处。
6.6小结
从上述结果可以看出，当静水船速小于等于水速时，不可能通过调整船头始终指向对岸的方法到达对岸。
只有当船速大于水速时，才可能用以上方法到达对岸，因为此时静水船速除了可以提供向对岸的分速度外，还能提供逆水流方向的分速度以防船向下漂移。
【题目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。
首先编写如下函数：
%---------------------题目9函数jingzh.m源程序-----------------------
functiondx=jingzh（t,x）
r1=1;
r2=1;
n1=100;
n2=100;
s1=0.5;
s2=2;
dx=[r1*x
（1）*（1-x
（1）/n1-s1*x
（2）/n2）;
r2*x
（2）*（1-s2*x
（1）/n1-x
（2）/n2）]
然后运行以下程序
%---------------------题目9jingzhv.m源程序-------------------------
0.2:
15;
x0=[10,10];
[t,x]=ode23s（@jingzh,ts,x0）;
[t,x];
t'
x,y'
得到以下图像及数据：
t
x
y
10.0000
12.8352
12.1993
16.2682
14.5058
20.3184
16.7489
24.9673
18.7198
30.1504
20.2070
35.7606
21.0409
2.1000
41.6580
21.1322
2.4000
47.6833
20.4906
2.7000
53.6716
19.2253
3.0000
59.4906
17.4859
3.3000
65.0146
15.4468
3.6000
70.1313
13.3133
3.9000
74.7922
11.2052
4.2000
78.9486
9.2409
4.5000
82.5937
7.4878
4.8000
85.7407
5.9731
5.1000
88.4168
4.7016
5.4000
90.6623
3.6593
5.7000
92.5242
2.8212
6.0000
94.0516
2.1579
6.3000
95.
升级会员 