贡珀茨模型

贡珀茨模型（Gompertz model）是一种描述逐渐饱和增长的微分方程模型。它让每单位现有规模的增长率，随着当前规模接近一个参考上限而下降。与洛吉斯蒂模型相比，两者都可产生趋向平台的曲线，但增长减速的方式和拐点位置不同。

从固定相对增长率到逐渐减慢

用 $N (t)$ 表示一个培养系统中的总生物量，单位为 mg，时间以日计。假设培养条件在观察期内不变，生物量可用连续正数近似，没有单独的收获或外来输入项。以下参数、曲线均为虚构教学模型，不是培养实验的实测数据。

指数增长与衰减中， $N^{'} = r N$ 使相对增长率 $N^{'} / N$ 保持恒定。要描述后期逐渐饱和，需要让它依赖当前规模。设 $K > 0$ 是平台规模，采用 $\frac{N^{'}}{N} = r \ln \frac{K}{N} ⟺ N^{'} = r N \ln \frac{K}{N}, N (0) = N_{0} > 0,$ 其中 $r > 0$ ，单位为每日； $N_{0}$ 、 $N$ 、 $K$ 使用相同质量单位，因而对数内的比值无量纲。这就是本文采用的贡珀茨方程。Peter Howard《Modeling with ODE》§4，p.11

若 $0 < N < K$ ，则 $\ln (K / N) > 0$ ，生物量增长；若 $N > K$ ，对数为负，生物量下降。 $N = K$ 是平衡点。参数 $K$ 表示这个模型在固定条件下的长期平台，不是从一个很小的初始样本就能确定的物理常量。

主例取 $K = 1000 m g, r = 0.4 d^{- 1}, N_{0} = 1000 e^{- 2} \approx 135.335 m g .$ 初始时 $\ln (K / N_{0}) = 2$ ，相对增长率为 $0.8 d^{- 1}$ ，绝对增长速率为 $0.8 N_{0} \approx 108.268$ mg/d。这里 $r = 0.4$ 不是初始相对增长率；它还要乘上当前的对数差距。

一个换元，把非线性方程变成衰减方程

令 $u (t) = \ln \frac{K}{N (t)} .$ 它度量的是当前生物量与平台之间的对数差距。由于 $K$ 固定，利用链式法则有 $u^{'} = - \frac{N^{'}}{N} = - r u .$ 这正是指数衰减方程。初值为 $u_{0} = \ln (K / N_{0})$ ，所以 $u (t) = u_{0} e^{- r t} .$ 再由 $N = K e^{- u}$ 换回原变量，得到 $N (t) = K \exp [- \ln (\frac{K}{N_{0}}) e^{- r t}] .$ 把 $t = 0$ 代入，指数变为 $- \ln (K / N_{0})$ ，恰好恢复 $N_{0}$ ；再求导，得到 $N^{'} = r N u$ ，恢复原方程。这样同时核验了初值与变化规律。

上图对数差距u从二指数降向零；下图相应生物量从约一百三十五逐渐增向一千，u等于一时生物量等于一千除e — 衰减的是对数差距，增长的是生物量；同一竖虚线表示同一时刻。

在主例中，解简化为 $N (t) = 1000 e^{- 2 e^{- 0.4 t}}$ 。1 d 和 3 d 时只需依次算内层指数、乘以 2、再算外层指数，得到约 $261.678$ mg 和 $547.502$ mg。

随着 $t \to \infty$ ， $e^{- r t} \to 0$ ，所以 $N (t) \to K$ 。若初值低于 $K$ ， $u_{0} > 0$ ，有限时刻始终有 $u (t) > 0$ ，因此曲线单调上升但不会越过 $K$ 。初值高于 $K$ 时， $u_{0} < 0$ ，同一个公式给出从上方单调下降到 $K$ 的解。平衡附近令 $N = K + ε$ ，由 $\ln (1 + ε / K) \approx ε / K$ ，可见 $ε^{'} \approx - r ε$ ：参数 $r$ 控制接近平台时的小偏差衰减速度。

生物量何时增长得最快

“生物量还在增长”与“增长速率还在增加”是两件事。令 $F (N) = r N \ln (K / N),$ 则 $F^{'} (N) = r [\ln (K / N) - 1], N^{″} = F^{'} (N) N^{'} = r^{2} N \ln (K / N) [\ln (K / N) - 1] .$ 在 $0 < N < K$ 内，前面的 $r^{2} N \ln (K / N)$ 为正，因此加速度的符号由最后一项决定：

当 $N < K / e$ ，有 $\ln (K / N) > 1$ ，增长速率继续增加；
当 $N > K / e$ 且仍小于 $K$ ，增长速率开始下降；
在 $N = K / e$ ，增长速率达到最大值 $r K / e$ ，时间曲线改变凹凸性。

主例中，令 $u (t) = 2 e^{- 0.4 t} = 1$ ，就得到拐点时刻 $t_{*} = \frac{\ln 2}{0.4} \approx 1.73287 d, N (t_{*}) = \frac{1000}{e} \approx 367.879 m g .$ 最大增长速率为 $400 / e \approx 147.152$ mg/d。此时相对增长率为 $r$ ，已经低于初始的 $2 r$ ，但总规模变大，使绝对增量达到了峰值。

以当前生物量为横轴的绝对增长速率，贡珀茨曲线峰值在K除e，洛吉斯蒂抛物线峰值在K除二，两端对应零增长 — 图横轴是生物量而非时间；峰值位置决定增长时间曲线的拐点。贡珀茨曲线在零处的值按右极限延伸。

一般地，增长解的拐点时刻是 $t_{*} = \ln u_{0} / r$ 。只有 $N_{0} < K / e$ ，即 $u_{0} > 1$ ，拐点才发生在观察起点之后。若初值已经处于 $K / e$ 与 $K$ 之间，随后全程减速增长；不能要求每一段观察曲线都呈现完整的“S”形。

达到平台的九成，需要多久

设目标质量为 $N_{*} \in (N_{0}, K)$ 。解方程前先把质量转换为对数差距： $\ln \frac{K}{N_{*}} = u_{0} e^{- r t} .$ 两边均为正数，可以再取对数，得到 $t = \frac{1}{r} \ln (\frac{u_{0}}{\ln (K / N_{*})}) .$ 要求达到 $N_{*} = 900$ mg 时， $\ln (K / N_{*}) = \ln (10 / 9) \approx 0.1053605$ ，故 $t = \frac{1}{0.4} \ln (\frac{2}{\ln (10 / 9)}) \approx 7.35879 d .$ 900 mg 可以在有限时间达到；目标若换成恰好 1000 mg，分母中的对数为零，对应无限等待时间。平台的 90%、99% 与平台本身，应分别处理。

若 $K$ 已由外部条件固定，还可把一系列严格位于 $0 < N < K$ 的观测值变换为 $z (t) = \ln \ln \frac{K}{N (t)} = \ln u_{0} - r t .$ 这是斜率为 $- r$ 的直线。在主例中， $z (0) = \ln 2$ 、 $z (2) = \ln 2 - 0.8$ ，两点的差除以 2，确实给出 $- 0.4$ 。

纵轴为log log K除N，主例在零日为log二、二日为log二减零点八，直线斜率负零点四 — 变换要求K已指定且0小于N小于K；平台附近的微小测量变化会被放大。

如果 $K$ 也未知，不能先用一个任意平台值把数据拉直，再把结果当作已知参数。尤其在 $N$ 接近 $K$ 时，变换对质量的导数为 $- 1 / [N \ln (K / N)]$ ，绝对值很大，会放大测量误差。选择模型和估计参数需要观察范围、误差机制与独立检验，不能只凭变换后的直线外观。

与洛吉斯蒂增长怎样公平比较

经典洛吉斯蒂方程为 $N^{'} = r_{L} N (1 - N / K)$ ，相对增长率随规模线性下降；贡珀茨方程则使用 $r_{G} \ln (K / N)$ 。对同样的 $K, N_{0}$ ，若令两个参数数值相同， $r_{G} = r_{L} = r$ ，它们在平台附近都满足 $(N - K)^{'} \approx - r (N - K)$ ，具有相同的局部松弛速率；但初始增长速率并不相同。

例如主例初始占比为 $x_{0} = e^{- 2}$ 。贡珀茨初始相对增长率为 $2 r = 0.8$ 每日，洛吉斯蒂初始相对增长率只有 $r (1 - e^{- 2}) \approx 0.345866$ 每日。图中的较快上升来源于不同的增长律和这个参数比较约定，不是实际数据已经证明一种模型优于另一种。

平台同为一千、初值同为一千乘e的负二次方、r同为零点四的两条增长曲线，贡珀茨先到平台附近，洛吉斯蒂较慢 — 本图固定的是同一平台、初值与平台附近松弛参数，没有固定初始斜率。

对任意 $0 < x < 1$ ，定义 $h (x) = - \ln x - (1 - x)$ ，则 $h (1) = 0$ 、 $h^{'} (x) = 1 - 1 / x < 0$ 。因此向左离开 1 后 $h (x) > 0$ ，即 $\ln (1 / x) > 1 - x$ 。所以在相同的中间规模、相同 $r, K$ 下，贡珀茨的增长速率更大。到达任意给定中间规模的时间为 $\int_{N_{0}}^{N_{*}} 1 / F (N) d N$ ；被积函数较小，便说明它更早达到目标。

若想改为“相同初始斜率”的比较，必须选择 $r_{L} = r_{G} \frac{\ln (K / N_{0})}{1 - N_{0} / K} .$ 此时两个参数不再相等，平台附近的衰减速度也不同。比较之前先指定保持哪些性质相同，才有清楚的结论。另一个不依赖这些数值选择的结构区别是：洛吉斯蒂增长速率的最大值发生在 $K / 2$ ，贡珀茨发生在 $K / e$ 。

名称、适用域与模型边界

模型以 Benjamin Gompertz 命名。他在 1825 年发表的工作研究人类死亡规律和年金计算，讨论了死亡强度随年龄增加的指数规律；不能把本文的培养量增长方程直接当成那篇论文中的实验结论。Gompertz，1825，Philosophical Transactions 115，513–583；Thomas B. L. Kirkwood 对原文的历史评述，2015

本文方程首先定义在 $N > 0$ 。当 $N ↓ 0$ ，乘积 $N \ln (K / N)$ 趋于零，因此可以把右侧连续延伸为零；但相对增长率 $r \ln (K / N)$ 却趋于无穷。极小规模时，这种“每单位规模可增长得任意快”的性质可能不合理。它不表示真实培养系统可以从零质量自行出现正质量，也提醒使用者应检查低密度阶段是否属于模型的适用范围。

如果资源随时间补充、平台随环境改变，或系统存在明显时滞、随机灭绝、收获与外来输入，就需要修改机制。解出一条光滑饱和曲线只是所给假设下的数学结果；在真实培养系统中使用它，还须有数据支持这些假设。

参考资料

Peter Howard，Modeling with ODE，2009，§4，pp.10–12：指数、洛吉斯蒂、贡珀茨种群方程的并列定义。
Vladimir Dobrushkin，Brown University，Maple Tutorial, Part 1.2: Populations，Gompertz growth：对数换元、显式解与平衡结构。
Benjamin Gompertz，On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies，1825，115：513–583：名称所关联的原始研究。
Thomas B. L. Kirkwood，Deciphering death: a commentary on Gompertz (1825)，2015，370：20140379：原始工作的历史背景。
相关条目：指数增长与衰减、洛吉斯蒂模型、微分方程、平衡点与稳定性、数学建模。