实验4常微分方程数值解
实验目的：
1. 练习数值积分的计算；
2. 掌握用MATLAB软件求微分方程初值问题数值解的方法； 3. 通过实例学习用微分方程模型解决简化的实际问题；
4. 了解欧拉方法和龙格——库塔方法的基本思想和计算公式，及稳定性等概念。 实验内容：
3.小型火箭初始质量为1400kg，其中包括1080kg燃料，火箭竖直向上发射是燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升是空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度，速度，加速度，及火箭到达最高点是的高度，速度和加速度，并画出高度，速度，加速度随时间变化的图形。 解答如下：
这是一个典型的牛顿第二定律问题，分析火箭受力情况； 先规定向上受力为正数
建立数学模型：
A燃料未燃尽前，在任意时刻（t<60s）
火箭受到向上的-F=32000N， 向下的重力G=mg，g=9.8，
向下的阻力f=kv^2, 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 [dx]=rocket[t,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; %绘出火箭高度与时间的关系曲线 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/m^2/s'),pause 火箭高度随时间变化的曲线
火箭速度随时间变化的曲线
* *
火箭加速度随时间变化的曲线
数据过多，故截取部分如下
* *
第一列为时间，第二列为火箭高度，第三列为火箭速度
515253545556575859609823.74810083.1810343.5910604.9410867.2111130.3911394.4911659.4911925.4212192.25258.9677259.9233260.8613261.794262.7221263.6462264.5671265.482266.3898267.2908 由此可以，在t=60s时，即火箭燃料燃尽瞬间，引擎关闭瞬间，火箭将到达12912m的高度，速度为267,29m，加速度a=0.9m/s^2
B燃料燃尽之后，与A 类似，分析受力如下
火箭受到向上的F=0
向下的重力G=mg，g=9.8，
向下的阻力f=kv^2, k=0.4, v表示此时火箭速度； 此时火箭收到的合力为F2=（-mg-f）； 火箭的初始质量为320kg，恒定
根据牛顿第二定律,加速度a=F2/m=-g-0.4v^2/320; 程序实现如下
function [ dx ] = rocket2( t,x ) %建立以rocket2为名的函数 dx=[x(2);-9.8-0.4.*x(2).^2/320]; %以向量的形式建立方程
ts=60:120; %给出初始时刻，估计终点时刻
* *
x0=[12190,267.26]; %给出x初始值 [t,x]=ode15s(@rocket2,ts,x0); %调用ode15s计算 [t,x]
plot(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.^2/320; %给出加速度的具体表达式 plot(t,a),grid; %绘出火箭加速度随时间变化的曲线 title('a-t');
xlabel('t/s'),ylabel('a/m^2/s'),pause 得到的曲线图形如下 火箭高度随时间的变化曲线
* *
从图中可以大致看出，最高点在13km左右，
火箭速度随时间的变化曲线
* *
加速度随时间变化曲线如下
* *
数据表格大致如下
0.0064 1.2822 0.0093 -0.0021 0.0065 1.2905 0.0074 -0.0017 0.0066 1.2971 0.0059 -0.0014 0.0067 1.3023 0.0046 -0.0012 0.0068 1.3063 0.0034 -0.0011 0.0069 1.3092 0.0023 -0.0010 0.0070 1.3110 0.0013 -0.0010 0.0071 1.3117 0.0003 -0.0010 0.0072 1.3116 -0.0007 -0.0010 0.0073 1.3104 -0.0017 -0.0010 从图表中可以看出，在71s左右速度到达0，即此时到达最高处，高度为13117m加速度为-9.8m/m/s^2;
本题总结：
这道题是典型的物理牛顿力学的题目，通过受力的正确分析，可以知道，以[h,v]为向量建立微分方程即可求解，h的微分是速度v，速度v的微分是加速度a
解题过程中存在的难点是：取值步长不太容易确定，而且是哪种算法不确定， 先用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，根据物理中路程与速度的关系，知路程的微分为速度v，由此中，小船在x,y方向上的速度分别为： X：dt=v1-v2*Y:dt=-v2*dy
ydx
x√x2+y^2
√x2+y^2 * *
初始条件为
dxdtdydt
(t0)=v1 *（1） (t0)=v2 *（2）
现求其解析解,利用微积分知识 *（1）÷*（2）得dy=
dx
v1∗√x2+y^2−v2y
+y;*(3) x
将*（3）右端带根号部分分子分母均÷y得 =dy
dx
v1∗√(x/y)2+1−v2
+y;
x
令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+y^2 +x/y; 根据隐函数与显函数的关系可得到
x2解析解：x=-2[(−d)^(-k)-(d)^(1+k)]
已知d=100m,v1=1m/s,v2=2m/s;
先利用龙格库塔方法求解渡河时间，及任意时刻小船的位置（x，y），及航行曲线，与解析解比较
此时，k=0.5，d=100 用MATLB编写源程序如下：
dyy
* *
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; %大致估算，确定终点值，给定步长为”1” x0=[0,-100]; [t,x]=ode23s(@boat,ts,x0); [t,x];
plot(t,x),grid, gtext('x(t)');gtext('y(t)'); pause;
plot(t,x(:,1)), grid; xlabel('t/s'); ylabel('x/m'); pause;
plot(t,x(:,2)),grid, title('y-t');
xlabel('t/s'),ylabel('y/m'); 图21
%给出x初始值 %调用ode23s计算 %绘出x(t),y(t)的函数曲线（图21） %绘出x随时间变化的曲线（图22） %绘出y随时间变化的曲线（图23） * *
图22
图23
* *
得到数据如下
01234567800.9900371.9597412.9088283.8368224.7429025.6273376.4893587.327888-100-98-96.0002-94.0008-92.002-90.0042-88.0073-86.0119-84.0184
* *
606162636465666768696.5375025.5960414.6374223.6645442.6803231.6878170.690109########2.15E-081.63E-08-1.74536-1.26513-0.86066-0.53305-0.28315-0.11154-0.0185000
从表格中数据可知，在大约67s时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); end plot(x,y); grid; gtext('x'); gtext('y');
* *
由此可以看出，由数值解和解析解得到的x-y曲线相差不多，所以可以认为解析解正确 改变水流速度v1，只要在原有程序基础上重新复制给v1=0,0.5,1.5,2,同时适当改变终点值即可现实现程序如下
A:v1=0，
d=-100;v1=0;v2=2; s=(x(1).^2+x(2).^2).^0.5; dx=[v1-v2.*x(1)./s;-v2.*x(2)./s] ts=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'),ylabel('y/m');
图形如下
* *
从此图形中我们可以看到，船并未偏离x=0的点，我们也可以从直观想象中的得到，当水速v1=0时，只要出发时，船头对准目标点，船将一直朝着直线向目的地行进
* *
40414243444546474849500-200-180-160-140-120-100-80-60-40-20######## 从表格中数据我们也可以很清楚地看到路程与时间是成明显的线性关系的，这是与我们水速为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)./s] ts=0:60; x0=[0,-100];
[t,x]=ode15s(@boat22,ts,x0); [t,x];
plot(t,x),grid,
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,
* *
* *
48495051525354555.0447154.350443.5713872.6861641.6693410.474746-7.15145-5.54958-4.01218-2.57055-1.2705-0.21628 实现程序过程中发现，终点值设为60，而在53s之后不再有出现，所以可以认为在54s之前就已经达到对岸 现在重新改变水速
* *
C:v1=1.5
d=-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; 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
* *
* *
1051061071081091101111121131141154.6506994.1507013.6507023.1507022.6507022.1507031.6507031.1507030.6507030.150703-0.00752-0.00477-0.00285-0.00158-0.0008-0.00034-0.00012################1.53E-07 与v1=0.5s类似，在114s之后不再给出数据，而我们设定的终点值是120，所以可以大致在
* *
114s时到达B点 现改变水速
D:v1=2
程序实现如下 d=-100;v1=2;v2=2; s=(x(1).^2+x(2).^2).^0.5; dx=[v1-v2.*x(1)./s;-v2.*x(2)./s] ts=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');
xlabel('t/s'),ylabel('y/m') 曲线如下
我们可以从此图中看出，当y=0时，x=50，也就是说，船根本没本法到达正对着的B点，而只能到达对岸，我们可以直观地理解假设我们的船在开始时与水速反向，这时船与水的合速度是0，故船无法前进，而根据方程知道，在船尚未到达对岸之前，只有当船的x轴向速度在模式刻超过水速，才有可能克服因水速而在x轴偏离B点的距离，而重新返回，到达B点，而在此给定的速度中，我们可以看到，船速的分量不可能超过水速，因此不可能到达B
* *
* *
本题总结：
这道题跟课本例题缉私船有较大的相似处，大体是模仿例题而做，所以在画图上显得有点繁琐，如果采用subplot作图，将会节约一些篇幅，使得作业版面看起来更整洁一点
3.两种群的竞争模型如下：
dxxyr1x(1s1)dtn1n2dyyxr2y(1s2)dtn2n1
其中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不变（或保持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] = 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=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的变化图形如下：
* *
读图可知，当t<10时，种群x数目随时间增大,种群y的数量随时间时减小，当t>10后，两种群数量基本不变，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变化的图像如下：
由上图我们可以推断出，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值得大小并不会改变种群的竞争能力。而随着时间的推移，竞争优势者的数量逐渐趋于环境对其的最大容量，而竞争劣势者最终依旧会灭亡。
值得注意的现象是，当初始值大于或等于最大容量时，无论竞争是否具有优势，种群数量都会在最开始的时间内有所下降。这在生态学上也是可以理解的。
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,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变化的图像如下：
从上图可以看出，初始值的差别只会影响进入稳定区域前的变化过程，与最终的稳定状态无关。与之前的情况是一样的。 (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变化的图像如下：
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=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=10时，x,y随t的变化及y随x变化的图像如下：
* *
r1=1,r2=1,n1=100,n2=10,x0=y0=10时，x,y随t的变化及y随x变化的图像如下：
由以上两图可以看出，最终的稳定态的情况还与最大容量有关，最大容量不同最终的竞争结果也不同。
d.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变化的图像如下：
由以上两图可知：初始值的差别可能导致两物种最终稳定种群数量的差别，当y0远大于x0时，最终在竞争中获胜的就可能改变为y，而不再是x。
我们还可以发现，当s均大于1时，不会再像s均小于1时那样，可以出现二者 平衡共存的现象，而是必有一方灭绝。
本题总结：
综上所述，s1,s2小于1时消耗生存资源的严重程度较轻，所以甲乙物种可以共存，但两者都大不到最大值；当其中之一大于1时，对应作用的物种就会由于生存资源的过度消耗而灭
* *
绝；当s1,s2都大于1时，两物种竞争激烈，最后s1,s2中更大者对应作用的物种灭绝。所谓物尽天择，自然资源是有限的，需要更少资源就能生存的物种在竞争中将占有优势。
还有一种情况，当s1,s2都大于1但相等时，由于方程的对称性，甲乙两物种都能生存下来，但都不能达到最大值。这种情况在自然界里几乎时不可能的。 s1=s2=1.5的结果图：
最后甲乙两物种数量都稳定在40上。
* *
实验总结：
此次实验花费时间较长，但是学到了不少东西，比如数值求解的具体实现，不同算法的速度，快速做多张图像的Subplot.
1.在应用知识方面，不仅用到MATLAB有关数值解的程序设计，还有在微积分三中学到的微分方程求解析式
2.在此次实验中基本用到的是刚性方程求解方法ode23s,ode15s，实验过程中使用过ode23,ode45,但速度太慢。
3.在作图上，在做第三题时需要做的曲线实在太多，才知道subplot的强大，第三题的程序和解答都要比过河问题更简洁
4.在解题思路上，大致是在理解题意的基础上，建立起物理模型，，此次主要是物理受力分析问题，具体用到的是牛顿第二定律，加上部分几何知识，列出方程，再用数学方法求解。
因篇幅问题不能全部显示，请点此查看更多更全内容