数学建模第4章

1、4.1 欧拉方法和龙格库塔方法求微分方程数值解，画出解的图形，对结果进行比较分析。以下方程供选择：d) 幂级数解解：原方程化为记，并用t代替x，则原方程化为：；且；于是可以用龙格库塔法求解。【模型求解】：用Matlab作龙格库塔法求解：% chapter 41.d%此函数是微分方程组function Xdot=ch41dfun(t,x)Xdot=x(2),-cos(t)*x(1)'%ch41d.mfunction I=ch41d(a)x0=1,0't,x=ode45('ch41dfun',0,a,x0);y=x(:,1);plot(t,y,'r'

2、);hold on;gtext('y');y1=1-1/jiecheng(2)*t.2+2/jiecheng(4)*t.4-9/jiecheng(6)*t.6+55/jiecheng(8)*t.8;plot(t,y1,'b');gtext('y1');t,y,y1hold off;运行程序可以得到 图1>> ch41d(10) 图2>> ch41d(15) 图3>> ch41d(20) 图4【结果分析】：由图1得：y(龙格库塔方法)和y1(级数近似解)在0到大约1.5的区间内是完全吻合的，从x=1.5之后，两条

3、曲线开始分离。之后y的变化趋势可见图2图4，呈幅度越来越大的上下震荡。4.2 小型火箭初始重量为1 400kg，其中包括1 080kg燃料。火箭竖直向上发射时燃料的燃烧率为18kg/s，由此产生32 000N的推力，火箭引擎在燃料用尽时关闭。设火箭上升时的空气阻力正比于速度的平方，比例系数为0.4kg/m。求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。解：解：先解出速度的微分方程如下：clearsyms a m t vm = 1400-18*t;a = (32000-9.8*m-0.4*v2)/m a =(18280+882

4、/5*t-2/5*v2)/(1400-18*t) 利用龙哥库塔算法画出从0时刻到60秒（燃料用尽）的图像odefun=(t,v)(-(882*t)/5 - (2*v2)/5 + 18280)/(18*t - 1400);t,v=ode45(odefun,0:0.1:60,0);subplot(2,2,1)plot(t,v);grid ona=diff(v)/0.1;t2 = 0:0.1:59.9;subplot(2,2,2)plot(t2,a);grid ons = cumsum(v).*0.1;subplot(2,2,3)plot(t,s);grid on 所以v(60) = 267.3m/

5、s a(60) = 0.9127m/s2 s(60) = 12200m4.3某容器盛满水后，底端直径为d0的小孔开启（如图1），根据水力学知识，当水面高度为h时，谁从小孔中流出的速度为v=0.6*（g*h）0.5（其中g为重力加速度，0.6问哦小孔收缩系数）1） 若容器为倒圆锥形（如图1），现测得容器高和上底直径都为1.2m，小孔直径d为3cm，为水从小孔中流完需要多少时间；2min时水面高度是多少。2） 若容器为倒葫芦形（如图2），现测得容器高1.2m，小孔直径d为3cm，由底端（记x=0）向上每隔0.1m测出容器的直径D（m）如表1所示，问水从小孔中流完需要多少时间；2min时水面高度是多

6、少。 图1X/m00.10.20.30.40.50.60.70.80.91.01.11.2D/m0.030.050.080.140.190.330.450.680.981.101.201.131.00 表1【模型分析】 由题知，水从小孔中流出，不仅与容器有关，还与水流速度v=0.6*（2*g*h）0.5有关。 第一小题容器是圆锥形，比较规则，但是由于水不断从小孔流出，容器中水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决。由（1）知，水面直径等于水深。水深为h时，流量为0.6(/4)d2*(gh)0.5，0.6*(g*h)(0.5)*(d0/2)2*dt=/4*h2*dh则水

7、深下降dh所需时间 ：dt=-(/4)h2*dh/0.6(/4)d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5水深由1.2m至0定积分得水从小孔流完的时间：T（其中已知d=0.03m，g=9.8m*s(-2） 对于第二问：设两分钟（120S）后水深为X m ，由 dt=-(/4)h2*dh/0.6*(/4)*d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5 则263.93-120 =X2.5/1.5*d2*(g)0.5 以d=0.03m，g=9.8m*s(-2代入上式得 水深：X 第二小题容器为倒葫芦形，比较不规则，比较复杂，不仅要考虑水不断从小孔流出，容器中

8、水的高度是不断变化的，水流速度没有一定的公式，所以要用到微积分解决，还要注意表1的倒葫芦形的不断变化，水深的高度变化是不规则的但仍可以用微积分。 由（2）知容器高1.2m，水深为h时，流量为0.6(/4)d2*(gh)0.5，由于不同高度，倒葫芦形半径不同，用欧拉方程和龙格库塔方法则水深下降dh所需时间 ：dt=t(k+1)-t(k)=-(/4)h2*dh/0.6(/4)d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5然后利用循环for k=1：length(L)，t(k)=(h(k+1)-h(k)*(/4)*d(k)2)/(0.6*(/4)*d2*(g(1.2-h(k)0.5

9、)，T=sum(t).可以求得水从小孔流完的总时间。 对于第二问：设两分钟（120S）后水深为X m ，由 S=0,利用条件，当120-s<0.0001时s=s+t(k)，x(k) 把d=0.03m，g=9.8m*s(-2)代入上式得 水深：X 【模型求解】 开始条件h=1.2m，d0=1.2m，g=9.8m*s(-2)，d=0.03m 水深为h时，流量Q为0.6(/4)d2*(gh)0.5，则水深下降dh所需时间 ：dt=-(/4)h2*dh/0.6(/4)d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5水深由1.2m至0定积分得水从小孔流完的时间：T（其中已知d=0.

10、03m，g-9.8m/s） 设两分钟（120S）后水深为X m ，由 dt=-(/4)h2*dh/0.6*(/4)*d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5 则263.93-120 =X2.5/1.5*d2*(g)0.5 以d=0.03，g=9.8代入上式得 水深：X 用欧拉方程和龙格库塔方法则水深下降dh所需时间 ：dt=t(k+1)-t(k)=-(/4)h2*dh/0.6(/4)d2*(gh)0.5=-h1.5*dh/0.6d2*(g)0.5 k1=0.15*sqrt(g*(x(n)*d2/(-43.6359*x(n)8+213.0457*x(n)7-414.873

11、*x(n)6+410.2075*x(n)5-218.8936*x(n)4+62.553*x(n)3-8.3215*x(n)2+0.49619*x(n)+.014892)2; k2=0.15*sqrt(g*(x(n)-h*k1)*d2/(-43.6359*(x(n)-h*k1)8+213.0457*(x(n)-h*k1)7-414.873*(x(n)-h*k1)6+410.2075*(x(n)-h*k1)5-218.8936*(x(n)-h*k1)4+62.553*(x(n)-h*k1)3-8.3215*(x(n)-h*k1)2+0.49619*(x(n)-h*k1)+.014892)2;x(n

12、+1)=x(n)-h*(k1+k2)/2;【附录】(1)g=9.8;d=0.03;syms hy=-h(1.5)/(0.6*d2*sqrt(2*g);T=int(y,h,1.2,0);eval(T) x=(T-120)*(1.5*d2*sqrt(2*g)(0.4);eval(x)运行结果如下>> sy1ans = 263.9316ans =0.9416(2)clear;g=9.8;d=0.03;k1=0;k2=0;h=0.4;x(1)=1.2;for n=1:1000 k1=0.15*sqrt(g*(x(n)*d2/(-43.6359*x(n)8+213.0457*x(n)7-41

13、4.873*x(n)6+410.2075*x(n)5-218.8936*x(n)4+62.553*x(n)3-8.3215*x(n)2+0.49619*x(n)+.014892)2; k2=0.15*sqrt(g*(x(n)-h*k1)*d2/(-43.6359*(x(n)-h*k1)8+213.0457*(x(n)-h*k1)7-414.873*(x(n)-h*k1)6+410.2075*(x(n)-h*k1)5-218.8936*(x(n)-h*k1)4+62.553*(x(n)-h*k1)3-8.3215*(x(n)-h*k1)2+0.49619*(x(n)-h*k1)+.014892)

14、2; x(n+1)=x(n)-h*(k1+k2)/2;end x(300) t=0:h:1000*h; plot(t,x); axis(0,400,0,1.21);运行结果如下>> sy2ans =1.0278水从小孔流完需要263s，2min时水面高度是0.94m点头绪4.4 放射性废物的处理：将装满放射性废物的圆桶扔到水深300ft的海底，圆桶体积55gal，装满放射性废物时的圆桶重量为527.436lbf，在海中浮力为470.327lbf。此外，下沉时圆桶还要受到海水的阻力，阻力与下沉速度成正比，比例系数为0.08lbf s/ft。大量破坏性实验发现当圆桶速度超过40ft/s

15、时，就会因与海底冲撞而发生破裂。（1）建立解决上述问题的微分方程模型。（2）用数值和解析两种方法求解微分方程，并回答谁赢得了官司解：【模型假设】假设装满放射性废物的圆桶在水中仅受重力、浮力和阻力的作用。设圆桶下沉距离为S，速度为v，由于阻力和下沉速度成正比，比例系数为k=0.08*0.4536*9.8/0.3848。设圆桶质量为m，重力为G，浮力为F，阻力为R，则有R=kv=0.08*0.4536*9.8*v/0.3848。根据牛顿第二定律可得微分方程组对此微分方程组，利用龙格-库塔方法能给出数值解。（一）用数值法求解微分方程【算法设计】有两个思路：思路一是利用龙格-库塔方法直接解上面关于下沉

16、距离S和下沉速度v的微分方程组。思路二是利用龙格-库塔方法解关于速度v和时间t的一阶微分方程，即，再利用辛普森公式对速度求积分已得到下沉的距离S；分别求出S和v与时间的关系后，对比同一时间下的S和v，通过判断下沉距离S=300ft时速度是否大于40ft/s或判断速度达到40ft/s时下沉的距离是否大于300ft（实际上下沉距离不可能大于海深，这里是只是一种数学上的处理），则可以判断圆桶与海底冲撞时是否会破裂，从而判断官司的胜负。思路一：【程序】function dx=waste(t,x) %建立名为waste的函数M文件m=527.436*0.4536;G=527.436*0.4536*9.8

17、; %圆桶所受重力F=470.327*0.4536*9.8; %圆桶所受浮力R=0.08*0.4536*9.8*x(2)/0.3848; %圆桶所受阻力dx=x(2);(G-F-R)/m; %以向量形式表示方程clear allts=0:0.1:25; %终点试探x0=0,0; %初始值opt=odeset; %默认误差限t,x=ode45(waste,ts,x0); %使用5级4阶龙格库塔公式计算t,x %输出t,S(t),v(t)plot(t,x(:,2),grid,title('v(t)-t曲线');%作出v(t)图形，以便考察速度变化趋势xlabel('t

18、9;);ylabel('v(t)');figure;%在新窗口中作图plot(t,x),axis(0 25 0 350),grid, %按照数值输出作S(t),v(t)的图形gtext('S(t)'),gtext('v(t)')hold on,Smax=300*0.3048;plot(0,25,Smax,Smax,'r-'),hold off%作出S=Smax的图形，以便比对hold on,Vmax=40*0.3048;plot(0，25,Vmax,Vmax,'m-'),hold off%作出v=Vmax的图形，以

19、便比对gtext('Smax'),gtext('Vmax')【计算结果】MATLAB输出图像如下（相关分析见后）：图1速度-时间图像图2思路二：【程序】function dx=waste1(t,x,G,F,m,b)dx=(G-F-b*x)/m;%列写微分方程clear all;h=0.1; %步长ts=0:h:2000;x0=0; %初始值G=527.436*0.4536*9.8; %重力F=470.327*0.4536*9.8; %浮力m=527.436*0.4536; %质量b=0.08*0.4536*9.8/0.3048; %比例系数opt=odeset;

20、t,x=ode45(waste1,ts,x0,opt,G,F,m,b);%调用5级4阶龙格-库塔方程plot(t,x),grid %作v(t)曲线xlabel('t');ylabel('v(t)'); %用辛普森公式对速度积分求下沉深度T=25; %估计时间25s沉到海底for i=0:2:10*T y=x(1:(i+1); k=length(y); a1=y(2:2:k-1);s1=sum(a1); a2=y(3:2:k-1);s2=sum(a2); z4(i+2)/2)=(y(1)+y(k)+4*s1+2*s2)*h/3; %辛普森公式endi=0:2:10

21、*T;figure;Smax=300*0.3048;Vmax=40.*0.3048;plot(x(i+1),z4'),grid;hold on,plot(0,25,Smax,Smax,'r:'),grid %作Smaxhold on,plot(Vmax,Vmax,0,350,'m:'),grid,hold off %作Vmaxlegend('S(v)','Smax','Vmax');xlabel('v');ylabel('S(v)'); figure;plot(i,z4

22、9;); %作S（t）曲线grid;xlabel('t');ylabel('S(t)');【输出结果】速度-时间图像图3下沉深度-速度图像：图4从上面的S(v)图像可以清楚地看到S=Smax（海深）时速度已超过了Vmax。可见圆桶下沉到海底时会破裂。下沉深度-时间图像：图5（二）用解析法求解微分方程【算法设计】对于上面的方程组令a=1.06111b=4.87676*10-3则原方程组等价于 （*）v(0)=0一阶线性常微分方程(*)式的特征方程为+b=0解得特征根为=-b=-4.87676*10-3所以（*）式齐次解的形式为v=A*exp(-bt),其中A为待定

23、的常系数（*）式的一个特解v=a/b=217.59所以（*）的通解可写成v=A*exp(-bt)+217.59，其中A为待定的常系数代入初始条件v(0)=0，得A=-217.59所以，v(t)=-217.59*exp(-4.87676*10-3t)+217.59=44618*exp(-4.87676*10-3t)+217.59*t-44618【程序】用MATLAB编写程序如下：clear allt=0:0.1:25;Smax=300*0.3048;Vmax=40*0.3048;v=-217.59*exp(-4.87676*t*0.001)+217.59;s=44618*exp(-4.87676

24、*t*0.001)+217.59*t-44618;plot(t,v,'b-',t,s,'g-'),gridlegend('速度-时间曲线','下沉深度-时间曲线')hold on,plot(0,25,Smax,Smax,'r:'),gridhold on,plot(Vmax,Vmax,0,350,'m:'),grid,hold offgtext('Smax'),gtext('Vmax')【计算结果】输出图形与上面用数值法得到的结果（图2）相同图5【结果的数学分析】从图

25、1速度-时间曲线可以看到，开始时圆桶的速度随着时间逐渐增大，加速度（曲线的斜率）逐渐减小，可见圆桶做的是加速度逐渐减小的变加速运动。可见若没有题中的约束条件（海深S=300ft），则时间充分大后（t>1500s）后，圆桶的加速度减至0，速度趋于一个定值，约为275m/s。根据上面的v(t)-t图像能大致判断圆桶速度达到40ft/s(约为12.2m/s)所用的时间。结合题目中的约束条件重新作图，判断下沉距离S=300ft时速度是否大于40ft/s或判断速度达到40ft/s时下沉的距离是否大于300ft，从而判断圆桶与海底冲撞时是否会破裂。图2和图5中Smax=300*0.3048m,Vma

26、x=40*0.3048m/s。从图中可见，圆桶还没沉到海底时，其速度已达到规定的圆桶不破裂的最大速度。由此可见，当圆桶下沉到海底时，必然会与海底碰撞破裂。显然是工程师们赢了官司！【结果的实际意义】从这道题目可以看到，利用Matlab构造数学模型可以分析和模拟实际生活问题。比如像这道题，利用Matlab可以模拟出圆桶下沉的速度和深度关系与下沉时间的关系曲线，从而得出圆桶在下沉到海底时是否会因与海底碰撞而破裂的结论。可见，适当运用数学模型可以在环境保护上等实际生活的很多地方作出贡献。4.5一只小船渡过宽为d的河流（见下图），目标是起点A正对着的另一岸B点。已知河水流速v1与船在静水中的速度v2之比

27、为k。a)建立小船航线的方程，求其解析解。b)设d=100m,v1=1m/s,v2=2m/s,用数值解法求渡河所需时间、任意时刻小船的位置及航行曲线，作图，并与解析解比较。c)若流速v1为0,0.5,1.5,2(m/s)，结果将如何？BxxydVV2V1CAy【模型建立】：如图，以B为原点，沿河岸向右为x轴正向，垂直河岸向下为y轴正向，建立坐标系。设在t时刻，船在x方向上的位移是x(t)，在y方向上的位移是y(t),则：在t时刻，船在x方向上的速度是，在y方向上的速度是，将船的速度V和静水速度V1在x,y轴方向上分解，可得：，又因为船头始终指向B点，所以，所以：，即这就是本题的微分方程。初始条

28、件为：x(0)=0,y(0)=-d【模型求解】：1.解析解:令，将直角坐标化成极坐标，则原微分方程化成:;化简后为：消去t，得：于是解得：，其中2.数值解法：根据上述微分方程，用龙格库塔方法求解：程序如下：%ch45fun.m%第四章第5题%此函数是微分方程组function Xdot=ch45fun(t,x,V1,V2)%x(1)代表x;x(2)代表y%加入限制条件，防止无限循环if (norm(x)>1e-5) Xdot=V1-V2*x(1)/sqrt(x(1)2+x(2)2),-V2*x(2)/sqrt(x(1)2+x(2)2)'else Xdot=0,0'end%

29、ch45.m%第四章第5题function I=ch45(V1,V2,d,a)x0=0,-d'%初始条件t,x=ode45(ch45fun,0,a,x0,V1,V2);plot(x(:,1),x(:,2),'r');hold on;%作图t,x(:,1),x(:,2)%打印数据%下面是作出精确解的图象seta=linspace(-pi/2,0,100);rou=d*(abs(tan(seta/2).(V2/V1)./sin(seta);xp=-rou.*cos(seta);yp=-rou.*sin(seta);plot(xp,yp,'b');hold o

30、ff;需要说明：在ch45fun.m中，加入了(norm(x)>1e-5)的限制条件，以保证在船离B点足够近时中止运算。否则无法得出正确结果。运行结果：（中间数据省略）可知渡河时航线如图，所用时间为：66.7秒。当v1=0,0.5,1.5,2时，结果分别为（图像）：v1=0说明在静水中，船沿直线到达B点。这与直观经验相符合渡河时间为50秒；v1=0.5渡河时间为：53.3秒v1=1.5渡河时间为：114.3秒v1=2渡河时间:从图上看出，到t=1000秒时，船已到达对岸，但是并没有到达B点，而是在B点下游50米处。由于船头指向B点，即船头指向逆流方向，且船速（静水）等于水速，可知船将保持

31、原地不动。也就是说，船永远到达不了B点。所以渡河时间为无穷大。v1=2.5渡河时间：与v1=2时情况类似，船能到达对岸，但是是在B点下游。由于船速（静水）小于水速，船将被水冲得顺流而下，同样永远到达不了B点。【结果分析】：1.以上图中红线是龙格库塔方法的结果，蓝线是解析解的结果。从图中可以看出，两种方法的结果基本上是相符的。但在接近B点时，解析解将无法得到正确解（因为会出现分母为零的情况，而Matlab不会计算极限）。当V1=0时，解析解也得不到正确解（同样是因为出现了分母为零）。2.随着V1的变化，船的航线也在变化。当V10时，航线为直线；0<V1<2时，航线呈一个类似于抛物线的

32、曲线，V1越大，该“类抛物线”的顶点的横坐标越大，纵坐标也越大。即水速越大，航线的顶点顺流而下的距离越大，离对岸越近；当V1=2时，船会静止在B点下游d/2的地方；V1>2时，船在接近对岸后会顺流而下。这一切都与我们的直观感觉相符合。4.6海防某部缉私艇上的雷达发现正东方向c海里处有一艘走私船正以一定速度行驶，缉私艇立即以最大速度前往拦截。用雷达进行跟踪时，可保持缉私艇的速度方向始终指向走私船。建立任意时刻缉私艇的位置和缉私艇航线的数学模型，确定缉私艇追上走私船的位置，求出追上的时间。解：（1） 只需应用余弦公式代替即可：建立直角坐标系如图，设在时刻缉私艇发现走私船，此时缉私艇的位置在，

33、走私船的位置在。走私船以速度与轴成方向行驶，缉私艇以速度按指向走私船的方向行驶。在任意时刻缉私艇位于点，而走私船到达点，直线PQ与缉私艇航线（图中曲线）相切，切线与轴正向夹角为。缉私艇在方向的速度分别为，由三角形PQR写出表达式，可以同样得到微分方程。   （2）此时可以得到走私船的速度函数。进行分段求解。4.7 用x(t)表示时刻t食饵（如食用鱼）的密度，即一定区域内的数量，y(t)表示捕食者（如鲨鱼）的密度。假设食饵独立生存时的（相对）增长率为常数r>0，即，而捕食者的存在使食饵的增长率减小，设减小量与捕食者密度成正比，比例系数为a>0，则。捕食者

34、离开食饵无法生存，设它独自存在时死亡率为常数d>0，即，而食饵的存在为捕食者提供了食物，使捕食者的死亡率减小，设减小量与食饵密度成正比，比例系数为b>0，则。实际上，当bx>d时捕食者密度将增长。                          给定食饵和捕食者密度的初始值x0, y0，由上可知x(t), y(t)满足以下方程：

35、                                    解x(t), y(t)描述了食饵和捕食者密度随时间的演变过程。由于得不到x(t), y(t)的解析解，需要用数值算法求解。    

36、0;  【 模型分析】1）相轨线y(x)的解析解：由方程组消去后得到，是可分离变量方程                                    两边积分得到y(x) 的通解   

37、0;                                         其中常数c由初始条件确定。可以证明，c在一定范围内时上式表示的曲线是封闭曲线，这等价于是周期函数。2）平均密度：记

38、的周期为，为了求x(t)在一个周期的平均值，将模型中的改写作       两边在周期内积分，容易算出x(t)平均值为                    类似地可得y(t)在一个周期的平均值              

39、             【MATLAB计算实验 】选取模型参数r=1，d=0.5， a=0.1，b=0.02，x0=25，y0=2，时间跨度tf=15。 编写如下M文件：function xdot=shier(t,x,r,d,a,b)           xdot=diag(r-a*x(2),-d+b*x(1)*x;   &#

40、160;           % 以向量形式表示方程     运行主程序：r=1;d=0.5;a=0.1;b=0.02;x0=25,2;tf=15;ts=0:0.1:tf;                       &

41、#160;         t,x=ode45(shier,ts,x0,r,d,a,b);plot(t,x),grid,%gtext('fontsize12x(t)'),gtext('fontsize12y(t)'),  % 将标记x(t) , y(t)的字体放大plot(x(:,1),x(:,2),grid,xlabel('x'),ylabel('y')   输出图形：      4.8两种群相互竞争模型如下：，其中x(t),y(t)分别为甲乙两种群的数量，为它们的固有增长率，为它们的最大容量。的含义是，对于供养甲的资源来说，单位数量的乙（相对）的消耗为单位数量甲（相对）消耗的倍，对可以作相应解释。该模型无解析解，试用数值方法研究以下问题：（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