数学实验第二次作业——常微分方程数值求解[实验相关]

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

2、如下：这是一个典型的牛顿第二定律问题，分析火箭受力情况；先规定向上受力为正数建立数学模型：A燃料未燃尽前，在任意时刻（t60s）火箭受到向上的-F=32000N，向下的重力G=mg，g=9.8，向下的阻力f=kv2, k=0.4, v表示此时火箭速度；此时火箭收到的合力为F1=（F-mg-f）；火箭的初始质量为1400kg，燃料燃烧率为-18kg/s;此刻火箭质量为m=1400-18*t根据牛顿第二定律知，加速度a=F1/m=(F-mg-f)/(m-r*t)= (32000-(0.4.*v.2)-9.8.*(1400-18.*t)由此可利用龙格-库塔方法来实现，程序实现如下Function d

3、x=rockett,x %建立名为rocket的方程m=1400;k=0.4;r=-18;g=9.8; %给出题目提供的常数值dx=x(2);(32000-(k*x(2)2)-g*(m+r*t)/(m+r*t);%以向量的形式建立方程a=(32000-(k*x(2)2)-g*(m+r*t)/(m+r*t); %给出a的表达式End;ts=0:60; %根据题目给定燃烧率计算出燃料燃尽的时间，确定终点x0=0,0; %输入x的初始值t,x=ode15s(rocket,ts,x0); %调用ode15s计算t,x;h=x(:,1);v=x(:,2);plot(t,x(:,1),grid; %绘出火

4、箭高度与时间的关系曲线title(h-t);xlabel(t/s);ylabel(h/m),pause;plot(t,x(:,2),grid ; %绘出火箭速度与时间的曲线关系title(v-t);xlabel(t/s);ylabel(v/m/s),pause;a=(32000-(0.4.*v.2)-9.8.*(1400-18.*t)/(1400-18.*t);plot(t,a),grid; %绘出火箭加速度与时间的曲线关系title(a-t);xlabel(t/s),ylabel(a/m2/s),pause 火箭高度随时间变化的曲线火箭速度随时间变化的曲线火箭加速度随时间变化的曲线数据过多，

5、故截取部分如下第一列为时间，第二列为火箭高度，第三列为火箭速度由此可以，在t=60s时，即火箭燃料燃尽瞬间，引擎关闭瞬间，火箭将到达12912m的高度，速度为267,29m，加速度a=0.9m/s2B燃料燃尽之后，与A 类似，分析受力如下火箭受到向上的F=0向下的重力G=mg，g=9.8，向下的阻力f=kv2, k=0.4, v表示此时火箭速度；此时火箭收到的合力为F2=（-mg-f）；火箭的初始质量为320kg，恒定根据牛顿第二定律,加速度a=F2/m=-g-0.4v2/320;程序实现如下function dx = rocket2( t,x ) %建立以rocket2为名的函数dx=x(2

6、);-9.8-0.4.*x(2).2/320; %以向量的形式建立方程ts=60:120; %给出初始时刻，估计终点时刻x0=12190,267.26; %给出x初始值t,x=ode15s(rocket2,ts,x0); %调用ode15s计算t,xplot(t,x(:,1),grid; %绘出火箭高度随时间变化的曲线title(h-t);xlabel(t/s),ylabel(h/m),pause;plot(t,x(:,2),grid; %绘出火箭速度随时间的变化曲线title(v-t);xlabel(t/s),ylabel(v/m/s),pause;v=x(:,2);a=-9.8-0.4*v

7、.2/320; %给出加速度的具体表达式plot(t,a),grid; %绘出火箭加速度随时间变化的曲线title(a-t);xlabel(t/s),ylabel(a/m2/s),pause得到的曲线图形如下火箭高度随时间的变化曲线从图中可以大致看出，最高点在13km左右，火箭速度随时间的变化曲线加速度随时间变化曲线如下数据表格大致如下从图表中可以看出，在71s左右速度到达0，即此时到达最高处，高度为13117m加速度为-9.8m/m/s2;本题总结：这道题是典型的物理牛顿力学的题目，通过受力的正确分析，可以知道，以h,v为向量建立微分方程即可求解，h的微分是速度v，速度v的微分是加速度a解题

8、过程中存在的难点是：取值步长不太容易确定，而且是哪种算法不确定， 先用ode15s速度较快，ode23s速度差不太多，其他两种速度较慢，等待时间较长5一只小船渡过宽为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）,情况又如何建立数学模型：在任意时刻t，小船位于（x，y），此时速度为v，根据物理中路程与速度的关

9、系，知路程的微分为速度v，由此中，小船在x,y方向上的速度分别为：X：dxdt=v1-v2*xx2+y2Y:dydt=-v2*yx2+y2初始条件为dxdtt0=v1 *（1）dydt t0=v2 *（2）现求其解析解,利用微积分知识*（1）*（2）得dxdy=v1*x2+y2-v2y+xy;*(3)将*（3）右端带根号部分分子分母均y得dxdy=v1*(x/y)2+1-v2+xy;令x/y=p，得到dx/dy=dp/dy*y+p;dx/dy=p;分离变量有y(-k)=c(1+p2+p);代入初值可确定当t=0时，y=-d,x=0,p=x/y=0,C=(-d)(-k)y(-k)/d(-k)=1

10、+x2y2 +x/y;根据隐函数与显函数的关系可得到解析解：x=-d2(y-d)(-k)-(yd)(1+k)已知d=100m,v1=1m/s,v2=2m/s;先利用龙格库塔方法求解渡河时间，及任意时刻小船的位置（x，y），及航行曲线，与解析解比较此时，k=0.5，d=100用MATLB编写源程序如下：function dx = boat(t,x) %建立名为boat的函数 d=-100;v1=1;v2=2; %给定常数值s=(x(1).2+x(2).2).0.5; dx=v1-v2.*x(1)./s;-v2.*x(2)./s; %用向量形式建立方程ts=0:70; %大致估算，确定终点值，给定

11、步长为”1”x0=0,-100; %给出x初始值t,x=ode23s(boat,ts,x0); %调用ode23s计算t,x;plot(t,x),grid, %绘出x(t),y(t)的函数曲线（图21）gtext(x(t);gtext(y(t);pause;plot(t,x(:,1), %绘出x随时间变化的曲线（图22）grid;xlabel(t/s);ylabel(x/m);pause;plot(t,x(:,2),grid, %绘出y随时间变化的曲线（图23）title(y-t);xlabel(t/s),ylabel(y/m);图21 图22图23得到数据如下从表格中数据可知，在大约67s时

12、y=0即船到达对岸目的地，为比较，先进行解析解的求解设计程序如下：function x=f(y)k=0.5;x=-0.5.*(-0.01).k.*y.(k+1)+0.5.*(-0.01).(-k).*y.(-k+1);y=0:-0.1:-100;for i=0:1:1000;x(:,i+1)=xy(-i/10);endplot(x,y);grid;gtext(x);gtext(y);由此可以看出，由数值解和解析解得到的x-y曲线相差不多，所以可以认为解析解正确改变水流速度v1，只要在原有程序基础上重新复制给v1=0,0.5,1.5,2,同时适当改变终点值即可现实现程序如下A:v1=0，d=-1

13、00;v1=0;v2=2;s=(x(1).2+x(2).2).0.5;dx=v1-v2.*x(1)./s;-v2.*x(2)./sts=0:60;x0=0,-100;t,x=ode15s(boat21,ts,x0);t,x;plot(x(:,1),x(:,2),grid,title(y-x);pause;plot(t,x(:,1),grid;xlabel(t/s);ylabel(x/m);plot(t,x(:,1),grid;title(x-t);xlabel(t/s),ylabel(x/m),pause;plot(t,x(:,2),grid,title(y-t);xlabel(t/s),yl

14、abel(y/m);图形如下从此图形中我们可以看到，船并未偏离x=0的点，我们也可以从直观想象中的得到，当水速v1=0时，只要出发时，船头对准目标点，船将一直朝着直线向目的地行进从表格中数据我们也可以很清楚地看到路程与时间是成明显的线性关系的，这是与我们水速为0的必然结果，由此也可以验证我们模型基本正确现改变水速B:v1=0.5程序实现如下d=-100;v1=0.5;v2=2;s=(x(1).2+x(2).2).0.5;dx=v1-v2.*x(1)./s;-v2.*x(2)./sts=0:60;x0=0,-100;t,x=ode15s(boat22,ts,x0);t,x;plot(t,x),g

15、rid,gtext(x(t);gtext(y(t);pause;plot(t,x(:,1),grid;xlabel(t/s);ylabel(x/m);plot(t,x(:,1),grid;title(x-t);xlabel(t/s),ylabel(x/m),pause;plot(t,x(:,2),grid,实现程序过程中发现，终点值设为60，而在53s之后不再有出现，所以可以认为在54s之前就已经达到对岸现在重新改变水速C:v1=1.5d=-100;v1=1.5;v2=2;s=(x(1).2+x(2).2).0.5;dx=v1-v2.*x(1)./s;-v2.*x(2)./s ts=0:120

16、;x0=0,-100;t,x=ode15s(boat23,ts,x0);t,x;plot(x(:,1),x(:,2),grid,title(y-x);pause;plot(t,x(:,1),grid;xlabel(t/s);ylabel(x/m);plot(t,x(:,1),grid;title(x-t);xlabel(t/s),ylabel(x/m),pause;plot(t,x(:,2),grid,title(y-t);xlabel(t/s),ylabel(y与v1=0.5s类似，在114s之后不再给出数据，而我们设定的终点值是120，所以可以大致在114s时到达B点现改变水速D:v1=2

17、程序实现如下d=-100;v1=2;v2=2;s=(x(1).2+x(2).2).0.5;dx=v1-v2.*x(1)./s;-v2.*x(2)./sts=0:300;x0=0,-100;t,x=ode15s(boat24,ts,x0);t,x;plot(x(:,1),x(:,2),grid,title(y-x);pause;plot(t,x(:,1),grid;xlabel(t/s);ylabel(x/m);plot(t,x(:,1),grid;title(x-t);xlabel(t/s),ylabel(x/m),pause;plot(t,x(:,2),grid,title(y-t);xla

18、bel(t/s),ylabel(y/m)曲线如下我们可以从此图中看出，当y=0时，x=50，也就是说，船根本没本法到达正对着的B点，而只能到达对岸，我们可以直观地理解假设我们的船在开始时与水速反向，这时船与水的合速度是0，故船无法前进，而根据方程知道，在船尚未到达对岸之前，只有当船的x轴向速度在模式刻超过水速，才有可能克服因水速而在x轴偏离B点的距离，而重新返回，到达B点，而在此给定的速度中，我们可以看到，船速的分量不可能超过水速，因此不可能到达B本题总结：这道题跟课本例题缉私船有较大的相似处，大体是模仿例题而做，所以在画图上显得有点繁琐，如果采用subplot作图，将会节约一些篇幅，使得作业

19、版面看起来更整洁一点3.两种群的竞争模型如下：其中x(t)，y(t)分别为甲乙两种群的数量；r1，r2为他们的固有增涨率；n1，n2为他们的最大容量。s1的含义是，对于供养甲的资源而言，单位数量乙（相对n2）的消耗为单位数量甲的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不变（或保持s11），计算并分析所得

20、结果；若s1=1.5，s2=0.7,再分析结果，由此你能得到什么结论，请用各参数生态学上的含义做出解释。 （3） 试验当s1=0.8，s2=0.7时会有什么结果；当s1=1.5，s2=1.7时又会有什么结果。你能解释这些结果吗？（1）程序如下：函数：function dx = species1(t,x ); %建立以species1为名的函数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

21、=1:0.01:100; %给定终点值，确定步长x0=10,10; %给定x初始值t,x=ode45(species1,ts,x0); %调用ode45计算A=t,x;plot(t,x(:,1),t,x(:,2),grid; %绘出种群x,y随时间变化的曲线 gtext(x(t);gtext(y(t);xlabel(t);ylabel(x,y);title(x,y-t图像);pause;plot(x(:,1),x(:,2),grid; %绘出x y种群数量随时间变化的曲线xlabel(x);ylabel(y);title(y-x图像);种群x,y随t的变化图形如下：读图可知，当t10后，两种群

22、数量基本不变，x=100,y=0灭绝。从生态学角度分析，以上反映了优胜劣汰的自然规律。很显然，种群x更具有竞争优势，随着自然的长期演变，x得以繁衍，而y最终被灭绝。种群y与x的数量关系如下图：上图为相图(x,y)，由上图也可以得出，x与y自初始点(10，10)点演变至(100,0)，种群y随x的增长，虽然在开始时也有一定的增长，但在种群x增长至一定限度以后，种群y的竞争劣势开始显现，随着x的增长，y逐渐减少，最终灭绝。（2）当s1,s2保持不变：a.其他赋值不变，r1=1,r2=10时，x,y随t的变化及y随x变化的图像如下：其他赋值不变，r1=10,r2=1时，x,y随t的变化及y随x变化的

23、图像如下：由上图我们可以推断出，r在生态学中代表种群的增长率，而从图中我们也可以知道，r对种群的竞争能力并没有起到决定性的作用。虽然在初始过程中，r对种群的数量有显著的影响，但随着时间的推移，达到稳定态时，种群的数量并没有因为r的变化而改变。 b.其他赋值不变，n1=10,n2=100时，x,y随t的变化及y随x变化的图像如下：其他赋值不变，n1=100,n2=10时，x,y随t的变化及y随x变化的图像如下：n在生态学中的意义为，环境对该种群的最大容量，我们可以从图中知道，n值得大小并不会改变种群的竞争能力。而随着时间的推移，竞争优势者的数量逐渐趋于环境对其的最大容量，而竞争劣势者最终依旧会灭

24、亡。 值得注意的现象是，当初始值大于或等于最大容量时，无论竞争是否具有优势，种群数量都会在最开始的时间内有所下降。这在生态学上也是可以理解的。c.其他赋值不变，x0=5,y0=50时，x,y随t的变化及y随x变化的图像如下：其他赋值不变，x0=50,y0=5时，x,y随t的变化及y随x变化的图像如下： 有以上两图，我们可以初步判定，初始值的大小只能决定开始阶段的竞争优势，但随着时间的推移，这种优势不断减少，而最终达到稳定状态时与初始值并没有关系。综合以上几种情况，我们可以推断，s才是决定种群在生态中的竞争能力的根本因素。当s1=1.5,s2=0.7时a.其他赋值不变，r1=1,r2=10时，x

25、,y随t的变化及y随x变化的图像如下：其他赋值不变，r1=10,r2=1时，x,y随t的变化及y随x变化的图像如下：同中情形相同，改变r值只能改变初始暂态过程的斜率，对最终的稳态没有作用。b.其他赋值不变，n1=10,n2=100时，x,y随t的变化及y随x变化的图像如下：其他赋值不变，n1=100,n2=10时，x,y随t的变化及y随x变化的图像如下：上图可知，同情形相同，最大容量n不会对稳定态的情况产生影响c.其他赋值不变，x0=5,y0=50时，x,y随t的变化及y随x变化的图像如下：其他赋值不变，x0=50,y0=5时，x,y随t的变化及y随x变化的图像如下：从上图可以看出，初始值的差

26、别只会影响进入稳定区域前的变化过程，与最终的稳定状态无关。与之前的情况是一样的。(3) 当s1=0.8,s2=0.7时：a.当r1=r2=1,n1=n1=100,x0=y0=10时，x,y随t的变化及y随x变化的图像如下：可以看出这种情况与上文中情况有明显的差别： 当参数s1和s2均小于1时，两种种群的数量均在初始时快速上升。两条曲线上升至某一点后x增长率下降，直至种群数量降低至一稳定值；而y则持续上升，但也无法达到最大容量；之后两种种群稳定保持在一定数量，二者平衡共存，生态系统达到稳定。b.r1=10,r2=1,n1=100,n2=100,x0=y0=10时，x,y随t的变化及y随x变化的图

27、像如下：r1=1,r2=10,n1=100,n2=100,x0=y0=10时，x,y随t的变化及y随x变化的图像如下：r1=1,r2=1,n1=100,n2=100,x0=5,y0=50时，x,y随t的变化及y随x变化的图像如下：r1=1,r2=1,n1=100,n2=100,x0=50,y0=5时，x,y随t的变化及y随x变化的图像如下：由图像可知，与之前的情形相同，改变增长率r和初始值同样只能改变达到稳定前的初始过程，并不能影响稳定态的情况。c.r1=1,r2=1,n1=10,n2=100,x0=y0=10时，x,y随t的变化及y随x变化的图像如下： r1=1,r2=1,n1=100,n2

28、=10,x0=y0=10时，x,y随t的变化及y随x变化的图像如下：由上图可知，改变最大容量会影响到稳定态的情况，由此可得出，当s均小于1时，环境对种群的最大容量n和s共同决定稳定态的情况。当s1=1.5,s2=1.7时a.当r1=r2=1,n1=n1=100,x0=y0=10时，x,y随t的变化及y随x变化的图像如下：x最终达到最大容量，y最终完全灭绝。b.r1=10,r2=1,n1=100,n2=100,x0=y0=10时，x,y随t的变化及y随x变化的图像如下：r1=1,r2=10,n1=100,n2=100,x0=y0=10时，x,y随t的变化及y随x变化的图像如下：由以上两图可以看出，固有增长率的差别也会影响到最终竞争的结果：当y0的固有增长率超过x0很多时，最终y就有可能达到最大容量而x可能会灭绝。c.r1=1,r2=1,n1=10,n2=100,x0=y0=