数学建模章绍辉版第四章作业

第四章作业

第二题：

针对严重的交通情况，国家质量监督检验检疫局发布的国家标准，车辆驾驶人员血液中的酒精含量大于或等于20mg/100ml,小于80mg/100ml 为饮酒驾车，血液中的酒精含量大于或等于80mg/100ml 的为醉酒驾车。

下面分别考虑大李在很短时间内和较长时间内（如2个小时）喝了三瓶啤酒，多长时间内驾车就会违反新的国家标准。

1、 问题假设

大李在短时间内喝下三瓶啤酒后，酒精先从吸收室（肠胃）吸收进中心室（血液和体液），然后从中心室向体外排除，忽略喝酒的时间，根据生理学知识，假设 （1） 吸收室在初始时刻t=0时，酒精量立即为

32

D ;在任意时刻，酒精从吸收室吸收进中心室的速率（吸收室在单位时间内酒精含量的减少量）与吸收室的酒精含量成正比，比例系数为1k ；

（2） 中心室的容积V 保持不变；在初始时刻t=0时，中心室的酒精含量为0；在任意时

刻，酒精从中心室向体外排除的速率（中心室在单位时间内酒精含量的减少量）与

中心室的酒精含量成正比，比例系数为2k ；

（3） 在大李适度饮酒没有酒精中毒的前提下，假设1k 和2k 都是常量，与饮酒量无关。

2、 |

3、 符号说明

酒精量是指纯酒精的质量，单位是毫克；

酒精含量是指纯酒精的浓度，单位是毫克/百毫升；

~t 时刻（小时）

； ()~x t 在时刻t 吸收室（肠胃）内的酒精量（毫克）

； 0~D 两瓶酒的酒精量（毫克）；

(t)~c 在时刻t 吸收室（血液和体液）的酒精含量（毫克/百毫升）

； 2()~c t 在时刻t 中心室（血液和体液）的酒精含量（毫克/百毫升）；

~V 中心室的容积（百毫升）

； 1~k 酒精从吸收室吸收进中心室的速率系数（假设其为常数）；

@

2~k 酒精从中心室向体外排除的速率系数（假设其为常数）；

3~k 在短时间喝下三瓶酒的假设下是指短时间喝下的三瓶酒的酒精总量除以中心室体积，

即03/2D V ；而在较长时间内（2小时内）喝下三瓶酒的假设下就特指03/4D V .

4、 模型建立和求解

（1） 酒是在很短时间内喝的：

记喝酒时刻为0t =（小时），设(0)0c =，可用()2113

212

()k t k t k k c t e e k k --=

--来计算血液中的酒精含量，此时12k k 、为假设中所示的常数，而033155.792D k V ??

==

???

. 下面用MATLAB 程序画图展示血液中酒精含量随时间变化并且利用fzero 函数和fminbnd 函数来得到饮酒驾车醉酒驾车对应的时间段，以及血液中酒精含量最高的时刻。

MATLAB 程序如下：

k1=;k2=;k3=; *

c=@(t)(k1.*k3)./(k1-k2).*(exp(-k2.*t)-exp(-k1.*t)); f=@(t)c(t)-20; g=@(t)c(t)-80; h=@(t)-c(t);

t1(1)=fzero(f,1);t1(2)=fzero(f,12), t2(1)=fzero(g,1);t2(2)=fzero(g,12) [t3,c3]=fminbnd(h,0,24) fplot(c,[0,20],'k') hold on

plot([0,20],[20,20],'k',[0,20],[80,80],'k') \

hold off

xlabel('时刻t （小时），从开始喝酒算起') ylabel('血液中的酒精含量（mg/100ml ）')

title('短时间喝下三瓶酒时，血液中酒精含量随时间的变化过程') gtext(',20)') gtext(',20)') gtext(',80)') gtext(',80)') gtext(',')

）

运行结果如下：

t1 =

t2 =

c3 =

-

所绘图形如下：

024

68101214161820

2040

60

80100

120

140

时刻t （小时），从开始喝酒算起

血液中的酒精含量（m g /100m l ）

短时间喝下三瓶酒时，血液中酒精含量随时间的变化过程

结果分析：

所以，当[0.06891,0.38052)(4.1125,11.589]t ∈?时，20()80c t ≤<，属饮酒驾车。当[0.38052,4.1125]t ∈时，属醉酒驾驶；当 1.307t =时，血液中的酒精含量最高为毫克/百毫升。

（2） 酒是在2小时内喝的：

可假设三瓶啤酒是在2小时内匀速喝的. 同样记喝酒时刻为0t =（小时），设(0)0c =，则吸收室的酒精量1()x t 满足分段的初值问题

101

11120111213, (0)0,0243, (2)(1),2

4k

D dx k x x t dt

D dx k x x e t dt

k -?=-+=≤≤????=-=-≥??

11101120

3(1) 0t 24

()3(1) t 24

k t

k k t D e k x t D e e --?-≤≤??=??-≥??

】

于是中心室内的酒精含量2()c t 满足分段的初值问题

1122232222728

(1) (0)0,02 (2),2

k t

k t dc k c k e c t dt dc k c k e c k t dt

--?=-+-=≤≤????=-+=≥?? 解得

1221456210

9 0t 2() 2k t k t

k t k t

k e k e k c t k e k e t ----?-+≤≤?=?-≥?? 其中

334D k V

=

，3412k k k k =

-，14

52

k k k k =，36

2= k k k ，1273(1)k k k e =- 1

2

228456k k k k e k e

k --=-+，7

912

k k k k =

-，22122()10

89k k k k k e k e -=+

因为1 2.0079k =，20.1855k =以及0/103.86D V =，所以

377.896k =，442.743k =，5462.66k =，6419.92k = 74243.1k =，8101.43k =，92328.3k =，10207.82k =

@

下面用MATLAB 程序画图展示血液中酒精含量随时间变化并且利用fzero 函数和fminbnd 函数来得到饮酒驾车醉酒驾车对应的时间段，以及血液中酒精含量最高的时刻。

MATLAB 程序如下：

k1=;k2=;k3=;

k4=;k5=;k6=;k9=;k10=;

c1=@(t)(k4.* exp(-k1.*t)-k5.*exp(-k2.*t)+k6).*(t>=0&t<=2)+... ( k10.* exp(-k2.*t)-k9.*exp(-k1.*t)).*(t>2); f1=@(t)c1(t)-20; g1=@(t)c1(t)-80; h1=@(t)-c1(t); ，

t1(1)=fzero(f1,1);t1(2)=fzero(f1,12),

t2(1)=fzero(g1,1);t2(2)=fzero(g1,12),

[t3,c3]=fminbnd(h1,0,20)

fplot(c1,[0,20],'k')

hold on

plot([0,20],[20,20],'k',[0,20],[80,80],'k')

hold off

xlabel('时刻t（小时），从开始喝酒算起')

ylabel('血液中的酒精含量（mg/100ml）')

title('短时间喝下三瓶酒时，血液中酒精含量随时间的变化过程') ~

gtext(',20)')

gtext(',20)')

gtext(',80)')

gtext(',80)')

gtext(',')

运行结果如下：

t1 =

t2 =

·

t3 =

c3 =

所绘图形如下：

024

68101214161820

20

40

60

80

100

120

时刻t （小时），从开始喝酒算起

血液中的酒精含量（m g /100m l ）

短时间喝下三瓶酒时，血液中酒精含量随时间的变化过程

结果分析：

所以，当[0.62321,1.6366)(5.1412,12.62]t ∈?时，20()80c t ≤<，属饮酒驾车。当

[1.6366,5.1412]t ∈时，属醉酒驾驶；当 2.6328t =时，血液中的酒精含量最高，为毫克/

百毫升. @

下面用图形比较两种不同假设下血液中酒精含量的变化过程。

MATLAB 程序如下：

k1=;k2=;k3=;

k4=;k5=;k6=;k9=;k10=;

c=@(t)(k1.*k3)./(k1-k2).*(exp(-k2.*t)-exp(-k1.*t));

c1=@(t)(k4.* exp(-k1.*t)-k5.*exp(-k2.*t)+k6).*(t>=0&t<=2)+... ( k10.* exp(-k2.*t)-k9.*exp(-k1.*t)).*(t>2);

plot(0::20,c(0::20),'--k',0::20,c1(0::20),'k',2,c1(2),'.k') xlabel('时刻t （小时），从开始喝酒算起') ylabel('血液中的酒精含量（mg/100ml ）') $

title('短时间喝下三瓶酒时，血液中酒精含量随时间的变化过程')

legend('很短时间内喝三瓶啤酒','两小时内匀速喝下三瓶啤酒','函数的分段点')

所绘图形如下：

020

40

60

80

100

120

140

时刻t （小时），从开始喝酒算起

血液中的酒精含量（m g /100m l ）

短时间喝下三瓶酒时，血液中酒精含量随时间的变化过程

第四题：

研究将鹿群放入草场后，草和鹿两个种群的相互作用，草的生长服从Logistic 规律，年固有增长率，最大密度为3000个密度单位，在草最茂盛时，每只鹿每年吃掉个密度单位的草，若没有草，鹿群的年死亡率高达,而在草最茂盛的时候草对鹿的死亡的补偿率为。

1、建立差分方程组模型，比较将100只鹿放入密度为1000和密度为3000的两种草场的情况下，草和鹿两个种群的数量演变过程。

`

（1）符号说明：

x ~k 第k 年草场的密度单位 ~k y 第k 年草场上鹿的数量

r ~草场上草的年固有增长率

a ~由于捕食导致的草的密度单位减少的速度大小

d ~如果没有草，鹿群的年死亡率 b ~草对鹿群死亡的补偿率 N ~草的最大密度单位

（2）模型的建立与求解： —

基于以上假设，由于草的生长服从Logistic 模型，建立差分方程组模型如下所示：

11(1)k k k k k k k k

k k k x ax y x x rx N N

bx y y y dy N

++=+-

-=-+

令1k k x x x +==，1k k y y y +==，与上述方程组联立得到平衡点为0(0,0)P 、1(,0)P N 、

2()(

,)dN rN b d P b ab

- 以下用MATLAB 实现差分方程组模型。

MATLAB 程序如下： n=20;r=;a=;b=;d=; N=3000;t=1:n+1;

x1(1)=1000;y1(1)=100; for k=1:n

x1(k+1)=x1(k)+r*x1(k)*(1-x1(k)/N)-a*x1(k)*y1(k)/N y1(k+1)=y1(k)-d*y1(k)+b*x1(k)*y1(k)/N end

subplot(2,2,1),plot(t,x1,'k^',t,y1,'kv'),axis([-1,21,0,3000]), xlabel('第k 年'),ylabel('数量'),gtext('草密度'),gtext('鹿数量'), title('草和鹿随时间的演变');

subplot(2,2,2),plot(x1,y1,'ko'),axis([1000,3000,0,1000]), xlabel('草密度'),ylabel('鹿数量'),title('相平面图'); x2(1)=3000;y2(1)=100; for k=1:n

x2(k+1)=x2(k)+r*x2(k)*(1-x2(k)/N)-a*x2(k)*y2(k)/N y2(k+1)=y2(k)-d*y2(k)+b*x2(k)*y2(k)/N end

subplot(2,2,3),plot(t,x2,'k^',t,y2,'kv'),axis([-1,21,0,3000]), xlabel('第k 年'),ylabel('数量'),gtext('草密度'),gtext('鹿数量'), title('草和鹿随时间的演变');

subplot(2,2,4),plot(x2,y2,'ko'),axis([1000,3000,0,1000]), xlabel('草密度'),ylabel('鹿数量'),title('相平面图');

运行结果如下：

x1 = ~ Columns 1 through 8

1000 1480

Columns 9 through 16 1550

Columns 17 through 21

y1 =

Columns 1 through 8

100 60

Columns 9 through 16

\

Columns 17 through 21

x2 =

Columns 1 through 8

3000 2840 2570 2149 Columns 9 through 16

1698 1907

Columns 17 through 21

）

y2 =

Columns 1 through 8

100 160

Columns 9 through 16

Columns 17 through 21

{

所绘图像如下：

1000

2000

3000

第k 年

数量

草和鹿随

时间的演变

1000

1500

200025003000

草密度鹿数量

相平面图

5

1015

20

1000

2000

3000第k 年

数量

1000

1500

200025003000

0500

1000

草密度

鹿数量

相平面图

由图像中可以看出2()(

,)dN rN b d P b ab

-为渐进稳定的平衡点。

|

2、建立常微分方程组模型，重做以上问题。

记草和鹿在第t 年的数量分别记为()x x t =和()y y t =，其余的符号假设与（1）相同，建立常微分方程组模型为：

(1)dx x axy rx dt N N

dy bxy dy dt N

=--=-+

令(1)0x axy rx N N -

-=且0

k bxy

dy N -+=，得到常微分方程组的临界点为0(0,0)P 、1(,0)P N 、2

()

(,)dN rN b d P b ab

-，以下运用数值解，方向场和特征线等技巧，用matlab 绘制常微分方程组的解曲线图，并加以分析。

为了消去方程中的参数t ，将两式相除，得到：

()dy Ndy bxy

dx rx N x axy

-+=-- 当将100头鹿放入密度为1000与密度为3000的草场中，初始值分别为(0)1000x =，

(0)100y =与(0)3000x =，(0)100y =.

利用matlab实现上述过程的源程序与运行结果如下：函数m文件：

function dx=fun(t,x)

dx=zeros(2,1);

\

dx(1)=x(1)**(1-x(1)/3000)*x(1)*x(2)/3000;

dx(2)=*x(2)+*x(1)*x(2)/3000;

主程序：

[t1,x1]=ode45('fun',0:20,[1000 100]);

[t1,x1]

subplot(2,2,1)

plot(t1,x1(:,1),'.-k',t1,x1(:,2),'.-k'),grid on

axis([0 20 0 3000])

title('草和鹿随时间的演变（x_0=1000）')

gtext('草密度'),gtext('鹿数量')

:

xlabel('第k年'),ylabel('数量')

subplot(2,2,2)

plot(x1(:,1),x1(:,2),'.-k'),grid on

axis([1000 3000 0 800])

title('相平面图（x_0=1000）')

xlabel('草密度'),ylabel('鹿数量')

[t2,x2]=ode45('fun',0:20,[3000 100]);

[t2,x2]

subplot(2,2,3)

plot(t2,x2(:,1),'.-k',t2,x2(:,2),'.-k'),grid on

?

axis([0 20 0 3000])

title('草和鹿随时间的演变（x_0=3000）')

gtext('草密度'),gtext('鹿数量')

xlabel('第k年'),ylabel('数量')

subplot(2,2,4)

plot(x2(:,1),x2(:,2),'.-k'),grid on

axis([1000 3000 0 800])

title('相平面图（x_0=3000）')

xlabel('草密度'),ylabel('鹿数量')

》

@

所绘图像如下：

5

10

1520

草和鹿随时间的演变（x 0=1000）第k 年

数量

1000

1500

200025003000相平面图（x 0=1000）

草密度

鹿数量

05

101520

草和鹿随时间的演变（x 0=3000）第k 年

数量

1000

1500

200025003000

020*******

800相平面图（x 0=3000）草密度

鹿数量

通过图像发现该常微分方程组的平衡点2()(

,)dN rN b d P b ab

-是渐进稳定的，而0(0,0)P 、1(,0)P N 不是渐进稳定的平衡点。

草与鹿群的数量演变过程如下表所示：

h 时用欧拉方法计算微分方程组的数值解。差分方程相当于取步长1