章绍辉数学建模第四章

第四章

2.

符号说明：

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

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

t ～时刻（小时）；

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

D0～在段时间内喝下2瓶啤酒后吸收室内的酒精量（毫克）；

c1(t) ～在时刻t吸收室（肠胃）的酒精含量（毫克/百毫升）；

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

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

k1～酒精从吸收室吸收进入中心室的速率系数；

k2～酒精从中心室向体外排除的速率系数；

k3～假如所喝入的酒精完全吸收进中心室却没有排除出体外，中心室的酒精含量将达到的最大值（毫克/百毫升）；

模型假设：

假设一：酒是在很短时间内喝的

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

（1）吸收室在初始时刻t=0时，酒精量立即为3D0/2；在任意时刻，酒精从吸收室吸收进入中心室的速率（吸收室在单位时间内酒精含量的减少量）与吸收室的酒精含量成正比，比例系数为k1；

（2）中心室的容积V保持不变；在初始时刻t=0时，中心室的酒精含量为0；在任意时刻，酒精从中心室向体外排除的速率（中心室在单位时间内酒精含量的减少量）与中心室的酒精含量成正比，比例系数为k2；

（3）在大李适度饮酒没有酒精中毒的前提下，假设k1和k2都是常数，与饮酒量无关；

假设二：酒是在较长一段时间（如2小时）内喝的

大李在较长一段时间（如2小时）内喝下三瓶啤酒后，酒精先从吸收室（肠胃）吸收进中心室（血液和体液），然后从中心室向体外排除，忽略喝酒的时间，根据生理学知识，假设

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

（2）中心室的容积V 保持不变；在初始时刻t=0时，中心室的酒精含量为0；在任意时刻，酒精从中心室向体外排除的速率（中心室在单位时间内酒精含量的减少量）与中心室的酒精含量成正比，比例系数为k 2；

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

模型建立和求解：

假设一：酒是在很短时间内喝的 根据4.1.7小节饮酒驾车案例，假设

k 1=2.0079，k 2=0.1855

则在很短时间内喝下三瓶酒时

0333

103.86155.7922

D k V =

=?= 记喝酒时间为t=0时，c 2(t)=0，则可根据

(

)

21

--13212

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

来计算t 时刻血液中的酒精含量c 2(t)，Matlab 代码如下： M 文件fun2_1.m

k1=2.0079;k2=0.1855;k3=155.79;

c2=@(t)(k1.*k3)./(k1-k2).*(exp(-k2.*t)-exp(-k1.*t)); hold on ;

plot(0:0.01:20,c2(0:0.01:20)); plot([0,20],[20,20],[0,20],[80,80]); xlabel('时刻t （小时）');

ylabel('血液中酒精含量c2（毫克/百毫升）');

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

输出图像

由图像可知c 2(t)的函数图像与y=20和y=80分别都有2个交点，故可继续添加代码求出这4个交点坐标，代码如下： M 文件fun2_1.m 续

f=@(t)c2(t)-20; g=@(t)c2(t)-80;

ft=[fzero(f,1),fzero(f,20)] gt=[fzero(g,1),fzero(g,20)]

输出结果

ft =

0.0689 11.5887 gt =

0.3805 4.1125

即

（1）当[0.0689,0.3805]

[4.1125,11.5887]t ∈时，220()80c t ≤<，属于饮酒驾车；

（2）当[0.3805,4.1125]t ∈时，2()80c t ≥，属于醉酒驾驶； 假设二：酒是在较长一段时间（如2小时）内喝的 根据4.1.7小节饮酒驾车案例，假设

k 1=2.0079，k 2=0.1855

则在较长一段时间（如2小时）内喝下三瓶酒时

0333

103.8677.89544

D k V =

=?= 记喝酒时间为t=0时，c 1(t)=0，c 2(t)=0，则吸收室的酒精量x 1(t)满足分段初值问题

110

11

1210

11

11

3,(0)0,0t 2 43

,(2)(1),2

4k dx D k x x dt

dx D k x x e k dt

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

?=-=-≥??

解得

11101120

3

(1) ,024

()3(e 1)e , 2 4k t k k t D e t k x t D t --?-≤≤??

=??-≥?? 故中心室内的酒精含量c 2(t)满足分段初值问题

1112

22322222321(1) ,(0)0,02(1),(2), 2 k t k k t dc k c k e c t dt

dc k c k e e c p t dt

--?=-+-=≤≤???

?=-+-=≥?? 解得

121213133212

2122223212

,0t 2

()(1) , 2

k t k t

k k t k t k k k k e e k k k k k k c t k e p e e t k k ----?-+≤≤?

--?=?

-?-≥?-?

其中

1

2

22313312

12

2

122k k

k k k k p e e k k k k k k --=

-

+-- 1

2

21222()

32112

(1)k k k

k k e p p e

e

k k --=+-

通过Matlab 画出c 2(t)的图像，代码如下： M 文件fun2_2.m

k1=2.0079;k2=0.1855;k31=155.79;k3=k31/2; c2=@(t)(k1.*k31)./(k1-k2).*(exp(-k2.*t)-exp(-k1.*t)); hold on ;

plot(0:0.01:20,c2(0:0.01:20),'--');

plot([0,20],[20,20],[0,20],[80,80]);

xlabel('时刻t（小时）');

ylabel('血液中酒精含量c2（毫克/百毫升）');

p1=k3*exp(-2*k1)/(k1-k2)-k1*k3*exp(-2*k2)/(k1*k2-k2*k2)+k3/k2;

p2=p1*exp(2*k2)+k3*(exp(2*k1)-1)*exp(2*(k2-k1))/(k1-k2);

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

c22=@(t)p2.*exp(-k2.*t)-k3.*(exp(2.*k1)-1).*exp(-k1.*t)./(k1-k2);

plot(0:0.01:2,c21(0:0.01:2));

plot(2:0.01:20,c22(2:0.01:20));

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

legend('很短时间内喝的','较长一段时间（如2小时）内喝的');

输出图像

由图像可知c2(t)的函数图像与y=20和y=80分别都有2个交点，故可继续添加代码求出这4个交点坐标，代码如下：

M文件fun2_2.m续

f=@(t)c2(t)-20;

g=@(t)c2(t)-80;

ft=[fzero(f,1),fzero(f,20)]

gt=[fzero(g,1),fzero(g,20)]

输出结果

ft =

0.6233 12.6196 gt =

1.6366 5.1412

即

（1）当[0.6233,1.6366]

[5.1412,12.6196]t ∈时，220()80c t ≤<，属于饮酒驾车；

（2）当[1.6366,5.1412]t ∈时，2()80c t ≥，属于醉酒驾驶。 3.

考虑3.4.2小节“酵母培养物的增长”案例，建立微分方程模型，模拟酵母培养物的增长，故假设t （小时）时刻的酵母生物量为x(t)（克），且t=0时x(t)=x 0，则x(t)满足微分方程初值问题

()dx

t x dt

ρ=，0(0)x x = 其中()t ρ为酵母培养物的增长量，根据3.4.2小节可知

()(1()/N)t r x t ρ=-

其中r>0称为酵母培养物固有增长率，N>0称为酵母培养物的最大容量，则得到阻滞增长方程

(1)dx x

rx dt N

=- 用分离变量法得

00()()rt

Nx x t x N x e

-=

+- 则根据表3.2通过Matlab 拟合出x(t)的参数r 、N 和x =，代码如下 M 文件fun3_1.m

t=0:18;

x=[9.6,18.3,29,47.2,71.1,119.1,174.6,257.3,350.7,441,513.3,559.7,... 594.8,629.4,640.8,651.1,655.9,659.6,661.8]; f=@(b,t)b(2).*b(3)./(b(3)+(b(2)-b(3)).*exp(-b(1).*t)); [b1,r1]=nlinfit(t(1:19),x(1:19),f,[0.5,660,9.6]) s=sum(r1.^2) figure(1) hold on ; plot(t,x,'*');

plot(0:0.01:18,f(b1,0:0.01:18));

xlabel('时间t（小时）'); ylabel('生物量x(t)（克）'); legend('观测值','模拟值'); figure(2)

plot(t,r1,'*',[0,20],[0,0]); axis([0,20,-40,40]); xlabel('时间t（小时）'); ylabel('模拟误差');

输出图像

输出结果

b1 =

0.5470 663.0220 9.1355 r1 =

Columns 1 through 8

0.4645 2.6701 2.4465 2.6143 -2.3504 1.6469 -5.1784 -2.1484 Columns 9 through 16

1.7686 5.0607 3.8475 -4.8434 -7.4295

2.9714 -0.5418 0.7993 Columns 17 through 19

0.2996 0.8932 1.2820 s =

194.3254

即求得固有增长率r=0.5470，酵母培养物最大容量N=663.0220，初始值x 0=9.1355，误差平方和s=194.3254，酵母培养物生物量增长的经验公式

0.54700.5470663.02209.13556057.037

()9.1355(663.02209.1355)9.1355653.8865t t

x t e e

--?=

=+-+ 且模拟效果较好。 4.

符号说明

x k ～第k 年草场的密度单位； y k ～第k 年操场上鹿的数量； r ～草场上草的年固有增长率； c ～草场密度单位减少的速度大小； d ～没有草，鹿群的年死亡率； b ～草对鹿群死亡的补偿率； N ～草场的最大密度单位；

模型建立与求解

（1）比较将100只鹿放入密度为1000和密度为3000的两种草场的情况下，草和鹿两个总群的数量演变过程，因为草的生长服从Logistic 模型，建立差分方程组模型为

11(1) k k k k k k k k k k k x cx y x x rx N N bx y y y dy N

++=+--=??-?

+???? 平衡点为0(0,0)P 、1(,0)P N 和2()

(,)dN rN b d P b cb

-，通过Matlab 计算代码如下： M 文件fun4_1.m

r=0.8;c=1.6;d=0.9;b=1.5;N=3000; hold on ; figure(1)

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

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

plot(1:21,x1(1:21),'*',1:21,y1(1:21),'o'); xlabel('第k 年'); ylabel('种群密度'); gtext('草场密度'); gtext('鹿群密度'); figure(2) plot(x1,y1,'*'); xlabel('草场密度'); ylabel('鹿群密度'); figure(3)

x2(1)=3000;y2(1)=100; for k=1:20

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

plot(1:21,x2(1:21),'*',1:21,y2(1:21),'o'); xlabel('第k 年'); ylabel('种群密度'); gtext('草场密度'); gtext('鹿群密度'); figure(4) plot(x2,y2,'*'); xlabel('草场密度'); ylabel('鹿群密度');

输出图像

figure （1）

figure（2）

figure（3）

figure（4）

输出结果

x1 =

Columns 1 through 8

1000.0 1480.0 2032.5 2502.3 2759.3 2824.5 2787.3 2692.5

Columns 9 through 16

2548.5 2355.9 2126.1 1889.9 1692.9 1574.9 1550.0 1604.8

Columns 17 through 21

1708.6 1823.1 1912.7 1955.1 1946.9

y1 =

Columns 1 through 8

100.0000 60.0000 50.4000 56.2598 76.0148 112.4757 170.0941 254.0580

Columns 9 through 16

367.4356 504.9419 645.2937 750.5099 784.2288 742.2214 658.6837 576.3405

Columns 17 through 21

520.0917 496.3208 502.0420 530.3204 571.4523

x2 =

Columns 1 through 8

3000.0 2840.0 2718.8 2570.0 2378.2 2149.0 1910.4 1707.1

Columns 9 through 16

1580.6 1547.6 1596.6 1698.0 1813.6 1907.0 1954.2 1950.1

Columns 17 through 21

1905.9 1841.9 1779.7 1736.4 1720.3

y2 =

Columns 1 through 8

100.0000 160.0000 243.2000 354.9293 491.5830 633.7025 744.2929 785.3710

Columns 9 through 16

748.9037 666.7620 582.6070 523.3679 496.6857 500.0647 526.8075 567.4258

Columns 17 through 21

610.0197 642.3314 655.7913 649.1482 628.5048

（2）比较将100只鹿放入密度为1000和密度为3000的两种草场的情况下，草和鹿两个总群的数量演

变过程，因为草的生长服从Logistic 模型，建立常微分方程组模型为

(1) dx x cxy

rx dt N N dy bxy dy dt N

????

=--=-+??? 平衡点为0(0,0)P 、1(,0)P N 和2()(

,)dN rN b d P b cb

-，并通过两式相除消去dt 得 ()dy Ndy bxy

dx rx N x cxy

-+=-- 通过Matlab 计算代码如下： 函数fun4_2 M 文件fun4_2.m

function dx=fun4_2(t,x)

r=0.8;c=1.6;d=0.9;b=1.5;N=3000; dx=zeros(2,1);

dx(1)=x(1)*r*(1-x(1)/N)-c*x(1)*x(2)/N; dx(2)=-d*x(2)+b*x(1)*x(2)/N;

主程序 M 文件fun4_3.m

hold on ; figure(1)

[t1,x1]=ode45('fun4_2',0:20,[1000,100]); plot(t1,x1(:,1),'-',t1,x1(:,2),'-'); xlabel('第k 年'); ylabel('种群密度'); gtext('草场密度'); gtext('鹿群密度'); figure(2)

plot(x1(:,1),x1(:,2),'-'); xlabel('2Y3??ü?è'); ylabel('?1èo?ü?è'); figure(3)

[t2,x2]=ode45('fun4_2',0:20,[3000,100]); plot(t2,x2(:,1),'-',t2,x2(:,2),'-'); xlabel('第k 年'); ylabel('种群密度'); gtext('草场密度'); gtext('鹿群密度'); figure(4)

plot(x2(:,1),x2(:,2),'-'); xlabel('草场密度'); ylabel('鹿群密度');

输出图像

figure（1）

figure（2）

figure（3）

figure（4）

输出结果x1 =

1000.0 100.0 1521.3 76.2 2024.0 75.4 2385.5 92.9 2574.7 131.6 2622.1 197.0 2562.8 294.1 2424.9 417.1 2239.2 545.0 2045.3 646.5 1880.8 700.4 1767.4 707.1 1708.7 683.9 1694.0 649.9 1709.0 618.3 1739.0 595.0 1771.3 582.0 1798.1 577.7 1815.3 579.9 1822.6 585.7 1822.1 592.4

x2 =

3000.0 100.0 2848.8 175.4 2684.1 284.8 2483.9 422.1 2257.2 561.9 2038.1 667.9 1860.9 718.1 1746.1 717.3 1691.1 686.9 1682.6 648.1 1704.2 614.0 1739.2 590.4 1774.8 578.0 1802.8 575.0 1819.7 578.5 1825.9 585.4 1824.0 592.9 1817.5 599.2 1809.6 603.2 1802.6 605.0 1797.7 605.0

模型分析

（1）相同点：两个模型给出的草场密度和鹿群密度的变化方向一致，都会趋向于稳定状态；

（2）不同点：差分方程模型相对于是取步长为1时的欧拉方法计算微分方程的数值解，而常微分方程模型使用Matlab中的ode45函数步长比较小，使得常微分方程模型的精度比差分方程模型的精度大，也使得差分方程模型的变化趋势比微分方程模型大。