撰文缘由

笔者在撰写本文前不得不说一下，B站里面的确有许许多多的宝藏UP主，本人也是在听入迷了相关饮酒驾车建模视频的前提下心中油然而生写作热情。因此本文只是相关视频的代码复现，并会在相应基础上进行一定程度的拓展。
笔者简介：CCNU计算机科学与技术jdb大二一员，也喜欢日漫唱歌梅老板和弹钢琴。
原视频名称:数学建模案例 6 饮酒驾驶的问题
(https://www.bilibili.com/video/BV1MR4y1g7ns?spm_id_from=333.999.0.0)

模型建立

在建立模型前我们作出合理假设:
1）人体体液中酒精浓度与血液中酒精浓度相同。
2）人体自身产生的酒精忽略不计，即正常不饮酒情况下，人体体液中的酒精浓度看作0。
3）人体体液对酒精吸收速度与当前肠胃中酒精含量成正比，比例系数为a。
4）酒精代谢速度与当前体液中酒精含量成正比，比例系数为b。
5）整个过程中人没有摄入任何影响代谢的药物和作剧烈运动。
6）人体的吸收率和代谢率是恒定的。
7）一瓶啤酒中约含酒精21700mg。
8）体重为70kg的人体液体积约为420（其中每单位100ml）。

“一口气”饮酒模型

下面我们对模型的参数进行说明:
1.Q:饮酒酒桶所含酒精量大小。
2.y(t):肠胃中所含酒精量(随时间变量t变化)。
3.x(t):体液中所含酒精量(随时间变量t变化)。
4.c(t):体液中所含酒精浓度(由x(t)/v0得到)。
5.v0:体液体积(其中每单位100ml)。
6.a*y(t):吸收速率(假设与含酒精量成正比关系)。
7.b*y(t):代谢速率(假设与酒精浓度成正比关系)。

所谓“一口气”喝酒的意思就是酒后驾车嫌疑人以迅雷不及掩耳之势抢铃儿响叮当之态喝下了一桶酒，这样的话酒会全部直接进入肠胃中呆着。那这样的话我们就可以开始建立微分方程模型。

我们分别对y(t)和c(t)进行分对象的讨论,对肠胃中的y(t)而言，单位时间中的变化量就是-a*y。于是乎我们可以列出微分方程组进行表示如下:

dy/dt=-a*t
y(0)=Q0 (其中Q0为酒桶中酒量大小，作为初始条件)

对于体液中的酒精浓度c(t)而言,在单位时间里面，你既有从肠胃中输送来的酒精，又有通过新陈代谢消失掉的酒精,那么我们可以得到下面的微分方程组:

dc/dt=a * Q0 * e^(-a*t) / v0 - b * c
c(0)=c0 (其中c0为未饮酒时体液中的酒精浓度，作为初始条件)

你咋一看这a * Q0 * e^(-at)有点陌生捏？
其实a * Q0 * e^(-at)就是上一个方程的y(t)的求解表达式。
准确的说就是dc/dt=a * y(t) / v0-b * c表达式的偷懒但是很精确的表示形式。

对于求解微分方程组你可以手算也可以心算还可以用除了不会生孩子其他啥都会的MATLAB求解，本人用的是MATLAB R2020b版本。

代码段如下表示:

%%一口气饮酒，肠胃与体液中酒精浓度随时间的变化 
dsolve('Dy=-a*y','y(0)=Q','t')%肠胃中酒精浓度的变化
dsolve('Dc=a*Q*exp(-a*t)/v0-b*c','c(0)=c0','t')%体液中酒精浓度的变化

得到的体液中的酒精浓度表达式就是:
c(t)=e ^ (-b * t)* (c0 + (Q * a) / (v0 * (a - b))) - (Q * a * e(b * t - a * t) * e^(-b * t))/(v0*(a - b))

不过这样的形式也太复杂了吧，但是我们可以利用案例中的数据进行曲线拟合，并且利用lsqcurvefit()函数求出正比例系数a和b。下面是案例的数据展示:

监测时间(h)	体液中酒精浓度(mg/100ml)
0.25	30
0.5	68
0.75	75
1	82
1.5	82
2	77
2.5	68
3	68
3.5	58
4	51
4.5	50
5	41
6	38
7	35
8	28
9	25
10	18
11	15
12	12
13	10
14	7
15	7
16	4

我们不妨取Q0=3*21700（一口气3瓶啤酒的量！），v0=420，c0=0的初始参数条件值，经过拟合参数a,b之后作出监测时间与体液中酒精浓度的曲线图如下:

如果以20mg/(100ml)的量为界线，那我们知道该同志饮酒一次后至少需要12小时左右才能够恢复到正常状态。

可见饮酒驾车情形下的危险系数之大!

第一个模型我们就简单介绍到这里。

“匀速”饮酒模型

这时有人就不服气了，他说我要是慢慢地喝酒，久到一个天长地久沧海桑田之后酒精浓度也就会趋于正常不是吗?

当然我们不讨论这种情况，但是我们将讨论一个人在2个小时里面将酒喝完后他的体液中酒精浓度变化情况。

其实这个模型和“一口气”饮酒模型唯一的不同点就是酒精进入肠胃是有速度的啦！
我们知道酒精以匀速进入肠胃，那么速度f1(t)可以表示为:

f1(t)=Q0/T ,t<=T
f1(t)=0 ,t>T (T=2的条件)

由此不难得出肠胃的酒精浓度表达式为:

dy/dt=Q0/T-a * y , t<=2
dy/dt=-a * y , t>2

可以解出y(t)的表达式，最后由a * y(t)得出吸收速率的分段关系式：

f(t)=(Q0 * (1-e^(-a * t)) ) / 2 ,t<=2
f(t)=-(Q0 * (1-e^(2 * a))) / 2 ,t > 2

求出吸收速率后，我们体液中酒精浓度变化的微分方程组也就是水到渠成了。由此体液中酒精浓度表达式分段如下:

dc/dt=f(t) / v0 - b * c(t)
f(t)=Q0/T-a * y(t) ,t<=2
f(t)=-a * y , t>2
c(0)=0
c(T+)=c(T-) (左右极限要连续嘛)

那这样的话我们就可以得到c(t)的表达式，但是这个太复杂了，没有意义展示，详细部分会在代码区里面全部展示。c(t)关于监测时间的相应的曲线图也可以绘制如下:

其实2小时喝完酒也就比一口气喝完酒好拿磨一丢丢而已，大哥不说二哥，酒后驾车是真滴危险呐！那么第二个模型我们也就简单介绍到这里了。

“匀加速”饮酒模型（稍微拓展）

那你又有话要说了，凭什么？为什么喝酒的速度会是匀速的？确实，我们可以让模型更加精确一点点而已，根据高中学过的加速度，我们可以建立“匀加速”饮酒模型?试试看吧。

对第二问稍作改动，我们让那个人在2个小时里面以匀加速方式喝酒并且将酒喝完。然后我们不难得出相应的计算加速度的公式。

由 1/2 * a *T ^2=Q0 得到
a=Q0 / 2

那么我们就可以修改一下酒精进入肠胃的速度表达式f2(t)

dy/dt=Q0 * t^2 / 2 - a * y , t<=2
dy/dt=-a * y , t>2
y(T+)=y(T-) （左右极限连续嘛）

然后我们就可以得到匀加速喝酒方式下肠胃中酒精浓度的表达式:

y(t)=Q0/(a^3) - (Q0 * e ^ (-a * t))/(a ^ 3)+(Q0 * t ^2) / (2 * a)-(Q0 * t)/(a ^2), t<=2
y(t)=-(e ^(-a * t) * (Q0-Q0 * e ^ (2 * a)+2 * Q0 * a * e^(2 * a)-2 * Q0 * a ^ 2 *e^(2 * a)))/(a ^3), t>2

然后同理求得体液中的酒精浓度c(t)与时间的关系式，当然这种模型的计算结果相当复杂，符号表达式过于冗余，不建议使用。详细的符号表达式部分会在代码区里面全部展示。

代码区

%%一口气饮酒，肠胃与体液中酒精浓度随时间的变化 
dsolve('Dy=-a*y','y(0)=Q','t')%肠胃中酒精浓度的变化
 dsolve('Dc=a*Q*exp(-a*t)/v0-b*c','c(0)=c0','t')%体液中酒精浓度的变化
 %%两小时均匀饮酒，肠胃与体液中酒精浓度随时间的变化
 dsolve('Dy=Q0/2-a*y','y(0)=0','t') %两小时内肠胃中酒精浓度变化
 dsolve('Dy=-a*y','y(2)=(Q0 - Q0*exp(-a*2))/(2*a)','t')%两小时后肠胃中酒精浓度变化
 %进而可以得到吸收速率f(t)=a*y(t)==(Q0 - Q0*exp(-a*t))/2  (t<=2)   
 %-(exp(-a*t)*(Q0 -Q0*exp(2*a)))/2   t>2
 dsolve('Dc=(Q0 - Q0*exp(-a*t))/(2*v0)-b*c','c(0)=0','t')%两小时内体液中酒精浓度变化
 dsolve('Dc=(-a*-(exp(-a*t)*(Q0 - Q0*exp(2*a)))/(2*a))/v0-b*c','c(2)=exp(-b*2)*((Q0*exp(-a*2)*exp(b*2))/(2*(a*v0 - b*v0))+ (Q0*exp(b*2))/(2*b*v0))+ (Q0*a*exp(-b*2))/(2*v0*b^2 - 2*a*v0*b)','t')
 %两小时以后体液中酒精浓度变化
 
  dsolve('Dy=Q0*t^2/2-a*y','y(0)=0','t') %匀加速两小时内肠胃中酒精浓度变化
  dsolve('Dy=-a*y','y(2)=Q0/a^3 - (Q0*exp(-a*2))/a^3 + (Q0*2^2)/(2*a) - (Q0*2)/a^2','t')%匀加速两小时后肠胃中酒精浓度变化
  
  dsolve('Dc=(Q0/a^3 - (Q0*exp(-a*t))/a^3 + (Q0*t^2)/(2*a) - (Q0*t)/a^2)*a/v0-b*c','c(0)=0','t')%匀加速两小时内体液中酒精浓度变化
  dsolve('Dc=(a*(exp(-a*t)*(Q0 - Q0*exp(2*a) + 2*Q0*a*exp(2*a) - 2*Q0*a^2*exp(2*a)))/a^3)/v0-b*c','c(2)=(exp(-b*2)*((2*Q0*exp(b*2))/b + (2*Q0*exp(b*2 - a*2))/(a - b) + (2*Q0*a*exp(b*2)*(a + b))/b^3 + (Q0*a^2*2^2*exp(b*2))/b - (2*Q0*a*2*exp(b*2)*(a + b))/b^2))/(2*a^2*v0) - (Q0*a*exp(-b*2))/(v0* b^4 - a*v0*b^3)','t')
  %匀加速两小时以后体液中酒精浓度变化
 

 %%利用模型一拟合参数a(吸收率),b(代谢率)
 t=[0.25,0.5,0.75,1,1.5,2,2.5,3,3.5,4,4.5,5,6,7,8,9,10,11,12,13,14,15,16];
 c=[30,68,75,82,82,77,68,68,58,51,50,41,38,35,28,25,18,15,12,10,7,7,4];
 x0=[2,0.2];
 x=lsqcurvefit(@jiujiafun,x0,t,c)
 %%绘图观察
t=linspace(0,10,1000);
c1=170.5*(exp(-0.1842.*t)-exp(-2.0261.*t));
plot(t,c1,'--r')
hold on
tt=linspace(0,2,200);
c2=228.411*(1.8419+0.1842*exp(-2.0261.*tt)-2.0261*exp(-0.1842.*tt));
plot(tt,c2,'-b');
ttt=linspace(2,20,800);
c3=228.411*(-10.4117*exp(-2.0261.*ttt)+0.9025*exp(-0.1842.*ttt));
hold on
plot(ttt,c3,'--b');
legend('一口气','两小时内','两小时后');
%%可以找找找体液酒精浓度降为20mg/100ml的大致时间

%各曲线上体液酒精浓度最大的时刻和体液酒精浓度值
[a,b]=max(c1);
a
b=b/100
[e,f]=max(c2);
e
f=f/100
[r,s]=max(c3);
r
s=s/100+2