章绍辉数学建模第四章.docx
《章绍辉数学建模第四章.docx》由会员分享，可在线阅读，更多相关《章绍辉数学建模第四章.docx（22页珍藏版）》请在冰豆网上搜索。
章绍辉数学建模第四章
第四章
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；在任意时刻，酒精从吸收室吸收进入中心室的速率（吸收室在单位时间内酒精含量的减少量）与吸收室的酒精含量成正比，比例系数为k1；
（2）中心室的容积V保持不变；在初始时刻t=0时，中心室的酒精含量为0；在任意时刻，酒精从中心室向体外排除的速率（中心室在单位时间内酒精含量的减少量）与中心室的酒精含量成正比，比例系数为k2；
（3）在大李适度饮酒没有酒精中毒的前提下，假设k1和k2都是常数，与饮酒量无关；
模型建立和求解：
假设一：
酒是在很短时间内喝的
根据4.1.7小节饮酒驾车案例，假设
k1=2.0079，k2=0.1855
则在很短时间内喝下三瓶酒时
记喝酒时间为t=0时，c2（t）=0，则可根据
来计算t时刻血液中的酒精含量c2（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））;
holdon;
plot（0:
0.01:
20,c2（0:
0.01:
20））;
plot（[0,20],[20,20],[0,20],[80,80]）;
xlabel（'时刻t（小时）'）;
ylabel（'血液中酒精含量c2（毫克/百毫升）'）;
title（'短时间内喝下三瓶酒时，血液中酒精含量随时间的变化过程'）;
输出图像
由图像可知c2（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.068911.5887
gt=
0.38054.1125
即
（1）当
时，
，属于饮酒驾车；
（2）当
时，
，属于醉酒驾驶；
假设二：
酒是在较长一段时间（如2小时）内喝的
根据4.1.7小节饮酒驾车案例，假设
k1=2.0079，k2=0.1855
则在较长一段时间（如2小时）内喝下三瓶酒时
记喝酒时间为t=0时，c1（t）=0，c2（t）=0，则吸收室的酒精量x1（t）满足分段初值问题
解得
故中心室内的酒精含量c2（t）满足分段初值问题
解得
其中
通过Matlab画出c2（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））;
holdon;
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.623312.6196
gt=
1.63665.1412
即
（1）当
时，
，属于饮酒驾车；
（2）当
时，
，属于醉酒驾驶。
3.
考虑3.4.2小节“酵母培养物的增长”案例，建立微分方程模型，模拟酵母培养物的增长，故假设t（小时）时刻的酵母生物量为x（t）（克），且t=0时x（t）=x0，则x（t）满足微分方程初值问题
，
其中
为酵母培养物的增长量，根据3.4.2小节可知
其中r>0称为酵母培养物固有增长率，N>0称为酵母培养物的最大容量，则得到阻滞增长方程
用分离变量法得
则根据表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）
holdon;
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.5470663.02209.1355
r1=
Columns1through8
0.46452.67012.44652.6143-2.35041.6469-5.1784-2.1484
Columns9through16
1.76865.06073.8475-4.8434-7.42952.9714-0.54180.7993
Columns17through19
0.29960.89321.2820
s=
194.3254
即求得固有增长率r=0.5470，酵母培养物最大容量N=663.0220，初始值x0=9.1355，误差平方和s=194.3254，酵母培养物生物量增长的经验公式
且模拟效果较好。
4.
符号说明
xk～第k年草场的密度单位；
yk～第k年操场上鹿的数量；
r～草场上草的年固有增长率；
c～草场密度单位减少的速度大小；
d～没有草，鹿群的年死亡率；
b～草对鹿群死亡的补偿率；
N～草场的最大密度单位；
模型建立与求解
（1）比较将100只鹿放入密度为1000和密度为3000的两种草场的情况下，草和鹿两个总群的数量演变过程，因为草的生长服从Logistic模型，建立差分方程组模型为
平衡点为
、
和
，通过Matlab计算代码如下：
M文件fun4_1.m
r=0.8;c=1.6;d=0.9;b=1.5;N=3000;
holdon;
figure
（1）
x1
（1）=1000;y1
（1）=100;
fork=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;
fork=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=
Columns1through8
1000.01480.02032.52502.32759.32824.52787.32692.5
Columns9through16
2548.52355.92126.11889.91692.91574.91550.01604.8
Columns17through21
1708.61823.11912.71955.11946.9
y1=
Columns1through8
100.000060.000050.400056.259876.0148112.4757170.0941254.0580
Columns9through16
367.4356504.9419645.2937750.5099784.2288742.2214658.6837576.3405
Columns17through21
520.0917496.3208502.0420530.3204571.4523
x2=
Columns1through8
3000.02840.02718.82570.02378.22149.01910.41707.1
Columns9through16
1580.61547.61596.61698.01813.61907.01954.21950.1
Columns17through21
1905.91841.91779.71736.41720.3
y2=
Columns1through8
100.0000160.0000243.2000354.9293491.5830633.7025744.2929785.3710
Columns9through16
748.9037666.7620582.6070523.3679496.6857500.0647526.8075567.4258
Columns17through21
610.0197642.3314655.7913649.1482628.5048
（2）比较将100只鹿放入密度为1000和密度为3000的两种草场的情况下，草和鹿两个总群的数量演变过程，因为草的生长服从Logistic模型，建立常微分方程组模型为
平衡点为
、
和
，并通过两式相除消去dt得
通过Matlab计算代码如下：
函数fun4_2M文件fun4_2.m
functiondx=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
holdon;
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（'²Ý³¡ÃܶÈ'）;
ylabel（'¹ȺÃܶÈ'）;
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.0100.0
1521.376.2
2024.075.4
2385.592.9
2574.7131.6
2622.1197.0
2562.8294.1
2424.9417.1
2239.2545.0
2045.3646.5
1880.8700.4
1767.4707.1
1708.7683.9
1694.0649.9
1709.0618.3
1739.0595.0
1771.3582.0
1798.1577.7
1815.3579.9
1822.6585.7
1822.1592.4
x2=
3000.0100.0
2848.8175.4
2684.1284.8
2483.9422.1
2257.2561.9
2038.1667.9
1860.9718.1
1746.1717.3
1691.1686.9
1682.6648.1
1704.2614.0
1739.2590.4
1774.8578.0
1802.8575.0
1819.7578.5
1825.9585.4
1824.0592.9
1817.5599.2
1809.6603.2
1802.6605.0
1797.7605.0
模型分析
（1）相同点：
两个模型给出的草场密度和鹿群密度的变化方向一致，都会趋向于稳定状态；
（2）不同点：
差分方程模型相对于是取步长为1时的欧拉方法计算微分方程的数值解，而常微分方程模型使用Matlab中的ode45函数步长比较小，使得常微分方程模型的精度比差分方程模型的精度大，也使得差分方程模型的变化趋势比微分方程模型大。
升级会员 