章绍辉数学建模第四章

数学建模第四章2.符号说明：酒精量是指纯酒精的质量，单位是毫克；酒精含量是指纯酒精的浓度，单位是毫克/百毫升；t～时刻（小时）；x(t)～在时刻t吸收室（肠胃）内的酒精量（毫克）；1D～在段时间内喝下2瓶啤酒后吸收室内的酒精量（毫克）；0c(t)～在时刻t吸收室（肠胃）的酒精含量（毫克/百毫升）；1c(t)～在时刻t中心室（血液和体液）的酒精含量（毫克/百毫升）；2V～中心室的容积（百毫升）；k～酒精从吸收室吸收进入中心室的速率系数；1k～酒精从中心室向体外排除的速率系数；2k～假如所喝入的酒精完全吸收进中心室却没有排除出体外，中心室的酒精含量将达到的最大值（毫克3/百毫升）；模型假设：假设一：酒是在很短时间内喝的大李在短时间内喝下三瓶啤酒后，酒精先从吸收室（肠胃）吸收进中心室（血液和体液），然后从中心室向体外排除，忽略喝酒的时间，根据生理学知识，假设（1）吸收室在初始时刻t=0时，酒精量立即为3D/2；在任意时刻，酒精从吸收室吸收进入中心室的速0率（吸收室在单位时间内酒精含量的减少量）与吸收室的酒精含量成正比，比例系数为k；1（2）中心室的容积V保持不变；在初始时刻t=0时，中心室的酒精含量为0；在任意时刻，酒精从中心室向体外排除的速率（中心室在单位时间内酒精含量的减少量）与中心室的酒精含量成正比，比例系数为k；2（3）在大李适度饮酒没有酒精中毒的前提下，假设k和k都是常数，与饮酒量无；关12假设二：酒是在较长一段时间（如2小时）内喝的大李在较长一段时间（如2小时）内喝下三瓶啤酒后，酒精先从吸收室（肠胃）吸收进中心室（血液和体液），然后从中心室向体外排除，忽略喝酒的时间，根据生理学知识，假设第1页数学建模（1）吸收室在初始时刻t=0时，酒精量为0；在任意时刻，酒精从吸收室吸收进入中心室的速率（吸收室在单位时间内酒精含量的减少量）与吸收室的酒精含量成正比，比例系数为k；1（2）中心室的容积V保持不变；在初始时刻t=0时，中心室的酒精含量为0；在任意时刻，酒精从中心室向体外排除的速率（中心室在单位时间内酒精含量的减少量）与中心室的酒精含量成正比，比例系数为k；2（3）在大李适度饮酒没有酒精中毒的前提下，假设k和k都是常数，与饮酒量无关；12模型建立和求解：假设一：酒是在很短时间内喝的根据4.1.7小节饮酒驾车案例，假设k=2.0079，k=0.185512则在很短时间内喝下三瓶酒时D33k30103.86155.79V22记喝酒时间为t=0时，c(t)=0，则可根据2kkctee1()--kt-kt13kk2-212来计算t时刻血液中的酒精含量c(t)，Matlab代码如下：2M文件fun2_1.mk1=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('短时间内喝下三瓶酒时，血液中酒精含量随时间的变化过程');输出图像第2页数学建模由图像可知c(t)的函数图像与y=20和y=80分别都有2个交点，故可继续添加代码求出这4个交点坐2标，代码如下：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.58874.1125gt=0.3805即t[0.0689,0.3805]ct2[4.1125,11.5887]20()80（1）当时，，属于饮酒驾车；t（2）当ct()802[0.3805,4.1125]时，，属于醉酒驾驶；假设二：酒是在较长一段时间（如2小时）内喝的根据4.1.7小节饮酒驾车案例，假设k=2.0079，k=0.185512第3页数学建模则在较长一段时间（如2小时）内喝下三瓶酒时D33k30103.8677.895V44记喝酒时间为t=0时，c(t)=0，c(t)=0，则吸收室的酒精量x(t)满足分段初值问题112dxD3,(0)0,0t2kxx041111dtdx1D0k43kxxek),2,(2)(112kdt1111解得3D(1et,02)kt1430kxt()11Dt(e1)e,22kkt4110故中心室内的酒精含量c(t)满足分段初值问题2dc2kckect(1),(0)0,02kt1dt2232dckckeecpt(1),(2),222k1ktdt122321解得kkkkee,0t2ktkt23133k2kk121kkk2122ct()2ke(1)2k1peet,2ktkt32kk1212其中kkkkpee22k2k3133k2kk121kkk21122ke(1)2k1ppee2k2(kk)32kk122121通过Matlab画出c(t)的图像，代码如下：2M文件fun2_2.mk1=2.0079;k2=0.1855;k31=155.79;k3=k31/2;c2=@(t)(k1.*k31)./(k1-k2).*(exp(-k2.*t)-exp(-k1.*t));holdon;第4页数学建模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小时）内喝的');输出图像由图像可知c(t)的函数图像与y=20和y=80分别都有2个交点，故可继续添加代码求出这4个交点坐2标，代码如下：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=第5页数学建模0.623312.6196gt=1.63665.1412即t[0.6233,1.6366]ct20()80，属于饮酒驾车；[5.1412,12.6196]2（1）当时，t（2）当ct[1.6366,5.1412]()80，属于醉酒驾驶。时，23.考虑3.4.2小节“酵母培养物的增长”案例，建立微分方程模型，模拟酵母培养物的增长，故假设t（小时）时刻的酵母生物量为x(t)（克），且t=0时x(t)=x，则x(t)满足微分方程初值问题0dxtxxx0()(0)，dtt()其中为酵母培养物的增长量，根据3.4.2小节可知trxt()(1()/N)其中r>0称为酵母培养物固有增长率，N>0称为酵母培养物的最大容量，则得到阻滞增长方程dxxrx(1)dtN用分离变量法得Nxxt()0xNxe()rt00则根据表3.2通过Matlab拟合出x(t)的参数r、N和x，代码如下=M文件fun3_1.mt=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));第6页数学建模xlabel('时间t（小时）');ylabel('生物量x(t)（克）');legend('观测值','模拟值');figure(2)plot(t,r1,'*',[0,20],[0,0]);axis([0,20,-40,40]);xlabel('时间t（小时）');ylabel('模拟误差');输出图像第7页数学建模输出结果b1=0.5470663.0220r1=9.1355Columns1through80.46452.67012.44653.84751.28202.6143-2.3504-7.42951.6469-5.1784-0.5418-2.14840.7993Columns9through161.76865.0607-4.84342.9714Columns17through190.29960.8932s=194.3254即求得固有增长率r=0.5470，酵母培养物最大容量N=663.0220，初始值x=9.1355，误差平方和0s=194.3254，酵母培养物生物量增长的经验公式663.02209.13556057.037e9.1355(663.02209.1355)9.1355653.8865xt()e0.5470t0.5470t且模拟效果较好。4.符号说明x～第k年草场的密度单位；ky～第k年操场上鹿的数量；kr～草场上草的年固有增长率；c～草场密度单位减少的速度大小；d～没有草，鹿群的年死亡率；b～草对鹿群死亡的补偿率；N～草场的最大密度单位；模型建立与求解（1）比较将100只鹿放入密度为1000和密度为3000的两种草场的情况下，草和鹿两个总群的数量演变过程，因为草的生长服从Logistic模型，建立差分方程组模型为第8页数学建模xcxyxxrxkkk(1)NNk1kkbxykkNyydyk1kkdNrNbd()(0,0)(,0)(,)，通过Matlab计算代码如下：PPNP平衡点为、0和bcb21M文件fun4_1.mr=0.8;c=1.6;d=0.9;b=1.5;N=3000;holdon;figure(1)x1(1)=1000;y1(1)=100;fork=1:20x1(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;endplot(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:20x2(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;endplot(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）第9页数学建模figure（2）第10页数学建模figure（3）figure（4）第11页数学建模输出结果x1=Columns1through81000.01480.02032.52126.12502.31889.91955.12759.31692.91946.92824.52787.31550.02692.51604.8Columns9through162548.52355.91574.9Columns17through211708.61823.11912.7y1=Columns1through8100.000060.000050.400056.259876.0148112.4757170.0941254.0580Columns9through16367.4356504.9419645.2937750.5099784.2288742.2214658.6837576.3405Columns17through21520.0917496.3208502.0420530.3204571.4523x2=Columns1through83000.02840.02718.81596.62570.01698.01736.42378.21813.61720.32149.01907.01910.41954.21707.11950.1Columns9through161580.61547.6Columns17through211905.91841.91779.7y2=Columns1through8100.0000160.0000243.2000354.9293491.5830633.7025744.2929785.3710Columns9through16748.9037666.7620582.6070523.3679496.6857500.0647526.8075567.4258Columns17through21610.0197642.3314655.7913649.1482628.5048（2）比较将100只鹿放入密度为1000和密度为3000的两种草场的情况下，草和鹿两个总群的数量演第12页数学建模变过程，因为草的生长服从Logistic模型，建立常微分方程组模型为dxxcxyrx(1)dtdyNbxyNdtdyNdNrNb(dPNP和b)P0(0,0)(,0)(,)，并通过两式相除消去dt得平衡点为、cb12dyNdybxydxrxNx)cxy(通过Matlab计算代码如下：函数fun4_2M文件fun4_2.mfunctiondx=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.mholdon;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年');yl