数学试验第二次作业

1、3.问题：小型火箭初始质量为1400kg,其中包括1080kg燃料，火箭竖直向上发射时燃料燃烧率为18kg/s,由此产生32000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比速度的平方，比例系数为0.4kg/m,求引擎关闭瞬间火箭的高度、速度、加速度及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。模型：设速度为v,根据牛顿第二定律，可得微分方程在0<t<60时，dv/dt=(32000-0.4*vA2-9.8*(-18*t+1400)/(-18*t+1400)在引擎关闭，火箭上升至最高点之前dv/dt=-(18*vA2+9.8*320)/

2、320;计算方法：编写dv/dt的函数functiondv=rocket1(t,v)dv=(32000-0.4*vA2-9.8*(-18*t+1400)/(-18*t+1400);functiondv=rocket2(t,v)dv=-(18*vA2+9.8*320)/320;用龙格-库塔方法求这两个常微分方程，再利用梯形公式求出v对t的积分，得到火箭上升的高度。并输出60s时的速度、加速度和高度，以及总的上升高度。ts1=0:0.1:60;v0=0;t1,v1=ode45(rocket1,ts1,v0);H1=trapz(t1,v1)dv=(32000-0.4*v1(601)A2-9.8*(-

3、18*t1(601)+1400)/(-18*t1(601)+1400)ts2=60:0.1:71.3;v0=v1(601)t2,v2=ode45(rocket2,ts2,v0);H2=trapz(t2,v2)H=H1+H2dv1=(32000-0.4*v1.A2-9.8.*(-18.*t1+1400)./(-18.*t1+1400);dv2=-(18.*v2.A2+9.8.*320)./320;dv=dv1;dv2;plot(t,dv),grid,t=t1;t2;v=v1;v2;pause,plot(t,v),grid,h(1)=0;fori=2:length(t);h(i)=trapz(t(

4、1:i),v(1:i);endpause,plot(t,h'),grid,实验结果：t=60s时，v=267.2612407732609m/sa=0.914984734734975m/s2H=12.18976913272247km后一段上升高度为0.9255331614101061km总上升高度为13111530229413258km显然此时的加速度为-9.8m/s2高度图像如下：速度图像如下:加速度图像如下:实验结果分析与讨论：从图上可知，60s之后，速度、加速度发生突变，速度在较短的时间里变为0。6.题目：一只小船渡过宽为d的河流，目标是起点A正对着的另一岸B点。已知河水流速v1与

5、船在静水中的速度v2之比为ko(图略)(1) 建立描述小船航线的数学模型，求其解析解；(2) 设d=100m,v1=1m/s,v2=2m/s,用数值解求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。(3) 若流速v1=0,0.5,1.5,2(m/s),结果将如何。模型：以水流方向为x轴，A到B的方向为y轴，小船在任意时刻的坐标为(x,y,t);由图容易得到速度与位移之间的关系：dydtv2*y初始条件为(0,.一dx两式相除得一dy令p0,"k则之=包*丫p,上式化为dp*yp=kjp_1pyv2dydydydp推出：ln(、p21+p)=kln(cy)推出：p0

6、.5(cy)k0.5(cy)k故有-0.5(cy)k0.5(cy)kyk1kk1k有x0.5cy0.5cy代入初值得c计算方法：首先编写船过河的微分方程，存入boat.mfunctiondx=boat(t,x)d=100;v1=1;v2=2;s=sqrt(x(1)A2+(d-x(2)A2);dx=-v2*x(1)/s+v1;v2*(d-x(2)/s;再编写解析解方程：functionx=boat2(y,k)x=0.5*(-0.01)A(-k)*y.A(1-k)-0.5*(-0.01)Ak*y.A(1+k);ts=0:0.01:100;x0=0,-100;option=odeset('r

7、eltol',1e-6,'abstol',1e-9);t,x=ode15s(boat,ts,x0,option,1,2);%这里需使用解刚性方程得龙格一库塔公式计算plot(t,x),gridgtext('x(t)')gtext('y(t)')pauseplot(x(:,1),x(:,2),grid,y=0:-0.01:-100;k=0.5;x=boat2(y,0.5);pause,plot(x,y),grid;计算结果：数值解x,y-t图像：200-20-40-60XQ10。1一1111II11i|1i/111|1n1ii|l>I

8、I/NII1jT|l1IIici*ndH111|1114II1P1yl!PII|iiiJ1111iii4ii/-JIII1l!111!1nJj1ii110W203040506070解析解x,y-t图像:x-y图像:当t=66.66s时，x=0.0069734m,y=-1.9451e-006m,可认为小船已经达到目的地B点。结果分析与讨论：可以发祥数值解、解析解的差距不是很大。将v1改为0,0.5,1.5,2m/s,v1=0,t=50s;10090:h70bi50403020100L_J1-1ri-1J1JJL-0.8-06-C4-02()0.204D6OSv1=0.5,t=53.5s;v1=1

9、.5,t=114.5s;120100SOGO40201H1111-11-1,111MMM-:1''1；11。152030v1=2,小船不可能到达B点。9.问题：两种群相互竞争模型如下：?y、x(t)r1x(1x/n1s1一),y(t)3(1S2x/ni),“其中x,y分别为甲乙两种群的数量，r1,r2为它们的固有增长率，n1,n2为它们的最大容量。si的含义是，对于供养甲的资源而言，单位数量乙的消耗为单位数量甲消耗的s1倍，S2对可作相应的解释。该模型无解析解，试用数值解法研究以下问题：(1)设r1=r2=1,n1=n2=100,s1=0.5,s2=2,x0=y0=10计算x(

10、t),y(t),画出它们的图形及相图，说明时间t充分大以后x和y的变化趋势(2)改变r1,r2,n1,n2,x0,y0,但s1,s2不变，计算并分析所得结果；若s1=1.5(>1),s2=0.7(<1),再分析结果。由此你能得到什么结论，请用各参数生态学上的含义作出解释。(3)试验当s1=0.8(<1),s2=0.7(<1)时会出现什么结果；当s1=1.5(>1),s2=1.7(>1)时又会有什么结果。能解释这些结果吗？模型及实验结果：function微分方程题目已给出，故直接列方程：dx=grow(t,x,r1,r2,n1,n2,s1,s2)r1=1;r2

11、=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)额量时间12010080604020ppipii11姻1_1Lyif/i/_11111111L11111111-L11111111JI1f.1/1IIfa1111ii1-1r-一1ihh11111iJLJJi-2005101620253035M4660甲乙数量关系252015105Q巧102030405060703090100110当t充分大时，x趋于100,y趋于0.实验分析与讨论：(2)当改变除了si和s2

12、以外的参数时，考虑最终结果的变化。functiondx=grow(t,x,r1,r2,n1,n2,s1,s2)r1=1;r2=1;n1=50;n2=100;s1=0.5;s2=5;dx=r1*x(1)*(1-x(1)/n1-s1*x(2)/n2);r2*x(2)*(1-s2*x(1)/n1-x(2)/n2)图形如下:6050403020100-10051016202S30354650冢I里时间甲乙数量关系10864.20j1-I21i1015202530354045SOS5可见，只要s1<1,s2>1,甲的数量最终趋于n1,而乙的数量趋于0.改变si和s2,如下：functiond

13、x=grow(t,x,r1,r2,n1,n2,s1,s2)r1=1;r2=1;n1=100;n2=100;s1=1.5;s2=0.7;dx=r1*x(1)*(1-x(1)/n1-s1*x(2)/n2);r2*x(2)*(1-s2*x(1)/n1-x(2)/n2)120100ao6040200-2005101620253035404650甲乙数量关系对比得出：只要si或者s2大于1,最终数量就会趋于n1或者n2.对比以上数据，我们可以得出结论：只要乙的资源消耗比不过甲，最终甲就会达到最大容量,而乙趋于0.反之亦然。(3)当si和s2均小于1时，结果如下：functiondx=grow(t,x,r

14、1,r2,n1,n2,s1,s2)r1=1;r2=1;n1=100;n2=100;s1=0.8;s2=0.7;dx=r1*x(1)*(1-x(1)/n1-s1*x(2)/n2);r2*x(2)*(1-s2*x(1)/n1-x(2)/n2)TOGO5040302010051016202S30英404550冢I里时间甲乙数量关系TO605040302010101520253035404550S5当s1,s2都大于1时:functiondx=grow(t,x,r1,r2,n1,n2,s1,s2)r1=1;r2=1;n1=100;n2=100;s1=1.5;s2=1.7;dx=r1*x(1)*(1-x(1)/n1-s1*x(2)/n2);r2*x(2)*(1-s2*x(1)/n1-x(2