实验4 常微分方程的数值解【实验目的】1．掌握用MATLAB软件求微分方程初值问题数值解的方法；2．通过实例用微分方程模型解决简化的实际问题；3．了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
【实验内容】题3小型火箭初始重量为1400kg，其中包括1080kg燃料。
火箭竖直向上发射时燃料燃烧率为18kg/s，由此产生32000N的推力，火箭引擎在燃料用尽时关闭。
设火箭上升时空气阻力正比于速度的平方，比例系数为m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
模型及其求解火箭在上升的过程可分为两个阶段，在全过程中假设重力加速度始终保持不变，g=s2。
在第一个过程中，火箭通过燃烧燃料产生向上的推力，同时它还受到自身重力（包括自重和该时刻剩余燃料的重量）以及与速度平方成正比的空气阻力的作用，根据牛顿第二定律，三个力的合力产生加速度，方向竖直向上。
因此有如下二式：a=dv/dt=/m=/(1400-18t)dh/dt=v又知初始时刻t=0，v=0，h=0。
记x(1)=h，x(2)=v，根据MATLAB 可以求出0到60秒内火箭的速度、高度、加速度随时间的变化情况。
程序如下：function [ dx ] = rocket( t,x )a=[*x(2)^2)/(1400-18*t)];dx=[x(2);a];endts=0:1:60;x0=[0,0];[t,x]=ode45(@rocket,ts,x0);h=x(:,1);v=x(:,2);a=[*(v.^2))./(1400-18*t)];[t,h,v,a];数据如下：t h v a 000因此，在引擎关闭的瞬间，火箭的速度为s，高度为，加速度为s2。
（2）在第二个阶段，火箭的重量保持不变，没有向上的推力，只收到重力和空气阻力。
因此有如下关系式：a=dv/dt=/m=dh/dt=v假设在80秒之前达到最高点，以60秒时的速度、高度、加速度为初始值进行计算，程序如下：function [ dx ] = rocket2( t,x )dx=[x(2);*x(2)^2/];endts2=60:1:80;x1=[,];[t2,x2]=ode45(@rocket2,ts2,x1);h2=x2(:,1);v2=x2(:,2);a2=*v2.^2./;[t2,h2,v2,a2];数据如下：t2 h2 v2 a2可以看到：在第60秒瞬间，加速度发生了突变，从s2突变为s2；而在第71秒至第72秒之间，速度从正变为负，即速度反向，说明在第71秒中的某个时刻速度为零，火箭达到了最高点。
因此需要对这个时间段进行分析，并且找到速度减小到0的时刻和此时的高度。
以为步长，在71s到72s中重新求解微分方程的数值解。
可见在t=时，速度为，可视为速度为零点，此时最大高度为，加速度。
综合（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点。
已知河水流速v 1与船在静水中的速度v 2之比为k 。
（1）建立描述小船航线的数学模型，求其解析解；（2）设d=100m ，v 1=1m/s ，v 2=2m/s ，用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较；（3）若流速v 1=0，，，2m/s ，结果将如何。
模型及其求解（1）.假设在航行过程中，人们不知道水流的速度，小船的方向始终指向目标B ，因此若以B 为原点，我们可以得到如下方程组：解初值为（x, y ）=( 0, -d) 的常微分方程组，得到解析解为：其中c=–1/d=–,故有事实上，若用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,v 1=1m/s,v 2=2m/s 。
求数值解时，由解析解可以看出，此为刚性方程。
为避免用ode45s 求解时间过长。
采用ode23s 进行求解。
假定100s 可以到达对岸。
ts=0::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/m1dx v dt dy dt ⎧=-⎪⎪⎨⎪=⎪⎩1111022k k k kx c y c y +--+-=1111()(0.01)(0.01)22k k k k x y y y +--=--+-可知在t=时，船到达对岸B点。
做x,y与时间t的关系图：从曲线上可以看出，0到30s这段时间内，y的增长几乎呈线性关系，即小船几乎研直线匀速前进。
现在求解析解并将之与数值解对比：function [ x] = jiexi(y,k)x=*^k*y.^(k+1)+*^(-k).*y.^(1-k);endy=-100:1:0;x2=jiexi(y,;plot(x2,y,'ro',x(:,1),x(:,2));legend('解析解','数值解',-1);从轨迹曲线中也可以看到，用数值解得到结果与解析解几乎重合，可信度很高。
（3）当改变v1的值时，解析式中的k值将发生变化，此时将对结果产生影响。
利用MATLAB计算和绘图也可以发现，渡河时间及航行曲线都发生了变化。
v1=0时，k=0。
说明河水静止不动，这种情况下，小船的总速度就是它在静水中的速度，于是沿着AB直线便可到达对岸，计算结果表示，t=50s时，小船到达B点。
v1=s时，k=，得到的曲线如下：这种情况下，t=时，小船到达B点，比起v1=1m/s时，小船在x方向上走过的距离缩短了大约一半，总路程缩短了许多。
v1=s时，k=，得到的曲线如下：这种情况下， t=，小船到达对岸，渡河时间明显长了许多。
而且由数值解的疏密程度也可以看出，小船花费较少时间久到达x的最大位移处，但是却花了相当长的时间回到x=0的目的地，可见s的河水流速给小船到达终点造成了巨大阻碍。
v1=2m/s时，k=1，得到的曲线如下：这种情况下，小船无论如何都无法到达B点，只能到达B点下游50米处。
从解析式中也可以看到，k=1时，有x(y)=+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 航行轨迹相同，但时间缩短为。
而且继续实验当v1=4m/s，v2=8m/s时，船的航行轨迹也与前两种情况相同，但过河时间为，说明过河时间与船速成倒数关系，这是从航行轨迹相同、路程相等可以很自然导出的关系。
在实际问题中，人们通常会对水的流速做初步判断，以调整船行驶的方向，而不是单纯的每个时刻都将方向对准目的地。
同时由于水的流速也不是始终不变的，在不同的流域和不同的时刻流速都可能不同，本题中只是采取了最为简单的数学模型进行近似计算。
题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=，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=，s2=，再分析结果。
由此你能得到什么结论，请用各参数生态学上的含义做出解释。
（3）实验s1=，s2=和s1=，s2=时会有什么结果，并作出解释。
模型及其求解（1）编写程序如下：function [ dx ] = compt( t,x )r1=1;r2=1;n1=100;n2=100;s1=;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)];endts=0::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)由统计数据可以看出，t=17时，乙种群的数量已经为0，之后甲种群的数量达到饱和。