随着 COMSOL Multiphysics® 软件 6.3 版本的发布，我们现在可以在当前时步中使用瞬态研究中的上一个解，即可以使用较早时间的解影响当前时间的解。这一功能有很多应用，包括求解延迟微分方程。这篇博客，让我们通过一个与众不同的示例来了解这项新功能。
山间草地中的思索
前不久，我在山间徒步旅行时发现了一只旱獭（又名土拨鼠）。随手拍了几张照片后，我开始思考这种可爱的食草动物的生命周期。我的思绪很快就飘到了如何构建一个方程来描述在山间草地中觅食的旱獭种群数量的变化。
这只毛茸茸的旱獭是这篇博客文章的灵感来源。
回想大学时读过的一些关于数值方法的教科书，我想起了各种可能派上用场的方程，包括 Lotka–Volterra 方程和逻辑斯谛微分方程（logistic differential equation）。受这些方程的启发，我们假设应该写两个常微分方程（ODE）：一个是草地上旱獭数量的变化率 M(t)，另一个是可食用生物量的变化率 B(t)。
第一个方程描述旱獭种群变化率，其中包含与繁殖 R(t) 和死亡 D(t) 有关的项，稍后我们将详细讨论。第二个方程包含三个常数：增长率 G_0 ，草地上可生长植被的最大质量 B_{max} ，旱獭每天的消耗率 E_0。
现在，我们来看旱獭种群的方程。首先考虑与当前种群数量线性相关的繁殖项：
另一方面，旱獭的死亡率应该与可用资源有关，因此作为中间步骤，我们将定义每只旱獭每天的觅食面积：
式中，A_0 为草地面积。我们假设草场面积越大，遇到捕食者的几率就越大。同时假设每个区域的 捕食者数量固定为 P(t) = P_0。因此，旱獭死亡率与当前的旱獭数量 M(t)、每只旱獭的日觅食面积以及捕食者密度成正比。于是，我们得到一个取决于两个未知数的死亡率方程。
进一步假设整个旱獭种群都会立即从冬眠中苏醒，我们将通过求解这些方程来预测夏季种群数量的变化，之后旱獭种群会重新回到洞穴中过冬。我们还将设置一个停止条件，以便在旱獭数量或植被数量下降到临界值以下时停止求解，这意味着种群或环境崩溃。
图片显示了如何模拟一组常微分方程。不同的常微分方程可以使用不同的单位。瞬态求解器 中还添加了 停止 条件功能。
如上图所示，在 COMSOL Multiphysics® 中可以使用 全局常微分方程和微分代数方程 接口求解这类方程。该接口定义了我们要求解的两个变量 M(t) 和 B(t)，以及每个变量的时间导数方程和初始值 M_\text{init}=M(t=0) 和 B_\text{init}=B(t=0)。我们可以通过求解这些方程来预测旱獭种群在夏季的变化情况。
植被生物量（绿色）和旱獭数量（棕色）的动态方程的解，其中标注了夏季末期的旱獭数量。夏季刚开始时，旱獭数量急剧上升，随后由于捕食增加而下降。
结果表明，刚开始旱獭的数量急剧上升，并迅速消耗掉可用的食物。这导致了一个调整过程，即由于捕食的增加而导致种群数量下降。再过一段时间，旱獭数量达到一个平衡状态。(顺便给感兴趣的读者提个醒：要确定这是否是一个 稳定的 平衡状态是相当有挑战性的）。
延迟微分方程解释附加因素
目前，我们得到了一些看似合理的结果，但也会有一些担忧。至少有两个重要因素需要考虑：
- 旱獭从冬眠中醒来后就开始交配，但并不是立即就能繁殖，而是有一个妊娠期。
- 如果旱獭的数量长期增加，那么草地就会吸引更多的捕食者。
让我们来看看如何用数学方法来描述这两个因素，以及如何在模型中实现它们。针对第一个因素，我们希望将繁殖率作为受孕时种群数量的函数，并考虑妊娠期 T_g，在此期间繁殖率为零。因此，繁殖率方程变为:
式中，H(t-T_g) 为 Heaviside 阶跃函数。繁殖率还会受分式 \frac{M\left( T_g \right)}{M_\text{init}} 的影响而下降，即妊娠期内种群数量下降的总和。我们还假设，一旦年轻旱獭开始觅食，成年旱獭的捕食就会停止。因此，在 t=T_g 之后，这一比例保持不变。
现在，我们得到一个延迟微分方程，可以结合使用 if() 语句与 at(time,variable) 算子将其输入软件。繁殖项可以写为
R0*if(t > Tg,at(t-Tg,M)*at(Tg,M)/Minit,0)
at() 算子中的第一个参数即为需要计算的第二个参数的上一时步。为了在求解过程中使用参考上一时步的功能，我们需要修改求解器的设置，如下图所示。通过使用 严格的时间步进 ，并且输出时步不大于延迟时间，来确保我们表征到分娩季节的开始。否则，模型的求解方法与之前相同。
图片显示了如何在解中启用时间算子。
让我们来看看将非恒定繁殖率纳入考量时的种群数量。首先，我们观察到，在没有旱獭出生的时期，种群数量会下降。之后，种群数量开始上升并逐渐趋于平稳。
在出生率中引入延迟项时旱獭的数量（棕色）和植被的生物量（绿色）。
最后，对死亡率进行修改，以考虑如果旱獭的平均数量在足够长的时间跨度内增加，捕食者的数量也会增加的情况。也就是说，我们需要追踪过去几天内旱獭数量的时间平均值。
与其对 T_p 期间的所有时步进行积分，不如应用微积分基本定理，并写出过去几天内平均种群数量变化率的方程：
我们可以通过在模型中添加另一个全局方程来简单地求解这个方程，但将从 t=0 开始计算。因此，我们使用另一个 if() 语句来定义 M(t<0) = M_\text{init}。使用移动平均数扩大捕食者密度 P(t) = \left( P_0 /M_\text{init}\right) M_\text{ave}(t)，从而影响死亡率。再次求解，我们可以看到夏季结束时对旱獭种群的影响。
好消息！旱獭的数量增加了，它们可以在漫长的冬眠期钻进洞穴了。
仿真结果显示旱獭的数量（棕色）和生物量（绿色），其中繁殖率和死亡率都取决于上一个解。
实验观测证实旱獭数量有所增加。
关于群体平衡模型的简要说明
需要强调的是，文中所举的示例只是为了在一个简单易懂的模型中展示一些令人兴奋的新功能。实际上，真实的种群动态模型比这要复杂得多，通常使用分区或基于区间的方法。即，使用不同的变量跟踪总种群的不同部分。博客使用 COMSOL Multiphysics® 模拟 COVID-19 的传播中讨论的预测流行病的 SEIR（易感者、暴露者、感染者、恢复者和免疫者)模型就是这种方法的一个入门示例。
这种基于区间的方法在化学工程中尤为重要，通常用于跟踪液滴和沉淀物的尺寸。关于这方面的模拟，以下案例模型做了详细介绍：
注意事项和结束语
这篇博客，我们展示了在 COMSOL Multiphysics® 6.3 版本中模拟和求解延迟微分方程的基本原理。虽然这个简单的示例模型可以产生一些非常有趣的动力学，但请记住，此示例仅用于学习。只要理解了这个简单的示例，你就可以处理更复杂的建模任务了。例如，将这些时间延迟微分方程与空间微分方程结合使用。软件中的每个接口都支持这些新的算子。无论您想在多物理场仿真中使用哪种接口组合，参考上一个解的技巧都与我们在文中展示的完全相同。
想尝试自己动手模拟文中讨论的模型吗？请在 COMSOL 案例库中下载相关的 MPH 文件：
评论 (0)