实验4常微分方程数值解.docx

1、实验4常微分方程数值解实验4 常微分方程数值解化工系 毕啸天 2010011811【实验目的】1. 练习数值微分的计算；2. 掌握用MATLAB 软件求微分方程初值问题数值解的方法；3. 通过实例学习用微分方程模型解决简化的实际问题；4. 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。【实验内容】题目3小型火箭初始重量为1400kg，其中包括1080kg 燃料。火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N 的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最

2、高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。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-18t4. 5.根据牛顿第二运动定律，有。代入数据有解出由以上5条分析，我们得到了一个常微分方程组： 初值条件为：v0=0,h(0)=0,t60s.3.2 程序代

3、码根据常微分方程组的初值问题，在MATLAB中计算数值解。记x(1) = h，x(2)= v，x =(x(1), x(2)T首先编写M文件function dx = Rocket(t,x)dx=x(2);(32000-0.4*x(2)2)/(1400-18*t)-9.8;%以向量形式表示微分方程endts=0:60 %终点时间为60s，步长定义为1即可x0=0,0;t,x=ode45(Rocket,ts,x0);t,xplot(t,x(:,1),grid,title(图1.高度-时间) xlabel(t/s) ylabel(h/m) pauseplot(t,x(:,2),grid, title

4、(图2.速度-时间) xlabel(t/s)ylabel(v/(m/s)pausea=(32000-0.4*x(:,2).2)./(1400-18*t)-9.8plot(t,a),grid,title(图3.加速度-时间)xlabel(t/s)ylabel(a/(m/s2)t,x(:,1),x(:,2),a由MATLAB计算得到从开始上升到关闭引擎瞬间的情况如下表：thva00013.05716.573713.18913.305226.44426.57713.453359.76240.06213.4974106.5753.53513.4335166.7966.8913.2616240.2780

5、.02112.9857326.7292.82912.6128425.79105.2212.1529536.99117.1111.61710659.8128.4311.02111793.63139.1410.3812937.85149.189.7083131091.8158.559.0209141254.7167.238.3309151425.9175.227.6502161604.8182.556.9901171790.8189.226.3593181983.1195.275.7646192181.2200.755.2095202384.5205.74.6946212592.4210.184.

6、222222804.5214.193.7943233020.6217.793.412243240.1221.013.073253462.7223.922.7726263687.9226.562.5044273915.6228.972.2677284145.6231.142.0633294377.8233.111.8898304611.9234.911.7433314847.7236.571.6178325085238.141.5062335323.8239.611.4095345564.1240.991.3293355805.8242.281.265366048.7243.51.2139376

7、292.9244.681.1708386538.1245.831.1303396784.5246.961.0947407032248.051.0663417280.5249.11.0456427530.2250.121.0308437780.9251.141.0178448032.5252.151.0024458285.1253.160.98757468538.8254.150.97632478793.4255.120.96964489049256.070.96629499305.6257.030.96239509563.1257.990.95272519821.5258.950.941185

8、210081259.90.933735310341260.830.932785410603261.750.936335510865262.670.936955611128263.610.925855711392264.540.913795811657265.460.910635911923266.350.916126012190267.260.91701由数据可以看出，引擎关闭瞬间，火箭的高度为h=12190m，速度v=267.26m/s，加速度a=0.91701m/s2。3.3燃烧耗尽后物理模型分析当引擎关闭后，火箭由重力作用上升，燃料耗尽后火箭质量为320kg。由牛顿第二定律可再次列出微分

9、方程如下:初值条件即为(1)中的末态条件,t60s.3.4程序代码记y(1) = h，y(2)= v，y =(y(1), y(2)T首先编写M文件function dy = Rocket2(t,y)dy=y(2);-9.8-0.4*y(2).2/320;%以向量形式表示微分方程endts=0:60 %终点时间为60s，步长定义为1即可x0=0,0;t,x=ode45(Rocket,ts,x0);t,xfor n=1:2000 %此处利用一个循环语句找出时间终点T=80-0.01*n; %估计出时间不会晚于80tss=60:0.02:T;y0=x(61,1),x(61,2);option=ode

10、set(reltol,1e-3,abstol,1e-6);t2,y=ode45(Rocket2,tss,y0,option);t2,y;if y(:,2)=0breakendendplot(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)pausea=(32000-0.4*x(:,2).2)./(1400-18*t)-9.8;a2=-9.8-0.4

11、*y(:,2).2/320;plot(t,a,b,t2,a2,r),grid,title(图3.加速度-时间)xlabel(t/s)ylabel(a/(m/s2)由MATLAB数据计算可知，当t=71.31s 时，火箭上升到最大高度h= 13115m，此时火箭的速度v=0.019874，几乎可认为已经停止，加速度a=-9.8m/s2。【小结】由以上三图的曲线可以看出，起初火箭迅速加速，而且速度值不大时，燃料动力的贡献要超过空气阻力，而质量的减少又会使加速度变大，故而两个因素共同作用使加速度较为稳定。此后速度增加，阻力变大，加速度渐渐变为0.当燃料耗尽时，推动力消失，加速度突变负，而此时速度由增

12、大转为减少。此时速度迅速下降，但是由于它还是正值，故高度上升。速度为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 建立数学模型及分析此种模型的前提是船并不事先知道水速，否则

13、只要调整一个合适的角度，直接沿直线通过即可。现小船不知道水速，则它的策略应为始终使船头瞄准B点。对速度进行xy两个方向的分解，可列出常微分方程如下：其初值条件为x(0)=0,y(0)=-d.6.2 常微分方程解析解两式相除，得变量代换，令, 则有dxd代入有分离变量积分可得代入初值，可解出整理得解析解最终表达式为6.3 MATLAB数值解先建立M文件函数function dx=Guohe(t,x,v1,v2)s=(x(1)2+x(2)2)0.5;dx=v1-x(1)/s*v2;-x(2)/s*v2;endts=0:0.01:80;x0=0,-100;option=odeset(reltol,1

14、e-6,abstol,1e-9);t,x=ode23s(Guohe,ts,x0,option,1,2);plot(t,x),gridtitle(图1.分位移-时间)xlabel(t/s)ylabel(x,y/m)gtext(x(t)gtext(y(t)pauseplot(x(:,1),x(:,2),grid, title(图2.过河轨迹) xlabel(x/m) ylabel(y/m) t,x(:,1),x(:,2)节选部分数据如下:txy0000.010.0099990.0099990.020.0199960.0199960.030.0299910.0299910.040.0399840.0

15、399840.050.0499750.0499750.060.0599640.0599640.070.0699510.0699510.080.0799360.0799360.090.0899190.0899190.10.09990.09990.110.109880.109880.120.119860.119860.130.129830.129830.140.13980.13980.150.149770.149770.160.159740.159740.170.169710.169710.180.179680.179680.190.189640.189640.20.19960.199666.5

16、0.16685-0.001113366.510.15685-0.0009838466.520.14685-0.0008623966.530.13685-0.00074866.540.12685-0.0006434966.550.11685-0.0005460366.560.10685-0.0004565766.570.096853-0.0003751166.580.086853-0.0003016466.590.076853-0.0002361866.60.066853-0.0001787166.610.056853-0.0001292466.620.046853-8.7773e-00566.

17、630.036853-5.4301e-00566.640.026853-2.8827e-00566.650.016853-1.1351e-00566.660.0068532-1.8739e-00666.676.7389e-0110【小结】小船起初速度较快，相当于在顺水被冲下去。后来快到达终点时较慢，相当于逆水而行。6.4解析解的图像先建立函数function x=Guohe2(y,k) x=0.5*(-0.01)(-k)*y.(1-k)-0.5*(-0.01)k*y.(1+k); endy=0:-0.01:-100; x=Guohe2(y,0.5); plot(x,y)，grid; title

18、(图4.过河轨迹，解析解) xlabel(x/m) ylabel(y/m) 可以看出，由数值解法得到的数据绘成的图与解析解的结果几乎一样。6.5改变流速v1=0ts=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),gridtitle(图1.分位移-时间)xlabel(t/s)ylabel(x,y/m)gtext(x(t)gtext(y(t)pauseplot(x(:,1),x(:,2),grid, title(图2.过河轨迹) xlab

19、el(x/m) ylabel(y/m) 6.6 改变流速v1=0.5m/s程序基本相同，在此略去，只给出图像。6.7 改变流速v1=1.5m/s6.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)，画出它们的

20、图形及相图(x,y),说明时间t充分大以后x(t)，y(t)的变化趋势（人们今天看到的已经是自然界长期演变的结局）。（2） 改变r1,r2,n1,n2,x0,y0,但s1,s2不变（或保持s11），计算并分析所得结果；若s1=1.5（1）,s2=0.7(1),再分析结果。由此你得到什么结论，请用各参数生态学上的含义作出解释。（3） 试验当s1=0.8（1）, s2=0.7(1）, s2=1.7(1)时又会有什么结果。能解释这些结果吗？9.1 模型分析此题的模型很类似教材中的“弱肉强食”模型，只是多考虑了竞争资源的问题。微分方程描述了竞争关系和捕食关系。通过对微分方程组进行数值求解，即得出竞争关

21、系的直观表示。9.2 程序代码设x(1)=x,x(2)=y编写M文件如下:function dx=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);endts=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,xplot(t,x),grid,title(图1.时间-种群数量)xlabel(t)ylabe

22、l(x,y)gtext(x(t)gtext(y(t)pauseplot(x(:,1),x(:,2),gridtitle(图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，而乙彻底灭

23、绝。与之前变化不同的是，达到稳定需要更长的时间，约25年。可以看出，增长率减少，种群数量变化的速度会减慢。(2). r1=5,r2=0.4数据分析：最终结果仍然没有改变！但不同的是，由于此时甲物种的固有增长率远高于乙，即竞争优势十分巨大，因此乙物种在初期甚至几乎没有增长就迅速衰减到0.二者约10年时间就达到了平衡。(2). r1=0.4,r2=5数据分析：至此我们基本可以下结论，改变增长率的值不会改变它们最终的演化结果。图中数据表明，甲物种增长率小，增长速度很慢，相当于减缓了乙物种的灭绝进程。乙初期数量迅速增长，但是后来随着甲的增多渐渐灭绝。不同的是，乙大约在11年时灭绝，但是甲的数量持续增长

24、到约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个参数，都不会改变

25、甲乙最终的宿命。甲总会优胜，达到环境容量值，而乙永远逃脱不了灭绝的命运。9.4改变s1,s21.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，都在过量地消耗着对方的资源，但是更大的一

26、方会被消耗更多的资源，最终导致灭亡。【本题总结】s1,s2小于1时，互相消耗程度较轻，因此可以达到平衡共存的状态。但两者都无法达到容量值，因为互相在制约着对方的生长。当其中之一大于1时，它就会因为被消耗资源而灭亡；当s1,s2都大于1时，两物种竞争激烈，最后s1,s2中更小者在竞争中获胜，另一物种灭绝。正所谓物竞天择，适者生存。自然中能供养物种的资源是十分有限的，谁能更好地适应环境，谁就能在竞争中取得最终的胜利。这对我们人类也是一个警示。【作业总结】 如果说上次作业让我简单地了解了数学建模，那么这次作业是让我更加深刻地理解了数学建模的重要意义。当我编完程序，图像显示了种群数量的变化时，我震惊了！一个简单的计算就可以反映出生活