酒后血液中酒精含量的数学模型1
酒后血液中酒精含量的数学模型
摘要
针对酒后驾车普遍存在并致交通肇事居高不下的现实 ,掌握饮酒后不同时刻血液中酒精的浓度非常必要。本文根据药物动力学知识,首先用微分方程建立了基本模型并推导出在长时间、 瞬时间和分段瞬时内饮酒的数学模型 ,从理论上完整的描述了人体血液中的酒精含量的变化过程。其次,根据所给数据 ,利用数学软件Matlab 对基本的模型进行了拟合 ,得出基本模型中的待定系数,并得出了人在不同情况下饮酒后的酒精含量与时间的关系图 从图中可以很好的反映出人体血液中的酒精含量的变化规律,它们的变化规律与实际变化相吻合 ,从而证明了所建的模型基本符合要求,进而可以根据关系图讨论题中的问题。运用微积分理论 ,建立微分方程并推导出在长时间、 瞬时间和分段瞬时内饮酒的数学模型 ,检验结果表明模型正确 ,理论数据与实际相吻合。从数学理论上解决了不同体重、不同时间饮用不同量的酒后在不同时刻血液中的酒精含量。并得出了人在不同情况下饮酒后的酒精含量与时间的关系图,从图中可以很好的反映出人体血液中的酒精含量的变化规律,它们的变化规律与实际变化相吻合 ,从而证明了所建的模型基本符合要求,进而可以根据关系图讨论题中的问题。
关键词:吸收速率 消除速率 数学模型 非线性数据拟合 Matlab 微分方程
1 问题的提出
据报载,2010年,全国共接报道路交通事故3906164起,同比上升35.9%。其中,
涉及人员伤亡的道路交通事故219521起,造成65225人死亡、254075人受伤,直接财产损失9.3亿。而2003年全国道路交通事故死亡人数仅仅为10.4372万,其中因饮酒驾车造成的占有相当的比例。
针对这种严重的道路交通情况,国家质量监督检验检疫局2004年5月31日发布了新的《车辆驾驶人员血液、呼气酒精含量阈值与检验》国家标准,新标准规定,车辆驾驶人员血液中的酒精含量大于或等于20毫克/百毫升,小于80毫克/百毫升为酒后驾车(原标准是小于100毫克/百毫升),血液中的酒精含量大于或等于80毫克/百毫升为醉酒驾车(原标准是大于或等于100毫克/百毫升)。
李强在中午12点喝了一瓶啤酒,下午6点检查时符合新的驾车标准,紧接着他在吃晚饭时又喝了一瓶啤酒,为了保险起见他呆到凌晨2点才驾车回家,又一次遭遇检查时却被定为饮酒驾车,这让他既懊恼又困惑,为什么喝同样多的酒,两次检查结果会不一样呢?
我们将参考下面给出的数据并自己收集资料建立饮酒后血液中酒精含量的数学模型,并讨论以下问题:
1. 对李强碰到的情况做出解释,是意外还是必然?和什么因素相关?
2. 王东在短时间内喝了三瓶啤酒,试问多长时间内驾车就会违反上述标准?(分酒后和醉酒两种情况讨论),如果王东在半个小时内喝了半斤37度白酒,情况又会如何?
3. 王东几乎一口气喝了三瓶啤酒,休息一小时后,又在两个小时内喝了一斤37度白酒,试问多长时间内驾车就会违反上述标准?(分酒后和醉酒两种情况讨论)
4. 怎样估计血液中的酒精含量在什么时间最高,试着对第2问或者第3问做出回答。
5. 根据你的模型论证:如果天天喝酒,是否还能开车?能否对符合要求的喝酒方式和喝酒的量给出直观的说明。
6. 根据你做的模型并结合新的国家标准写一篇短文,给想喝一点酒的司机如何驾车提出忠告。
参考数据
1. 人的体液占人的体重的65%至70%,其中血液只占体重的7%左右;而药物(包括酒精)在血液中的含量与在体液中的含量大体是一样的。
2. 体重约70kg的某人在短时间内喝下3瓶啤酒后,隔一定时间测量他的血液中酒精含量(毫克/百毫升),得到数据如下:
2 问题的分析
考虑饮酒后酒精在人体内的变化情况,酒精被饮入体内首先是进入胃中,然后再随着血液循环进入体液,然后再由体液分解排出体外。所以我们可以对问题进行如下化简:在酒精吸收和分解的过程中,我们考虑酒精在进入胃的过程中没有损失,而胃内的酒精只是在向体液中渗透,并不考虑体液中的酒精反向渗透回胃内。酒精进入体液中后在进入血液中,我们把血液和体液分开来考虑,虽然血液只占了体液的一小部分,对结果也会有一定的误差。所以我们建立模型时就把胃看成了一个空间,把血液和体液分开看成
两个空间,而这三个空间的关系是酒精从胃渗透向体液再进入血液,而血液液中的酒精只是通过分解排出。
(1)对问题1的分析
我们运用模型3酒是分段瞬时喝完的(间歇饮酒方式)的数学模型,将首次饮酒后经过了六个小时后再次饮酒,这个时候血液中的酒精浓度计算应该是,首次饮酒在血液中的残留继续分解,而第二次饮酒还要经过一个吸收和分解的过程,所以再过八小时测出的血液中的酒精浓度和首次饮酒也有关系。
(2)对问题2的分析
我们运用模型2酒是瞬时喝完(单次饮酒方式)的数学模型,在短时间内喝一定量的酒,经过模型可直接求解出血液中酒精浓度等于20毫克/百毫升所对应的时间和血液中酒精浓度为等于80毫克/百毫升所对应的时间。而对于长时间饮酒,我们可以认为酒是匀速饮入,我们对时间进行分割,运用模型1酒在T小时内喝完(连续均匀饮酒方式)的数学模型,同样可以求解出血液中酒精浓度等于20毫克/百毫升所对应的时间和血液中酒精浓度为等于80毫克/百毫升所对应的时间。
(3)对问题3的分析
我们运用模型4酒先瞬时喝D1,t1小时后又在T小时内喝了D2的数学模型,经过此模型可直接求解出血液中酒精浓度等于20毫克/百毫升所对应的时间和血液中酒精浓度为等于80毫克/百毫升所对应的时间。
(4)对问题4的分析
无论是短时间饮酒,还是长时间饮酒,我们都可以根据模型1酒在T小时内喝完(连续均匀饮酒方式)和模型2酒是瞬时喝完的(单次饮酒方式)很容易求出血液中酒精的含量在何时最大。
(5)对问题5的分析
根据前面几个问题可以得出每天不论短时间内还是长时间内饮过量的酒,都不能在很长一个时间中恢复标准,即使长时间内均匀时间段中喝少量的酒,人体血液中的酒精的含量也会积少成多,超过标准。
(6)对问题6的分析
对广大爱喝酒的朋友提一些建议和意见,希望他们要爱惜生命,不要饮酒上路。 本文从数学理论上较详细的讨论各种不同情况下不同时刻血液中酒精含量,理论与实际相吻合。通过讨论 ,可将科学诉诸于大众 ,打消驾车者的侥幸心理 ,防患于未然。
3 模型的假设
假设1:酒精从胃部向体液的转移速率,及向外排除的速率分别与胃部和体液中的
酒精浓度成正比。
假设2:体液总体积保持不变。
假设3:进入胃里的酒精全部扩散到体液里,。
假设4:酒精并不会从体液反向渗入到胃部。
假设5:酒精只会通过血液排出体外。
假设6:在较短时间内喝酒的情况下,酒精量是瞬时进入胃里的。
假设7:在较长一段时间内喝酒的情况下,酒精量是匀速进入胃里的。
假设8:酒精被正常吸收和排出,排除呕吐等一些非正常的排出情况
假设9:忽略人对酒精的敏感度以及对酒精的分解能力存在的个性化差异
假设10:啤酒瓶的容量为600毫升,酒精浓度为5%,其他规格不考虑。
4 符号的说明
5 模型建立与求解
5.1基本模型
酒精在人体内的变化速率 ,即动力学方程为:
dx
dtr1(t)r2(t)
因为酒精消除速率r(t)正比于胃内的酒精含量x(t)[1]。
1)(
于是有 r2(t)kx(t) (2) 式(1)可变为 dxr1(t)kx(t) (3) dt
用q(t)表示t时刻被血液吸收的酒精量
则 dqr1(t) dt
5.1.1酒在T小时内喝完(连续均匀饮酒方式)
1、当0#tT时,此时假设酒精是连续均匀摄入体内 ,胃内的酒精量为:
tFDq(t) T
根据药物学知识 ,人体吸收酒精的速率与胃内的酒精量成正比;人体排出酒精的x0(t)速率与体液中的酒精量成正比[2],将有
tr1(t)k1FDq(t) T
于是有
dqtk1FDq(t), q(0)0 dtT
解之得
q(t)
特别地 FD1k1tt(e1) Tk1
FD1k1TT(e1) Tk1q(T)
从而有
r1(t)(1ek1t) dx=r1(t)-kx(t),x(0)=0(初始时胃内酒精含量) dt
可得 再由
x(t)FD11k1tktkt(ee)(1e) (4) Tkk1k
FD11k1tktkt(ee)(1e) (5) Tkk1k c(t)
2、当t³T时,此时酒已喝完。胃内的酒精量为: FD-q(t),此时
dq
dtk1FDq(t)
由q(T)的取值可得:q(t)ek1t(FDek1tN) 其中常数NFD(1ek1T
k)
1T
从而,rk
1(t)k1Ne1t代入方程(1), 并由式(4)求出,可以解得:
x(t)k1N
kkek1tMekt
1
其中常数MFDk1(ekT1)
TK(kk
1)
c(t)1
Vk1Nek1tMekt kk1
当0#tT时,dc(t)
dt>0,c单调递增,当t³T时,令dc(t)
1(t)dt=0,得lnek0t01
t10
2k
0k1
此时体液中酒精浓度c(t) 达到最大值,记为 Tm1。
5.1.2酒是瞬时喝完的(单次饮酒方式)
由第一部分的讨论,当T0时,N?FD,M?FDk1
(k-k
1)
由式(6)、(7)可得
x(t)k1FD
k(ektek1t)
1k
c(t)k1FD
V(k(ektek1t)
1k)
令dc(t)
dt=0得 t1k
kln1
1kk
此时体液中酒精浓度c(t) 达到最大值,记为 Tm2。
6) 7) (8)(9)(10) ( (
当时间充分大时 , r1=0 ,方程 (1)变为
dxkx(t) dt
于是 c(t)Cekt (11)
5.1.3酒是分段瞬时喝完的(间歇饮酒方式)
由于在日常生活中喝酒是分段瞬时的 ,可假设初瞬时摄入酒精量为 D1,t0时瞬
时摄入酒精量为D2,则:
1、当0#tt¢时
q(t)FD1(1ek1t)
x(t)k1FD1ktk
kk(ee1t)
1
c(t)k1FD1
V(kk)(ektek1t)
1
2、当t³t0时
由于又摄入酒精量D2,胃内酒精量为D1+D2-q¢(t),则:
dq(t)
dtk1F(D1D2)q(t)
且q¢(t)=FDkt
1(1-e10)
由此解出
q(t)F(D1D2)(FD1FDktk
2e10)e1t
rt
1k1(FD1FDk
2e10)ek1t
由(1)式及
x(t)k1FD1(ekt0ek1t0
k)
1k
可得 x(t)k1kt0k1t
kk(FD1FD2e)ekt(FD1FD2e0)ek1t (12)
1
c(t)k1kt0k1t0
V(k(FDFDe)ekt(FDFDek1t (13)
1k)1212)e
对于分多时段瞬时摄入酒精的情况可类似讨论。
c1(t)1c(t) 2
5.1.4酒是先瞬时喝D1,t1小时后又在T小时内喝了D2。
1、 当0tt1时,由模型二可得
q(t)FD1(1ek1t)
x(t)kFD1ktk1t(ee) k1k
kFD1(ektek1t) (k1k)Vc(t)
2、 当t1tT时 可知胃内的酒精量为:tt1FD2FD1q(t) Tt1
tt1则有 r1(t)k1FD2FD1q(t)
Tt1
tt1dq于是有 k1FD2FD1q(t),q(t1)FD1(1ek1t1) dtTt1
解得 q(t)FD2FD2FD2FD2 t(ek1t1FD1)ek1tFD1t1Tt1k1(Tt1)Tt1Tt1
令OFD2FD2FD2FD2化简得 ,Pek1t1FD1,QFD1t1Tt1k1(Tt1)Tt1Tt1
q(t)OtPek1tQ
从而有 r1(t)Ok1Pek1t
再由 kFD1kt1k1t1dxr1(t)kx(t),x(t1)(ee) dtk1k
可得 x(t)A(1ek(tt1))B(ek1tet1(kk1)kt)C(ektet(kk1)kt);
其中 AFD2kFD2FD2kFD,B1(FD1t1),C11. k(Tt1)kk1Tt1k1(Tt1)k1k
c(t)ABC(1ek(tt1))(ek1tet1(kk1)kt)(ektet(kk1)kt). VVV
3、 当tT时
体内未被吸收的酒精量为:FD2FD1q(t)
有 dqk1FD2FD1q(t),q(T)OTPek1TQ dt
解得 q(t)(FD1FD2)(1ek1(tT))(OTPek1TQ)ek1(tT)
进而解得 x(t)(k1(FD1FD2q(T))k1(tT)k(tT))(ee)x(T)ek(tT); kk1
q(T)OTPek1TQ
k(Tt1)其中 )B(ek1Tet1(kk1)kT)C(ekTeT(kk1)kT) x(t)A(1e
.
所以 c(t)(k1(FD1FD2q(T))k1(tT)k(tT)x(T)k(tT))(ee)e; (kk1)VV
5.2模型求解
由假设知啤酒的酒精浓度为5%,规格为600ml,密度为0.9g/ml,可计算得一瓶啤酒的酒精含量m=5%*600ml*0.9g/ml=27g,那么两瓶啤酒的酒精含量D为54g。运用题目所给的体重约70kg的某人在短时间内喝下2瓶啤酒后,隔一定时间测量的血液
FD中酒精含量参考数据,假设取初始值k=0,k1=1,=110,将这些值根据非线性V
拟合的最小二乘法,使用Matlab软件可以求得:
FDk=0.1949,k1=1.8887,=107.6807 V
并得出数据建立的拟合曲线如下图1:
90
80
70
60
50
40
30
20
10
[1**********]416
图1 样本数据酒精含量拟合图
那么得出模型2中的c(t)的函数表达式:c(t)120.0712(e0.1949te1.8887t)
5.2.1问题1的解答
根据样本数据已经计算得出短时间内喝完2瓶啤酒的酒浓度
c(t)120.0712(e0.1949te1.8887t),那么可以得出短时间内喝完1瓶啤酒的酒精浓度c1(t)60.0356(e0.1949te1.8887t),我们假设大李的体重仍为70kg,并且大李在中午12点是瞬时喝下,中午12点至下午6点经过6个小时可以计算得出大李血液中的酒精浓度c(6)=18.6435毫克/百毫升
精浓度c(t)= 16.5467毫克/百毫升
c(t)=20.0295 >20毫克/百毫升的国家新标准,不能通过检查。所以得出如果李强在晚
上7:15之前吃晚饭,则将能通过检查,而在7:15及之后吃晚饭,则将不能通过检查。而题目中是说李强没有通过检查,那么李强应该是在7:15之后吃的晚饭。 5.2.2问题2的解答
我们假设王东的体重仍为70kg,他是在短时间内喝下3瓶啤酒的,根据模型2中c(t)的函数表达式ct
k1FD
(ektek1t),只需将系数扩大1.5倍即得该情况下
V(k1k)
的图象,由图2所示。从图2可看出当c(t)=20时该值在时间坐标轴上有2个交点,利用Matlab软件解出2个交点的值为t1=0.0705,t2=11.2766。当c(t)=80时该值在时间坐标轴上也有2个交点,同理解出2个交点的值为t3=0.3847,t4=4.1593。很明显,在t1~t2(0.0705~11.2766)这段时间内,c(t)³20,故在这段时间内驾车属于饮酒驾车, 在t3~t4(0.3847~4.1593)这段时间内,c(t)³80,故在这段时间内驾车属于醉酒驾车。
140
120
100
80
60
40
20
[1**********]16
图2 瞬时喝下3瓶啤酒的曲线拟合图
如果王东在半个小时内喝了半斤37度白酒,根据模型1中2个时间段内c(t)的函数表达式c(t)
FD11k1tktkt(ee)(1e)即可得该情况下的函数图象,由
Tkk1k
图3所示。
从图3可看出当c(t)=20时该值在时间坐标轴上有2个交点,利用Matlab软件解出2个交点的值为t1=0.2813,t2=11.5287。当c(t)=80时该值在时间坐标轴上也有2个交点,同理解出2个交点的值为t3=0.6585,t4=4.4112。很明显,在
t1~t2(0.2813~11.5287)这段时间内,c(t)³20,故在这段时间内驾车属于饮酒驾车,在t3~t4(0.6585~4.4112)这段时间内,c(t)³80,故在这段时间内驾车属于醉酒驾车。
140
120
100
80
60
40
20
[1**********]16
图3 半小时内喝下半斤37度白酒的曲线拟合图 5.2.3问题3的解答
王东几乎一口气喝了三瓶啤酒,休息一小时后,又在两个小时内喝了一斤37度白酒,根据模型4中3个时间段内c(t)的各个函数表达式即可得该情况下的函数图象,由图4所示。从图4可看出当c(t)=20时该值在时间坐标轴上有2个交点,利用Matlab软件解出2个交点的值为t1=0.0705,t2=19.2036。当c(t)=80时该值在时间坐标轴上也有2个交点,同理解出2个交点的值为t3=0.3847,t4=12.0907。很明显,
c(t)³20,在t1~t2(0.0705~19.2036)这段时间内,故在这段时间内驾车属于饮酒驾车,
在t3~t4(0.3847~12.0907)这段时间内,故在这段时间内驾车属于醉酒驾车。 c(t)³80,
400350
[***********]
[1**********]0
图4 问题3的曲线拟合图 5.2.4问题4的解答
ekT1lnk1T
由模型1知问题2中在王东半小时内喝下半斤白酒的情况下,当t2时,
kk1
体液中酒精浓度c(t) 达到最大值,记为 Tm1,计算得Tm1=1.6121。
由模型2知问题2中在王东瞬时喝下三瓶啤酒的情况下,当t1
k1
ln1时 k1kk
此时体液中酒精浓度c(t) 达到最大值,记为 Tm2,计算得Tm2=1.3409。
5.2.5问题5的解答
根据对上述模型求解的分析,得出以下结论:(1)每天不论短时间内还是长时间内饮过量的酒,都不能在很长一个时间中恢复标准,问题2中的两个数据:在饮酒11.2766小时内血液中酒精浓度大于20mg/dml,11.5287小时内血液中酒精浓度大于20mg/dml,这样就不能驾驶。(2)即使长时间内均匀时间段中喝少量的酒,人体血液中的酒精的含量也会积少成多,超过标准。因此,我们要做到适量、适时的喝酒。 5.2.6问题6的解答
给想喝酒的驾驶员的忠告
酒后驾车是导致交通事故的重要危险因素之一。目前我国在饮酒与交通安全方面的
形式十分严峻,据我国公安部统计,近10年来因饮酒所导致的道路交通事故,人员伤亡及经济损失仍逐年增加,2002年因饮酒所导致的道路事故数、死亡人数、受伤人数和经济损失分别达到1996年的262.3%,184.4%,325.3%和144.3%。
针对这种严重的道路交通情况,国家质量监督检验检疫局2004年5月31日发布了新的车辆驾驶人员血液、呼气酒精含量阈值与检验》国家标准,新标准对车驾驶人员血液中的酒精含量加强了限制,要求更为严格。
从生物医药学角度看,少量的酒精能增加血液里的胆固醇而使得冠状动脉不易形成血栓,这是酒精的抗血栓药性。然而酒精是一种抑制性药物,会影响人的整个神经系统,影响到人的大脑思维和行动。饮酒驾车十分危险!在现实生活中,由于种种原因人们难免有时饮酒,然而饮酒驾车新标准的发布对人们的酒后驾车提出了更为严格的要求,这让不少人困惑不已,那么有没有两全其美的办法呢?我们通过对血液中酒精浓度的变化情况进行科学的分析和研究,告诉大家如何饮酒而不会导致饮酒驾车,同时也告诉大家关于人体内酒精浓度的一些数学奥秘。
1) 在摄入同样酒精量的前提下,瞬时喝酒比慢速平稳喝酒更能让人体内的酒精在短时间内降到较低的水平。
2) 一个体重为70kg的驾车人员在瞬时喝下1瓶640ml,酒精度为3.5%的啤酒后,必须在喝完酒的6小时之后才能通过交通安检,安全地开车。
3) 饮酒后,血液系统中达到酒精含量最大值的时间与酒精摄入量无关人体血液中的酒精含量曲线为单峰曲线,即只存在一个极大值。
4) 对一个正常的驾驶员而言,是容许天天喝酒的,因为它存在一个稳态,但是还是受时间和量的限制。
5) 不管是慢速喝酒,还是快速喝酒,人体内达到酒精浓度最大的时间与酒精摄入量无关,只与酒精在体内的扩散速度有关,即只与人的体质有关。
除此之外,在各种应酬中还应掌握几个原则:饮酒要适量,以舒服为原则,不能逞英雄;身体不适,情绪不良尽量避免饮酒,以免加重症状。
司机一滴泪,亲人两行泪!司机朋友们,为了你的安全,他人的安全,更为了让你的亲人放心,请你们健康的饮酒!安全的驾车!
6 模型的改进与推广
本文对人在饮酒后血液中的酒精含量的变化过程建立了基本模型。这个模型考虑很全面 ,在理论上较好的描述了饮酒后血液中酒精含量的计算公式 ,并运用微积分思想解决了人在不同的情况下饮酒后血液中的酒精含量的计算,模型都可以进行合理的分析,由于没有收集资料, 条件所限 ,未进行大量数据验证。所以没有很好的了解人饮酒后酒精在人体中变化的情况 ,只是做了简单的假设,从而引起误差。这个模型还可以得出人在喝一定量的酒后血液中酒精含量和时间的关系,从而司机可根据这个关系来判断饮酒后安全驾车的时间。如果考虑人在一天的不同时刻 ,其体内的新陈代谢速度不同 ,或考虑饮酒者的个体差异 ,可以对模型做进一步的优化和改进,可以得出更好的结果。比如人在饮酒时及饮酒后的呼吸过程和人的各种器官,如肝脏对酒精的降解过程以及排泄过程都可以减少血液中的酒精含量 。 参考文献:
【1】 刘来福,曾文艺 数学模型与数学建模 北京 北京师范大学出版社,2002 【2】 陈瀚珠,等 实用内科学 北京 人民卫生出版社 ,1997 【3】 http://finance sina com.cn/6/20030123/17263060.shtml 【4】 陈荣江,孙用明,张万琴等 饮酒后血液中酒精含量模型的建立与初步验证 2004
附件:
血液酒精浓度检测
Blood alcohol concentration measurement
以下是maple对建立的模型进行求解的代码 > #自定义函数
BAC:=proc(eqn1,eqn2,init1,init2,result) local sols,c1,c2,t0;
sols:=dsolve({eqn1,eqn2,init1,init2},result); c1:=simplify(op(2,sols))/V; c2:=c(t)=rhs(c1);
t0:=solve(rhs(diff(c2,t))=0,t); print (simplify(op(1,sols)));#x0 print (simplify(op(2,sols)));#x1 print (c2);#C
if t00 then print (t=t0) end if;#t end proc:
> #1.1在T小时内喝完,t
A:=k1*(t/T*F*D-q(t)): B:=k*x(t):
eqn1:=diff(x(t),t)=A-B: eqn2:=diff(q(t),t)=A: #定义初始值 init1:=x(0)=0: init2:=q(0)=0:
#调用自定义函数,得出解析解
BAC(eqn1,eqn2,init1,init2,{q(t),x(t)});
> #1.2在T小时内喝完,t>=T时 #定义方程
A:=k1*(F*D-q(t)): B:=k*x(t):
eqn1:=diff(x(t),t)=A-B: eqn2:=diff(q(t),t)=A: #定义初始值
init1:=q(T)=D*F*(-1+k1*T+exp(-k1*T))/(k1*T):
init2:=x(T)=D*F*(-k*exp(k*T)+exp(k*T)*k1+exp(-T*(-k+k1))*k-k1)*exp(-k*T)/(T*k*(-k+k1)):
#调用自定义函数,得出解析解
BAC(eqn1,eqn2,init1,init2,{q(t),x(t)});
> #2.1瞬间喝完 #定义方程
A:=k1*(F*D-q(t)): B:=k*x(t):
eqn1:=diff(x(t),t)=A-B: eqn2:=diff(q(t),t)=A: #定义初始值 init1:=q(0)=0: init2:=x(0)=0:
#调用自定义函数,得出解析解
BAC(eqn1,eqn2,init1,init2,{q(t),x(t)});
> #t->infinity
eqn3:=diff(x(t),t)=-k*x(t);
sol:=dsolve({eqn3,x(0)=C},x(t));
> #3.1分段喝完,t
A:=k1*(F*D[1]-q(t)): B:=k*x(t):
eqn1:=diff(x(t),t)=A-B: eqn2:=diff(q(t),t)=A: init1:=q(0)=0: init2:=x(0)=0:
BAC(eqn1,eqn2,init1,init2,{q(t),x(t)});
> #3.2分段喝完,t>=T[0]
A:=k1*(F*(D[1]+D[2])-q(t)): B:=k*x(t):
eqn1:=diff(x(t),t)=A-B: eqn2:=diff(q(t),t)=A:
init1:=q(T[0])=F*D[1]*(1-exp(-k1*T[0])):
init2:=x(T[0])=-k1*F*D[1]*(exp(-T[0]*(-k+k1))-1)*exp(-k*T[0])/(-k+k1): BAC(eqn1,eqn2,init1,init2,{q(t),x(t)});
二、模型求解编码: 1.最小二乘法拟合 Alcohol.m 文件代码
tdata=[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]; cdata=[30 68 75 82 82 77 68 68 58 51 50 41 38 35 28 25 18 15 12 10 7 7 4]; k0=[1,0,110];
options=optimset('maxfunevals',1e4); k=lsqcurvefit(@myfun,k0,tdata,cdata) t=0:0.05:16; c=myfun(k,t);
plot(tdata,cdata,'+r',t,c);
myfun.m 函数代码
function c=myfun(k,t)
c=k(3)*k(1)/(k(1)-k(2))*(exp(-k(2)*t)-exp(-k(1)*t));
2.问题1的代码
Alcohol1.m 文件代码 t=6; t1=14; k=0.1949; k1=1.8887; a=107.6807/2;
c=a*k1/(k1-k)*(exp(-k*t)-exp(-k1*t))
c1=a*k1/(k-k1)*((exp(k1*t)+1)*(exp(-k1*t1)-exp(t*(k-k1))*exp(-k*t1))-exp(-k*t1)+...
exp((k-k1)*t)*exp(-k*t1))
说明:在令t=7.25可求7点十五分的解
3.问题2的代码 (1)短时间
Alcohol2.m 文件代码
t=fsolve('myfun2',10,optimset('Display','off'))
myfun2.m 函数代码 function c=myfun2(t) k=0.1949; k1=1.8887; a=107.6807/2;
c=3*a*k1/(k1-k)*(exp(-k*t)-exp(-k1*t))-20;
说明:可通过alcohol2.m文件中的初值t=10或1,和myfun2.m文件c中20变为80求出全部结果。
(2)半小时
Alcohol2.m 文件代码
t=fsolve('myfun1',10,optimset('Display','off'))
myfun1.m 函数代码 function c=myfun1(t) T=0.5; k=0.1949; k1=1.8887; a=107.6807/2;
N=3*a/(k1*T)*(1-exp(k1*T));
M=-3*a*k1*(exp(k*T)-1)/(T*k*(k-k1));
c=-k1*N/(k-k1)*exp(-k1*t)+M*exp(-k*t)-20; 说明:思想和上面相同
4.问题3代码
Alcohol2.m 文件代码
t=fsolve('myfun3',10,optimset('Display','off'))
myfun3.m 函数代码 function c=myfun3(t) k=0.1949; k1=1.8887; a=107.6807/2; t1=1; T=3;
O=6*a/(T-t1);
P=6*a/(T-t1)/k1*exp(k1*t1)-3*a; Q=3*a-6*a/(T-t1)*t1-6*a/k1/(T-t1); qtv=O*T+P*exp(-k1*T)+Q; ct=myfun4(T);
c=(9*a-qtv)*k1/(k-k1)*(exp(-k1*(t-T))-exp(-k*(t-T)))+ct*exp(-k*(t-T));
myfun4.m 函数代码 function c=myfun4(t) k=0.1949; k1=1.8887; a=107.6807/2; t1=1; T=3;
c=6*a/k/(T-t1)*(1-exp(-k*(t-t1)))-k1/(k-k1)*(3*a-6*a/(T-t1)*t1-6*a/k1/(T-t1))*...
(exp(-k1*t)-exp(t1*(k-k1)-k*t))+3*a*k1/(k1-k)*(exp(-k*t)-exp(t*(k-k1)-k*t));