常微分方程的数值解教程文件.doc

常微分方程的数值解

精品文档
实验4 常微分方程的数值解

【实验目的】
1．掌握用MATLAB软件求微分方程初值问题数值解的方法；
2．通过实例用微分方程模型解决简化的实际问题；
3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。

【实验内容】
题3
小型火箭初始重量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。

模型及其求解
火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=9.8m/s2。
在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。因此有如下二式：
a=dv/dt=(F-mg-0.4v2)/m=(32000-0.4v2)/(1400-18t)-9.8
dh/dt=v
又知初始时刻t=0，v=0，h=0。记x(1)=h，x(2)=v，根据MATLAB可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。程序如下：
function [ dx ] = rocket( t,x )
a=[(32000-0.4*x(2)^2)/(1400-18*t)]-9.8;
dx=[x(2);a];
end

ts=0:1:60;
x0=[0,0];
[t,x]=ode45(@rocket,ts,x0);
h=x(:,1);
v=x(:,2);
a=[(32000-0.4*(v.^2))./(1400-18*t)]-9.8;
[t,h,v,a];
数据如下：
t
h
v
a
0
0
0
13.06
1.00
6.57
13.19
13.30
2.00
26.44
26.58
13.45
3.00
59.76
40.06
13.50
4.00
106.57
53.54
13.43
5.00
166.79
66.89
13.26
6.00
240.27
80.02
12.99
7.00
326.72
92.83
12.61
8.00
425.79
105.22
12.15
9.00
536.99
117.11
11.62
10.00
659.80
128.43
11.02
11.00
793.63
139.14
10.38
12.00
937.85
149.18
9.71
13.00
1091.79
158.55
9.02
14.00
1254.71
167.23
8.33
15.00
1425.93
175.22
7.65
16.00
1604.83
182.55
6.99
17.00
1790.78
189.22
6.36
18.00
1983.13
195.27
5.76
19.00
2181.24
200.75
5.21
20.00
2384.47
205.70
4.69
21.00
2592.36
210.18
4.22
22.00
2804.52
214.19
3.79
23.00
3020.56
217.79
3.41
24.00
3240.08
221.01
3.07
25.00
3462.65
223.92
2.77
26.00
3687.88
226.56
2.50
27.00
3915.58
228.97
2.27
28.00
4145.60
231.14
2.06
29.00
4377.76
233.11
1.89
30.00
4611.86
234.91
1.74
31.00
4847.68
236.57
1.62
32.00
5085.02
238.14
1.51
33.00
5323.85
239.61
1.41
34.00
5564.11
240.99
1.33
35.00
5805.77
242.28
1.27
36.00
6048.72
243.50
1.21
37.00
6292.87
244.68
1.17
38.00
6538.11
245.83
1.13
39.00
6784.48
246.96
1.09
40.00
7031.96
248.05
1.07
41.00
7280.54
249.10
1.05
42.00
7530.19
250.12
1.03
43.00
7780.85
251.14
1.02
44.00
8032.49
252.15
1.00
45.00
8285.12
253.16
0.99
46.00
8538.75
254.15
0.98
47.00
8793.39
255.12
0.97
48.00
9049.01
256.07
0.97
49.00
9305.58
257.03
0.96
50.00
9563.08
257.99
0.95
51.00
9821.52
258.95
0.94
52.00
10080.93
259.90
0.93
53.00
10341.30
260.83
0.93
54.00
10602.62
261.75
0.94
55.00
10864.86
262.67
0.94
56.00
11127.98
263.61
0.93
57.00
11392.04
264.54
0.91
58.00
11657.03
265.46
0.91
59.00
11922.96
266.35
0.92
60.00
12189.78
267.26
0.92

因此，在引擎关闭的瞬间，火箭的速度为267.26m/s，高度为12189.78m，加速度为0.92m/s2。
（2）在第二个阶段，火箭的重量保持不变，没有向上的推力，只收到重力和空气阻力。因此有如下关系式：
a=dv/dt=(-mg-0.4v2)/m=-0.4v2/320-9.8
dh/dt=v
假设在80秒之前达到最高点，以60秒时的速度、高度、加速度为初始值进行计算，程序如下：

function [ dx ] = rocket2( t,x )
dx=[x(2);-0.4*x(2)^2/320-9.8];
end

ts2=60:1:80;
x1=[12189.78,267.26];
[t2,x2]=ode45(@rocket2,ts2,x1);
h2=x2(:,1);
v2=x2(:,2);
a2=-0.4*v2.^2./320-9.8;
[t2,h2,v2,a2];
数据如下：
t2 h2 v2 a2
60.00 12189.78 267.26 -99.08
61.00 12416.32 192.70 -56.22
62.00 12584.73 147.43 -36.97
63.00 12715.67 116.11 -26.65
64.00 12819.59 92.79 -20.56
65.00 12902.81 74.29 -16.70
66.00 12969.22 58.96 -14.15
67.00 13021.42 45.73 -12.41
68.00 13061.17 33.95 -11.24
69.00 13089.63 23.12 -10.47
70.00 13107.61 12.91 -10.01
71.00 13115.55 3.02 -9.81
72.00 13113.66 -6.80 -9.86
73.00 13101.90 -16.78 -10.15
74.00 13079.96 -27.19 -10.72
75.00 13047.26 -38.34 -11.64
76.00 13002.90 -50.62 -13.00
77.00 12945.47 -64.56 -15.01
78.00 12872.96 -80.96 -17.99
79.00 12782.31 -101.08 -22.57
80.00 12668.88 -127.03 -29.97

可以看到：在第60秒瞬间，加速度发生了突变，从0.92m/s2突变为-99.08m/s2；而在第71秒至第72秒之间，速度从正变为负，即速度反向，说明在第71秒中的某个时刻速度为零，火箭达到了最高点。因此需要对这个时间段进行分析，并且找到速度减小到0的时刻和此时的高度。
以0.1为步长，在71s到72s中重新求解微分方程的数值解。
71.00 13115.16 2.98 -9.81
71.10 13115.41 2.00 -9.81
71.20 13115.56 1.02 -9.80
71.30 13115.61 0.04 -9.80
71.40 13115.57 -0.94 -9.80
71.50 13115.42 -1.92 -9.80
可见在t=71.3时，速度为0.04，可视为速度为零点，此时最大高度为13115.61，加速度-9.80。
综合（1），（2），可以绘出高度，速度和加速度随时间的变化曲线。
plot(t,h,t2,h2),xlabel('t/s'),ylabel('h/m'),title('高度随时间变化曲线');
plot(t,v,t2,v2),xlabel('t/s'),ylabel('v/(m/s)'),title('速度随时间变化曲线');
aa=[a',a2']'；tt=[t',t2']'；
plot(tt,aa),xlabel('t/s'),ylabel('a/(m/s2)'),title('加速度随时间的变化曲线');

题6
一只小船渡过宽为d的河流，目标是起点A正对着的另一岸B点。已知河水流速v1与船在静水中的速度v2之比为k。
（1）建立描述小船航线的数学模型，求其解析解；
（2）设d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较；
（3）若流速v1=0，0.5，1.5，2m/s，结果将如何。

模型及其求解
（1）.假设在航行过程中，人们不知道水流的速度，小船的方向始终指向目标B，因此若以B为原点，我们可以得到如下方程组：

解初值为（x, y）=( 0, -d) 的常微分方程组，得到解析解为：

其中 c=–1/d=–0.01,故有
事实上，若用matlab中计算微分方程的语句：
[x,y]=dsolve('Dx=v1-v2*x/sqrt(x^2+y^2)','Dy=-v2*y/sqrt(x^2+y^2)','x(0)=0','y(0)=-d','t');
则显示“Warning: Explicit solution could not be found.”即无法得到x,y关于t的分立解。
（2）d=100m,v1=1m/s,v2=2m/s。求数值解时，由解析解可以看出，此为刚性方程。为避免用ode45s求解时间过长。采用ode23s进行求解。假定100s可以到达对岸。
ts=0:0.1:100;
x0=[0,-100];
[t,x]=ode23s(@boat,ts,x0);
[t,x];
plot(t,x);gtext('x');gtext('y');xlabel('时间t/s');ylabel('距离/m');title('x,y与时间t的关系');
可以得到数据如下（部分）：
t/s
x/m
y/m
10.00
8.93
-80.04
20.00
15.41
-60.36
30.00
18.87
-41.49
40.00
18.71
-24.31
50.00
14.48
-10.34
50.10
14.42
-10.23
66.50
0.19
-0.00
66.60
0.09
-0.00
66.70
0.00
0
可知在t=66.7s时，船到达对岸B点。做 x,y与时间t的关系图：从曲线上可以看出，0到30s这段时间内，y的增长几乎呈线性关系，即小船几乎研直线匀速前进。
现在求解析解并将之与数值解对比：
function [ x] = jiexi(y,k)
x=-0.5*(-0.01)^k*y.^(k+1)+0.5*(-0.01)^(-k).*y.^(1-k);
end

y=-100:1:0;
x2=jiexi(y,0.5);
plot(x2,y,'ro',x(:,1),x(:,2));legend('解析解','数值解',-1);

从轨迹曲线中也可以看到，用数值解得到结果与解析解几乎重合，可信度很高。
（3）当改变v1的值时，解析式中的k值将发生变化，此时将对结果产生影响。利用MATLAB计算和绘图也可以发现，渡河时间及航行曲线都发生了变化。
v1=0时，k=0。说明河水静止不动，这种情况下，小船的总速度就是它在静水中的速度，于是沿着AB直线便可到达对岸，计算结果表示，t=50s时，小船到达B点。

v1=0.5m/s时，k=0.25，得到的曲线如下：
这种情况下，t=53.33s时，小船到达B点，比起v1=1m/s时，小船在x方向上走过的距离缩短了大约一半，总路程缩短了许多。
v1=1.5m/s时，k=0.75，得到的曲线如下：

这种情况下， t=114.28s，小船到达对岸，渡河时间明显长了许多。而且由数值解的疏密程度也可以看出，小船花费较少时间久到达x的最大位移处，但是却花了相当长的时间回到x=0的目的地，可见1.5m/s的河水流速给小船到达终点造成了巨大阻碍。
v1=2m/s时，k=1，得到的曲线如下：

这种情况下，小船无论如何都无法到达B点，只能到达B点下游50米处。从解析式中也可以看到，k=1时，有x(y)=-0.005y2+50。曲线呈开口朝左的抛物线状。从速度合成的三角形上来看，由于v1和v2长度相等，v1的方向也已确定，它们的合速度的方向与v1的夹角不可能大于90°，也就是说在x的分位移始终在增加，不可能减小。即使v2沿着水流逆向，也只能使合速度为零，此时正好是小船停在B点下游50米处的情况。类似的问题在高中物理出现过。

综合上述曲线，有下图：

仔细观察解析式可以发现影响船过河轨迹仅是k值，即船与河水的相对速度，我们不妨作此假设。
如果我们用v1=2m/s，v2=4m/s与前面的结果做下对比（我们仅做数值解的对比，因为解析解相同）。结果证明当v1=2m/s，v2=4m/s 时船的航行轨迹与v1=1m/s，v2=2m/s 航行轨迹相同，但时间缩短为33.33s。而且继续实验当v1=4m/s，v2=8m/s时，船的航行轨迹也与前两种情况相同，但过河时间为16.66s，说明过河时间与船速成倒数关系，这是从航行轨迹相同、路程相等可以很自然导出的关系。
在实际问题中，人们通常会对水的流速做初步判断，以调整船行驶的方向，而不是单纯的每个时刻都将方向对准目的地。同时由于水的流速也不是始终不变的，在不同的流域和不同的时刻流速都可能不同，本题中只是采取了最为简单的数学模型进行近似计算。

题9两种群相互竞争的模型如下：
X’(t)=r1*x*(1-x/n1-s1*y/n2);
Y’(t)=r2*y*(1-s2*x/n1-y/n2);
其中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，s2=0.7，再分析结果。由此你能得到什么结论，请用各参数生态学上的含义做出解释。
（3）实验s1=0.8，s2=0.7和s1=1.5，s2=1.7时会有什么结果，并作出解释。

模型及其求解
（1）编写程序如下：
function [ dx ] = compt( 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)];
end
ts=0:0.5:20;
x0=[10,10];
[t,x]=ode23s(@compt,ts,x0);
[t,x];
plot(t,x);title('t充分大后x(t),y(t)的变化趋势');xlabel('时间t');ylabel('种群数量');gtext('x(t)');gtext('y(t)');
plot(x(:,1),x(:,2));xlabel('x');ylabel('y');title('相图x,y');

设定t从0到20，得出t,x(t),y(t)的数值关系
t
X(t)
Y(t)
0
10.0000
10.0000
0.5000
15.0556
13.7363
1.0000
21.8041
17.4479
1.5000
30.1504
20.2070
2.0000
39.6707
21.1843
2.5000
49.6905
20.1332
3.0000
59.4906
17.4859
3.5000
68.4731
14.0317
4.0000
76.2336
10.5329
4.5000
82.5937
7.4878
5.0000
87.5743
5.0991
5.5000
91.3224
3.3592
6.0000
94.0516
2.1579
6.5000
95.9832
1.3614
7.0000
97.3207
0.8479
7.5000
98.2306
0.5233
8.0000
98.8409
0.3209
8.5000
99.2457
0.1960
9.0000
99.5118
0.1194
9.5000
99.6855
0.0726
10.0000
99.7982
0.0441
10.5000
99.8709
0.0267
11.0000
99.9177
0.0162
11.5000
99.9477
0.0098
12.0000
99.9668
0.0060
12.5000
99.9790
0.0036
13.0000
99.9867
0.0022
13.5000
99.9916
0.0013
14.0000
99.9947
0.0008
14.5000
99.9967
0.0005
15.0000
99.9979
0.0003
15.5000
99.9987
0.0002
16.0000
99.9992
0.0001
16.5000
99.9995
0.0001
17.0000
99.9997
0.0000
17.5000
99.9998
0.0000
18.0000
99.9999
0.0000
18.5000
99.9999
0.0000
19.0000
100.0000
0.0000
19.5000
100.0000
0.0000
20.0000
100.0000
0.0000

由统计数据可以看出，t=17时，乙种群的数量已经为0，之后甲种群的数量达到饱和。由图线更能很好的分析出甲乙两个种群的变化趋势。刚开始时两种群的数量都在增长，但环境的容纳量有限，甲乙两种群不可避免的竞争，于是甲的数量持续增长而乙数量在达到一个峰值后逐渐减少，两种群数量相差悬殊。最后乙灭绝而甲繁荣。
究其原因，s1<s2，表明甲消耗供养对方的食物的能力要远强于乙。甲不但可以放心吃供养自己的食物，还能大量消耗供养乙的食物。由于能够获取更多的食物，这样表现为甲的生存能力——种群增长率要大于乙。所以甲种群走向繁荣，乙种群濒于灭绝。

（2）探究r1，r2，n1，n2，x0，y0对种群数量的影响。除探究的因素以外，要让其他因素保持（1）中数值不变。
A. 首先实验自然增长率r。对于r1>r2，显然甲增加的会更快，结果与（1）相同。而对于r2>r1，不妨设r1=1,r2=4,作图。

可见开始的时候，凭借着高自然增长率，乙种群数量超过了甲。但之后甲种群数量开始稳步上升，乙种群数量则不可避免的下降，最后甲种群再次繁荣，乙种群灭绝。且达到最后稳定的时间比（1）要短，也许是食物消耗过快的结果。
由此可见，乙自然增长率r的提高并不能让乙种群避免灭绝的命运。

B. 其次探究最大容量n的影响。显而易见，如果n1>n2，则趋势与（1）相同。仅仅是甲的最大数量更大了。
对于n2>n1，不妨设n2=300,n1=100。作图：

可见乙种群最大容量的增大，也不能避免他灭绝的命运。从生态学上来看，这样的结果符合常理。
C. 之后探究种群数量初值。对于x0>y0，显然趋势和最后结果与（1）相同，甲的增长会更加迅速。对于x0<y0，设x0=10,y0=40，作图：

可见，乙种群初始数量的增大，虽然在开始时有过一个高峰，但是并不能避免最后灭绝的命运。从生态学上解释，当一个种群生存能力很弱时，即使初始数量很多，最后也会灭绝。
下面探究s1，s2。设s1=1.5。s2=0.7。其他因素与（1）相同。

可见结果发生了逆转：甲种群最终走向了衰落，几乎灭绝，而乙种群走向了繁荣。
由此看出，能否消耗更多供养自己和供养对手的食物和资源，是一个种群能否繁衍生存的关键因素。如果自己的食物和对手的食物都不能有效地“为我所用”，那么一个种群就难逃灭绝的命运。
数据如下：
t
X(t)
Y(t)
0
10.0000
10.0000
0.5000
14.1649
14.8747
1.0000
18.8098
21.1830
1.5000
23.1789
28.6762
2.0000
26.4025
36.8127
2.5000
27.9436
44.9840
3.0000
27.7660
52.6886
3.5000
26.2263
59.6679
4.0000
23.8150
65.8441
4.5000
20.9810
71.2424
5.0000
18.0622
75.9298
5.5000
15.2627
79.9823
6.0000
12.7023
83.4531
6.5000
10.4474
86.4135
7.0000
8.5052
88.9142
7.5000
6.8680
91.0077
8.0000
5.5082
92.7472
8.5000
4.3924
94.1805
9.0000
3.4860
95.3523
9.5000
2.7558
96.3033
10.0000
2.1717
97.0702
10.5000
1.7070
97.6851
11.0000
1.3390
98.1756
11.5000
1.0487
98.5653
12.0000
0.8202
98.8738
12.5000
0.6408
99.1173
13.0000
0.5003
99.3091
13.5000
0.3903
99.4597
14.0000
0.3043
99.5779
14.5000
0.2372
99.6705
15.0000
0.1848
99.7429
15.5000
0.1440
99.7995
16.0000
0.1121
99.8437
16.5000
0.0873
99.8782
17.0000
0.0680
99.9051
17.5000
0.0529
99.9261
18.0000
0.0412
99.9424
18.5000
0.0321
99.9552
19.0000
0.0250
99.9651
19.5000
0.0194
99.9728
20.0000
0.0151
99.9788

(3)当s1=0.8，s2=0.7时，两者都小于1。预计到竞争会很激烈，以ts=0:0.5:100;作图：

可见，甲乙两个种群竞争十分激烈。开始时增长率几乎相同，到后来乙种群数量超过了甲种群，但是甲种群也并没有因此走向不可逆转的衰亡，最后两者种群数量达到稳定。
部分数据如下：
t
X(t)
Y(t)
0
10.0000
10.0000
0.5000
14.7714
14.8619
1.0000
20.7694
21.0796
1.5000
27.5290
28.2659
2.0000
34.2315
35.6596
2.5000
40.0812
42.4607
3.0000
44.5901
48.1110
3.5000
47.7058
52.4606
4.0000
49.6508
55.6472
4.5000
50.7366
57.9272
5.0000
51.2511
59.5615
5.5000
51.3929
60.7411
6.0000
51.3135
61.6167
6.5000
51.1212
62.3048
7.0000
50.8523
62.8450
7.5000
50.5623
63.3042
8.0000
50.2601
63.6927
8.5000
49.9595
64.0316
9.0000
49.6719
64.3400
9.5000
49.3972
64.6178
10.0000
49.1372
64.8662
20.0000
46.3654
67.3807
30.0000
45.6776
67.9864
40.0000
45.5066
68.1363
50.0000
45.4661
68.1717
60.0000
45.4571
68.1796
70.0000
45.4551
68.1813
80.0000
45.4547
68.1817
83.0000
45.4546
68.1817
83.5000
45.4546
68.1818
84.0000
45.4546
68.1818
当t=30以后，两个种群的数量基本保持不变，已经可以视为到达稳定状态。
分析：s1=0.8，s2=0.7时，两者都小于1,表示两种群食用供养自己资源的能力要强于食用供养对手资源的能力。在这种情况下，两种群竞争激烈，而没有一方占据绝对优势，所以两种群最后的数量都能够达到稳定值。还可以看出，s微小的差异，会造成最后种群数量较大的不同。
若s1=1.5，s2=1.7，两者均大于1。作图如下：

很惊奇的，最后甲走向了繁荣而乙走向了灭绝。
分析：s1=1.5，s2=1.7时，两者都大于1,表示两种群食用供养对手资源的能力要强于食用供养自己资源的能力。相当于“攻强守弱”。而这种竞争情形下，更强消耗能力成为生存的决定性因素，此种能力更强的一方将会竞争中胜出。与s都小于1的情况截然不同。

收集于网络，如有侵权请联系管理员删除