数学建模笔记CH2:微积分方法建模

第二章:微积分方法建模

overview:

涉及连续的变量,可以用微积分求解,求得解析式便于下一步分析。有些离散的变量也可以演变成连续变量进行求解。当我们描述实际对象的某些特性随时间(或空间)而演变的过程,分析它的变化规律,预测它的未来性态时,通常要建立对象的动态模型。建模时首先要根据建模目的和对问题的具体分析作出简化假设,然后按照对象内在的或可以类比的其它对象的规律列出微分方程,求出方程的解并将结果翻译回实际对象,就可以进行描述、分析或预测了。

2.1飞机降落曲线

飞机降落时,其曲线是一条三次抛物线,如下图:

水平速度为常数u,出于安全考虑,飞机的垂直加速度不得超过g/10,已知飞行高度h,要求在坐标原点降落,求开始下降点x0允许的最小值。

求曲线解析式

设飞机降落曲线为:

y = a x 3 + b x 2 + c x + d y = ax^3+bx^2+cx+d y=ax3+bx2+cx+d

依据假设可知:
y ( 0 ) = 0 , y ( x 0 ) = h , y ′ ( 0 ) = 0 , y ′ ( x 0 ) = 0 y(0)=0,y(x_0)=h,y'(0)=0,y'(x_0)=0 y(0)=0,y(x0)=h,y(0)=0,y(x0)=0,带入上式可得:

{ y ( 0 ) = d = 0 , y ′ ( 0 ) = c = 0 , y ( x 0 ) = a x 3 + b x 2 + c x + d = h , y ′ ( x 0 ) = 3 a x 2 + 2 b x 0 + c = 0 x = 0 \small \begin{cases} y(0)=d=0,\\ y'(0)=c=0,\\ y(x_0)=ax^3+bx^2+cx+d=h,\\ y'(x_0)=3ax^2+2bx_0+c=0x=0\\ \end{cases} y(0)=d=0,y(0)=c=0,y(x0)=ax3+bx2+cx+d=h,y(x0)=3ax2+2bx0+c=0x=0
所以:
a = − 2 h x 0 3 , b = 3 h x 0 2 , c = 0 , d = 0 a = -\frac{2h}{{x_0}^3},b=\frac{3h}{{x_0}^2},c=0,d=0 a=x032h,b=x023h,c=0,d=0

其曲线为:
y = − h x 0 2 ( 2 x 0 x 3 − 3 x 2 ) y=-\frac{h}{{x_0}^2}(\frac{2}{x_0}x^3-3x^2) y=x02h(x02x33x2)

求最佳着陆点

飞机垂直速度是y关于时间t的导数:
d y d t = − h x 0 2 ( 6 x 0 x 2 − 6 x ) d x d t \frac{dy}{dt}=-\frac{h}{{x_0}^2}(\frac{6}{x_0}x^2-6x)\frac{dx}{dt} dtdy=x02h(x06x26x)dtdx

其中, d x d t \frac{dx}{dt} dtdx是飞机水平速度,所以 d x d t = u \frac{dx}{dt}=u dtdx=u,则:
d y d t = − 6 h u x 0 2 ( x 2 x 0 − x ) \frac{dy}{dt}=-\frac{6hu}{{x_0}^2}(\frac{x^2}{x_0}-x) dtdy=x026hu(x0x2x)

垂直加速度为: d 2 y d t 2 = − 6 h u x 0 2 ( 2 x x 0 − 1 ) d x d t = − 6 h u 2 x 0 2 ( 2 x x 0 − 1 ) \frac{d^2y}{dt^2}=-\frac{6hu}{{x_0}^2}(\frac{2x}{x_0}-1)\frac{dx}{dt}=-\frac{6hu^2}{x_0^2}(\frac{2x}{x_0}-1) dt2d2y=x026hu(x02x1)dtdx=x026hu2(x02x1)

a ( x ) = d 2 y d t 2 a(x)=\frac{d^2y}{dt^2} a(x)=dt2d2y,则:

∣ a ( x ) ∣ = 6 h u 2 x 0 2 ∣ 2 x x 0 − 1 ∣ , x ∈ [ 0 , x 0 ] |a(x)|=\frac{6hu^2}{{x_0}^2}|\frac{2x}{x_0}-1|,x\in[0,x_0] a(x)=x026hu2x02x1,x[0,x0]

所以,垂直加速度最大值为:

m a x ∣ a ( x ) ∣ = 6 h u 2 x 0 2 , x ∈ [ 0 , x 0 ] max|a(x)|=\frac{6hu^2}{{x_0}^2},x\in[0,x_0] maxa(x)=x026hu2,x[0,x0]

由于题目限制: m a x ∣ a ( x ) ∣ ≤ g 10 max|a(x)|\leq\frac{g}{10} maxa(x)10g,所以:

x 0 ≥ u 60 h g x_0\geq u\sqrt{\frac{60h}{g}} x0ug60h ,此为 x 0 x_0 x0的最小值。

2.2经济增长模型

发展经济、增加生产有两个重要因素,一是增加投资(扩大厂房、购买设备、技术革新等),二是增加劳动力。恰当调节投资增长和劳动力增长的关系,使增加的产量不致被劳动力的增长抵消,劳动生产率才能不断提高,从而真正起到发展经济的作用。为此,需要分析产量、劳动力和投资之间变化规律,从而保证经济正常发展。
假设:
Q ( t ) Q(t) Q(t)——地区、部门、企业在某时间t的产量
L ( t ) L(t) L(t)——地区、部门、企业在某时间t的劳动力
K ( t ) K(t) K(t)——地区、部门、企业在某时间t的资金
Z ( t ) Z(t) Z(t)——每个劳动力在某时间t所占有的产量(劳动力生产率)

道格拉斯(Douglas)生产函数

由于我们只关心产量、劳动力、投资的相对增长,而不是绝对量,所以定义:
产量指数: i Q ( t ) = Q ( t ) Q ( 0 ) i_Q(t) = \frac{Q(t)}{Q(0)} iQ(t)=Q(0)Q(t),劳动力指数: i L ( t ) = L ( t ) L ( 0 ) i_L(t) = \frac{L(t)}{L(0)} iL(t)=L(0)L(t),投资指数: i K ( t ) = K ( t ) K ( 0 ) i_K(t) = \frac{K(t)}{K(0)} iK(t)=K(0)K(t) (1)
令: ξ ( t ) = ln ⁡ i L ( t ) i K ( t ) \xi(t)=\ln\frac{i_L(t)}{i_K(t)} ξ(t)=lniK(t)iL(t), ψ ( t ) = ln ⁡ i Q ( t ) i K ( t ) \psi(t)=\ln\frac{i_Q(t)}{i_K(t)} ψ(t)=lniK(t)iQ(t) (2)
根据统计数据,散点 ( ξ , ψ ) (\xi,\psi) (ξ,ψ)在直角坐标系下的图像大致为:
在这里插入图片描述
由图可知,大多数的点处在一条经过原点的直线附近,故 ξ \xi ξ ψ \psi ψ有如下关系:

ψ = γ ξ ( 0 < γ < 1 ) \psi=\gamma\xi (0<\gamma<1) ψ=γξ(0<γ<1) (3)
带入上式可得:
i Q ( t ) = i L γ ( t ) i K 1 − γ ( t ) i_Q(t)={i_L}^\gamma(t){i_K}^{1-\gamma}(t) iQ(t)=iLγ(t)iK1γ(t) (4)
a = Q ( 0 ) L − γ ( 0 ) K 1 − γ ( 0 ) a=Q(0)L^{-\gamma}(0)K^{1-\gamma}(0) a=Q(0)Lγ(0)K1γ(0),则由(1)和(4)可知:
Q ( t ) = a L γ ( t ) K 1 − γ ( t ) Q(t)=aL^\gamma(t)K^{1-\gamma}(t) Q(t)=aLγ(t)K1γ(t) (5)
式(5)就是经济学中著名的Douglas生产函数,它表明产量余劳动力、投资之间的关系,即:
Q ˙ Q = γ L ˙ L + ( 1 − γ ) K ˙ K \frac{\dot{Q}}{Q}=\gamma\frac{\dot{L}}{L}+(1-\gamma)\frac{\dot{K}}{K} QQ˙=γLL˙+(1γ)KK˙ (6)
其中 Q ˙ \dot{Q} Q˙ L ˙ \dot{L} L˙ K ˙ \dot{K} K˙表示Q、L、K关于t的导数。
这个函数表明,相对增长量, Q ˙ Q \frac{\dot{Q}}{Q} QQ˙ L ˙ L \frac{\dot{L}}{L} LL˙ K ˙ K \frac{\dot{K}}{K} KK˙之间呈线性关系,当 γ → 1 \gamma\rightarrow1 γ1时产量增长主要依靠劳动力增长,当 γ → 0 \gamma\rightarrow0 γ0时产量增长主要依靠投资。亦称 γ \gamma γ产量对劳动力的弹性系数

劳动生产率增长条件

定义劳动生产率 Z ( t ) = Q ( t ) L ( t ) Z(t)=\frac{Q(t)}{L(t)} Z(t)=L(t)Q(t),则 Z ˙ Z = Q ˙ Q − L ˙ L \frac{\dot{Z}}{Z}=\frac{\dot{Q}}{Q}-\frac{\dot{L}}{L} ZZ˙=QQ˙LL˙,将(6)代入 Z ( t ) Z(t) Z(t)可得:
Z ˙ Z = ( 1 − γ ) ( K ˙ K − L ˙ L ) \frac{\dot{Z}}{Z}=(1-\gamma)(\frac{\dot{K}}{K}-\frac{\dot{L}}{L}) ZZ˙=(1γ)(KK˙LL˙) (7)
由(7)可见,只要 K ˙ K > L ˙ L \frac{\dot{K}}{K}>\frac{\dot{L}}{L} KK˙>LL˙就能保证 Z ˙ Z > 0 \frac{\dot{Z}}{Z}>0 ZZ˙>0,即劳动生产率的提高需要由投资的相对增长大于劳动力的相对增长为前提条件

2.3存贮模型

原料、商品的存贮问题,存的太多,存贮费用过高;存的太少,无法及时满足需求。目的:制定最有存贮策略,即:多长时间顶一次货,每次顶多少货,才能使总费用最小。

模型一

模型假设
  • 每次订货费用为 C 1 C_1 C1,每天每吨货物贮存费用为 C 2 C_2 C2(已知)
  • 每天货物的需求量r吨为已知
  • 订货周期为T天,每次订货Q吨,当贮存量降到0时订货立即到达。
模型建立

订货周期、订货量余每天的需求量存在关系:
Q = r T Q=rT Q=rT
订货以后贮存量 q ( t ) q(t) q(t)均匀地下降,即 q ( t ) = Q − r t q(t)=Q-rt q(t)=Qrt,如下图:
在这里插入图片描述
则一个订货周期的费用:
{ 订 货 费 : C 1 , 存 贮 费 : ∫ 0 T q ( t )   d t = 1 2 C 2 Q T = 1 2 C 2 r T 2 \small \begin{cases} 订货费:C_1,\\ 存贮费:\int_{0}^{T}q(t)\,dt=\frac{1}{2}C_2QT=\frac{1}{2}C_2rT^2\\ \end{cases} {:C1,:0Tq(t)dt=21C2QT=21C2rT2
即: C ( T ) = C 1 + 1 2 C 2 r T 2 C(T)=C_1+\frac{1}{2}C_2rT^2 C(T)=C1+21C2rT2
则一个订货周期每天的费用为: C ‾ ( T ) = C 1 T + 1 2 C 2 r T \overline{C}(T)=\frac{C_1}{T}+\frac{1}{2}C_2rT C(T)=TC1+21C2rT

模型求解

d C ‾ d T = 0 \frac{d\overline{C}}{dT}=0 dTdC=0,可求得: T = 2 c 1 R C 2 T=\sqrt{\frac{2c_1}{RC_2}} T=RC22c1
进而: Q = 2 C 1 r C 2 Q=\sqrt{\frac{2C_1r}{C_2}} Q=C22C1r
此模型称为经济订货批量公式,简称EOQ公式

模型二 允许缺货的存贮模型

模型假设
  • 每次订货费用为 C 1 C_1 C1,每天每吨货物贮存费用 C 2 C_2 C2(已知)
  • 每天的货物需求量r吨(已知)
  • 订货周期为T天,订货量Q,允许缺货,每天每吨货物缺货费 C 3 C_3 C3(已知)
模型建立

缺货时的存储量q视为负值,则 q ( t ) q(t) q(t)的图像如下,货物在 t = T 1 t=T_1 t=T1时用完,于时 Q = r T 1 Q=rT_1 Q=rT1
在这里插入图片描述
则一个订货周期内总费用为:
{ 订 货 费 : C 1 , 存 贮 费 : ∫ 0 T 1 q ( t )   d t = 1 2 C 2 Q T 1 = 1 2 C 2 Q 2 r , 缺 货 费 : ∫ T 1 T ∣ q ( t ) ∣   d t = C 3 2 r ( T − T 1 ) 2 = C 3 2 r ( r T − Q ) 2 \small \begin{cases} 订货费:C_1,\\ 存贮费:\int_{0}^{T_1}q(t)\,dt=\frac{1}{2}C_2QT_1=\frac{1}{2}C_2\frac{Q^2}{r},\\ 缺货费:\int_{T_1}^{T}|q(t)|\,dt=\frac{C_3}{2}r(T-T_1)^2=\frac{C_3}{2r}(rT-Q)^2 \end{cases} :C1,:0T1q(t)dt=21C2QT1=21C2rQ2,:T1Tq(t)dt=2C3r(TT1)2=2rC3(rTQ)2
即: C ( T , Q ) = C 1 + 1 2 C 2 Q 2 1 r + 1 2 r C 3 ( r T − Q ) 2 C(T,Q)=C_1+\frac{1}{2}C_2Q^2\frac{1}{r}+\frac{1}{2r}C_3(rT-Q)^2 C(T,Q)=C1+21C2Q2r1+2r1C3(rTQ)2
则平均每天的费用为:
C ‾ ( T , Q ) = C ( T , Q ) T = C 1 T + C 2 Q 2 2 r T + C 3 ( r T − Q ) 2 2 r T \overline{C}(T,Q)=\frac{C(T,Q)}{T}=\frac{C_1}{T}+\frac{C_2Q^2}{2rT}+\frac{C_3(rT-Q)^2}{2rT} C(T,Q)=TC(T,Q)=TC1+2rTC2Q2+2rTC3(rTQ)2

模型求解

{ ∂ C ‾ ∂ T = 0 ∂ C ‾ ∂ Q = 0 \small \begin{cases} \frac{\partial\overline{C}}{\partial T}=0\\ \frac{\partial\overline{C}}{\partial Q}=0 \end{cases} {TC=0QC=0
求出T、Q的最优解,分别记为 T ′ 、 Q ′ T'、Q' TQ
T ′ = 2 C 1 r C 2 C 2 + C 3 C 3 , Q ′ = 2 C 1 r C 2 C 3 C 2 + C 3 T'=\sqrt{\frac{2C_1}{rC_2}\frac{C_2+C_3}{C_3}},Q'=\sqrt{\frac{2C_1r}{C_2}\frac{C_3}{C_2+C_3}} T=rC22C1C3C2+C3 ,Q=C22C1rC2+C3C3

分析

μ = C 2 + C 3 C 3 \mu=\sqrt{\frac{C_2+C_3}{C_3}} μ=C3C2+C3 ,与模型一相比,有:
T ′ = μ T , Q ′ = Q μ T'=\mu T,Q'=\frac{Q}{\mu} T=μT,Q=μQ
显然,T’ > T,Q’ < Q,即在允许缺货时应增大订货周期,减少订货批次;当缺货非 C 3 C_3 C3相对于存贮费 C 2 C_2 C2而言越大时, μ \mu μ越小,T’和Q’越接近T和Q。

2.4城市人口统计模型

模型一(估算城市现有人口)

城市人口密度常用 P ( r ) = b r 2 + a P(r)=\frac{b}{r^2+a} P(r)=r2+ab或者 P ( r ) = a e − b r ( a , b > 0 ) P(r)=ae^{-br}(a,b>0) P(r)=aebr(a,b>0)来近似表示,其中r是距城中心的距离。则计算距离市中心C区域内的人口数N可以这样:从城市中心画一条射线,把这条线上从0到C之间n等分,每小区间长度为 Δ r \Delta r Δr,每小区间确定一个圆环,第j个圆环面积为从城市中心为:
π r j 2 − π r j − 1 2 = π r j 2 − π ( r j − Δ r ) 2 = 2 π r j Δ r − π ( Δ r ) 2 ≈ 2 π r j Δ r \pi{r_j}^2-\pi{r_{j-1}}^2=\pi{r_j}^2-\pi(r_j-\Delta r)^2=2\pi r_j\Delta r-\pi(\Delta r)^2\approx2\pi r_j \Delta r πrj2πrj12=πrj2π(rjΔr)2=2πrjΔrπ(Δr)22πrjΔr ( Δ r 很 小 ) (\Delta r很小) (Δr)

第j个圆环上的人口数近似为 P ( r j ) 2 π r j Δ r P(r_j)2\pi r_j \Delta r P(rj)2πrjΔr,所以:

N ≈ ∑ j = 1 n P ( r j ) 2 π r j Δ r N\approx \sum_{j=1}^{n}P(r_j) 2\pi r_j \Delta r Nj=1nP(rj)2πrjΔr

n → ∞ n\rightarrow\infty n时,

N = ∫ 0 C P ( r ) 2 π r   d r N=\int_{0}^{C}P(r)2\pi r\,dr N=0CP(r)2πrdr

模型二(预测城市未来人口)

P(t)表示t时刻的城市人口数量,人口变化手下面规则的影响:

  • t时刻的净增长人口以每年r(t)的比率增加
  • 在一段世界内,由于自然死亡和人口迁移, T 1 T_1 T1时刻的人口数 P ( T 1 ) P(T_1) P(T1)的一部分在 T 2 T_2 T2时刻仍然存在,用 h ( T 2 − T 1 ) P ( T 1 ) h(T_2-T_1)P(T_1) h(T2T1)P(T1)来表示,这里 0 < h ( T 2 − T 1 ) < 1 0<h(T_2-T_1)<1 0<h(T2T1)<1 T 2 − T 1 T_2-T_1 T2T1是这段时间的长度。

把(0,T]时间划分为n等分,每个小区间长度为 Δ t \Delta t Δt
假设初始时刻人口为P(0),到时刻T将只剩h(T)P(0)。当 Δ t \Delta t Δt很小的时候,从 t j − 1 t_{j-1} tj1 t j t_j tj,净增长的人口比率近似为常数 r ( t j ) r(t_j) r(tj)。这段时间净增的人口数近似为 r ( t j ) Δ t r(t_j)\Delta t r(tj)Δt,在 t j t_j tj时刻的人口到T时刻只剩 h ( T − t j ) r ( t j ) Δ t h(T-t_j)r(t_j)\Delta t h(Ttj)r(tj)Δt。所以在T时刻的总人口数近似为:

P ( T ) ≈ h ( T ) P ( 0 ) + ∑ j = 1 n h ( T − t j ) r ( t j ) Δ t P(T)\approx h(T)P(0)+\sum_{j=1}^{n}h(T-t_j)r(t_j)\Delta t P(T)h(T)P(0)+j=1nh(Ttj)r(tj)Δt

n → ∞ n\rightarrow\infty n,得:

P ( T ) = h ( T ) P ( 0 ) + ∫ 0 T h ( T − t ) r ( t )   d t P(T)=h(T)P(0)+\int_{0}^{T}h(T-t)r(t)\,dt P(T)=h(T)P(0)+0Th(Tt)r(t)dt

2.7古生物年代确定

主要思路: C 14 C^{14} C14的半衰期为5568年,根据 C 14 C^{14} C14的蜕变减少量的变化来判断生物的死亡时间。

模型假设

  • 地球周围大气的 C 14 C^{14} C14可以认为不变,现代生物和古代生物体内的 C 14 C^{14} C14蜕变速度一致。
  • 由原子物理学, C 14 C^{14} C14的蜕变速度与该时刻的 C 14 C^{14} C14含量成正比。

模型建立

记, x ( t ) x(t) x(t)代表t时刻生物体内 C 14 C^{14} C14的含量。由假设可知: d x d t = − k x , k > 0 ( 1 ) \frac{dx}{dt}=-kx,k>0 (1) dtdx=kx,k>01
(1)的通解为 x ( t ) = C e − k t x(t)=Ce^{-kt} x(t)=Cekt,设生物死亡时间为 t 0 t_0 t0,代入可知, C = x 0 C=x_0 C=x0,于时:
x ( t ) = x 0 e − k t ( 2 ) x(t)=x_0e^{-kt} (2) x(t)=x0ekt(2)
C 14 C^{14} C14的半衰期为T,则有:
x ( T ) = x 0 2 ( 3 ) x(T)=\frac{x_0}{2} (3) x(T)=2x03
将(3)代入(2)可得 k = l n 2 T k=\frac{ln2}{T} k=Tln2,故可解得:
t = T l n 2 l n x 0 x ( t ) ( 4 ) t=\frac{T}{ln2}ln\frac{x_0}{x(t)} (4) t=ln2Tlnx(t)x04
由于 x 0 和 x ( t ) x_0和x(t) x0x(t)难以测量,故使用另一种方法求t:
对(2)两边求导:
x ′ ( t ) = − x 0 k e − k t = − k x ( t ) x'(t)=-x_0ke^{-kt}=-kx(t) x(t)=x0kekt=kx(t)
x ′ ( 0 ) = − k x ( 0 ) = − k x 0 x'(0)=-kx(0)=-kx_0 x(0)=kx(0)=kx0
两式相除,可得: x ′ ( t ) x ′ ( 0 ) = x 0 x ( t ) \frac{x'(t)}{x'(0)}=\frac{x_0}{x(t)} x(0)x(t)=x(t)x0,代入(4)可得:
t = T l n 2 l n x ′ ( 0 ) x ′ ( t ) t=\frac{T}{ln2}ln\frac{x'(0)}{x'(t)} t=ln2Tlnx(t)x(0)
所以只需要测出标本 C 1 4 C^14 C14的平均蜕变速度(单位:次/分钟),即 x ′ ( t ) x'(t) x(t),和现在的 C 1 4 C^14 C14的蜕变速度 x ′ ( 0 ) x'(0) x(0),即可求出t。

2.8预测人口的增长

模型一:Malthhus指数增长模型

模型假设

1.某国/地区再t时刻的人口数x(t)为连续可微函数。
2.人口的增长率r是常数,即,单位时间人口的增长量与当时的人口成正比

模型建立

x 0 x_0 x0为初始时刻的人口,即 x ( 0 ) = x 0 x(0)=x_0 x(0)=x0
则从 t t t t + Δ t t+\Delta t t+Δt内人口的增长量为:
x ( t + Δ t ) − x ( t ) = r x ( t ) Δ t x(t+\Delta t)-x(t)=rx(t)\Delta t x(t+Δt)x(t)=rx(t)Δt
可导出下面的微分方程:
{ d x d t = r x x ( 0 ) = x 0 \small \begin{cases} \frac{dx}{dt}=rx\\ x(0)=x_0 \end{cases} {dtdx=rxx(0)=x0
解得: x ( t ) = x 0 e r t , r > 0 x(t)=x_0e^{rt} ,r>0 x(t)=x0ert,r>0

模型二:Logistic阻滞增长模型

模型假设

1.同模型一
2.当人口增加到一定数量后,增长率随着人口的继续增长而逐渐减少,且r(x)为x的线性函数 r ( x ) = r − s x r(x)=r-sx r(x)=rsx,其中r为x=0时的增长率,成为固有增长率。
3.自然资源和环境条件所能容纳的最大人口数量为 x m x_m xm,称作最大人口容量

模型建立

x = x m x=x_m x=xm时,增长率为0,即 r ( x m ) = 0 r(x_m)=0 r(xm)=0,进而: S = r x m S=\frac{r}{x_m} S=xmr,所以: r ( x ) = r ( 1 − x x m ) r(x)=r(1-\frac{x}{x_m}) r(x)=r(1xmx)
其中的 r , x m r,x_m r,xm是根据人口统计数据确定的常数, x m x_m xm常由经验决定(模型缺点: x m x_m xm不易轻易地得到)。
模仿模型一可得:
{ d x d t = r ( 1 − x x m ) x x ( 0 ) = x 0 \small \begin{cases} \frac{dx}{dt}=r(1-\frac{x}{x_m})x\\ x(0)=x_0 \end{cases} {dtdx=r(1xmx)xx(0)=x0
解得: x ( t ) = x m 1 + ( x m x 0 ) e − r t x(t)=\frac{x_m}{1+(\frac{x_m}{x_0})e^{-rt}} x(t)=1+(x0xm)ertxm
在这里插入图片描述

2.9药物在体内地分布与排除

模型背景

药物进入机体后,在随血液输送到各器官和组织的过程中,不断地被吸收、分布、代谢,最终排出体外。药物在血液中的浓度( μ g / m v \mu g/mv μg/mv)称血药浓度。血药浓度的大小直接影响到药物的疗效,浓度太低不能达到预期的治疗效果,浓度太高又可能导致中毒、副作用太强或造成浪费。
因此,研究药物在体内吸收、分布和排除的动态过程,对于新药研制时剂量的确定、给药方案设计等药理学和临床医学的发展具有重要的指导意义和实用价值。
为了研究目的,将一个机体划分成若干个房室,每个房室是机体的一部分,比如中心室和周边室。在一个房室内药物呈均匀分布,而在不同的房室之间按一定规律进行转移。

模型假设

1.药物进入机体后,全部进入中心室(血液较丰富的心、肺、肾等器官和组织),中心室的容积在给药过程中保持不变;
2.药物从中心室排出体外,与排除的数量相比,药物的吸收可以忽略;
3.药物排除的速率与中心室的血药浓度成正比。

模型建立

f 0 ( t ) : 给 药 速 度 f_0(t):给药速度 f0(t)
c ( t ) : 中 心 室 血 药 浓 度 c(t):中心室血药浓度 c(t)
x ( t ) : 中 心 室 药 量 x(t):中心室药量 x(t)
V : 中 心 室 容 积 V:中心室容积 V
k : 排 除 速 率 系 数 k:排除速率系数 k
在这里插入图片描述

求各种给药方式下血药浓度变化情况

上述变量有如下关系:
x ˙ = f 0 ( t ) − k x \dot{x}=f_0(t)-kx x˙=f0(t)kx
x ˙ + k x = f 0 ( x ) \dot{x}+kx=f_0(x) x˙+kx=f0(x)
x ( t ) = V c ( t ) x(t)=Vc(t) x(t)=Vc(t)
可得: c ˙ ( t ) + k c ( t ) = f 0 ( t ) V ( 1 ) \dot{c}(t)+kc(t)=\frac{f_0(t)}{V} (1) c˙(t)+kc(t)=Vf0(t)1

  • 1.快速静脉注射(指数模型)
    给药量为D,则初始条件 c ( 0 ) = D V , f 0 ( t ) = 0 c(0)=\frac{D}{V},f_0(t)=0 c(0)=VD,f0(t)=0
    解得: c ( t ) = D V e − k t ( 2 ) c(t)=\frac{D}{V}e^{-kt} (2) c(t)=VDekt2
  • 2.恒速静脉注射
    设持续时间为 τ \tau τ,注射速率为 k 0 k_0 k0,则有:
    ( x ≤ t ≤ τ ) (x\le t\le\tau) (xtτ)时, f 0 ( t ) = k 0 f_0(t)=k_0 f0(t)=k0,初始条件为c(0)=0,
    ( t ≥ τ ) (t\ge\tau) (tτ)时, f 0 ( t ) = 0 f_0(t)=0 f0(t)=0,初始条件 c ( τ ) = k 0 V k ( 1 − e − k τ ) c(\tau)=\frac{k_0}{Vk}(1-e^{-k\tau}) c(τ)=Vkk0(1ekτ),所以(1)的解为:
    c ( t ) = { k 0 V k ( 1 − e − k t ) , 0 ≤ t ≤ τ k 0 V k ( 1 − e − k τ ) e − k ( t − τ ) , t ≥ τ ( 3 ) \small c(t)= \begin{cases} \frac{k_0}{Vk}(1-e^{-kt}) ,0\le t\le\tau \\ \frac{k_0}{Vk}(1-e^{-k\tau})e^{-k(t-\tau)},t\ge\tau (3) \end{cases} c(t)={Vkk0(1ekt),0tτVkk0(1ekτ)ek(tτ),tτ3
    在这里插入图片描述
  • 3.口服或肌肉注射
    在药物输入中心室之前先有一个将药物吸入血液的过程,可以看作有一个吸收室,药物由吸收室进入中心室,药物由吸收室进入中心室额转移速率系数记为 k 1 k_1 k1,给药量D,吸收室药量 x 0 ( t ) x_0(t) x0(t)。则有:
    在这里插入图片描述
    { x 0 ˙ = − k 1 x 0 x 0 ( t ) = D \small \begin{cases} \dot{x_0}=-k_1x_0\\ x_0(t)=D \end{cases} {x0˙=k1x0x0(t)=D
    上式可推出: x 0 ( t ) = D e − k 1 t x_0(t)=De^{-k_1t} x0(t)=Dek1t
    于时 f 0 ( t ) = k 1 D e − k 1 t f_0(t)=k_1De^{-k_1t} f0(t)=k1Dek1t初始条件c(0)=0,(1)的解为:
    c ( t ) = k 1 D V ( k 1 − k ) ( e − k t − e − k 1 t ) , k 1 > k ( 4 ) c(t)=\frac{k_1D}{V(k_1-k)}(e^{-kt}-e^{-k_1t}) , k_1>k (4) c(t)=V(k1k)k1D(ektek1t),k1>k4
    在这里插入图片描述

2.10导弹跟踪

背景

在发射导弹时刻(t=0),导弹位于坐标原点(0,0),飞机位于(a,b),飞机研平行于x轴的方向以常速 v 0 v_0 v0飞行。导弹在时刻t的位置为点(x,y),其速度为常值 v 1 v_1 v1,导弹在飞行过程中,按照制导系统时钟指向飞机。请确定导弹的飞行轨迹以及击中飞机所需的时间T。
在这里插入图片描述

模型建立与求解

首先建立导弹的运动方程。导弹飞行曲线在点M(x,y)处的切线方程为:
Y − y = d y d x ( X − x ) = d y d t d x d t ( X − x ) Y-y=\frac{dy}{dx}(X-x)=\frac{\frac{dy}{dt}}{\frac{dx}{dt}}(X-x) Yy=dxdy(Xx)=dtdxdtdy(Xx)
其中(x,y)为切线上动点的坐标。由于点 A ( x A , b ) A(x_A,b) A(xA,b)应位于切线上,且 x A = a + v 0 t x_A=a+v_0t xA=a+v0t,所以:
b − y = d y d t d x d t ( a + v 0 t − x ) b-y=\frac{\frac{dy}{dt}}{\frac{dx}{dt}}(a+v_0t-x) by=dtdxdtdy(a+v0tx)
从而,导弹的飞行轨迹为:
{ d x d t ( b − y ) = d y d t ( a + v 0 − x ) ( 1 ) ( d x d t ) 2 + ( d y d t ) 2 = v 1 2 ( 2 ) \small \begin{cases} \frac{dx}{dt}(b-y)=\frac{dy}{dt}(a+v_0-x) (1)\\ (\frac{dx}{dt})^2+(\frac{dy}{dt})^2=v_1^2 (2) \end{cases} {dtdx(by)=dtdy(a+v0x)1(dtdx)2+(dtdy)2=v122
由(1)可得: d x d y ( b − y ) = a + v 0 − x \frac{dx}{dy}(b-y)=a+v_0-x dydx(by)=a+v0x
两边对t求导,得: d 2 x d 2 y d y d t ( b − y ) − d x d y d y d t = v 0 − d x d t \frac{d^2x}{d^2y}\frac{dy}{dt}(b-y)-\frac{dx}{dy}\frac{dy}{dt}=v_0-\frac{dx}{dt} d2yd2xdtdy(by)dydxdtdy=v0dtdx
即: d 2 x d 2 y d y d t ( b − y ) = v 0 ( 3 ) \frac{d^2x}{d^2y}\frac{dy}{dt}(b-y)=v_0 (3) d2yd2xdtdy(by)=v03
由(2)得: ( d y d t ) 2 [ 1 + ( d x d t d y d t ) 2 ] = v 1 2 (\frac{dy}{dt})^2[1+(\frac{\frac{dx}{dt}}{\frac{dy}{dt}})^2]=v_1^2 (dtdy)2[1+(dtdydtdx)2]=v12
即: d y d t = v 1 [ 1 + ( d x d y ) 2 ] 1 2 \frac{dy}{dt}=\frac{v_1}{[1+(\frac{dx}{dy})^2]^{\frac{1}{2}}} dtdy=[1+(dydx)2]21v1
代入(3)可得导弹得运动方程:
d 2 x d 2 y ( b − y ) = λ [ 1 + ( d x d y ) 2 ] 1 2 ( 4 ) \frac{d^2x}{d^2y}(b-y)=\lambda[1+(\frac{dx}{dy})^2]^{\frac{1}{2}} (4) d2yd2x(by)=λ[1+(dydx)2]214
其中, λ = v 0 v 1 \lambda=\frac{v_0}{v_1} λ=v1v0
x ( 0 ) = 0 , x ( b ) = a + v 0 T x(0)=0,x(b)=a+v_0T x(0)=0,x(b)=a+v0T(在T时刻击中目标) (5)
接下来求(4)满足(5)的解:
p = d x d y p=\frac{dx}{dy} p=dydx,则 d p d y = d 2 x d y 2 \frac{dp}{dy}=\frac{d^2x}{dy^2} dydp=dy2d2x
(4)可化为: d p d y ( b − y ) = λ ( 1 + p 2 ) 1 2 ( 6 ) \frac{dp}{dy}(b-y)=\lambda(1+p^2)^{\frac{1}{2}} (6) dydp(by)=λ(1+p2)216
l n [ p + ( 1 + p 2 ) 1 2 ] = − λ l n ( b − y ) + c 1 ( 7 ) ln[p+(1+p^2)^{\frac{1}{2}}]=-\lambda ln(b-y)+c_1 (7) ln[p+(1+p2)21]=λln(by)+c17
(6)的初始条件为 p ( 0 ) = a b p(0)=\frac{a}{b} p(0)=ba,令 c 1 = l n ( k b λ ) c_1=ln(kb^{\lambda}) c1=ln(kbλ),则:
l n ( p + 1 + p 2 ) = l n ( b − y ) − λ + l n ( k b λ ) = l n [ k b λ ( b − y ) λ ] ln(p+\sqrt{1+p^2})=ln(b-y)^{-\lambda}+ln(kb^\lambda)=ln[\frac{kb^\lambda}{(b-y)^\lambda}] ln(p+1+p2 )=ln(by)λ+ln(kbλ)=ln[(by)λkbλ]
p + 1 + p 2 = k b λ ( b − y ) λ p+\sqrt{1+p^2}=\frac{kb^\lambda}{(b-y)^\lambda} p+1+p2 =(by)λkbλ
于时可以得到降阶方程:
d x d y = 1 2 [ k b λ ( b − y ) λ − ( b − y ) λ k b λ ] \frac{dx}{dy}=\frac{1}{2}[\frac{kb^\lambda}{(b-y)^\lambda}-\frac{(b-y)^\lambda}{kb^\lambda}] dydx=21[(by)λkbλkbλ(by)λ]
其通解为:
x = 1 2 [ k b λ ( λ − 1 ) ( b − y ) ( λ − 1 ) + ( b − y ) λ + 1 ( λ + 1 ) k b λ ] + c ( 8 ) x=\frac{1}{2}[\frac{kb^\lambda}{(\lambda-1)(b-y)^{(\lambda-1)}}+\frac{(b-y)^{\lambda+1}}{(\lambda+1)kb^\lambda}]+c (8) x=21[(λ1)(by)(λ1)kbλ+(λ+1)kbλ(by)λ+1]+c8
根据初始条件x(0)=0,可得:
c = b [ ( 1 + k 2 ) λ + k 2 − 1 ] / 2 k ( 1 − λ 2 ) c=b[(1+k^2)\lambda+k^2-1]/2k(1-\lambda^2) c=b[(1+k2)λ+k21]/2k(1λ2)
所以导弹飞行轨迹方程为: x = 1 2 [ k b λ ( λ − 1 ) ( b − y ) ( λ − 1 ) + ( b − y ) λ + 1 ( λ + 1 ) k b λ ] + b [ ( 1 + k 2 ) λ + k 2 − 1 ] / 2 k ( 1 − λ 2 ) x=\frac{1}{2}[\frac{kb^\lambda}{(\lambda-1)(b-y)^{(\lambda-1)}}+\frac{(b-y)^{\lambda+1}}{(\lambda+1)kb^\lambda}]+b[(1+k^2)\lambda+k^2-1]/2k(1-\lambda^2) x=21[(λ1)(by)(λ1)kbλ+(λ+1)kbλ(by)λ+1]+b[(1+k2)λ+k21]/2k(1λ2)
又由 x ( b ) = a + v 0 T x(b)=a+v_0T x(b)=a+v0T得到导弹集中目标的时间为:
T = c − a v 0 = a 2 + b 2 − a λ v 1 ( 1 − λ 2 ) ( 10 ) T=\frac{c-a}{v_0}=\frac{\sqrt{a^2+b^2}-a\lambda}{v_1(1-\lambda^2)} (10) T=v0ca=v1(1λ2)a2+b2 aλ10

2.11食饵-捕食者系统

一个包含两个群体的系统,其中一个群体紧密地依赖于另一个群体,成为食饵-捕食者系统。假设:x(t):t时刻食饵的数量;y(t):t时刻捕食者的数量。
如果各自独立生活,则:
{ d x d t = λ x d y d t = − μ y ( λ , μ > 0 ) \small \begin{cases} \frac{dx}{dt}=\lambda x\\ \frac{dy}{dt}=-\mu y (\lambda,\mu>0) \end{cases} {dtdx=λxdtdy=μy(λ,μ>0)
如今两者生活在一起,则有:
{ d x d t = ( λ − α y ) x ( 1 ) d y d t = − ( μ − β x ) y ( 2 ) ( α , β > 0 ) \small \begin{cases} \frac{dx}{dt}=(\lambda-\alpha y)x (1)\\ \frac{dy}{dt}=-(\mu-\beta x)y (2)(\alpha,\beta>0) \end{cases} {dtdx=(λαy)x1dtdy=(μβx)y2(α,β>0)
上式称为Volterra-Lotka方程,初始条件为 x ( 0 ) = x 0 , y ( 0 ) = y 0 x(0)=x_0,y(0)=y_0 x(0)=x0,y(0)=y0(1)/(2)可得:
d y d x = ( β x − μ ) y ( λ − α y ) x \frac{dy}{dx}=\frac{(\beta x-\mu)y}{(\lambda-\alpha y)x} dxdy=(λαy)x(βxμ)y
可得通解:
− α y − β x + λ l n y + μ l n x = l n c -\alpha y-\beta x+\lambda lny+\mu lnx=lnc αyβx+λlny+μlnx=lnc y λ e α y x μ e β x = c \frac{y^\lambda}{e^{\alpha y}}\frac{x^\mu}{e^{\beta x}}=c eαyyλeβxxμ=c
将初始条件代入,可得特解,是xoy面上的一条闭轨线
在这里插入图片描述
当食饵较多时,捕食者增多因而食饵必定减少,使得捕食者也随之减少,从而食饵又会增多。两者的数量如此起伏,周而复始,维持着生态平衡

2.12传染病模型

背景

建立传染病模型的目的是描述传染过程、分析受感染人数的变化规律、预报高潮来到的时间。
为了简单起见,假设传播期间内所观察地区人数N不变,不计生死迁移,时间以天为单位。

模型一 SI模型

模型假设
  • 1.人群分为健康者和病人,在t时刻这两类人所占比例分别为s(t),i(t),即s(t)+i(t)=1;
  • 2.平均每个病人每天接触人数是常数 λ \lambda λ,即每个病人平均每天使得 λ s ( t ) \lambda s(t) λs(t)个健康者受感染变成病人, λ \lambda λ称为日接触率
模型建立

根据模型假设2,在T时刻,每个病人每天可以使得 λ s ( t ) \lambda s(t) λs(t)个健康者变成病人,病人人属为Ni(t),故每天新增 λ N s ( t ) i ( t ) \lambda Ns(t)i(t) λNs(t)i(t)个患者,即:
N d i d t = λ N s i N\frac{di}{dt}=\lambda Nsi Ndtdi=λNsi,假设t=0时患者比例 i 0 i_0 i0,可得模型:
{ d i d t = λ i ( 1 − i ) ( 1 ) i ( 0 ) = i 0 \small \begin{cases} \frac{di}{dt}=\lambda i(1-i) (1)\\ i(0)=i_0 \end{cases} {dtdi=λi(1i)1i(0)=i0
式(1)的解为: i ( t ) = 1 1 + ( 1 t 0 e − λ t ) ( 2 ) i(t)=\frac{1}{1+(\frac{1}{t_0}e^{-\lambda t})} (2) i(t)=1+(t01eλt)12
在这里插入图片描述

模型解释
  • 1.当 i = 1 2 i=\frac{1}{2} i=21时, d i d t \frac{di}{dt} dtdi达到最大值,此时 t = t m = λ − 1 l n ( 1 i 0 − 1 ) t=t_m=\lambda^{-1}ln(\frac{1}{i_0}-1) t=tm=λ1ln(i011),也就是说,高潮到来时, λ \lambda λ越大,则 t m t_m tm越小。
  • 2.当 t → ∞ t\rightarrow\infty t时, i → 1 i\rightarrow 1 i1此时所有的人都被感染,因为SI模型没有考虑治愈病人。

模型二 SIS模型

在SI模型的基础上引入治愈,对SI模型进行修正。

模型假设
  • 1.同SI模型假设1
  • 2.同SI模型假设2
  • 3.病人每天被治愈的占病人总数的比例为 μ \mu μ,称作日治愈率。
模型修正

SI模型可修正为,t时刻每天有 N i μ Ni\mu Niμ的病人转变为健康。
{ d i d t = λ i ( 1 − i ) − μ i ( 3 ) i ( 0 ) = i 0 \small \begin{cases} \frac{di}{dt}=\lambda i(1-i)-\mu i (3)\\ i(0)=i_0 \end{cases} {dtdi=λi(1i)μi3i(0)=i0
(3)的解为:
i ( t ) = { [ λ λ − μ + ( 1 i 0 − λ λ − μ ) e − ( λ − μ ) t ] − 1 , λ ≠ μ ( λ t + 1 i 0 ) − 1 , λ = μ ( 4 ) \small i(t)= \begin{cases} [\frac{\lambda}{\lambda-\mu}+(\frac{1}{i_0}-\frac{\lambda}{\lambda-\mu})e^{-(\lambda-\mu)t}]^{-1},\lambda\not ={\mu}\\ (\lambda t+\frac{1}{i_0})^{-1},\lambda=\mu (4) \end{cases} i(t)={[λμλ+(i01λμλ)e(λμ)t]1,λ=μ(λt+i01)1,λ=μ4
由(3)可以计算出使得 d i d t \frac{di}{dt} dtdi达到最大值的高潮时刻 t m t_m tm d i d t \frac{di}{dt} dtdi的最大值 ( d i d t ) m (\frac{di}{dt})_m (dtdi)m i = λ − μ 2 λ i=\frac{\lambda-\mu}{2\lambda} i=2λλμ时达到)
a = λ μ a=\frac{\lambda}{\mu} a=μλ,可知:
i ( ∞ ) = { 1 − 1 a , a > 1 0 , a ≤ 1 \small i(\infty)= \begin{cases} 1-\frac{1}{a},a>1\\ 0,a\le 1 \end{cases} i()={1a1,a>10,a1
在这里插入图片描述

SIR模型

模型假设
  • 1.人群分为健康者、病人、移出者(病愈免疫者),三类人在t时刻在总人数N中占比例分别为s(t)、i(t)、r(t),即s(t)+i(t)+r(t)=1
  • 2.病人日接触率为 λ \lambda λ,日治愈率为 μ \mu μ,传染期间接触数 σ = λ μ \sigma=\frac{\lambda}{\mu} σ=μλ
模型建立

i(t)随t的变化规律同模型二,对于r(t):
N d r d t = μ N i , 且 d s d t + d i d t + d r d t = 0 N\frac{dr}{dt}=\mu Ni,且\frac{ds}{dt}+\frac{di}{dt}+\frac{dr}{dt}=0 Ndtdr=μNi,dtds+dtdi+dtdr=0
于时可得模型:
{ d s d t = − λ s i d i d t = λ s i − μ i ( 5 ) s ( 0 ) = s 0 , i ( 0 ) = i 0 \small \begin{cases} \frac{ds}{dt}=-\lambda si\\ \frac{di}{dt}=\lambda si-\mu i (5)\\ s(0)=s_0,i(0)=i_0 \end{cases} dtds=λsidtdi=λsiμi5s(0)=s0,i(0)=i0
从(5)中消去dt,结合 σ \sigma σ的实际意义,可得:
{ d i d s = 1 σ s − 1 ( 6 ) i ∣ s = s 0 = i 0 \small \begin{cases} \frac{di}{ds}=\frac{1}{\sigma s}-1 (6)\\ i|_{s=s_0}=i_0 \end{cases} {dsdi=σs116is=s0=i0
(6)的解为:
i = ( s 0 + i 0 ) − s + 1 σ l n s s 0 ( 7 ) i=(s_0+i_0)-s+\frac{1}{\sigma}ln\frac{s}{s_0} (7) i=(s0+i0)s+σ1lns0s7
根据(5)(7)以及图像可分析s(t),i(t),r(t)的变化规律:

  • 1.无论 s 0 , i 0 s_0,i_0 s0,i0为多少, i ∞ = 0 i_\infty=0 i=0,即病人终将消失。
  • 2.最终未被感染的健康者比例 s ∞ s_\infty s时方程 s 0 + i 0 − s ∞ + 1 σ l n s ∞ s 0 = 0 ( 8 ) s_0+i_0-s_\infty+\frac{1}{\sigma}ln\frac{s_\infty}{s_0}=0 (8) s0+i0s+σ1lns0s=08 ( 0 , 1 σ ) (0,\frac{1}{\sigma}) (0,σ1)内的单根。
  • 3.若 s 0 > 1 σ s_0>\frac{1}{\sigma} s0>σ1,则当 s = 1 σ s=\frac{1}{\sigma} s=σ1时,i(t)达到最大值 i m = s 0 + i 0 − 1 σ ( 1 + l n σ s 0 ) i_m=s_0+i_0-\frac{1}{\sigma}(1+ln\sigma s_0) im=s0+i0σ1(1+lnσs0),i(t)先增后减至0.
  • 4.若 s 0 ≤ 1 σ s_0\le\frac{1}{\sigma} s0σ1,则 i ( t ) → 0 , s ( t ) → s ∞ i(t)\rightarrow0,s(t)\rightarrow s_\infty i(t)0,s(t)s

模型解释

  • 1. 1 σ \frac{1}{\sigma} σ1是一个阈值,当 s 0 > 1 σ s_0>\frac{1}{\sigma} s0>σ1时传染病会蔓延, s 0 ≤ 1 σ s_0\le\frac{1}{\sigma} s0σ1时就不会蔓延
  • 2. σ = λ μ \sigma=\frac{\lambda}{\mu} σ=μλ表示 λ \lambda λ越小, μ \mu μ越大, σ \sigma σ也越小,从而越有利。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值